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

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

tar xf v1.1.1.tar.gz

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

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

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

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

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

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