基因结构注释（4）：整合预测结果

参考链接

如何对基因组进行注释

从头预测，同源注释和转录组整合都会得到一个预测结果，相当于收集了大量证据，下一步就是通过这些证据定义出更加可靠的基因结构，这一步可以通过人工排查，也可以使用EVidenceModeler(EVM).

EVM对gene_prediction.gff3有特殊的要求，就是GFF文件需要反映出一个基因的结构，gene->(mRNA -> (exon->cds(?))(+))(+), 表示一个基因可以有多个mRNA，即基因的可变剪接, 一个mRNA都可以由一个或者多个exon(外显子), 外显子可以是非翻译区(UTR),也可以是编码区(CDS). 而GlimmerHMM, SNAP等

这三类根据人为经验来确定其可信度，从直觉上就是用PASA根据mRNA得到的结果高于从头预测。

软件下载：

wget -4 https://github.com/EVidenceModeler/EVidenceModeler/archive/v1.1.1.tar.gz
tar xf v1.1.1.tar.gz

权重文件创建

首先将EVM文件夹下的simple_example/weights.txt复制到自己的目录

第一列是来源类型(ABINITIO_PREDICTION, PROTEIN, TRANSCRIPT), 第二列对应着GFF3文件的第二列，第三列则是权重.

根据需要修改成自己的,用制表符分割

EVM运行

准备好权重文件后，可以运行EVM了

cat augustus.gff genemark.gff > denovo.gff
evidence_modeler.pl --genome genome.fa --weights weights.txt \
                    --gene_predictions denovo.gff \
                    --protein_alignments proteinprediction.gff \
                    --transcript_alignments \
                    transcripts.fasta.transdecoder.genome.gff3 \
                    >evm.out

最后生成的文件为evm.out,转为gff3格式

并行EVM

主要是为了让整合结果更快一点

分割原始数据, 用于后续并行.

EvmUtils/partition_EVM_inputs.pl --genome genome.fa \
        --gene_predictions denovo.gff3 \
        --protein_alignments proteinprediction.gff \
        --transcript_alignments transcripts.fasta.transdecoder.genome.gff3 \
        --segmentSize 100000 --overlapSize 10000 \
        --partition_listing partitions_list.out

创建并行运算命令并执行

EvmUtils/write_EVM_commands.pl --genome genome.fa --weights `pwd`/weights.txt \
      --gene_predictions denovo.gff3 \
      --protein_alignments proteinprediction.gff \
      --transcript_alignments transcripts.fasta.transdecoder.genome.gff3 \
      --output_file_name evm.out  \
      --partitions partitions_list.out >  commands.list
parallel --jobs 10 < commands.list

合并运行结果

EvmUtils/recombine_EVM_partial_outputs.pl --partitions partitions_list.out \
         --output_file_name evm.out

结果转成gff3

EvmUtils/convert_EVM_outputs_to_GFF3.pl  --partitions partitions_list.out 
         --output evm.out  --genome genome.fa
find . -regex ".*evm.out.gff3" -exec cat {} \; | bedtools sort -i - > EVM.all.gff

基因过滤与命名

注释过滤：对于初步预测得到的基因，还可以稍微优化一下，例如剔除编码少于50个AA的预测结果，将转座子单独放到一个文件中

gffread EVM.all.gff -g input/genome.fa -y tr_cds.fa
bioawk -c fastx '$seq < 50 {print $comment}' tr_cds.fa | cut -d '=' -f 2 > short_aa_gene_list.txt
grep -v -w -f short_aa_gene_list.txt EvM.all.gff > filter.gff

可以看到这个顺序还是不对，需要再排序，gene->mRNA->exon->CDS

vi sort_EVM.py

import sys
import re
fout = open(sys.argv[3],'w')

ref_dict={}
with open(sys.argv[1]) as gene_a:
        for line in gene_a:
                line_s = line.strip().split('\t')
                info = re.split('=|;',line_s[8])
                ID = info[1]
                ref_set = []
                for n in range(0,8):
                        ref_set.append(line_s[n])
                ref_dict.setdefault(ID,[]).append(ref_set)
                ref_dict.setdefault(ID,[]).append(info[3])

with open(sys.argv[2]) as mrna:
        for eachline in mrna:
                i = eachline.strip().split('\t')
                info1 = re.split('=|;',i[8])
                parent = info1[3]
                mrna_n = info1[1]
                ref_set1 = []
                for a in range(0,8):
                        ref_set1.append(i[a])
                mrna_h = '\t'.join(ref_set1)
                if parent in ref_dict:
                        vs = ref_dict[parent]
                        head = '\t'.join(vs[0])
                        fout.write('%s\tID=%s;Name=%s\n'%(head,parent,vs[1]))
                        fout.write('%s\tID=%s;Parent=%s;Name=%s\n'%(mrna_h,mrna_n,parent,vs[1]))

fout.close()


phthon3 sort_EVM.py filter.gff filter.gff EVM_sort.gff

接下来进行重命名

对每个基因实现编号，形如ABCD000010的效果，方便后续分析。如下代码是基于EVM_sort.gff

vi rename.py

#!/usr/bin/env python3
import re
import sys

if len(sys.argv) < 3:
    sys.exit()

gff = open(sys.argv[1])
prf = sys.argv[2]

count = 0
mRNA  = 0
cds   = 0
exon  = 0

print("##gff-version 3.2.1")
for line in gff:
    if not line.startswith("\n"):
        records = line.split("\t")
        records[1] = "."
    if re.search(r"\tgene\t", line):
        count = count + 10
        mRNA  = 0
        gene_id = prf + str(count).zfill(6)
        records[8] = "ID={}".format(gene_id)
    elif re.search(r"\tmRNA\t", line):
        cds   = 0
        exon  = 0
        mRNA  = mRNA + 1
        mRNA_id    = gene_id + "." + str(mRNA)
        records[8] = "ID={};Parent={}".format(mRNA_id, gene_id)
    elif re.search(r"\texon\t", line):
        exon     = exon + 1
        exon_id  = mRNA_id + "_exon_" + str(exon)
        records[8] = "ID={};Parent={}".format(exon_id, mRNA_id)
    elif re.search(r"\tCDS\t", line):
        cds     = cds + 1
        cds_id  = mRNA_id + "_cds_" + str(cds)
        records[8] = "ID={};Parent={}".format(cds_id, mRNA_id)
    else:
        continue

    print("\t".join(records))

gff.close()


rename.py EVM_sort.gff LH > renamed.gff

完美成功

结语

接下来我将进行功能注释，目前想到的是用Nr库、Swiss-Prot（也可以用TrEMBL，但sp要更可靠）、interproscan（包含了Pfam以及GO）、KEGG这四个基本的注释，对于真菌基因组，可能会再加上COG、CAZyme、病毒相关因子(PHI数据库)、次生代谢基因(AntiSMASH)。