总体思路是先通过Diamond将所有的转录本用blastp进行比对到UniprotKB蛋白序列非冗余库中。
将比对的结果获取了用于指导后期的ORF寻找中。
所用软件为：
Diamond
文件放置方式如下：
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
这个头文件可以用于统计读框长度分布等数据和可视化
