利用TransDecoder提取转录本CDS、UTR序列
前言：在对miRNA进行靶标预测分析时，需要特定物种的转录本作为靶向的数据库。通常对miRNA预测以3'UTR区域为主，这就需要对转录本数据库进行UTR区域的提取。
我在之前写过一个教程讲述了利用ExUTR包进行UTR序列的提取[1]，该脚本需要比对至UTR数据库，所以会丢失很大一部分没有比对上的序列。
在本文，利用TransDecoder包对转录本的ORF进行预测并比对至Uniprot数据库寻找同源支持，几乎所有的mRNA转录本都能获得ORF区域，并生成基于转录本的gff3注释文件。随后，利用python编程从TransDecoder的输出文件中，可以轻松提取出3'UTR、CDS和5'UTR序列，该方法是目前我认为最完整提取UTR的办法。
软件安装：
1.从网站下载
#创建文件夹并进入文件夹
mkdir -p ~/opt/biosoft && cd ~/opt/biosoft
#获取下载链接
wget https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.5.0.zip
#解压缩
unzip TransDecoder-v5.5.0.zip
#重命名
mv TransDecoder-TransDecoder-v5.5.0 TransDecoder-v5.5.0
2.利用conda安装
#激活虚拟环境
conda activate mirna
#在指定环境中安装包
conda install -n minran TransDecoder
#测试
TransDecoder help
运行TransDecoder
这里以我已经通过de novo组装获得的转录本或者自己收集到的基因fasta序列为起点开始，所以省略其他步骤。
输入的文件名为： transcripts.fasta
1. 预测转录本中长的开放阅读框, 默认是100个氨基酸，可以用-m修改。
TransDecoder.LongOrfs -t transcripts.fasta
输出的结果为两个文件夹，
transcripts.fasta.transdecoder_dir和transcripts.fasta.transdecoder_dir.__checkpoints。
2. 使用DIAMOND对上一步输出的transcripts.fasta.transdecoder.pep在蛋白数据库中进行搜索，寻找同源证据支持。
DIAMOND软件安装及使用详见其他教程[2]。
# BLASTP比对
diamond blastp -d uniprot_sprot.fasta -q transcripts.fasta.transdecoder_dir/longest_orfs.pep --evalue 1e-5 --max-target-seqs 1 > blastp.outfmt6
uniprot_sprot.fasta：为构建的蛋白数据库
--evalue： 设置比对的阈值1e-5
--max-target-seqs： 最多比对数为1
blastp.outfmt6 为输出的比对后序列
3. 预测可能的编码区
TransDecoder.Predict \
-t transcripts.fasta \
--retain_blastp_hits blastp.outfmt6
输出的文件包括：
transcripts.fasta.transdecoder.cds: 最终预测的CDS序列
transcripts.fasta.transdecoder.gff3: 最终ORF对应的GFF3
其中gff3文件内包含了CDS\five_primer_UTR\three_primer_UTR在转录本的信息，非常重要的文件。
--retain_blastp_hits <string> blastp output in '-outfmt 6' format. Any ORF with a blast match will be retained in the final output.
4.提取CDS、five_primer_UTR、three_primer_UTR序列
根据gff3文件的信息可以提取，这部分需要用到python编程：
- 对gff3文件进行预处理
#删除gff3文件中空白行,重定向为新的gff3文件
grep -v '^\s*$' 09_Bdo_endocytosis_trans_sequence.fasta.transdecoder.gff3 > 09_output.gff3
- 新建脚本并写入python代码
具体代码如下：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from Bio import SeqIO
# fasta = open("/data2/masw_data/seqdb/chr1A.fasta", "rU")
#转录本的fasta文件
record_dict = SeqIO.index("/home/user/input/09_Bdo_endocytosis_trans_sequence.fasta", "fasta")
five_UTR_sequence = open('five_UTR.fasta', 'w')
CDS_sequence = open('CDS.fasta','w')
three_UTR_sequence = open('three_UTR.fasta', 'w')
five_UTR = {}
CDS = {}
three_UTR = {}
with open('/home/user/input/09_output.gff3', 'r') as f:
for line in f:
line1 = line.strip().split()
chr = line1[0]
feature = line1[2]
start = line1[3]
end = line1[4]
direction = line1[6]
name = line1[0]
if feature == 'five_prime_UTR':
five_UTR[name] = (chr, start, end, direction)
if feature == 'CDS':
CDS[name] = (chr, start, end, direction)
if feature == 'three_prime_UTR':
three_UTR[name] = (chr, start, end, direction)
# get five_UTR_sequence
for key, value in five_UTR.items():
if value[3] == '+':
five_UTR_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1])-1:int(value[2])].seq))
if value[3] == '-':
five_UTR_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1]) - 1:int(value[2])].seq.reverse_complement()))
# get CDS
for key, value in CDS.items():
if value[3] == '+':
CDS_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1])-1:int(value[2])].seq))
if value[3] == '-':
CDS_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1]) - 1:int(value[2])].seq.reverse_complement()))
# get three_UTR_sequence
for key, value in three_UTR.items():
if value[3] == '+':
three_UTR_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1])-1:int(value[2])].seq))
if value[3] == '-':
three_UTR_sequence.write('>%s\n%s\n' % (key, record_dict[value[0]][int(value[1]) - 1:int(value[2])].seq.reverse_complement()))
five_UTR_sequence.close()
CDS_sequence.close()
three_UTR_sequence.close()
3）修改可执行权限
sudo chmod + x gff3_fasta.py
- 执行python脚本
./gff3_fasta.py
参考文献：
