面向大規(guī)模人群的基因組注釋系統(tǒng)

閆貞磊，國宏哲

(哈爾濱工業(yè)大學(xué) 計(jì)算學(xué)部，哈爾濱 150001)

1 材料與方法

1.1 數(shù)據(jù)下載

注釋數(shù)據(jù)來源于多個(gè)公共數(shù)據(jù)庫，下載的數(shù)據(jù)經(jīng)過格式及參考基因組版本轉(zhuǎn)換以備注釋使用。

1.1.1 Annovar整合數(shù)據(jù)下載

Annovar數(shù)據(jù)下載方式為：annotate_variation.pl -buildver hg38 -downdb -webfrom annovar dataname humandb/。從公共數(shù)據(jù)庫annovar[1]下載了包含人群變異頻率數(shù)據(jù)、疾病變異數(shù)據(jù)、祖源數(shù)據(jù)、SNP數(shù)據(jù)、全基因組關(guān)聯(lián)分析數(shù)據(jù)、保守性數(shù)據(jù)、變異位置及功能等數(shù)據(jù)。

1.1.2 CADD數(shù)據(jù)下載

CADD數(shù)據(jù)需要下載多版本數(shù)據(jù)，目前使用的是1.6(hg38/hg19)，注意下載的數(shù)據(jù)是拆分了snv和indel的，都要下載。

從公共數(shù)據(jù)庫CADD[2]下載了包含bscores、chromhmm、dbscSNV、encode、encodeMotifDB、exons、gerp、grantham、mirSVR、mirTargetScan_v71、mmsplice、mutation_density、phastCons、phyloP、reference、regulatoryBuilt、Remap、spliceai、transcripts、vep等數(shù)據(jù)。

1.1.3 UCSC數(shù)據(jù)庫數(shù)據(jù)下載

UCSC數(shù)據(jù)庫使用FTP下載方式在HG38 goldpath路徑下載數(shù)據(jù)，同時(shí)也通過Tablebrowser下載數(shù)據(jù)。

從UCSC[3]數(shù)據(jù)庫下載的數(shù)據(jù)分為phastCons20way placental、phastCons100way vertebrate、phyloP20way placental、phyloP100way vertebrate、Recombinate Rate、repeatmasker、TAD等。

1.1.4 Fantom 5數(shù)據(jù)庫

CAGE Promoters和CAGE Enhancers 下載于Fantom 5[4]數(shù)據(jù)庫，直接選擇hg38版本的F5.hg38.enhancers.bed.gz、hg38_fair+new_CAGE_peaks_phase1and2.bed.gz文件。

1.1.5 Ensembl Regulatory數(shù)據(jù)

Ensembl Regulatory Build 下載于Ensembl[5]數(shù)據(jù)庫，提供的鏈接可能直接復(fù)制下載會(huì)失效，可以按照鏈接提供的目錄自己逐級(jí)查找，下載的文件為homo_sapiens.GRCh38.Regulatory_Build.regulatory_features.20190329.gff.

1.1.6 ORegAnno3 下載(2016)

ORegAnno3[6]存儲(chǔ)為調(diào)控元件信息數(shù)據(jù)庫，數(shù)據(jù)處理后保留所有位置信息及調(diào)控原件信息內(nèi)容。

1.1.7 LINSIGHT得分

LINSIGHT[7]是一個(gè)統(tǒng)計(jì)模型，用于估計(jì)人類基因組中非編碼序列的陰性選擇。LINSIGHT得分衡量了非編碼位點(diǎn)上的陰性選擇概率，可用于對(duì)與遺傳疾病相關(guān)的SNVs進(jìn)行優(yōu)先排序，或量化調(diào)控序列(如增強(qiáng)子或啟動(dòng)子)的進(jìn)化約束。

1.1.8 miRNA數(shù)據(jù)

miRNA信息來源于miRBase r21[8]and snoRNABase v3[9]。

1.1.9 GenCode 數(shù)據(jù)

GenCode[10]數(shù)據(jù)包含基因位置、名稱、轉(zhuǎn)錄本、類型等信息。

1.2 數(shù)據(jù)格式及參考基因組版本轉(zhuǎn)換

所有已經(jīng)下載的數(shù)據(jù)參考基因組版本為hg38的數(shù)據(jù)不需要進(jìn)行參考基因組版本轉(zhuǎn)換，將其處理成注釋可用格式即可。參考基因組為HG19的數(shù)據(jù)需要進(jìn)行數(shù)據(jù)版本轉(zhuǎn)換，版本及格式轉(zhuǎn)換完成后的數(shù)據(jù)處理為varnote可用的知識(shí)庫，以供varnote工具對(duì)變異數(shù)據(jù)進(jìn)行注釋。

1.2.1 liftover等工具進(jìn)行版本轉(zhuǎn)換

liftover的獲取地址：wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver。同時(shí)，獲取坐標(biāo)注釋文件hg38Tohg19和hg19Tohg38，且獲取命令如下：

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz (下載得到的*over.chain.gz不需要解壓)其他版本的注釋文件可在UCSC數(shù)據(jù)庫查找。

Input文件 BED格式文件，BED格式文件只定義前三列：chr start end，無表頭(注：start不等于end)。坐標(biāo)轉(zhuǎn)換命令：liftover inputfile，over.chain.gz，outputfile，unmapfile。

1.2.2 其他版本轉(zhuǎn)換工具

部分文件轉(zhuǎn)換版本的時(shí)候用的是CrossMap.py[11]，尤其是hg18及更老版本轉(zhuǎn)換到hg38的版本的時(shí)候，可以實(shí)現(xiàn)一步轉(zhuǎn)換。CrossMap.py bed *_to_*.chain.gz input.bed out.bed。

1.2.3 bw文件處理

bw文件可以用于UCSC及IGV等可視化工具展示，但是注釋的時(shí)候不適合使用，要轉(zhuǎn)換問bed文件才能進(jìn)行注釋。bw文件轉(zhuǎn)換使用的是UCSC提供的bigWigToBedGraph[12]工具，運(yùn)行命令為：bigWigToBedGraph file.bw out.bedGraph(如果是HG19的數(shù)據(jù)按照上面版本轉(zhuǎn)換方法進(jìn)行版本轉(zhuǎn)換)。

1.2.4 VCF數(shù)據(jù)預(yù)處理

所有.vcf格式變異均處理為avinput(參見annovar avinput處理)格式，知識(shí)庫及得到的變異數(shù)據(jù)文件都要進(jìn)行相應(yīng)處理，vcf格式知識(shí)庫數(shù)據(jù)還包括GWAS、GRASP和CADD等文件。

1.3 數(shù)據(jù)注釋

1.3.1 流程結(jié)構(gòu)(見圖1)

圖1 基因組注釋流程圖Fig.1 Flow chart of genome annotation

1.3.2 注釋步驟

轉(zhuǎn)換avinput：avinput文件是annovar的標(biāo)準(zhǔn)輸入文件，其自帶數(shù)據(jù)庫均依據(jù)該文件中ref/alt的格式建立，當(dāng)前注釋流程的依賴數(shù)據(jù)庫中包含了很多annovar自帶的數(shù)據(jù)庫，且無法通過avinput的ref/alt格式逆推回原本的格式。同時(shí)，avinput文件定義的ref/alt格式可以有效地將整合后的vcf文件中的ref/alt標(biāo)準(zhǔn)化，如原本的ref/alt為C/T，但整合后有可能會(huì)變?yōu)锳C/AT，顯然整合后的ref/alt并不能與數(shù)據(jù)庫中的信息匹配，故標(biāo)準(zhǔn)化是保證注釋命中的關(guān)鍵步驟。由于annovar自帶的標(biāo)準(zhǔn)化過程限定條件太多，部分不滿足要求的行會(huì)被自動(dòng)過濾掉，且這一過程無法通過參數(shù)進(jìn)行調(diào)節(jié)，故這一步驟采用了改寫annovar代碼的方式，自行編寫出了python版本的轉(zhuǎn)換代碼。

拆分同位置上多allele：由于該注釋流程獲得的上游輸入文件是多個(gè)樣本整合后的vcf，而整個(gè)項(xiàng)目最終會(huì)需要將多次收到的輸入文件進(jìn)行整體整合。對(duì)于單個(gè)收到的vcf來說，其中樣本的基因型是依照當(dāng)前的vcf給出的，這一基因型會(huì)與最終整合的vcf不匹配。同時(shí)，轉(zhuǎn)換后的avinput每一個(gè)allele都是單獨(dú)一行，若不進(jìn)行拆分則不利于后續(xù)的整合步驟。

1.3.3 annovar注釋refGene

refGene是基于基因的注釋，varnote無法注釋該類型的數(shù)據(jù)，故采用annovar對(duì)該數(shù)據(jù)庫進(jìn)行注釋。

1.3.4 varnote注釋

生成varnote配置文件：基于注釋框架中指定的注釋庫、項(xiàng)目、輸出項(xiàng)目名自動(dòng)化生成varnote使用的配置文件。

通過avinput建立數(shù)據(jù)庫：varnote進(jìn)行注釋時(shí)，對(duì)于完全注釋不到的行，不會(huì)在結(jié)果文件中進(jìn)行輸出，這樣的輸出結(jié)果不符合實(shí)際需求以及后續(xù)整合的需要，且這一行為無法更改。故將轉(zhuǎn)換后的avinput文件建立為一個(gè)數(shù)據(jù)庫，通過自己注釋自己的方式，使得每一行信息都能得到注釋，從而保證結(jié)果文件中不會(huì)丟行。

注釋：通過讀取流程框架中配置文件的數(shù)據(jù)庫信息、參數(shù)信息拼接出varnote的實(shí)際執(zhí)行語句進(jìn)行調(diào)用注釋。

1.4 整合注釋結(jié)果

1.4.1 去除無用數(shù)據(jù)列

無用數(shù)據(jù)列包括原始vcf文件中包含的部分信息、avinput數(shù)據(jù)庫產(chǎn)生的額外信息等。在最后的整合過程中予以剔除。

1.4.2 重新對(duì)數(shù)據(jù)列排序

根據(jù)配置文件中指定的數(shù)據(jù)列排列順序，將注釋結(jié)果的信息列進(jìn)行重新排序。

1.4.3 整合refGene注釋結(jié)果

將annovar單獨(dú)注釋出的refGene信息整合到最終結(jié)果文件中。

1.4.4 整合樣本基因型

將拆分后的多allele樣本基因型合并到注釋結(jié)果中，形成最終結(jié)果文件。

1.5 AnnotSV注釋

1.5.1 注釋各類型SV變異

結(jié)構(gòu)變異注釋的輸入文件是根據(jù)SV類型產(chǎn)生的多個(gè)vcf文件，部分類型的文件可能不存在或內(nèi)容為空。注釋時(shí)給定的輸入文件參數(shù)并非是實(shí)際的文件名，而是前綴名，經(jīng)過拼接后自動(dòng)查找目錄中是否有相應(yīng)類型的文件，逐個(gè)調(diào)用annotSV進(jìn)行注釋。

1.5.2 整合SV注釋結(jié)果

由于輸入文件多個(gè)文件，產(chǎn)生的結(jié)果文件同樣是多個(gè)，且無法保證每種SV類型的結(jié)果文件都存在。對(duì)各個(gè)類型的結(jié)果文件進(jìn)行檢驗(yàn)，查看其是否存在，將存在的文件進(jìn)行合并拼接，得到最后的完整SV注釋結(jié)果。

2 結(jié)果與分析

2.1 注釋數(shù)據(jù)統(tǒng)計(jì)

2.1.1 人群變異頻率數(shù)據(jù)

為了進(jìn)行注釋收集整理了多類變異頻率數(shù)據(jù)，其中包含人群變異頻率數(shù)據(jù)，例如：the Exome Aggregation Consortium(外顯子組整合數(shù)據(jù)庫)、Genome Aggregation Database[13](簡稱gnomAD)是由各國研究者聯(lián)合發(fā)展起來的基因組突變頻率數(shù)據(jù)庫，匯集眾多大規(guī)模測序計(jì)劃的基因組數(shù)據(jù)、NHLBI 計(jì)劃匯集的外顯子組突變頻率數(shù)據(jù)、1000G 計(jì)劃匯集的全基因組突變頻率數(shù)據(jù)、KAVIAR[14]整合了HGP、1000G、CGI69、UK10K等共35個(gè)計(jì)劃的基因突變頻率數(shù)據(jù)、HRCR整合了32 k個(gè)樣本整合的基因組突變頻率數(shù)據(jù)、巴西人群外顯子組突變頻率數(shù)據(jù)[15]、中東地區(qū)人群外顯子組變異頻率數(shù)據(jù)、以及收集整理的中國人群變異頻率數(shù)據(jù)。我們對(duì)人群變異頻率數(shù)據(jù)的總數(shù)據(jù)量、SNP/INDEL數(shù)量、稀有突變、低頻突變、高頻突變的數(shù)量進(jìn)行了統(tǒng)計(jì)，各人群數(shù)據(jù)庫內(nèi)稀有變異比例最高，SNP數(shù)量大于INDEL數(shù)據(jù)量，位點(diǎn)數(shù)量最多的數(shù)據(jù)庫為gnomad，近3億個(gè)變異位點(diǎn)，其他全基因組數(shù)據(jù)庫包含變異位點(diǎn)數(shù)量大多小于1億個(gè)。統(tǒng)計(jì)結(jié)果(見圖2)。

圖2 人群變異頻率數(shù)據(jù)Fig.2 Illustration of frequency of population genome variation

2.1.2 臨床及表型注釋數(shù)據(jù)

為了了解變異是否是已知的致病變異我們整理了多個(gè)DNA變異和人類疾病關(guān)系知識(shí)庫，包括：60種人類腫瘤細(xì)胞系外顯子組變異頻率數(shù)據(jù)、預(yù)測內(nèi)含子單核苷酸變異的致病性影響數(shù)據(jù)庫、NCBI維護(hù)的疾病相關(guān)的人類基因組變異數(shù)據(jù)庫、EBI負(fù)責(zé)維護(hù)的一個(gè)收集已發(fā)表的GWAS研究的數(shù)據(jù)庫[16]、GRASP 2.0[17](SNP和表型關(guān)聯(lián)數(shù)據(jù)庫)、癌癥基因突變數(shù)據(jù)庫[18]，此外還有基于美國醫(yī)學(xué)遺傳學(xué)與基因組學(xué)學(xué)會(huì)(The American College of Medical Genetics and Genomics，ACMG)曾制定過序列變異解讀指南而分析變異數(shù)據(jù)的數(shù)據(jù)庫-InterVar[19]。我們統(tǒng)計(jì)了這些疾病及表型相關(guān)的數(shù)據(jù)庫內(nèi)的變異數(shù)量及分類情況，大多數(shù)變異數(shù)據(jù)和臨床關(guān)系 為“未明確”或者“可能致病/良性”，明確的“良性”及“致病性”位點(diǎn)比例較低，InterVar和clinvar統(tǒng)計(jì)結(jié)果(見圖3和圖4)。

圖3 InterVar臨床變異數(shù)據(jù)統(tǒng)計(jì)Fig.3 Statistics of clinical variation based on the InterVar database

圖4 Clinvar臨床變異數(shù)據(jù)統(tǒng)計(jì)Fig.4 Statistics of clinical variation based on the clinvar database

2.1.3 序列保守性數(shù)據(jù)庫

收集整理了多個(gè)評(píng)估序列保守性的數(shù)據(jù)庫以及對(duì)相應(yīng)保守性的得分信息，比如靈長類動(dòng)物基因組序列PhyloP[20]方法保守性打分?jǐn)?shù)據(jù)、哺乳動(dòng)物基因組序列PhyloP方法保守性打分?jǐn)?shù)據(jù)、脊椎動(dòng)物基因組序列PhyloP方法保守性打分?jǐn)?shù)據(jù)、靈長類動(dòng)物基因組序列PhastCons方法保守性打分?jǐn)?shù)據(jù)、哺乳動(dòng)物基因組序列PhastCons方法保守性打分?jǐn)?shù)據(jù)、脊椎動(dòng)物基因組序列PhastCons方法保守性打分?jǐn)?shù)據(jù)、由GERP++[21]方法計(jì)算得到的GerpN值等，同時(shí)還下載整理了祖源等位基因情況。以上的保守性打分?jǐn)?shù)據(jù)對(duì)基因組上所有的位點(diǎn)進(jìn)行打分，每個(gè)位點(diǎn)都有相應(yīng)的分值。

2.1.4 表觀注釋數(shù)據(jù)

為了對(duì)變異區(qū)域進(jìn)行表觀注釋，我們下載了多個(gè)相關(guān)的數(shù)據(jù)庫，其中包括識(shí)別promoter和enhancer區(qū)域的數(shù)據(jù)、細(xì)胞中染色質(zhì)標(biāo)記狀態(tài)比例數(shù)據(jù)、預(yù)測得到的轉(zhuǎn)錄因子motif數(shù)量的數(shù)據(jù)、與轉(zhuǎn)錄起始和終止位置的距離數(shù)據(jù)、表觀遺傳學(xué)標(biāo)記和檢測相應(yīng)轉(zhuǎn)錄調(diào)控原件的數(shù)據(jù)、染色質(zhì)拓?fù)浣Y(jié)構(gòu)域識(shí)別的數(shù)據(jù)等。

2.1.5 變異綜合得分

為確定每個(gè)變異的功能和致病性風(fēng)險(xiǎn)程度，多個(gè)數(shù)據(jù)庫基于多種類的注釋結(jié)果對(duì)變異數(shù)據(jù)進(jìn)行了綜合打分，這些數(shù)據(jù)庫包括了LINSIGHT得分，LINSIGHT用于估計(jì)人類基因組中非編碼序列的陰性選擇。LINSIGHT得分衡量了非編碼位點(diǎn)上的陰性選擇概率，可用于對(duì)與遺傳疾病相關(guān)的SNVs進(jìn)行優(yōu)先排序，或量化調(diào)控序列(如增強(qiáng)子或啟動(dòng)子)的進(jìn)化約束。還包含了在編碼區(qū)和非編碼區(qū)對(duì)遺傳變異進(jìn)行注釋的Eigen[22]得分、CADD得分、fitcons[23]得分等。各數(shù)據(jù)庫包含位點(diǎn)數(shù)量(見表1)。

表1 變異綜合得分總數(shù)Table 1 Total number of variants from multiple databases

2.1.6 變異基礎(chǔ)信息

變異區(qū)域的GC比例，以及變異所在的區(qū)域所在的功能區(qū)，變異在相應(yīng)的功能區(qū)域的相對(duì)位置，例如是否位于promoter區(qū)域，是否在編碼區(qū)域，在第幾個(gè)外顯子上，在第幾個(gè)內(nèi)含子上，這些變異的信息相關(guān)的數(shù)據(jù)都進(jìn)行了整理，以及用dbSNP(150版本)數(shù)據(jù)對(duì)每個(gè)位置的變異進(jìn)行了注釋。這些信息覆蓋整個(gè)基因組。

2.1.7 蛋白功能數(shù)據(jù)

基因變異會(huì)對(duì)蛋白功能造成相應(yīng)改變，基于氨基酸序列的同源性和物理性質(zhì)來預(yù)測氨基酸的替換對(duì)蛋白質(zhì)功能是否造成影響，用來評(píng)估基因變異的有害程度的數(shù)據(jù)有SIFT[26]得分及polyphen[27]得分，除了這兩種打分?jǐn)?shù)據(jù)我們還整理了MA[28]、LRT[29]等蛋白功能預(yù)測數(shù)據(jù)對(duì)基因組變異進(jìn)行打分。

2.1.8 miRNA數(shù)據(jù)

MicroRNA (miRNA) 是一類內(nèi)生的、長度約為20-24個(gè)核苷酸的小RNA，其在細(xì)胞內(nèi)具有多種重要的調(diào)節(jié)作用。每個(gè)miRNA可以有多個(gè)靶基因，而幾個(gè)miRNA也可以調(diào)節(jié)同一個(gè)基因。miRNA的基本信息及靶向關(guān)系的信息來源于miRBase和TragetScan[30]數(shù)據(jù)庫等數(shù)據(jù)庫，各庫包含的數(shù)據(jù)量(見表2)。

表2 miRNA數(shù)據(jù)量統(tǒng)計(jì)Table 2 Statistics of miRNA based on multiple databases

2.1.9 結(jié)構(gòu)變異注釋數(shù)據(jù)

結(jié)構(gòu)變異注釋使用的注釋工具室AnnotSV[31]，注釋數(shù)據(jù)使用的是軟件自帶的，注釋項(xiàng)目包括：變異的基本信息，比如變異所在的轉(zhuǎn)錄本，所在的CDS長度，變異長度等信息，還提供變異位置其他人群機(jī)構(gòu)變異頻率的注釋，比如DGV、gnomAD、1000G、IMH數(shù)據(jù)庫。變異致病性注釋使用的是dbVAR、ClinGen、clinvar等數(shù)據(jù)。結(jié)合以上所有信息對(duì)相應(yīng)的結(jié)構(gòu)變異進(jìn)行打分，得出致病性結(jié)論。

2.2 注釋流程涉及軟件及工具

2.2.1 LiftOver

將所有的數(shù)據(jù)轉(zhuǎn)換為同一個(gè)參考基因組之后進(jìn)行注釋是一個(gè)必經(jīng)的步驟，LiftOver是一個(gè)高效的位置信息轉(zhuǎn)換工具，在我們使用工具進(jìn)行轉(zhuǎn)換的時(shí)候，成功轉(zhuǎn)換的比例大于99.97%，2 000萬條變異的數(shù)據(jù)進(jìn)行基因組位置轉(zhuǎn)換在五分鐘內(nèi)完成，注釋消耗時(shí)間與目標(biāo)文件的大小成正比。

2.2.2 Annovar

Annovar[1]是一個(gè)高效的軟件工具，可以利用最新的信息對(duì)從不同基因組(包括人類基因組hg18、hg19、hg38，以及小鼠、蠕蟲、蒼蠅,酵母和許多其他基因組)檢測到的遺傳變異進(jìn)行功能注釋。給出帶有染色體、起始位置、結(jié)束位置、參考核苷酸和觀察到的核苷酸的變異列表。其功能包含基于基因的注釋、基于區(qū)域的注釋、基于過濾的注釋及其它輔助功能。

在工作流中，Annovar負(fù)責(zé)對(duì)樣本VCF進(jìn)行基于基因注釋，主要依賴于refGene數(shù)據(jù)庫。在測試過程中，針對(duì)500萬行的VCF文件，完成注釋工作需要耗時(shí)約8分鐘，針對(duì)2 500萬行的VCF文件，完成注釋工作需要耗時(shí)約40分鐘。可見ANNOVAR的工作效率是比較穩(wěn)定的，注釋消耗時(shí)間與目標(biāo)文件的大小成正比。

2.2.3 VARNOTE

VARNOTE[32]是ANNOVAR的作者最新開發(fā)的注釋工具，該工具主要針對(duì)基于過濾及基于區(qū)域的注釋，由于算法上的優(yōu)化，其注釋效率較ANNOVAR有極大的提升，且注釋過程中不再需要將依賴的數(shù)據(jù)庫完整讀入內(nèi)存中，有效地降低了對(duì)機(jī)群資源的占用。使用VARNOTE進(jìn)行注釋需要自行構(gòu)建依賴的數(shù)據(jù)庫，這一過程由于需要進(jìn)行排序、壓縮、建立索引等過程，相較于ANNOVAR比較繁瑣，但VARNOTE的效率提升主要便是依托于對(duì)依賴數(shù)據(jù)庫的處理，故這一過程必不可少。

在工作流中，VARNOTE負(fù)責(zé)對(duì)樣本VCF進(jìn)行所有的基于過濾以及區(qū)域的注釋工作，當(dāng)前包含的數(shù)據(jù)庫共56個(gè)。在測試過程中，針對(duì)500萬行的無INDEL位點(diǎn)VCF文件，VARNOTE對(duì)其進(jìn)行56個(gè)數(shù)據(jù)庫的注釋僅需要3分鐘。但在實(shí)際工作過程中，VARNOTE的注釋效率波動(dòng)較大。由于其算法原因，VARNOTE對(duì)INDEL位點(diǎn)的注釋效率較低，當(dāng)目標(biāo)VCF文件中INDEL位點(diǎn)的占比較大時(shí)，其注釋效率會(huì)有比較明顯的下降。但整體而言，VARNOTE的注釋效率仍然是比較理想的。

2.2.4 AnnotSV

AnnotSV是一款使用Tcl語言編寫的用于注釋以及評(píng)估結(jié)構(gòu)變異的軟件，其目的在于發(fā)現(xiàn)結(jié)構(gòu)變異的潛在致病性以及過濾出假陽性結(jié)構(gòu)變異。該軟件接受VCF格式以及BED格式的文件作為輸入，以tab分割的文本文檔作為輸出文件。AnnotSV在注釋時(shí)需要使用各種不同的參考數(shù)據(jù)庫來進(jìn)行支持，這部分不需要用戶自行處理，在注釋過程中若本地不存在相應(yīng)的數(shù)據(jù)庫，軟件會(huì)自動(dòng)進(jìn)行下載操作。

在工作流中，AnnotSV負(fù)責(zé)對(duì)結(jié)構(gòu)變異數(shù)據(jù)進(jìn)行注釋工作。在測試過程中，針對(duì)27 000行的結(jié)構(gòu)變異VCF文件，完成注釋工作需要13分鐘，針對(duì)49 000行的結(jié)構(gòu)變異VCF文件，完成注釋工作需要25分鐘。注釋消耗時(shí)間與目標(biāo)文件的大小大體上符合線性規(guī)律。

2.3 變異注釋結(jié)果統(tǒng)計(jì)

對(duì)注釋后的樣例數(shù)據(jù)結(jié)果進(jìn)行統(tǒng)計(jì)，發(fā)現(xiàn)在示例樣本中變異類型中大多數(shù)是單堿基替換變異，插入和缺失變異數(shù)量相近，并且遠(yuǎn)少于單堿基替換位點(diǎn)數(shù)量，樣例數(shù)據(jù)中位點(diǎn)變異數(shù)量超2 000萬個(gè)。多數(shù)變異位點(diǎn)位于基因間區(qū)及內(nèi)含子區(qū)域。因?yàn)闃永龜?shù)據(jù)中單堿基替換位點(diǎn)比例較高，蛋白結(jié)構(gòu)變異類型中同義突變及非同義突變數(shù)量最多，非同義突變數(shù)量大于同義突變。樣例數(shù)據(jù)涉及變異類型統(tǒng)計(jì)結(jié)果(見圖5)，樣例數(shù)據(jù)變異位點(diǎn)區(qū)域分布(見圖6)。

圖5 樣例數(shù)據(jù)變異類型數(shù)量Fig.5 Statistics of multiple types of variation

因?yàn)闃永龜?shù)據(jù)中單堿基替換位點(diǎn)比例較高，蛋白結(jié)構(gòu)變異類型中同義突變及非同義突變數(shù)量最多，遠(yuǎn)超移碼突變，非同義突變數(shù)量大于同義突變，而且占總突變數(shù)量的一半以上。移碼突變中，不是3或3的倍數(shù)堿基缺失或增加的移碼突變(編碼氨基酸序列改變)數(shù)量少于其他移碼突變數(shù)量。突變導(dǎo)致終止密碼子獲得或缺失在所有變異類型內(nèi)比例最低。編碼區(qū)蛋白結(jié)構(gòu)變異類型比例(見圖7)。

圖6 樣例數(shù)據(jù)變異所在區(qū)域分布比例Fig.6 Statistics of regional distribution of genomic variation

2.3.1 轉(zhuǎn)換和顛換

堿基顛換(transversion)是指在堿基置換中嘌呤與嘧啶之間的替代，而轉(zhuǎn)換(transition)則是一個(gè)嘌呤被另一個(gè)嘌呤，或者是一個(gè)嘧啶被另一個(gè)嘧啶替代。在樣例數(shù)據(jù)中，轉(zhuǎn)換變異數(shù)量是顛換變異數(shù)量兩倍有余(約2.1倍)，其中G>T和C>A堿基變異比例最高，均超過了全部變異數(shù)量的20%，T>G和A>C次之，均超過全部變異數(shù)量的10%，各種顛換變異的堿基比例均較低。各堿基變異情況統(tǒng)計(jì)結(jié)果(見圖8)。

圖7 樣例數(shù)據(jù)編碼區(qū)蛋白結(jié)構(gòu)變異類型比例Fig.7 Statistics of ratio of protein structural variationin encoding region

圖8 轉(zhuǎn)換及顛換變異比例Fig.8 Statistics of ratio of transition and transversion variations

2.3.2 樣本變異數(shù)量分布

樣例數(shù)據(jù)中每個(gè)樣本的變異數(shù)量均有所差異，經(jīng)過我們統(tǒng)計(jì)單樣本變異數(shù)量在470萬到500萬區(qū)間，樣本變異數(shù)量在485萬附近的樣本最多，樣本變異數(shù)量分布(見圖9)。

圖9 樣本變異數(shù)量分布Fig.9 Distribution of variation quantity

2.3.3 結(jié)構(gòu)變異數(shù)量分布

每個(gè)樣本中有多種類型的結(jié)構(gòu)變異，樣例數(shù)據(jù)的結(jié)構(gòu)變異數(shù)量在7 000到1.2萬之間，大部分樣本內(nèi)的結(jié)構(gòu)變異數(shù)量在9 000個(gè)左右，樣例數(shù)據(jù)結(jié)構(gòu)變異數(shù)量分布(見圖10)。

圖10 樣本結(jié)構(gòu)變異數(shù)量分布Fig.10 Distribution of the number of structural variations

3 結(jié) 論

1)基因組注釋是確定基因組序列中基因及編碼區(qū)位置的過程，其可以直接對(duì)序列賦予生物學(xué)意義，故基因組注釋是解讀基因組序列重要途徑?；虮磉_(dá)及變異對(duì)其生物學(xué)功能發(fā)揮具有決定性意義，通過基因組注釋，可以確定基因調(diào)控區(qū)域及編碼區(qū)域是否變異，以及變異的影響，從而幫助研究者從生物功能的角度理解基因組。個(gè)體表型的差異究其根本是由于遺傳物質(zhì)的差異造成的，通過基因組注釋可以幫助研究者發(fā)現(xiàn)個(gè)體差異的遺傳本質(zhì)，發(fā)現(xiàn)群體的遺傳特征。此外，通過對(duì)基因組進(jìn)行注釋，結(jié)合已有變異與表型的關(guān)聯(lián)知識(shí)，可以評(píng)估個(gè)體的健康情況并指導(dǎo)個(gè)性化用藥。

2)我們研究開發(fā)的面向大規(guī)模人群的自動(dòng)化基因組注釋系統(tǒng)利用國際具有高置信度的人類基因組變異數(shù)據(jù)庫，綜合運(yùn)用不同注釋方法和高效的索引技術(shù)，實(shí)現(xiàn)了對(duì)個(gè)人基因組和群體基因組的高效、快速，準(zhǔn)確和全面的基因組信息注釋。注釋系統(tǒng)對(duì)每一個(gè)變異位點(diǎn)實(shí)現(xiàn)了基因組區(qū)域類型、變異類型、基因類型等基礎(chǔ)基因組相關(guān)信息標(biāo)注。同時(shí)，系統(tǒng)對(duì)每一個(gè)變異位點(diǎn)進(jìn)行綜合功能性分析，例如：基因組轉(zhuǎn)錄、甲基化、功能損失以及各類疾病的風(fēng)險(xiǎn)預(yù)測等分析。該注釋系統(tǒng)計(jì)算并展示個(gè)人基因組的差異化特征，助力預(yù)測個(gè)體基因組變異的相關(guān)疾病的發(fā)病風(fēng)險(xiǎn)。同時(shí)，該系統(tǒng)可應(yīng)用大規(guī)模人群基因組數(shù)據(jù)分析中，提供全面且準(zhǔn)確的基因組變異數(shù)據(jù)的注釋服務(wù)，為不同地區(qū)和民族的人群遺傳學(xué)分析提供了堅(jiān)實(shí)基礎(chǔ)。