蛋白編碼基因的注釋

輸入文件的準(zhǔn)備

注釋要注意基因組數(shù)據(jù)完整性和連續(xù)性，評(píng)估組裝質(zhì)量礁哄，這決定著注釋結(jié)果的質(zhì)量和準(zhǔn)確性

需要以蛋白質(zhì)序列长酗、表達(dá)序列標(biāo)簽 (EST)、全長(zhǎng)轉(zhuǎn)錄本序列和 RNA-Seq 數(shù)據(jù)的形式等外部證據(jù)桐绒。這些與屏蔽重復(fù)序列的基因組比對(duì)夺脾，以幫助確定正確的基因結(jié)構(gòu)

方法學(xué)

計(jì)算階段

重復(fù)序列的識(shí)別

重復(fù)序列的識(shí)別可以通過RepeatModeler完成。最新版本茉继，即 RepeatModeler2咧叭，由多個(gè)不同的重復(fù)序列發(fā)現(xiàn)程序組成，即 RepeatScout烁竭、RECON菲茬、LTRharvest 和 LTR_retriever 。 RepeatScout 和 RECON 分別用于調(diào)查高度保守的 TE 和不太保守的 TE 的基因組派撕。 LTRHarvest 用于搜索長(zhǎng)末端重復(fù) (LTR) 轉(zhuǎn)座子婉弹，LTR_retriever 用于過濾來自 LTRHarvest 的假陽(yáng)性結(jié)果。

第一步：使用組裝的基因組構(gòu)建數(shù)據(jù)庫(kù)腥刹，如下所示：

~/RepeatModeler-2.0/BuildDatabase -name genome genome.fa

BuildDatabase命令采用以下參數(shù)：我們要?jiǎng)?chuàng)建的數(shù)據(jù)庫(kù)的名稱-name马胧，后跟輸入文件genome.fa

Repeat Modeler2并沒有用所有的scaffold/contig建數(shù)據(jù)庫(kù)，它會(huì)進(jìn)行采樣衔峰，因此每次運(yùn)行會(huì)有不同的結(jié)果佩脊。在這種情況下如果過濾掉短contig（如長(zhǎng)度小于1000bp的contig）蛙粘，就可以避免這些短contig被采樣。因此威彰，建議在建庫(kù)之前過濾掉短contig出牧。

創(chuàng)建數(shù)據(jù)庫(kù)后可以按如下命令運(yùn)行 RepeatModeler2

~/RepeatModeler-2.0/RepeatModeler -database genome -LTRStruct

選項(xiàng)-database表示要分析序列的數(shù)據(jù)庫(kù)名，即創(chuàng)建數(shù)據(jù)庫(kù)BuildDatabase提供的名稱歇盼。-LTRStruct讓RepeatModeler2運(yùn)行LTRHarvest 和 LTR_retreiver舔痕，并將結(jié)果與 RepeatScout 和 RECON 輸出相結(jié)合。

RepeatModeler2成功運(yùn)行后會(huì)產(chǎn)生兩個(gè)輸出文件：genome-families.fa和genome-families.stk豹缀。這兩個(gè)文件包含相同的信息伯复，即在基因組中從頭識(shí)別的共有重復(fù)序列。但是邢笙，genome-families.fa是 FASTA 格式啸如，而genome-families.stk是可以提交到 Dfam 數(shù)據(jù)庫(kù)的Stockholm-formatted文件

重復(fù)序列的屏蔽

識(shí)別出重復(fù)序列后，就可以在基因組數(shù)據(jù)上屏蔽他們氮惯。RepeatModeler2識(shí)別出的重復(fù)序列可以用RepBase數(shù)據(jù)庫(kù)補(bǔ)充叮雳。RepBase是最常用的DNA重復(fù)元件數(shù)據(jù)庫(kù)。以下代碼可以完成

perl ~/RepeatMasker/util/buildRMLibFromEMBL.pl RMRBSeqs.embl > RMRBSeqs.embl.fasta
cat genome-families.fa RMRBSeqs.embl.fasta > genome-families. RepBase.fa

然后就可以用RepeatMasker來屏蔽重復(fù)序列妇汗。以下是運(yùn)行RepeatMasker的命令

~/RepeatMasker -lib genome-families.RepBase.fa -pa 16 -xsmall genome.fa

RepeatMasker命令輸入以下文件：屏蔽重復(fù)序列所需的通過 RepeatModeler2 和 RepBase輸出的重復(fù)序列組合庫(kù) (lib) 以及需要被執(zhí)行屏蔽重復(fù)序列的組裝基因組(genome.fa)帘不。此外，-xsmall選項(xiàng)指示工具執(zhí)行softmask杨箭，-pa指定要使用的內(nèi)核數(shù)寞焙。完成后，RepeatMasker 輸出三個(gè)文件：一個(gè)屏蔽重復(fù)序列的基因組組裝文件 (genome.masked.fa)互婿， 一個(gè)包含屏蔽序列列表 的文件(genome.out)棺弊，以及一個(gè)總結(jié)序列重復(fù)內(nèi)容的文件 (genome.tbl)。

證據(jù)排列

Exonerate可以組裝蛋白序列和表達(dá)序列標(biāo)簽（EST)擒悬。如以下命令所示

~/exonerate --softmasktarget --model protein2genome --showvulgar no --showalignment no --showquerygff no --showtargetgff yes --percent 80 --ryo "AveragePercentIdentity: %pi\n" --query protein.fasta --target genome.masked.fa > genome.protein.exonerate
~/exonerate --softmasktarget --model est2genome --showvulgar no --showalignment no --showquerygff no --showtargetgff yes -percent 80 --ryo "AveragePercentIdentity: %pi\n" --query est. fasta --target genome.masked.fa > genome.est.exonerate

--softmasktarget選項(xiàng)允許對(duì)目標(biāo)序列 (genome.masked.fa) 進(jìn)行softmask。 --model選項(xiàng)的選擇取決于所使用的查詢序列的類型：protein2genome發(fā)出信號(hào) Exonerate 以查找與目標(biāo)匹配率至少 50% 的蛋白質(zhì)稻艰，而est2genome用于查找前10個(gè)匹配項(xiàng)到目標(biāo)的EST懂牧。根據(jù)所需的輸出格式設(shè)置選項(xiàng)--showvulgar、--showalignment尊勿、--showquerygff僧凤、--showtargetgff和--ryo。 --percent選項(xiàng)指的是百分比自我評(píng)分閾值元扔。

至于RNA-Seq reads躯保，可以通過三種方式處理：（1）它們可以與感興趣的基因組比對(duì)，比對(duì)后的reads作為證據(jù)提供澎语； (2) aligned reads可以組裝成轉(zhuǎn)錄本作為證據(jù)途事； (3) raw reads可以從頭組裝验懊，從頭組裝作為證據(jù)。第一種方法可以使用 Hisat2 執(zhí)行尸变，這是一種流行的拼接感知對(duì)齊器义图。下面顯示了如何使用 Hisat2 v. 2.1.0 的示例：

hisat2-build -p 16 genome.masked.fa genome.masked.fa
hisat2 -x genome.masked.fa -p 16 -1 sampleA.1.fq -2 sampleA.2. fq -S sampleA.sam

hisat2-build命令用于為屏蔽重復(fù)序列的基因組建立索引；它采用以下參數(shù)：輸入文件 genome.masked.fa召烂、要?jiǎng)?chuàng)建的索引文件的名稱（在本例中為 genome.masked.fa）和要使用的線程數(shù)（本例中為-p 16 ）碱工。然后可以使用hisat2命令執(zhí)行讀取比對(duì)，該命令采用以下參數(shù)：索引文件的基本名稱 -x奏夫、讀取對(duì)的第一個(gè)配對(duì) -1和讀取對(duì)的第二個(gè)配對(duì) -2怕篷。 -S選項(xiàng)指示 Hisat2 以 SAM 格式輸出讀取映射，然后可以使用 SAMtools將其轉(zhuǎn)換為排序的索引BAM文件酗昼。下面顯示了運(yùn)行 SAMtools v. 1.9 的示例命令

samtools view -Sb sampleA.sam | samtools sort - > sampleA. sorted.bam
samtools index sampleA.sorted.bam

已經(jīng)表明廊谓，第一種方法的基因預(yù)測(cè)在統(tǒng)計(jì)學(xué)上比涉及組裝 轉(zhuǎn)錄本組裝的方法（例如第二種和第三種方法）更穩(wěn)健。此外仔雷，第二種方法可能會(huì)導(dǎo)致在轉(zhuǎn)錄本組裝步驟中從相鄰的蹂析、間隔很近的基因中讀取 RNA-Seq 進(jìn)行合并，但denovo組裝可以避免這個(gè)問題碟婆。

組裝轉(zhuǎn)錄本與屏蔽重復(fù)序列的基因組的比對(duì)仍然可以用作注釋階段的證據(jù)來源电抚；這可以通過程序Program to Assemble Spliced Alignments (PASA)來實(shí)現(xiàn)。
首先竖共，通過“seqclean”命令對(duì)轉(zhuǎn)錄本序列（transcripts.fasta）進(jìn)行質(zhì)量修整（去除 Ns蝙叛、低質(zhì)量堿基和 polyA 尾），如下所示：

~/PASApipeline/bin/seqclean transcripts.fasta

上面產(chǎn)生了兩個(gè)文件公给，其中任何一個(gè)都可以用于下一步借帘；這些文件包括修剪后的轉(zhuǎn)錄序列 transcripts.fasta.clean)和一個(gè)制表符分隔的文件，其中包含執(zhí)行的修剪類型和序列中發(fā)生修剪的位置 transcripts.fasta.cln淌铐。然后可以按如下方式運(yùn)行PASA v.2.4.1 比對(duì)組裝流程

~/PASApipeline/Launch_PASA_pipeline.pl \
-c alignAssembly.config \
-C -R -g genome.masked.fasta \
-t all_transcripts.fasta.clean \
-T -u transcripts.fasta \
--ALIGNERS gmap

這里肺然，gmap 用于轉(zhuǎn)錄比對(duì)。 -c參數(shù)用于指定包含存儲(chǔ)輸出的 SQLite 數(shù)據(jù)庫(kù)或 MySQL 數(shù)據(jù)庫(kù)路徑的配置文件腿准。 “-C”和“-R”選項(xiàng)分別用于創(chuàng)建數(shù)據(jù)庫(kù)和運(yùn)行對(duì)齊/組裝流程际起。屏蔽重復(fù)序列程序集、干凈的轉(zhuǎn)錄本序列和原始轉(zhuǎn)錄本序列分別通過“-g”吐葱、“-t”和“-u”參數(shù)作為輸入文件提供街望。 “-T”參數(shù)向 PASA 發(fā)出信號(hào)，表明轉(zhuǎn)錄本序列已使用“seqclean”命令進(jìn)行了修剪弟跑。完成時(shí)會(huì)生成幾個(gè)文件灾前；感興趣的是 FASTA 格式 (PASA.assemblies.fasta) 和一般特征格式 3 (GFF3) 格式 (PASA.assemblies.gff3) 的 PASA 程序集，它們是注釋階段所需的孟辑。 GFF3 文件是一個(gè)制表符分隔的文件哎甲，如表 2中所述蔫敲。

基因預(yù)測(cè)

基因預(yù)測(cè)可以是denovo，也可以基于證據(jù)烧给。證據(jù)驅(qū)動(dòng)的基因預(yù)測(cè)是基因預(yù)測(cè)的首選和更準(zhǔn)確的方法燕偶。這種方法需要將外部證據(jù)來源納入基因預(yù)測(cè) 〈〉眨可用于此目的的流行工具是 BRAKER2指么。 BRAKER2通過RNA-Seq比對(duì)利用兩個(gè)單獨(dú)工具（即 GeneMark-EX 和 AUGUSTUS）的denovo基因預(yù)測(cè)功能進(jìn)行訓(xùn)練和基因結(jié)構(gòu)推斷。 GeneMark-EX 是一種強(qiáng)大的自訓(xùn)練算法榴鼎，可以學(xué)習(xí)物種特異性參數(shù)并識(shí)別真核基因組中的蛋白質(zhì)編碼基因伯诬；如果可用，它會(huì)使用外部證據(jù)來完善其模型巫财。 GeneMark-EX 識(shí)別的初始預(yù)測(cè)基因集隨后被 AUGUSTUS 用作訓(xùn)練集盗似，它通過整合外部證據(jù)產(chǎn)生最終的基因預(yù)測(cè)。當(dāng)具有良好轉(zhuǎn)錄組覆蓋率的 RNA-Seq 文庫(kù)可用于感興趣的物種時(shí)平项，可以使用以下 BRAKER2 命令：

perl ~/BRAKER/scripts/braker.pl \ --genome=genome.masked.fa \ --AUGUSTUS_BIN_PATH=~/Augustus/bin \ --AUGUSTUS_CONFIG_PATH=~/Augustus/config \ --GENEMARK_PATH=~/gm_et_linux_64 \ --BAMTOOLS_PATH=~/bamtools/bin \
--species=species1 \ --workingdir=Braker_output \ --bam=all.samples.bam \ --softmasking

“braker.pl”腳本采用以下參數(shù)：softmask程序集--genome赫舒、BAM 格式的單獨(dú)或合并的 Hisat2 映射--bam、AUGUSTUS 路徑-AUGUSTUS_BIN_PATH闽瓢、GeneMark-EX- -GENEMARK_PATH)和其他附屬程序接癌，如 Bamtools --BAMTOOLS_PATH，運(yùn)行 AUGUSTUS 所需的配置文件 --AUGUSTUS_CONFIG_PATH)和物種名稱 -species扣讼；這是 AUGUSTUS 所必需的缺猛。此外，--softmasking選項(xiàng)向 BRAKER2 發(fā)出信號(hào)椭符，表明輸入程序集已被softmask荔燎。 BRAKER2

完成后輸出三個(gè)文件：一個(gè)包含從 RNA-Seq 映射中提取的內(nèi)含子信息并分別提供給 GeneMark-EX 和 AUGUSTUS 用于基因訓(xùn)練和預(yù)測(cè)的文件 (hintsfile.gff)，一個(gè)包含 GeneMark-EX 預(yù)測(cè)的基因的文件（genemark.gtf ）和一個(gè)包含 AUGUSTUS 預(yù)測(cè)的基因的文件 (augustus.gtf)销钝。

注釋階段

一旦基因預(yù)測(cè)可用有咨，下一步就是將它們與對(duì)齊的外部證據(jù)結(jié)合起來（見注釋 9）≌艚。可以使用 EvidenceModeler (EVM)執(zhí)行此任務(wù)摔吏。 EVM 將基因預(yù)測(cè)與外部證據(jù)的比對(duì)結(jié)合到加權(quán)共識(shí)基因結(jié)構(gòu)中。 EVM v.1.1.1 需要以下輸入文件：FASTA 格式的重復(fù)掩蔽程序集纵装、GFF3 格式的證據(jù)（比對(duì)的蛋白質(zhì)和轉(zhuǎn)錄序列和基因預(yù)測(cè)）以及證據(jù)權(quán)重文件。至于權(quán)重文件据某，它包括三列橡娄，即所使用證據(jù)的類別、類型和權(quán)重癣籽； “class”表示所用證據(jù)的類型挽唉，“type”包含所用證據(jù)的來源滤祖，“weight”是應(yīng)用于每種證據(jù)類型的數(shù)值。圖 2 顯示了基于我們工作流計(jì)算階段產(chǎn)生的證據(jù)的權(quán)重文件示例瓶籽。 “ABINITIO_PREDICTION”指來自 BRAKER2 的 AUGUSTUS 和 GeneMark-EX 輸出文件匠童，“PROTEIN”指來自 Exonerate 的合并蛋白質(zhì)和 EST 比對(duì)，“TRANSCRIPT”指 PASA 組件塑顺。在這里汤求，我們使用 EVM 手冊(cè)中推薦的權(quán)重值。典型 EVM 運(yùn)行的第一步包括將輸入文件劃分為更小的數(shù)據(jù)集严拒，如下所示：

ABINITIO_PREDICTION AUGUSTUS    1
ABINITIO_PREDICTION GeneMark_hmm    1
PROTEIN exonerate   1
TRANSCRIPT  PASA    10

~/EvidenceModeler/EvmUtils/partition_EVM_inputs.pl \
--genome genome.masked.fasta \
--gene_predictions augustus.genemark.gff3 \
--protein_alignments exonerate.gff3 \
--transcript_alignments pasa.gff3 \
--segmentSize 100000 \
--overlapSize 10000 \
--partition_listing partitions_list.out

上述命令實(shí)質(zhì)上是根據(jù)相應(yīng)的 contigs 拆分 FASTA 和 GFF3 文件扬绪，較大的 contigs 進(jìn)一步拆分為較小的重疊片段。 --genome選項(xiàng)指定 FASTA 格式的重復(fù)掩蔽裝配裤唠。 --gene_predictions挤牛、--protein_alignments和--transcript_alignments文件分別指代 BRAKER2、Exonerate 和 PASA 輸出文件种蘸。 --segmentSize設(shè)置為小于 1 Mb 以減少內(nèi)存需求墓赴，--overlapSize設(shè)置為大于預(yù)期基因長(zhǎng)度的長(zhǎng)度以防止丟失完整的基因結(jié)構(gòu)。 --partition_listing參數(shù)允許 EVM 總結(jié)在此步驟中創(chuàng)建的分區(qū)列表航瞭。第二步包括創(chuàng)建一組 EVM 命令诫硕，這些命令可以在第一步中的各個(gè)分區(qū)上運(yùn)行，如下所示：

~/EvidenceModeler/EvmUtils/write_EVM_commands.pl \
--genome genome.masked.fasta \ --weights weights.txt \
--gene_predictions augustus.genemark.gff3 \
--protein_alignments exonerate.gff3 \
--transcript_alignments pasa.gff3 \
--output_file_name evm.out \
--partition partitions_list.out > commands.list

--output_file_name參數(shù)指定輸出文件名沧奴，commands.list文件包含可以運(yùn)行的命令痘括，如下所示：

~/EvidenceModeler/EvmUtils/execute_EVM_commands.pl commands. list | tee run.log

每個(gè)命令的退出值保存在“run.log”文件中；退出值為零表示 EVM 命令已成功運(yùn)行滔吠。最后纲菌，合并分區(qū)，解決冗余或不一致的預(yù)測(cè)疮绷。這可以使用以下命令完成：

~/EvidenceModeler/EvmUtils/recombine_EVM_partial_outputs.pl
--partitions partitions_list.out --output_file_name evm.out

原始evm.out文件以制表符分隔翰舌，詳細(xì)說明了完全支持每個(gè)外顯子結(jié)構(gòu)的證據(jù)集。也可以使用以下命令將此文件轉(zhuǎn)換為 GFF3 格式：

~/EvidenceModeler/EvmUtils/convert_EVM_outputs_to_GFF3.pl \
--partitions partitions_list.out \
--output evm.out \
--genome genome.masked.fasta

最后冬骚，可以用以下命令將每個(gè)分區(qū)的GFF3文件合并成一個(gè)文件椅贱。

find . -regex ".*evm.out.gff3" -exec cat {} \; > EVM.all.gff3

圖3顯示了所產(chǎn)生的GFF3文件的摘錄≈欢常可以提取FASTA格式的CDS和注釋基因的蛋白質(zhì)庇麦，如下所示。

~/EVidenceModeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl EVM.all.gff3 genome.fasta CDS > EVM.CDS.fasta ~/EVidenceModeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl EVM.all.gff3 genome.fasta prot > EVM.prot.fasta

質(zhì)控

一旦基因注釋完成喜德，檢查注釋的質(zhì)量是很重要的山橄。這個(gè)步驟不應(yīng)該被忽視，因?yàn)樗赡軐?duì)后續(xù)的分析產(chǎn)生不利的影響舍悯。例如航棱，習(xí)慣上使用一個(gè)物種的基因注釋來協(xié)助密切相關(guān)物種的注釋睡雇；如果存在不正確的注釋，它們可以從一個(gè)項(xiàng)目延續(xù)到另一個(gè)項(xiàng)目饮醇。幸運(yùn)的是它抱，有幾種方法可以檢查基因注釋的質(zhì)量。例如朴艰，可以用BUSCO v.4.0在基因組水平和基因水平上進(jìn)行兩種不同的分析观蓄，具體如下。

~/busco-master/src/busco/run_BUSCO.py -i genome.fa -m genome -o genome -l eudicotyledons_odb10 -c 16 
~/busco-master/src/busco/run_BUSCO.py -i EVM.prot.fasta -m protein -o EVM.prot -l eudicotyledons_odb10 -c 16

參數(shù)-i定義了要分析的輸入文件（基因組或蛋白質(zhì)FASTA）呵晚，參數(shù)-o定義了存儲(chǔ)輸出文件的文件夾蜘腌，參數(shù)-m指定了運(yùn)行BUSCO的模式（基因組、蛋白質(zhì)和轉(zhuǎn)錄組）饵隙。通過-l參數(shù)指定包含種群特定信息的世系文件（eudicotyledons_odb10）撮珠。使用的CPU數(shù)量是通過-c參數(shù)指定的。每次運(yùn)行的結(jié)果包括完整的金矛、重復(fù)的芯急、缺失的和片段的基因的數(shù)量和百分比∈豢。可以對(duì)兩次運(yùn)行的結(jié)果進(jìn)行比較娶耍，以確定所報(bào)告的BUSCOs的數(shù)量和百分比的任何差異。理想情況下饼酿，這種差異是不存在的或最小的榕酒。與基因組水平相比，在基因水平上有更多的缺失和破碎的BUSCO基因故俐，將表明基因注釋不完整想鹰，需要進(jìn)一步調(diào)查

另一種有用的方法是通過序列同源性對(duì)基因進(jìn)行功能注釋[1]。例如药版，鑒于同源基因在物種之間是保守的辑舷，我們期望它們具有高度相似的序列，并在不同的物種中發(fā)揮相同的作用槽片。序列同源性可以通過對(duì)注釋的基因進(jìn)行BLAST搜索來評(píng)估何缓，這些數(shù)據(jù)庫(kù)包括NCBI非冗余數(shù)據(jù)庫(kù)、SwissProt/UniProtKB和RefSeq还栓。使用BLAST v.2.9.0的例子命令顯示在下面碌廓。

blastp -db database -query EVM.prot.fasta -outfmt 6 -evalue 1e-3 -out EVM.prot.blastp

blastp命令需要以下參數(shù)：一個(gè)BLAST格式的數(shù)據(jù)庫(kù)-db，一個(gè)輸入文件-query剩盒，如何顯示結(jié)果-outfmt,在這種情況下選擇表格格式）谷婆，期望值e-value閾值-evalue和輸出文件的名稱-out。一個(gè)典型的BLAST輸出文件顯示在圖4中，表3中對(duì)各列進(jìn)行了解釋波材。

值得注意的是，有些基因在這些數(shù)據(jù)庫(kù)中可能沒有任何匹配的基因身隐；這種物種特有的基因可能是真正的生物變異的結(jié)果廷区，也可能是被錯(cuò)誤地注釋了。因此贾铝，進(jìn)行進(jìn)一步的檢查是很重要的隙轻。

一種方法是挖掘基因的蛋白域，因?yàn)檫@些蛋白域意味著基因的生化或調(diào)節(jié)功能垢揩；缺乏這種信息的基因?qū)⒈慌懦谙掠畏治鲋鈁25]玖绿。蛋白質(zhì)結(jié)構(gòu)域信息可以用InterProScan等工具推斷出來[26]。下面用InterProScan v.5.41-78.0舉例說明如何做到這一點(diǎn)叁巨。

~/interproscan-5.41-78.0/interproscan.sh -appl TIGRFAM, SUPER FAMILY,CDD,PRINTS,Pfam,PANTHER,ProDom,PrositeProfiles,PrositePatterns,Coils -i EVM.prot.fasta -o EVM.prot.interpro.out -f gff3

命令 "interproscan.sh "接受以下參數(shù)：以逗號(hào)分隔的要運(yùn)行的應(yīng)用程序列表-appl斑匪，輸入文件-i，輸出文件-o锋勺，以及輸出格式-f蚀瘸；本例為GFF3。一個(gè)典型的輸出文件顯示在圖5中庶橱。另一種方法是檢查CDS序列中是否存在起始（"ATG"）和終止（"TAG"贮勃、"TAA "和 "TGA"）密碼子；如果CDS序列中缺少這些密碼子苏章，相應(yīng)的基因可以被標(biāo)記為不完整寂嘉，并排除在下游分析之外。

本文翻譯自《Plant Bioinformatics：Methods and Protocols》 David Edwards 圖片和表格可前往查看