转录本预测编码区：在transcript中找CDS区

RNA-Seq数据组装得到的转录本具有编码蛋白质的能力。预测编码区对于确定转录产物在细胞中的分子作用至关重要!!!

通过Cufflinks或Stringtie得到的gtf文件，所以需要准备如下两个文件：

  1. transcripts.gtf: 你自己的RNA-seq数据经过Cufflinks，或Stringtie得到的gtf文件，记录预测转录本的GTF文件

  2. genome.fasta: 参考基因组序列

第一步：依据参考genome.fasta文件，从gtf文件中提取fasta序列。

conda activate python3 #进入环境
TransDecoder.Predict #启动TransDecoder

###基于mm10的 Tophat结果
示例：
gtf_genome_to_cdna_fasta.pl trans.gtf /xx/xx/Mus_musculus.GRCm38.dna.toplevel.fa > trans.fa
操作路径于：/home/u20111230014/workspace/tran_gtf_ZF 
PS: DPP-0_transcripts.gtf和 mm10.fa必须在同一个文件夹
实例：
gtf_genome_to_cdna_fasta.pl DPP-0_transcripts.gtf mm10.fa >  DPP-0_mm10_trans.fa
gtf_genome_to_cdna_fasta.pl DPP-1_transcripts.gtf mm10.fa >  DPP-1_mm10_trans.fa

###基于GRCm39的 Tophat结果
操作所在路径：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna

实例
将所有所需要的文件软连接到同一文件夹
ln -s /home/u20111230014/workspace/genome/GRCm39/Tophat_Cufflinks_DPP-0_clean.fq/transcripts.gtf DPP-0_GRCm39_transcripts.gtf 

gtf_genome_to_cdna_fasta.pl DPP-0_GRCm39_transcripts.gtf GRCm39.fa >  DPP-0_GRCm39_trans.fa
gtf_genome_to_cdna_fasta.pl DPP-1_GRCm39_transcripts.gtf GRCm39.fa >  DPP-1_GRCm39_trans.fa

第二步：将gtf文件转换为gff3格式(只是为了后续分析)。

示例
gtf_to_alignment_gff3.pl trans.gtf > trans.gff3

实例
###基于mm10的 Tophat结果
操作路径于：/home/u20111230014/workspace/tran_gtf_ZF 
gtf_to_alignment_gff3.pl DPP-0_transcripts.gtf >  DPP-0_mm10_trans.gff3
gtf_to_alignment_gff3.pl DPP-1_transcripts.gtf >  DPP-1_mm10_trans.gff3

###基于GRCm39的 Tophat结果
操作所在路径：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna
软连接脚本：

ln -s /home/u20111230014/workspace/tran_gtf_ZF/gtf_to_alignment_gff3.slurm tophat_gtf_to_alignment_gff3.slurm
实例
gtf_to_alignment_gff3.pl DPP-0_GRCm39_transcripts.gtf >  DPP-0_GRCm39_trans.gff3
gtf_to_alignment_gff3.pl DPP-1_GRCm39_transcripts.gtf >  DPP-1_GRCm39_trans.gff3

第三步：预测转录本中长的开放阅读框，-m参数设置为300个氨基酸。并进行blastp和hmmer比对。

示例
TransDecoder.LongOrfs -m 300 -t trans.fa
实例
TransDecoder.LongOrfs -m 300 -t DPP-0_mm10_trans.fa
TransDecoder.LongOrfs -m 300 -t DPP-1_mm10_trans.fa

第四步：为了提高预测的灵敏度，与已知的蛋白数据库进行比对，探索同源性。

1) blastp search

选用Swissprot数据库的蛋白序列，下载uniprot_sprot.fasta.gz，网址是ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/

解压缩

 gunzip uniprot_sprot.fasta.gz

建blast+库于路径：/home/u20111230014/workspace/genome  
#### 建数据库
makeblastdb -in uniprot_sprot.fasta -dbtype prot -parse_seqids -hash_index -out uniprot

## makeblastdb参数说明：

-in 输入数据库文件

-dbtype 蛋白库用prot,核酸用nucl

-out 所建数据库的名称

 -parse_seqids      

-hash_index 创建哈希序列

数据比对

###基于mm10.fa的比对
操作路径在：/home/u20111230014/workspace/tran_gtf_ZF 
index文件：uniprot的11个子文件，再加上一个uniprot_sprot.fasta，共12个
示例：
blastp -query trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6 &
实例：
blastp -query DPP-0_mm10_trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > DPP-0_mm10_blastp.outfmt6

！！！PS: 如果报错：BLAST Database error: No alias or index file found for nucleotide database [nt] in search path, 必须将uniprot的index，共12个文件，跟要跑的所有文件放在一起才行，如下：

(python3) [u20111230014@cpu10 tran_gtf_ZF]$ ll uniprot*
lrwxrwxrwx 1 u20111230014 u20111230014 43 Feb 23 16:23 uniprot -> /home/u20111230014/workspace/genome/uniprot
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pdb -> /home/u20111230014/workspace/genome/uniprot/uniprot.pdb
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phd -> /home/u20111230014/workspace/genome/uniprot/uniprot.phd
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phi -> /home/u20111230014/workspace/genome/uniprot/uniprot.phi
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.phr -> /home/u20111230014/workspace/genome/uniprot/uniprot.phr
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pin -> /home/u20111230014/workspace/genome/uniprot/uniprot.pin
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pog -> /home/u20111230014/workspace/genome/uniprot/uniprot.pog
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pos -> /home/u20111230014/workspace/genome/uniprot/uniprot.pos
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pot -> /home/u20111230014/workspace/genome/uniprot/uniprot.pot
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.psq -> /home/u20111230014/workspace/genome/uniprot/uniprot.psq
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.ptf -> /home/u20111230014/workspace/genome/uniprot/uniprot.ptf
lrwxrwxrwx 1 u20111230014 u20111230014 55 Feb 24 18:32 uniprot.pto -> /home/u20111230014/workspace/genome/uniprot/uniprot.pto
lrwxrwxrwx 1 u20111230014 u20111230014 63 Feb 23 16:28 uniprot_sprot.fasta -> /home/u20111230014/workspace/genome/uniprot/uniprot_sprot.fasta

如果还是报同样的错，就换一种方法：下载preformatted的数据库！

去到miniconda下 blast包里含update_blastdb.pl文件的路径，用perl 更新一下数据库
!!!blast软件所在的位置：/home/u20111230014/miniconda3/pkgs/blast-2.12.0-h3289130_3/bin

先查看数据库列表：
(python3) [u20111230014@cpu10 bin]$ perl update_blastdb.pl --showall
Connected to NCBI
16S_ribosomal_RNA
18S_fungal_sequences
28S_fungal_sequences
Betacoronavirus
ITS_RefSeq_Fungi
ITS_eukaryote_sequences
LSU_eukaryote_rRNA
LSU_prokaryote_rRNA
SSU_eukaryote_rRNA
env_nt
env_nr
human_genome
landmark
mito
mouse_genome
nr
nt
pataa
patnt
pdbaa
pdbnt
ref_euk_rep_genomes
ref_prok_rep_genomes
ref_viroids_rep_genomes
ref_viruses_rep_genomes
refseq_select_rna
refseq_select_prot
refseq_protein
refseq_rna
swissprot
tsa_nr
tsa_nt
taxdb
###再选择需要的数据库进行更新下载：
perl update_blastdb.pl nr #nr 是蛋白质， nt是核酸
!!!PS:nohup挂在后台下载，中途如果超时的话，需重新执行命令就行，之前下载好的不用重新执行！！！大概有十几个文件，下好批量解压命令：ls *.tar.gz | xargs -n1 tar xzvf

基于GRCm39.fa的比对

操作路径在：/home/u20111230014/workspace/genome/GRCm39/tophat_gtf_genome_to_cdna/
index文件：uniprot的11个子文件，再加上一个uniprot_sprot.fasta，共12个
实例
blastp -query DPP-0_mm10_trans.fa.transdecoder_dir/longest_orfs.pep -db uniprot -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > DPP-0_mm10_blastp.outfmt6

2) Pfam search

在建立索引数据库hmmpress Pfam-A.hmm成功后，使用 hmmscan 进行 Pfam 注释
conda activate python3 #进入环境
hmmbuild --help #启动
Pfam-A存放于：hmm=/home/u20111230014/workspace/genome/Pfam-A/Pfam-A.hmm
PS: 执行hmmscan 前输入上面的路径指示，或直接写入slurm脚本

进入存放转换格式好的样本tran.fa文件目录

cd /home/u20111230014/workspace/tran_gtf_ZF 
(python3) [u20111230014@cpu6 tran_gtf_ZF]$ ll *mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38711086 Feb 23 15:26 DPP-0_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38539068 Feb 23 15:28 DPP-1_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 45990886 Feb 23 15:31 DPP-2_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 39814002 Feb 23 15:34 DPP-3_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 33385081 Feb 23 15:37 DPP-4_mm10_trans.fa
-rw-rw-r-- 1 u20111230014 u20111230014 38240368 Feb 23 15:40 DPP-5_mm10_trans.fa

示例：
hmmscan --cpu 8 -o hmm_pfam.trans.txt --tblout hmm_pfam.trans.tbl --noali -E 1e-5 $hmm trans.fa

###参数说明
 -o FILE 将结果输出到指定的文件中。默认是输出到标准输出

 --tblout FILE  将蛋白质家族的结果以表格形式输出到指定的文件中。默认不输出该文件

--domtblout FILE  将蛋白结构域的比对结果以表格形式输出到指定的文件中。默认不输出该文件。该表格中包含query序列起始结束位点与目标序列起始结束位点的匹配信息

 --acc 在输出结果中包含 PF 的编号，默认是蛋白质家族的名称

--noali 在输出结果中不包含比对信息。输出文件的大小则会更小

-E FLOAT    default:10.0 设定 E_value 阈值，推荐设置为 1e-5

 -T FLOAT  设定 Score 阈值

--domE FLOAT    default:10.0  设定 E_value 阈值。该参数和 -E 参数类似，不过是 domain 比对设定的值

--cpu 线程
实例脚本：
hmm=/home/u20111230014/workspace/genome/Pfam-A/Pfam-A.hmm
hmmscan --cpu 25 -o DPP-0_mm10_hmm_pfam.trans.txt --tblout DPP-0_mm10_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-0_mm10_trans.fa
hmmscan --cpu 25 -o DPP-1_mm10_hmm_pfam.trans.txt --tblout DPP-1_mm10_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-1_mm10_trans.fa
####
hmmscan --cpu 25 -o DPP-0_GRCm39_hmm_pfam.trans.txt --tblout DPP-0_GRCm39_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-0_GRCm39_trans.fa

hmmscan --cpu 25 -o DPP-1_GRCm39_hmm_pfam.trans.txt --tblout DPP-1_GRCm39_hmm_pfam.trans.tbl --noali -E 1e-5 $hmm DPP-1_GRCm39_trans.fa

第五步：结合blastp和pfam的结果预测可能的编码区。

TransDecoder.Predict #启动
示例：
TransDecoder.Predict -t trans.fa --retain_pfam_hits hmm_pfam.trans.tbl --retain_blastp_hits blastp.outfmt6
实例：
TransDecoder.Predict -t DPP-0_GRCm39_trans.fa --retain_pfam_hits DPP-0_GRCm39_hmm_pfam.trans.tbl --retain_blastp_hits DPP-0_GRCm39_blastp.outfmt6

第六步: 生成基于参考基因组的编码区注释文件

示例：
cdna_alignment_orf_to_genome_orf.pl trans.fa.transdecoder.gff3 trans.gff3 trans.fa > trans.fa.transdecoder.genome.gff3

实例：
cdna_alignment_orf_to_genome_orf.pl DPP-0_mm10_trans.fa.transdecoder.gff3 DPP-0_mm10_trans.gff3 DPP-0_mm10_trans.fa > DPP-0_mm10_trans.fa.transdecoder.genome.gff3

最终输出结果文件如下:

  1. trans.fa.transdecoder.pep: 最终预测的CDS对应的蛋白序列

  2. trans.fa.transdecoder.cds: 最终预测的CDS序列

  3. trans.fa.transdecoder.gff3: 最终ORF对应的GFF3

  4. trans.fa.transdecoder.bed: 以BED格式存放ORF位置信息

  5. trans.fa.transdecoder.genome.gff3: 基于参考基因组的GFF3文件

BED和GFF3文件可以使用IGV可视化。

参考：（[小张聊科研]）https://mp.weixin.qq.com/s?src=11&timestamp=1645587577&ver=3637&signature=ydPEWFG06TSywWPWpAOTRgJNtSQ2BUEOlUYjQZV3OGiZKZxJmIp2mHe9bWlV3UDcnKn0p1K0SWstUxzPQyTmCKC7GdYDBRdHKhKYct*vh-mH-qsE4a4A-TBNazUb&new=1

生信分析师高端
http://dl2.kerwin.cn:8063/csdn/key/article-a_giant_pig-103065967/auth/1645608443-2022062317272318-0-2af5159d4f371312acc26ac982c48187

卡西莫多的礼物
https://blog.csdn.net/qq_35696312/article/details/104779347

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,362评论 5赞 477
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,330评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,247评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,560评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,580评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,569评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,929评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,587评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,840评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,596评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,678评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,366评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,945评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,929评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,165评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 43,271评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,403评论 2赞 342

推荐阅读更多精彩内容