本發(fā)明涉及生物信息學(xué)，尤其涉及一種m6a甲基化高通量測(cè)序數(shù)據(jù)的自動(dòng)化分析方法。

背景技術(shù)：

1、表觀遺傳學(xué)被認(rèn)為是在不改變dna序列的情況下研究遺傳變化的學(xué)科，它主要關(guān)注行為和環(huán)境如何引起和影響基因功能的變化，這些變化包括dna甲基化、組蛋白修飾、染色質(zhì)結(jié)構(gòu)和重塑等。作為表觀遺傳學(xué)的一個(gè)重要研究領(lǐng)域，迄今為止，已經(jīng)在編碼和非編碼rna(ncrnas)上發(fā)現(xiàn)了170多種不同類(lèi)型的rna修飾，其中大部分都是rna甲基化。作為發(fā)生在腺苷第六個(gè)碳原子的甲基化，n6-甲基腺苷(m6a)被認(rèn)為是最豐富的mrna內(nèi)部修飾，并且在lncrnas、mirna、snrna、snorna和rrna種都廣泛存在，幾乎覆蓋了整個(gè)轉(zhuǎn)錄組。

2、m6a修飾rna是由名為“m6a?writer”的多組分甲基化轉(zhuǎn)移酶復(fù)合物、去甲基化酶“eraser”和閱讀蛋白“reader”共同作用的動(dòng)態(tài)過(guò)程。m6a甲基化修飾過(guò)程的動(dòng)態(tài)性，能夠?qū)Σ煌飳W(xué)環(huán)境和生長(zhǎng)發(fā)育階段m6a的調(diào)控作用和細(xì)胞對(duì)外界環(huán)境和內(nèi)部信號(hào)的響應(yīng)情況進(jìn)行精準(zhǔn)呈現(xiàn)，體現(xiàn)了m6a甲基化的巨大研究潛力。隨著高通量測(cè)序技術(shù)的高速發(fā)展，merip-seq技術(shù)的出現(xiàn)能夠?qū)崿F(xiàn)不同條件下全轉(zhuǎn)錄組水平的m6a甲基化修飾圖譜繪制，極大的推進(jìn)了對(duì)m6a甲基化修飾在基因表達(dá)調(diào)控、細(xì)胞功能和疾病發(fā)生中作用的認(rèn)識(shí)。

3、merip-seq技術(shù)的巨大潛力促使計(jì)算生物學(xué)家開(kāi)發(fā)了一系列的分析工具，盡管表觀轉(zhuǎn)錄組領(lǐng)域?yàn)榇_保單個(gè)工具的可用性做出了巨大努力，但由于完整的merip-seq數(shù)據(jù)挖掘及分析過(guò)程涉及多個(gè)分析步驟及頻繁的數(shù)據(jù)轉(zhuǎn)換，且每項(xiàng)分析內(nèi)容可選軟件較多，不同的軟件原理存在差異，容易影響分析結(jié)果的準(zhǔn)確性和統(tǒng)一性。同時(shí)，當(dāng)項(xiàng)目多、數(shù)據(jù)量大時(shí)，人為手動(dòng)操作在耗費(fèi)大量時(shí)間的同時(shí)還可能導(dǎo)致人為錯(cuò)誤，造成分析結(jié)果與實(shí)際項(xiàng)目不匹配。因此，對(duì)merip-seq數(shù)據(jù)生物信息分析流程進(jìn)行自動(dòng)化改造，可以降低人力資源消耗，高效自動(dòng)可靠地得到所需分析結(jié)果。

4、目前生物信息分析流程常用的方式主要通過(guò)腳本組織和基于流程框架搭建，但前者不支持范式的依賴(lài)管理及故障修復(fù)，具有開(kāi)發(fā)效率低、難以維護(hù)等缺點(diǎn)。而目前基于后者開(kāi)發(fā)的merip-seq數(shù)據(jù)自動(dòng)化分析流程，例如meripseqpipe、merip-pe等，均存在分析模塊不全面且分析功能欠缺等問(wèn)題，這會(huì)導(dǎo)致對(duì)m6a修飾圖譜的理解不完整，從而影響后續(xù)功能研究或?qū)嶒?yàn)設(shè)計(jì)。

5、因此，需要一種在全轉(zhuǎn)錄組范圍內(nèi)定位rna分子上m6a甲基化位點(diǎn)的m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化分析方法，能夠在揭示不同條件(不同組織、不同病理狀態(tài)、不同時(shí)間點(diǎn)等)下rna甲基化模式的同時(shí)關(guān)聯(lián)甲基化修飾與基因表達(dá)調(diào)控，進(jìn)一步詮釋m6a甲基化修飾在基因表達(dá)調(diào)控中的功能。

技術(shù)實(shí)現(xiàn)思路

1、發(fā)明目的：本發(fā)明提供一種m6a甲基化高通量測(cè)序數(shù)據(jù)的自動(dòng)化分析方法，基于nextflow流程框架，具有項(xiàng)目可移植性，可進(jìn)行批量項(xiàng)目分析，便于控制分析內(nèi)容。

2、技術(shù)方案：本發(fā)明所述的一種m6a甲基化高通量測(cè)序數(shù)據(jù)的自動(dòng)化分析方法，包括如下步驟：

3、步驟1、獲取待分析下機(jī)原始數(shù)據(jù)，創(chuàng)建項(xiàng)目基本信息文件，每個(gè)項(xiàng)目包含多個(gè)樣本，用戶(hù)根據(jù)項(xiàng)目具體分析任務(wù)對(duì)在配置文件中相應(yīng)數(shù)據(jù)分析基礎(chǔ)模塊參數(shù)進(jìn)行設(shè)置；

4、步驟2、創(chuàng)建conda運(yùn)行環(huán)境變量配置文件environment.yml以及docker環(huán)境構(gòu)建腳本dockerfile，將兩者置入流程配置文件nextflow.config，并在配置文件nextflow.config中進(jìn)一步對(duì)m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析流程運(yùn)行時(shí)所有輸入文件、輸出文件、數(shù)據(jù)文件、日志記錄、腳本及應(yīng)用程序的相對(duì)路徑或絕對(duì)路徑及流程運(yùn)行情況配置；

5、步驟3、設(shè)置m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析基礎(chǔ)模塊的生物信息分析腳本，m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析基礎(chǔ)模塊包括檢查子模塊、上游預(yù)處理模塊、表觀轉(zhuǎn)錄組分析模塊、聯(lián)合分析模塊、結(jié)果報(bào)告模塊；

6、步驟4、按照m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析基礎(chǔ)模塊之間的依賴(lài)關(guān)系，配置測(cè)序數(shù)據(jù)分析模塊子模塊的執(zhí)行順序和并行條件；

7、步驟5、根據(jù)nextflow.config配置文件運(yùn)行參數(shù)設(shè)定，對(duì)測(cè)序數(shù)據(jù)分析模塊子模塊運(yùn)行隊(duì)列進(jìn)行檢查，判斷是否滿(mǎn)足依賴(lài)關(guān)系；若滿(mǎn)足則執(zhí)行；若不滿(mǎn)足，則返回錯(cuò)誤；

8、步驟6、查看結(jié)果輸出和運(yùn)行日志。

9、進(jìn)一步的，步驟1中，獲取待分析下機(jī)原始數(shù)據(jù)，數(shù)據(jù)類(lèi)型為fq、壓縮格式fq.gz、bam或bed。

10、進(jìn)一步的，步驟2中，創(chuàng)建conda運(yùn)行環(huán)境變量配置文件environment.yml以及docker環(huán)境構(gòu)建腳本dockerfile，將兩者置入流程配置文件nextflow.config具體包括如下步驟：

11、步驟21、定義全局參數(shù)，包括所有輸入文件、輸出文件、數(shù)據(jù)文件、日志記錄、腳本及應(yīng)用程序的相對(duì)路徑或絕對(duì)路徑以及進(jìn)程使用內(nèi)存、cpu及運(yùn)行時(shí)間限制；

12、步驟22、將environment.yml文件路徑提供給nextflow.config，通過(guò)nextflow根據(jù)每個(gè)分析基本單元指定的依賴(lài)關(guān)系創(chuàng)建并激活conda環(huán)境；

13、步驟23、指定執(zhí)行docker容器的方式；

14、步驟24、列舉所有m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析基本單元，按照基本單元的依賴(lài)關(guān)系對(duì)其執(zhí)行情況進(jìn)行定義。

15、進(jìn)一步的，步驟3中，m6a甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化數(shù)據(jù)分析基礎(chǔ)模塊包括檢查子模塊、上游數(shù)據(jù)預(yù)處理模塊、表觀轉(zhuǎn)錄組分析模塊、轉(zhuǎn)錄組分析模塊、聯(lián)合分析模塊、結(jié)果報(bào)告模塊；檢查子模塊用于檢查運(yùn)行環(huán)境和nextflow.config中參數(shù)是否滿(mǎn)足要求，上游數(shù)據(jù)預(yù)處理模塊用于預(yù)處理原始測(cè)序數(shù)據(jù)，表觀轉(zhuǎn)錄組分析模塊用于可視化m6a甲基化在全轉(zhuǎn)錄組水平上的修飾情況，轉(zhuǎn)錄組分析模塊用于提供與表觀轉(zhuǎn)錄組信息對(duì)應(yīng)的轉(zhuǎn)錄組信息，聯(lián)合分析模塊用于對(duì)條件特異性的m6a甲基化調(diào)控基因進(jìn)行預(yù)測(cè)，進(jìn)而提供揭示m6a甲基化調(diào)控功能機(jī)制的線索，結(jié)果報(bào)告模塊用于生成并整合分析結(jié)果。

16、進(jìn)一步的，搭建上游數(shù)據(jù)預(yù)處理模塊具體包括如下步驟：

17、(1)構(gòu)建數(shù)據(jù)過(guò)濾單元，將數(shù)據(jù)過(guò)濾軟件fastp和trim?galore對(duì)原始測(cè)序數(shù)據(jù)過(guò)濾代碼模塊化，默認(rèn)參數(shù)為自動(dòng)檢測(cè)測(cè)序接頭(adapter)序列并修剪，去除平均堿基質(zhì)量得分(q-score)≤20的序列、長(zhǎng)度小于36個(gè)堿基的讀段和低復(fù)雜度序列；

18、(2)構(gòu)建核糖體去除單元，將比對(duì)軟件hisat2過(guò)濾比對(duì)到核糖體rna上的數(shù)據(jù)的代碼模塊化；

19、(3)構(gòu)建參考基因組比對(duì)單元，將比對(duì)軟件hisat2/bowtie2/star將序列比對(duì)到參考基因組的代碼模塊化，流程配置文件nextflow.config中設(shè)置默認(rèn)工具為hisat2；

20、(4)構(gòu)建比對(duì)結(jié)果統(tǒng)計(jì)單元，將對(duì)于比對(duì)結(jié)果包括reads在基因組(轉(zhuǎn)錄本)上的覆蓋情況、基因飽和度分析、reads片段的長(zhǎng)度分布情況以及reads在基因組上一些區(qū)域如cdsexon、5’utr?exon等的分布數(shù)目統(tǒng)計(jì)的代碼及腳本模塊化。

21、進(jìn)一步的，搭建表觀轉(zhuǎn)錄組分析模塊具體包括如下步驟：

22、(1)構(gòu)建單樣本m6a修飾富集區(qū)域(peak)檢測(cè)單元，將使用macs2/exomepeak2/metpeak(流程配置文件nextflow.config設(shè)置默認(rèn)工具為macs2)對(duì)轉(zhuǎn)錄組范圍內(nèi)m6a發(fā)生甲基化修飾區(qū)域進(jìn)行識(shí)別的代碼模塊化；同時(shí)，為了確保peak?calling結(jié)果的高準(zhǔn)確性，將mspc(multiple?sample?peak?calling)/bedtools(流程配置文件nextflow.config設(shè)置默認(rèn)工具為bedtools)對(duì)不同方法(macs2/exomepeak2/metpeak)識(shí)別到的peak合并構(gòu)建形成“consensus?peak”的代碼進(jìn)行模塊化，用戶(hù)需要修改配置文件nextflow.config中的相關(guān)參數(shù)以保證至少有兩個(gè)及以上的工具進(jìn)行peak?calling；

23、(2)構(gòu)建組內(nèi)重復(fù)樣本間peaks合并單元，將在組內(nèi)存在兩個(gè)及以上生物學(xué)重復(fù)情況下，通過(guò)r包c(diǎn)hipseeker對(duì)富集倍數(shù)≥2，且在至少兩個(gè)生物學(xué)重復(fù)中出現(xiàn)的peaks以overlap≥50％進(jìn)行合并的代碼模塊化；

24、(3)構(gòu)建peak相關(guān)基因分析單元，將根據(jù)peak在基因組上的位置及注釋基因組上已知元件(如基因、轉(zhuǎn)錄起始位點(diǎn)、終止密碼子、3’utr等)信息找到其在基因組特征上的分布和在其周?chē)囟ǚ秶鷥?nèi)的基因等信息，包括peak在基因功能元件上的分布、peak在染色體上的分布以及peak富集程度分布分析的代碼進(jìn)行模塊化；

25、(4)構(gòu)建組間peak合并及多樣本聚類(lèi)分析單元，將對(duì)組間peak合并，并對(duì)不同樣本在合并peak集上的信號(hào)強(qiáng)度矩陣進(jìn)行層次聚類(lèi)分析并繪制熱圖的代碼模塊化；

26、(5)構(gòu)建組間peak差異分析單元，將組間共有或特有peak統(tǒng)計(jì)注釋及可視化代碼模塊化，同時(shí)將根據(jù)每個(gè)peak相對(duì)甲基化水平對(duì)比較組進(jìn)行差異分析代碼進(jìn)行模塊化；

27、(6)構(gòu)建motif分析單元，將使用homer識(shí)別peak中顯著motif序列的代碼模塊化；

28、(7)構(gòu)建基因功能富集分析單元：將得到的peak相關(guān)基因go/kegg功能富集分析代碼模塊化。

29、進(jìn)一步的，搭建轉(zhuǎn)錄組分析模塊，具體包括如下步驟：

30、(1)構(gòu)建表達(dá)水平定量單元，將定量軟件featurecounts對(duì)基因表達(dá)水平定量代碼模塊化；

31、(2)構(gòu)建組間表達(dá)差異分析單元，根據(jù)基因表達(dá)水平對(duì)比較組進(jìn)行差異分析代碼進(jìn)行模塊化；

32、(3)構(gòu)建基因功能富集分析單元，將得到的差異表達(dá)基因go/kegg功能富集分析代碼模塊化。

33、進(jìn)一步的，搭建結(jié)果報(bào)告模塊具體包括如下步驟：

34、(1)構(gòu)建上游數(shù)據(jù)預(yù)處理結(jié)果報(bào)告單元，包括數(shù)據(jù)質(zhì)控報(bào)告及比對(duì)結(jié)果統(tǒng)計(jì)報(bào)告；

35、(2)構(gòu)建表觀轉(zhuǎn)錄組結(jié)果報(bào)告單元，包括peak統(tǒng)計(jì)結(jié)果、peak相關(guān)基因分析情況、組間peak聚類(lèi)及差異情況，peak上motif序列富集情況、peak相關(guān)基因功能富集情況；

36、(3)構(gòu)建轉(zhuǎn)錄組結(jié)果報(bào)告單元，包括基因表達(dá)信息、組間基因表達(dá)情況差異及功能富集情況；

37、(4)構(gòu)建聯(lián)合分析報(bào)告單元，包括條件性m6a調(diào)控表達(dá)基因列表及相關(guān)功能富集情況。

38、進(jìn)一步的，步驟4中，根據(jù)項(xiàng)目信息提供數(shù)據(jù)類(lèi)型判斷數(shù)據(jù)能夠執(zhí)行的分析模塊，按照分析基礎(chǔ)模塊基本單元依賴(lài)關(guān)系，將其銜接起來(lái)，在保證系統(tǒng)穩(wěn)定性的前提下，實(shí)現(xiàn)流程的靈活性最大化。

39、進(jìn)一步的，步驟6中，檢查執(zhí)行數(shù)據(jù)分析過(guò)程中，程序出現(xiàn)問(wèn)題輸出保存在日志文件中的錯(cuò)誤信息，通過(guò)查看日志文件并改正錯(cuò)誤信息進(jìn)行故障恢復(fù)，以便再次運(yùn)行快速恢復(fù)到報(bào)錯(cuò)環(huán)節(jié)。

40、有益效果：與現(xiàn)有技術(shù)相比，本發(fā)明具有如下顯著優(yōu)點(diǎn)：(1)本發(fā)明能夠?qū)崿F(xiàn)rna甲基化高通量測(cè)序數(shù)據(jù)自動(dòng)化分析，且具有可移植性和故障修復(fù)等優(yōu)點(diǎn)，有助于減輕分析人員分析壓力，便于控制分析內(nèi)容；(2)按照依賴(lài)關(guān)系組織子模塊，能夠按照配置文件參數(shù)自動(dòng)銜接子模塊，在保證系統(tǒng)穩(wěn)定性的前提下，實(shí)現(xiàn)了流程實(shí)施的靈活性最大化；(3)能夠?qū)⒎治鼋Y(jié)果以網(wǎng)頁(yè)版報(bào)告形式進(jìn)行展現(xiàn)，并設(shè)置單個(gè)項(xiàng)目結(jié)果圖表超連接，易于用戶(hù)獲取分析結(jié)果。