GeneMarkGeorgia Institute of Technology開發(fā)的一系列基因預測工具莫瞬。真核生物基因組預測主要會用到GeneMark-ES/ET, 其中GeneMark-ES可用于無監(jiān)督自訓練咕缎，也就是只要提供一個基因組序列即可谴蔑，而GeneMark-ET則是在GeneMark-ES的基礎(chǔ)上整合了高通量的RNA-Seq轉(zhuǎn)錄本數(shù)據(jù)，工作流程如下
如果是學術(shù)代承、非盈利組織汁蝶，那么可以在http://exon.gatech.edu/GeneMark/license_download.cgi提交申請，之后就能得到軟件的下載鏈接。下載的工具只要解壓縮即可使用掖棉，不需要額外編譯墓律，但是需要安裝幾個perl模塊，
cpanm YAML Hash::Merge Logger::Simple Parallel::ForkManager
使用方法
在軟件目錄下的"README.GeneMark-ES-suite"介紹了如何使用軟件幔亥，對于GeneMark-ET而言耻讽，使用的腳本如下
gmes_petap.pl --sequence seq.fna --ET introns.gff --et_score 4 --cores 60
其中seq.fna是你的物種基因組序列，而introns.gff則是記錄內(nèi)含子的位置信息帕棉，--et_score是introns.gff中覆蓋該內(nèi)含子區(qū)域的read數(shù)针肥，默認是10，而--cores則是多線程參數(shù)香伴。
那么問題來了慰枕，如何獲取introns.gff？GeneMark-ET文章被NAR接收的時間是2013年9月即纲，目前我們一直推薦大家使用的HISAT2是2015年在Nature Methods發(fā)表具帮，而在HISAT2之前又快又好的比對工具是STAR是2013年1月。在這兩個工具出來之前崇裁，用的最多的RNA-seq比對工具就是目前連自己都宣布不維護居然還有很多培訓班在講的TopHat2匕坯，所以這篇文章的introns.gff文件就來自于TopHat2，還有兩個大家都未必聽過的工具(TrueSight, UnSplicer).
因此拔稳，在軟件目錄下下提供了bet_to_gff.pl用于將TopHat2的輸出結(jié)果中的junctions.bed轉(zhuǎn)成introns.gff
bed_to_gff.pl --bed tophat_out/junctions.bed --gff introns.gff --label tophat2
為了獲取junctions.bed，我特地去下載了TopHat2, 然后用了40線程進行比對锹雏，結(jié)果花了我將近3個小時的時間巴比，才將3G大小轉(zhuǎn)錄組數(shù)據(jù)回帖到基因組上。
這個速度實在是太慢了礁遵，我絕對忍受不了轻绞。突然間記憶中涌現(xiàn)出一位不愿透露姓名的朋友給別人講轉(zhuǎn)錄組時用到的STAR，它的輸出結(jié)果也有一個文件記錄著內(nèi)含子的信息佣耐。STAR的輸出文件"SJ.out.tab"存放著可信坍縮剪切位點政勃，每一列的格式說明如下，
- 第一列: 染色體
- 第二列: 內(nèi)含子起始（以1為基）
- 第三列: 內(nèi)含子結(jié)束（以1為基）
- 第四列：所在鏈兼砖，1(+)奸远，2(-)
- 第五列: 內(nèi)含子類型，0表示不是下面的任何一種功能讽挟，1表示GT/AG, 2表示:GT/AC,3表示GC/AG,4表示GT/GC,5表示AT/GC,6表示GT/AT
- 第六列: 是否是已知的注釋
- 第七列: 有多少唯一聯(lián)配支持
- 第八列: 有多少多重聯(lián)配支持
- 第九列: 最多可變剪切聯(lián)配
根據(jù)"README.GeneMark-ES-suite"里提到的introns.gff格式說明懒叛，
- 第一列：seqid
- 第二列: 來源
- 第三列是固定值，"intron"
- 第四列耽梅，第五列是內(nèi)含子起始第一個堿基位置和結(jié)束的最后一個堿基位置
- 第六列薛窥，覆蓋該外顯子的read數(shù)
- 第7列，+/-
- 第8列和第9列都為"."
其實很容易用一個awk命令實現(xiàn)格式轉(zhuǎn)換眼姐。如下是STAR所用到的命令
mkdir -p star_index
STAR \
--runThreadN 20 \
--runMode genomeGenerate \
--genomeDir star_index \
--genomeFastaFiles reference.fa
STAR
--runThreadN 20 \
--runMode alignReads \
--genomeDir star_index \
--readFilesIn read_1.fq.gz,read_2.fq.gz \
--readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --outWigType wiggle read2
awk 'BEGIN{OFS="\t"} $7 >= 2{if($4==1){st="+"}else{st="-"} print $1,"STAR","intron",$2,$3,$7,st,".","."}' SJ.out.tab > STAR.gff
后來我又發(fā)現(xiàn)诅迷，軟件目錄下還有一個star_to_gff.pl,
usage /opt/biosoft/gmes_petap/star_to_gff.pl --star [name] --gff [name] --label [label]
--star input name/s of junctions.bed from RNA-Seq alignment
--gff output intron coordinates in GFF format
--label [RNA_seq_junction] use this label in GFF to preserve name of the alignment tool
盡管它的使用說明--star居然寫著junctions.bed,但實際上--star接受的輸入就是SJ.out.tab（我看了一下源代碼）佩番，因此用法如下
star_to_gff.pl --star SJ.out.tab --gff SJ.gff --label STAR
將上述TopHat2，STAR(awk版本)罢杉，STAR(star_to_gff.pl版本)三者的GFF文件和BAME文件共同在IGV中展示時答捕，你會發(fā)現(xiàn)這三者幾乎是一模一樣
當然差異也是存在的，比如說下圖箭頭部分屑那。但是回溯到GFF文件時拱镐，我發(fā)現(xiàn)箭頭部分的intron僅有一條read支持，而tophat2和awk版本都是最低2條read支持持际。
最后使用genemark-et進行預測
gmes_petap.pl --sequence ../../../Assembly/5-final/ref.fa --ET STAR.gff --cores 60
最后會在當前文件下生成genemark.gtf沃琅。Genemark-et結(jié)果一定會有很高的假陽性，不過沒關(guān)系蜘欲，畢竟我們只需要它的模型輸出output/gmhmm.mod作為MAKER的輸入而已益眉。
