宏基因组测序分析（九）基因预测及去冗余

由于原核生物编码基因不包含内含子，由连续的编码区构成开放阅读框（ORF），编码区具有容易辨别的序列特征，基因预测准确性较高。

prokka 和 metaGeneMark 均为常用宏基因组注释软件，这里使用 prokka进行基因预测。

prokka是一个整合流程，调用prodigal进行ORF预测，并结合其他软件进行ncRNA等预测，输出多种格式结果

• 准备输入文件

B1.contigs.fa

基因预测

prokka \
--outdir B1_prokka \ # 输出目录
--prefix B1 \ # 输出文件前缀
--addgenes \ # 添加gene信息
--addmrna \ # 添加mrna信息
--locustag B1 \ # 基因ID添加tag
--kingdom Bacteria \ #物种分类Archaea|Bacteria|Mitochondria|Viruses
--metagenome \ # 宏基因组模式
--cpus 8 \ # 线程数
./B1.contigs.fa # 输入文件

• 主要结果：
  B1.faa : 预测的氨基酸序列文件
  B1.ffn ：预测出的转录本序列
  B1.fna ：输入的基因组文件
  B1.fsa ：添加Sequin tags的基因组文件
  B1.gbf ：Genbank格式注释结果
  B1.gff ：gff3格式注释结果
  B1.sqn ：ASN1 格式 Sequin文件，修改后可用于提交Genbank
  B1.tbl : Feature Table, 创建sequin文件使用
  B1.tsv : features: locus_tag,ftype,len_bp,gene,EC_number,COG,product
  B1.txt : 注释结果统计信息

• 对多个样品基因预测结果去冗余，得到 unigene:

## 将多样品转录本和氨基酸序列分别合并
cat ./*.ffn > all.trans.fa
cat ./*.faa > all.pep.fa

## 提取编码序列ID
sed -n "s/ˆ>\(\S\+\).*$/\1/p" all.pep.fa > all.cds.id

## 基于ID提取CDS序列
seqtk subseq all.trans.fa all.cds.id > all.cds.fa

## 对CDS进行聚类
cd-hit-est \
-i all.cds.fa \ # 输入数据
-o all.cds.cdhit \ # 输出数据
-c 0.95 \ # 相似性大于等于95%的为一类
-aS 0.9 \ # 控制短序列比对严格程度，表示比对区间要占到短序列的90%
-M 32000 \ #使用的内存32G
-T 8 # 使用的线程数

• 输出结果文件：
  all.cds.cdhit ： unigene 序列文件
  all.cds.cdhit.clstr：聚类cluster信息

• 基于 unigene ID提取其氨基酸序列

# 提取原始unigene ID
awk '$1 ~/ˆ>/{print $1}' all.cds.cdhit | sed 's/ˆ>//' > unigene.id

# 基于unigene ID 提取其氨基酸序列文件
seqtk subseq all.pep.fa unigene.id > unigene_pep.fasta

# 改名
mv all.cds.cdhit unigene_cds.fasta

• 输出结果：
  unigene_cds.fasta ：unigene CDS序列文件
  unigene_pep.fasta ：unigene 氨基酸序列文件

欢迎关注Bioinfor 生信云