一、注釋策略
進(jìn)行注釋的三種策略：
- 從頭注釋(de novo prediction)：通過已有的概率模型來預(yù)測基因結(jié)構(gòu)硼身，在預(yù)測剪切位點和UTR區(qū)準(zhǔn)確性較低
- 同源預(yù)測(homology-based prediction)：有一些基因蛋白在相近物種間的保守型搞硅急，所以可以使用已有的高質(zhì)量近緣物種注釋信息通過序列聯(lián)配的方式確定外顯子邊界和剪切位點
- 基于轉(zhuǎn)錄組預(yù)測(transcriptome-based prediction)：通過物種的RNA-seq數(shù)據(jù)輔助注釋，能夠較為準(zhǔn)確的確定剪切位點和外顯子區(qū)域佳遂。
每一種方法都有自己的優(yōu)缺點营袜，所以最后需要用EvidenceModeler(EVM)和GLEAN工具進(jìn)行整合，合并成完整的基因結(jié)構(gòu)丑罪〖园澹基于可靠的基因結(jié)構(gòu)凤壁，后續(xù)可才是功能注釋，蛋白功能域注釋跪另，基因本體論注釋拧抖，通路注釋等
基因結(jié)構(gòu)注釋應(yīng)是功能注釋的先決條件，完整的真核生物基因組注釋流程需要如下步驟：
- 必要的基因組重復(fù)序列屏蔽
- 從頭尋找基因, 可用工具為: GeneMarkHMM, FGENESH, Augustus, SNAP, GlimmerHMM, Genscan免绿；
- 同源蛋白預(yù)測, 內(nèi)含子分析: GeneWIse, Exonerate, GenomeThreader,TransDecoder唧席；
- 將EST序列，全長cDNA序列和Trinity/Cufflinks/Stringtie組裝的轉(zhuǎn)錄組和基因組聯(lián)配嘲驾；
- 如果第4步用到了多個數(shù)據(jù)來源淌哟，使用PASA基于重疊情況進(jìn)行聯(lián)配；
- 使用EvidenceModler根據(jù)上述結(jié)果進(jìn)行整合辽故；
- 使用PASA更新EVM的一致性預(yù)測绞绒，增加UTR注釋和可變剪切注釋；
- 必要的人工檢查榕暇。
二、從頭預(yù)測基因 ab initio
文件準(zhǔn)備：在水稻網(wǎng)站http://rice.hzau.edu.cn/rice_rs3/下載MH63RS3版本的水稻基因組作為練手?jǐn)?shù)據(jù)喻杈。
1.Augustus
安裝（使用conda安裝）
source activate genome_anno
conda install augustus=3.3
conda install -c bioconda augustus -y
使用
mkdir 01-augustsus && cd 01-augustsus
cp MH63RS3.fasta.gz ../01augustus/
mv MH63RS3.fasta genome.fa
seqkit split genome.fa #結(jié)果文件在genome.fa.split文件夾里彤枢。split是按照指定方式將染色體拆分成小部分。
find genome.fa.split/ -type f -name "*.fa" | parallel -j 30 augustus --species=rice --gff3=on >> test.gff #并行處理 ###注意選擇物種
join_aug_pred.pl < test.gff | grep -v '^#' > test.joined.gff
bedtools sort -i test.joined.gff > augustus.gff
2.GeneMark-ES/ET
安裝
#GeneMark-ES/ET是唯一一款支持無監(jiān)督訓(xùn)練模型, 軟件下載需要登記筒饰。
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_HbJ0A/gmes_linux_64.tar.gz
tar xf gm_et_linux_64.tar.gz
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_AoSqY/gm_key_64.gz
gzip -dc gm_key_64.gz > ~/.gm_key
cpan YAML Hash::Merge Logger::Simple Parallel::ForkManager
echo "export PATH=$PATH:$HOME/Tools/GeneMark/gmes_linux_64" >> ~/.bashrc
使用
gmes_petap.pl --ES --sequence genome.fa --cores 50
最后得到的是genemark.gtf缴啡，是標(biāo)準(zhǔn)的GTF格式，可以使用Sequence Ontology Project提供的gtf2gff3.pl進(jìn)行轉(zhuǎn)換瓷们。
兩種將gtf轉(zhuǎn)為gff格式的方法业栅。
wget http://genes.mit.edu/burgelab/miso/scripts/gtf2gff3.pl
chmod 755 gtf2gff3.pl
perl gtf2gff3.pl /database/genome_assembly/1_03_GeneMark/genemark.gtf | bedtools sort -i - > genemark.gff
gffread genemark.gtf -o- > genemark.gff
3.GlimmerHMM
安裝
wget https://ccb.jhu.edu/software/glimmerhmm/dl/GlimmerHMM-3.0.4.tar.gz
tar -xzf GlimmerHMM-3.0.4.tar.gz
使用
1.用訓(xùn)練好的水稻數(shù)據(jù)集。
glimmerhmm_linux_x86_64 genome.fa -d $HOME/Tools/GlimmerHMM/GlimmerHMM/trained_dir/rice -g -n 1 -o glimm.gff
#參數(shù)
USAGE: ./glimmerhmm_linux_x86_64 <genome1-file> <training-dir-for-genome1> [options]
Options:
-p file_name If protein domain searches are available, read them from file file_name
-d dir_name Training directory is specified by dir_name (introduced for compatibility with earlier versions)
-o file_name Print output in file_name; if n>1 for top best predictions, output is in file_name.1, file_name.2, ... , file_name.n f
-n n Print top n best predictions
-g Print output in gff format
-v Don't use svm splice site predictions
-f Don't make partial gene predictions
-h Display the options of the program
2.自行訓(xùn)練數(shù)據(jù)集
#1.手動訓(xùn)練谬晕，構(gòu)建數(shù)據(jù)集
$HOME/Tools/GlimmerHMM/GlimmerHMM/train/trainGlimmerHMM genome.part_Chr01.fa exon_file
#2.指定自己的訓(xùn)練結(jié)果文件夾
cut -f 1 -d ' ' genome.fa.tmp >genome.fa
glimmerhmm_linux_x86_64 genome.fa -d TrainGlimmM-MH63RS3 -g -n 1 -o glimm.gff
三碘裕、同源預(yù)測基因結(jié)構(gòu)
GenomeThreader
針對一個基因進(jìn)行同源預(yù)測
gth -genomic chr1_5k.fa -protein cer.fa -intermediate -gff3out
# 其中cer.fa就是AT1G02205.2的氨基酸序列
全基因組范圍預(yù)測流程如下：
準(zhǔn)備cDNA和或protein序列：在https://phytozome.jgi.doe.gov/p下載靠譜的物種的蛋白質(zhì)序列，如 Arabidopsis thaliana, Oryza sativa, Brassica rapa, 查找文獻(xiàn)尋找目前該物種的已有EST/cDNA序列攒钳，或者RNA-seq從頭組裝轉(zhuǎn)錄組帮孔。這里僅考慮用同源物種的蛋白序列進(jìn)行比對分析，轉(zhuǎn)錄組從頭組裝數(shù)據(jù)用于PASA整體比對到參考基因組和更新已有的基因解雇不撑。
分別測試不同物種的同源注釋結(jié)果文兢。
#run seperately
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Athaliana_167_TAIR10.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Athaliana.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/BrapaFPsc_277_v1.3.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Brapa.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Osativa_323_v7.0.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Osativa.gff3 &
當(dāng)然實際的同源注釋流程中不能是單個物種分別預(yù)測，應(yīng)該是將所有的蛋白序列進(jìn)行合并焕檬，然后用BLASTX找到最優(yōu)的聯(lián)配姆坚，之后用GenomeThreader進(jìn)行預(yù)測。PASA流程提到的UniRef90作為同源注釋的搜索數(shù)據(jù)庫可能是更好的選擇实愚，由于UniRef優(yōu)先選擇哪些人工審查兼呵、注釋質(zhì)量高兔辅、來源于模式動植物的蛋白，所以可靠性相對于直接使用同源物中可能更高萍程。
四幢妄、基于轉(zhuǎn)錄組的基因結(jié)構(gòu)預(yù)測
TransDecoder
安裝
wget https://github.com/TransDecoder/TransDecoder/archive/refs/tags/TransDecoder-v5.7.0.zip
unzip TransDecoder-v5.7.0.zip
$HOME/Tools/TransDecoder/TransDecoder-TransDecoder-v5.7.0/util
可以下載靠譜的轉(zhuǎn)錄組數(shù)據(jù)：https://phytozome.jgi.doe.gov/p
第一個步驟：用已有的注釋文件和參考基因組來進(jìn)行識別。
從有參考基因組的轉(zhuǎn)錄結(jié)果GTF文件預(yù)測編碼區(qū)域：
我們從cufflinks或stringtie輸出的gtf文件開始分析流程茫负，因此有兩個輸入文件：
? transcripts.gtf: 記錄預(yù)測轉(zhuǎn)錄本的GTF文件
? genome.fasta: 參考基因組序列
- 在含有轉(zhuǎn)錄本的gtf文件中提取fasta序列蕉鸳。
$HOME/Tools/TransDecoder/TransDecoder-TransDecoder-v5.7.0/util/gtf_genome_to_cdna_fasta.pl MH63.gtf genome.fa > transcripts.fasta
- 將gtf轉(zhuǎn)為gff
$HOME/Tools/TransDecoder/TransDecoder-TransDecoder-v5.7.0/util/gtf_to_alignment_gff3.pl MH63.gtf > MH63.gff3
- 提取最長的開放閱讀框，默認(rèn)是100個氨基酸忍法，可以用-m修改潮尝。默認(rèn)情況下，TransDecoder.LongOrfs將識別長度至少為100個氨基酸的開放閱讀框饿序∶闶В可以通過-m參數(shù)來降低這個值，但是要知道隨著最小長度的變短原探，ORF預(yù)測的假陽性率迅速增長乱凿。
$HOME/Tools/TransDecoder/TransDecoder-TransDecoder-v5.7.0/TransDecoder.LongOrfs -t transcripts.fasta
日志結(jié)果
Use file: /database/genome_assembly/03TransDecoder/transcripts.fasta.transdecoder_dir/longest_orfs.pep for Pfam and/or BlastP searches to enable homology-based coding region identification.
#可以用生成的longest_orfs.pep蛋白文件在Pfarm或者BlastP比對來識別編碼區(qū)域。識別他是否準(zhǔn)確咽弦。
Then, run TransDecoder.Predict for your final coding region predictions.
#最后的預(yù)測徒蟆。
- 最后生成一個基于有參基因組的編碼區(qū)域注釋文件:
cdna_alignment_orf_to_genome_orf.pl transcripts.fasta.transdecoder_dir/longest_orfs.gff3 transcripts.gff3 transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3
- 預(yù)測可能的編碼區(qū)
TransDecoder.Predict -t target_transcripts.fasta [ homology options ]
候選編碼區(qū)的最終集合可以在文件.transdecoder中找到。擴(kuò)展名包括.pep型型，.cds段审，.gff3和.bed。
$HOME/Tools/TransDecoder/TransDecoder-TransDecoder-v5.7.0/TransDecoder.Predict -t transcripts.fasta -retain_blastp_hits blastp.outfmt6
五闹蒜、三種方法的結(jié)果整合寺枉。
最后使用EvidenceModler根據(jù)上述結(jié)果進(jìn)行整合，使用PASA更新EVM的一致性預(yù)測绷落，增加UTR注釋和可變剪切注釋姥闪；在進(jìn)行適當(dāng)?shù)娜斯ぜm正。
https://blog.csdn.net/u012110870/article/details/82500684
https://yanzhongsino.github.io/2021/08/10/omics_genome.management/
未完待續(xù)
