TransDecoder寻找转录本编码区

  • 时间:
  • 来源:互联网
  • 文章标签:

TransDecoder寻找转录本编码区

  • 1. 软件安装
  • 2. 预测
  • 3. 可视化

背景
基础知识不是我的侧重点,感兴趣可以去进一步查阅(图来源于百度百科)。

在这里插入图片描述

物种:小鼠

1. 软件安装

1. 1 TransDecoder

conda的安装可以参考https://blog.csdn.net/qq_40794743/article/details/108112334

# 创建新的conda环境
conda create -y -n find_cds python=3
conda activate find_cds
conda install -y -c bioconda/label/cf201901 transdecoder

安装成功之后运行TransDecoder.Predict出现类似如下报错,提示没有安装URI/Escape模块。剩下的几行报错我也没仔细去看,应该是此包不存在的连锁效应。由于我是perl语言小白,所以通过biostars等平台得到了以下解决办法。
在这里插入图片描述
Launch Perl CPAN:如果使用服务器且没有root权限,选择local

perl -MCPAN -e shell

Install Perl URI/Escape Module:

install URI::Escape

在这里插入图片描述
中间省略,最后使用q退出。
在这里插入图片描述
安装完成之后,再次进行测试。

TransDecoder.Predict 

在这里插入图片描述
成功安装!

1.2 blast+

首先找到最新版本的blast,获取linux tar.gz编译版本的安装包。 ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/。之后使用wget命令进行下载,加-b参数将其挂到后台并生成日志文件。

# 下载
cd curreent_dir/
wget -bc ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.11.0/ncbi-blast-2.11.0+-x64-linux.tar.gz
mkdir blast
#解压缩
tar -zxvf ncbi-blast-2.11.0+-x64-linux.tar.gz -C blast
#配置路径
echo "export PATH=/db/home/shenwei/local/app/blast/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc
# 查看是否安装成功
blastp -version

在这里插入图片描述

1.3 hmmer

hmmer的安装很顺利,该软件与transdecoder安装在同一conda小环境下。

conda install -c bioconda/label/cf201901 hmmer
hmmbuild --help

在这里插入图片描述

2. 预测

2.1 文件准备

官网分别提供了从fasta文件开始和从stringtie或cufflinks等得到的gtf文件开始的process。我这里是通过stringtie得到的gtf文件,所以需要准备如下两个文件:

  • transcripts.gtf: 记录预测转录本的GTF文件
  • genome.fasta: 参考基因组序列

2.1 开始

第一步:依据reference fasta文件,从gtf文件中提取fasta序列。

gtf_genome_to_cdna_fasta.pl trans.gtf /xx/xx/Mus_musculus.GRCm38.dna.toplevel.fa > trans.fa

过程截图
在这里插入图片描述
第二步:将gtf文件转换为gff3格式(只是为了后续分析)。
文件格式可以参考https://blog.csdn.net/qq_40794743/article/details/108112334

gtf_to_alignment_gff3.pl trans.gtf > trans.gff3

第三步: 预测转录本中长的开放阅读框,-m参数设置为300个氨基酸。blastp和pfam比对。

TransDecoder.LongOrfs -m 300 -t trans.fa 

过程截图
在这里插入图片描述
第四步:为了提高ORF预测的灵敏度,与已知的蛋白数据库进行比对,探索同源性。

Blastp Search

*diamond也是一款蛋白序列比对软件,据说速度比blastp快500多倍。感兴趣可以去探索。 *

我这里使用Swissprot数据库的蛋白序列,网址是ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/。NCBI官网也提供了database,大家可以自行去下载。

# 下载
 wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
 
# 解压缩
 gunzip uniprot_sprot.fasta.gz
 
# 建数据库
makeblastdb -in uniprot_sprot.fasta -dbtype prot -parse_seqids -hash_index -out uniprot
## -in 输入数据库文件
## -dbtype 蛋白库用prot,核酸用nucl
## -out 所建数据库的名称
## -parse_seqids      
## -hash_index 创建哈希序列

# 数据比对
nohup blastp -query novel.coding.fa.transdecoder_dir/longest_orfs.pep \
-db uniprot -max_target_seqs 1 \
-outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6 &
## 参数说明:
### -query: 输入文件路径及文件名(.fasta格式);
### -db: 格式化了的数据库路径及数据库名(数据库可以从PDB/NCBI里下载所有的相关/整个库中的序列);
### -outfmt:输出的文件格式,6是tabular格式对应BLAST的m8格式;
### -evalue: 设置输出结果的e-value值;
### -max_target_seqs:找到最大的目标的数目,也可以用-num_descriptions,tabular格式输出结果的条数;
### -num_threads :线程数

blast建库结果文件
在这里插入图片描述

比对结果

# Query id, Subject id, % identity, alignment length, mismatches, gap openings, q.start, q.end, s.start, s.end, e-value, bit score
第一列为Query(递交序列),
第二列为数据库序列(目标序列subejct),
第三列为: identity
第四列为:比对长度
第五列为:错配数
第六列为:gap数
第七列和第八列为:Query开始碱基位置和结束碱基位置
第九列和第十列为:Subject开始碱基位置和结束碱基位置
第十一列为:期望值
第十二列为:比对得分

在这里插入图片描述

Pfam Search

Pfam 数据库中每个编号代表一个蛋白质家族。Pfam 分 A 和 B 两个数据库,其中 A 数据库是经过手工校正的高质量数据库, B 数据库虽然质量低些,依然可以用来寻找蛋白质家族的保守位点。Pfam v27.0 版本的数据库中, A 数据库包含 14,836 个蛋白质家族编号(以 PF 开头); B 数据库包含 20,000 个蛋白质家族编号 (以 PB 开头)。但是在V28之后,B数据库由于数据的不准,已经被抛弃了。

mkdir hmm_pfam;cd hmm_pfam
# 下载 PFAM 数据库
nohup wget -c ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam33.1/Pfam-A.hmm.gz &
# 解压缩
gunzip Pfam-A.hmm.gz 
# 得到 PFAM 数据库的 HMM 文件。 HMM 文件是文本文件,需要将其变成二进制格式,以加快运算速度,同时进行压缩

# 建立索引数据库
hmmpress Pfam-A.hmm

建库过程截图
在这里插入图片描述
建库结果文件
在这里插入图片描述

# 使用 hmmscan 进行 Pfam 注释
hmm=/mnt/data/hwb/hmm_pfam/Pfam-A.hmm
cd /xx/xx # 进入fa文件目录
nohup 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
##     线程

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

nohup TransDecoder.Predict -t trans.fa \
--retain_pfam_hits hmm_pfam.trans.tbl \
--retain_blastp_hits blastp.outfmt6 &

最终的编码区域预测结果将是TransDecoder.Predict,blastp,hmmer三者结果的交集。

最终结果文件:

  • base_freqs.dat: 碱基出现频率
  • trans.fa.transdecoder.pep: 最终预测的CDS对应的蛋白序列
  • trans.fa…transdecoder.cds: 最终预测的CDS序列
  • trans.fa.transdecoder.gff3: 最终ORF对应的GFF3
  • trans.fa.transdecoder.bed: 以BED格式存放ORF位置信息

当然,还有blastp和hmmer产生的结果文件。

3. 可视化

trans.fa.transdecoder.bed可以拖到GenomeView或igv进行可视化。

pfam和blastp的比对结果可以通过散点和density组合图的形式展示(感觉还不错)。
在这里插入图片描述

图来源于NC 文章

由于我这里的结果不便展示,所以我构造了一个数据集,来展示这个图的画法。

# 载入包,构造数据集
library(ggplot2)
library(gridExtra)
library(ggsci)
N<-200
x1 <- rnorm(mean=1.5, sd=0.5,N)
y1 <- rnorm(mean=2,sd=0.2, N)
x2 <- rnorm(mean=2.5,sd=0.5, N)
y2 <- rnorm(mean=2.5,sd=0.5, N)
x3 <- rnorm(mean=1, sd=0.3,N)
y3 <- rnorm(mean=1.5,sd=0.2, N)

data2 <- data.frame(x=c(x1,x2,x3),y=c(y1,y2,y3),class=rep(c("A","B","C"),each=200))

# 主图
scatter <- ggplot() + 
  geom_point(data=data2,aes(x=x,y=y,fill=class),shape=21,color="black",size=3)+
  scale_fill_aaas()+
  geom_vline(xintercept = 3,linetype = 2) +
  geom_hline(yintercept = 3,linetype = 2) +
  theme_test()+
  theme(legend.position=c(0.1,0.85))

# 顶部图
hist_top <- ggplot()+
  geom_density(data=data2,aes(x),colour='black',alpha=0.7,size = 1.2)+
  theme_classic()+
  theme(legend.position="none",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# 右部图
hist_right <- ggplot()+
  geom_density(data=data2,aes(y),colour='black',alpha=0.7,size = 1.2)+
  theme_classic()+
  coord_flip()+
  theme(legend.position="none",
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())
        
# 创建空白图层,可选
# empty <- ggplot() +
#   theme(panel.background=element_blank(),
#         axis.title.x=element_blank(), 
#         axis.title.y=element_blank(),
#         axis.text.x=element_blank(),
#         axis.text.y=element_blank(),
#         axis.ticks=element_blank())

# 合并
ggarrange(hist_top, empty, scatter, hist_right, ncol=2, nrow=2, 
          widths=c(4,1), heights=c(1,4),align = "hv")

在这里插入图片描述
感觉还不错,细节就交给AI吧。


##侵权联系删除!

本文链接http://www.taodudu.cc/news/show-1944915.html