利用Diamond-OrfPredictor获取可靠的转录本编码

总体思路是先通过Diamond将所有的转录本用blastp进行比对到UniprotKB蛋白序列非冗余库中。

将比对的结果获取了用于指导后期的ORF寻找中。

所用软件为：

Diamond

OrfPredictor

本地化的OrfPredictor下载地址

文件放置方式如下：

yeyt@ubuntu:~/biodata/diamondblast$ tree .
.
├── orfpredictor
│   └── OrfPredictor
│       ├── extractCDS.pl
│       ├── How2Run_OrfPredictor.pdf
│       ├── OrfPredictor_agreement.pdf
│       └── OrfPredictor.pl
├── stringtie75.fasta
└── uniprotKB.dmdb.dmnd

设置一个ORF.sh脚本

yeyt@ubuntu:~/biodata/diamondblast$ cat ORF.sh 
diamond blastx --db uniprotKB.dmdb.dmnd --query stringtie75.fasta -k 1 --outfmt 0 -o stringtie75.dmd.blastx0 --threads 18 #每次仅仅输出一条目标hit，并且输出结果格式为0，采用线程为18个
perl orfpredictor/OrfPredictor/OrfPredictor.pl stringtie75.fasta stringtie75.dmd.blastx0 1 both 0.001 stringtie75.orf.pep  # 然后用OrfPredictor进行预测，其中输入了Diamond结果和核酸序列fasta序列。

然后运行这个shell script 即可

最后得到的就是

yeyt@ubuntu:~/biodata/diamondblast$ ll -alt 
total 69839552
drwxrwxr-x 3 yeyt yeyt        4096 May 15 11:24 ./
-rw-rw-r-- 1 yeyt yeyt   545349352 May 15 01:25 ORF6frame.txt  #包含所有的转录本6个读框的翻译结果
-rw-rw-r-- 1 yeyt yeyt         890 May 15 01:25 noOrf.txt # 包含没有找到读框的序列名字
-rw-rw-r-- 1 yeyt yeyt    67096763 May 15 01:25 stringtie75.orf.pep #所有核酸序列所翻译得到的蛋白序列 
-rw-rw-r-- 1 yeyt yeyt   250819435 May 15 00:23 stringtie75.dmd.blastx0 #Diamond比对结果pairwise输出
-rw-rw-r-- 1 yeyt yeyt         244 May 14 22:51 ORF.sh #运行的脚本

然后可以看到序列的具体情况

yeyt@ubuntu:~/biodata/diamondblast$ head stringtie75.orf.pep 
>maker-Fvb1-1-snap-gene-0.15-mRNA-1 +1  1   1497
MVEPFNSIKKKMSVQVALPGSDMSRAFCKGASEIVLGMCDKVVNTDGEALPLSEEQRNKISDVINGFSCEALRTLCVAFK
DIESPSGAESIPEDGYTLIAVVGIKDPVRPAVREAVKTCLDAGITVRMVTGDNINTAKAIAKECGILTEDGLAIEGPDFR
KMSEQEMAEIIPKLQVMARSLPLDKHTLVKQLRNVHKDVVAVTGDGTNDAPALHEADIGLAMGIAGTEVAKENADVIIMD
DNFTSIVNVARWGRAVYINIQKFVQFQLTVNIVALMLNFVSACVSGSAPLTAVQMLWVNLIMDTLEPPHDGLMRRPPVGR
NTNFITGIMWRNIIGQSIYQITVLLVLNFCGMPLLNITVQKQVLLSRPSYSTVLYFARRGRRNRNYFFCFKIFSLNVFNE
INSRDMEKINVFGGLFQSYTFMVVMICTVVFQILIIEFLGTFAQTVPLSWELWLASVLIGAVSLLVAVVLKCIPVSIRKQ
DAKDHGDIHEPLLRGPELA
>maker-Fvb1-1-augustus-gene-0.13-mRNA-1 +3  432 701
MLEVKVHHAKESQTQLLPLFPTVNFLEVICPIPGMFQILLVGCWFHLSLVEGAMLSQKMSASYLSRSQLAVHHSKGYHLC
yeyt@ubuntu:~/biodata/diamondblast$ head -100 stringtie75.orf.pep 
>maker-Fvb1-1-snap-gene-0.15-mRNA-1 +1  1   1497
MVEPFNSIKKKMSVQVALPGSDMSRAFCKGASEIVLGMCDKVVNTDGEALPLSEEQRNKISDVINGFSCEALRTLCVAFK
DIESPSGAESIPEDGYTLIAVVGIKDPVRPAVREAVKTCLDAGITVRMVTGDNINTAKAIAKECGILTEDGLAIEGPDFR
KMSEQEMAEIIPKLQVMARSLPLDKHTLVKQLRNVHKDVVAVTGDGTNDAPALHEADIGLAMGIAGTEVAKENADVIIMD
DNFTSIVNVARWGRAVYINIQKFVQFQLTVNIVALMLNFVSACVSGSAPLTAVQMLWVNLIMDTLEPPHDGLMRRPPVGR
NTNFITGIMWRNIIGQSIYQITVLLVLNFCGMPLLNITVQKQVLLSRPSYSTVLYFARRGRRNRNYFFCFKIFSLNVFNE
INSRDMEKINVFGGLFQSYTFMVVMICTVVFQILIIEFLGTFAQTVPLSWELWLASVLIGAVSLLVAVVLKCIPVSIRKQ
DAKDHGDIHEPLLRGPELA
>maker-Fvb1-1-augustus-gene-0.13-mRNA-1 +3  432 701
MLEVKVHHAKESQTQLLPLFPTVNFLEVICPIPGMFQILLVGCWFHLSLVEGAMLSQKMSASYLSRSQLAVHHSKGYHLC
VRLRGSCSLR
>MSTRG.1.1  +2  248 517
MLEVKVHHAKESQTQLLPLFPTVNFLEVICPIPGMFQILLVGCWFHLSLVEGAMLSQKMSASYLSRSQLAVHHSKGYHLC
VRLRGSCSLR
>MSTRG.2.1  +2  251 382
MHLLDKRSTRALPVALMIHDNSTDRTAFVLATHHSNICPINFRW
>maker-Fvb1-1-snap-gene-0.18-mRNA-1 +1  1   210
MIKKLVVGPWVGSTGPPYGVHRSARPFYRRYAPGLNWPGRASGAVTLKKLECSKQAYALDTLAWDNIIGF

并且可以提出头文件信息

yeyt@ubuntu:~/biodata/diamondblast$ grep '>' stringtie75.orf.pep | head 
>maker-Fvb1-1-snap-gene-0.15-mRNA-1 +1  1   1497
>maker-Fvb1-1-augustus-gene-0.13-mRNA-1 +3  432 701
>MSTRG.1.1  +2  248 517
>MSTRG.2.1  +2  251 382
>maker-Fvb1-1-snap-gene-0.18-mRNA-1 +1  1   210
>MSTRG.3.2  -1  877 1398
>MSTRG.3.4  -1  919 1440
>MSTRG.3.3  -3  1332    1853
>MSTRG.3.1  -3  1998    2519
>snap_masked-Fvb1-1-processed-gene-0.6-mRNA-1   -3  3   284
yeyt@ubuntu:~/biodata/diamondblast$ grep '>' stringtie75.orf.pep > stringtie75.orf.pep.headinfo

这个头文件可以用于统计读框长度分布等数据和可视化