# 参考链接

如何对基因组进行注释

从头预测,同源注释和转录组整合都会得到一个预测结果,相当于收集了大量证据,下一步就是通过这些证据定义出更加可靠的基因结构,这一步可以通过人工排查,也可以使用 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 文件的第二列,第三列则是权重.

image-20201215134257484

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

image-20201215134654664

# 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

image-20201215140246554

可以看到这个顺序还是不对,需要再排序,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()
gffrename.py EVM_sort.gff LH > renamed.gff

image-20201215140853074

完美成功

# 结语

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

更新于 阅读次数

请我喝[茶]~( ̄▽ ̄)~*

amane 微信支付

微信支付

amane 支付宝

支付宝