gff3格式

基因結(jié)構(gòu)注釋文件一般為gff3的格式，一共是9列，依次為基因組序列id，注釋來源，類型，起始位置，終止位置，得分，正負(fù)鏈，相位，屬性。
基因結(jié)構(gòu)注釋文件中，基因包含mRNA，mRNA包含exon, CDS, UTR等信息，同時(shí)在注釋文件中除基因行外，其他行在第9列會通過Parent指明該行從屬的上一級ID，也就是一個(gè)基因的gene行、mRNA行、CDS行、exon行都會通過Parent層層關(guān)聯(lián)在一起。

EVM軟件整合出的基因結(jié)構(gòu)注釋文件一般如下所示，目前有一個(gè)需求，就是修改基因名稱，修改為類似AT01G000001這種形式

Contig1 EVM gene    6013    6784    .   -   .   ID=evm.TU.Contig1.1;Name=EVM%20prediction%20Contig1.1
Contig1 EVM mRNA    6013    6784    .   -   .   ID=evm.model.Contig1.1;Parent=evm.TU.Contig1.1;Name=EVM%20prediction%20Contig1.1
Contig1 EVM exon    6013    6447    .   -   .   ID=cds.evm.model.Contig1.1;Parent=evm.model.Contig1.1
Contig1 EVM exon    6530    6784    .   -   .   ID=cds.evm.model.Contig1.1;Parent=evm.model.Contig1.1
Contig1 EVM CDS 6013    6447    .   -   0   ID=cds.evm.model.Contig1.1;Parent=evm.model.Contig1.1
Contig1 EVM CDS 6530    6784    .   -   0   ID=cds.evm.model.Contig1.1;Parent=evm.model.Contig1.1

gffutils包

gffutils是一個(gè)用于解析gff3/gtf格式的python包，讀取gff3文件非常方便，利用這個(gè)包實(shí)現(xiàn)對基因結(jié)構(gòu)注釋文件重新命名python腳本如下

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import re
import sys
import argparse
import logging
import gffutils
from collections import defaultdict

def rename(args):
    seqid2name = dict()
    for line in open(args.change, 'r'):
        tem = line.strip().split()
        seqid2name[tem[0]] = tem[1]

    db = gffutils.create_db(args.gff,':memory:',force=True,keep_order=True,merge_strategy="create_unique", sort_attribute_values = True)
    mRNA_children=("five_prime_UTR","three_prime_UTR","CDS","exon")
    idmap={
        "CDS" : "cds",
        "exon" : "exon",
        "five_prime_UTR":"5utrp",
        "three_prime_UTR":"3utrp"
    }

    f_out = open('%s.rename.gff3' % args.prefix, 'w')
    seqid = None
    for gene in db.features_of_type("gene",order_by=('seqid','start','end')):
        if gene.seqid != seqid:
            genenum = 0
        genenum += args.addnum
        seqid = gene.seqid

        genename = '{0}G{1:06}'.format(seqid2name[seqid],genenum)
        f_out.write('{seqid}\t{source}\t{featuretype}\t{start}\t{end}\t{score}\t{strand}\t{frame}\tID={geneid}\n'.format(seqid=gene.seqid, source=gene.source, featuretype=gene.featuretype, start=gene.start, end=gene.end, score=gene.score, strand=gene.strand, frame=gene.frame, geneid=genename))

        for t,mRNA in enumerate(db.children(gene,featuretype="mRNA",order_by=('seqid','start','end'))):
            mrna_num = t + 1
            mrnaid = '{genename}.mRNA{num}'.format(genename=genename,num=mrna_num)
            f_out.write('{seqid}\t{source}\t{featuretype}\t{start}\t{end}\t{score}\t{strand}\t{frame}\tID={mrnaid};Parent={geneid}\n'.format(seqid=mRNA.seqid, source=mRNA.source, featuretype=mRNA.featuretype, start=mRNA.start, end=mRNA.end, score=mRNA.score, strand=mRNA.strand, frame=mRNA.frame, mrnaid=mrnaid, geneid=genename))
            
            numdict = defaultdict(int)
            for child in db.children(mRNA,featuretype=mRNA_children,order_by=("start",'end')):
                numdict[child.featuretype] += 1
                child_id = '{genename}.{childid}{num}'.format(genename=genename,childid=idmap[child.featuretype],num=numdict[child.featuretype])
                f_out.write('{seqid}\t{source}\t{featuretype}\t{start}\t{end}\t{score}\t{strand}\t{frame}\tID={child_id};Parent={mrnaid}\n'.format(seqid=child.seqid, source=child.source, featuretype=child.featuretype, start=child.start, end=child.end, score=child.score, strand=child.strand, frame=child.frame, child_id=child_id, mrnaid=mrnaid))

    f_out.close()
            
def main():
    logging.basicConfig(
        stream=sys.stderr,
        level=logging.INFO,
        format="[%(levelname)s] %(message)s"
    )

    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description="rename gff3 file")

    parser.add_argument('-g', '--gff', required=True, help='gff3 file')
    parser.add_argument('-c', '--change', required=True, help='a file, correspondence between sequence name and gene name prefix')
    parser.add_argument('-a', '--addnum', type=int,default=1, help='diff in gene number, such as if addnum = 10, xx1G000010, xx1G000020')
    parser.add_argument('-p', '--prefix', default='result', help='prefix of output')
    args = parser.parse_args()

    rename(args)


if __name__ == "__main__":
    main()

幫助文檔，用法

python rename_gff.py -h
usage: rename_gff.py [-h] -g GFF -c CHANGE [-a ADDNUM] [-p PREFIX]

rename gff3 file

optional arguments:
  -h, --help            show this help message and exit
  -g GFF, --gff GFF     gff3 file  #輸入注釋文件
  -c CHANGE, --change CHANGE    #序列id和更換前綴之間的對應(yīng)關(guān)系文件
                        a file, correspondence between sequence name and gene
                        name prefix
  -a ADDNUM, --addnum ADDNUM   #相鄰基因編號差值
                        diff in gene number, such as if addnum = 10,
                        xx1G000010, xx1G000020
  -p PREFIX, --prefix PREFIX   #輸出文件前綴
                        prefix of output

我輸入的change.bed文件內(nèi)容如下

Contig1 AT01
Contig2 AT02
Contig3 AT03
Contig4 AT04
Contig5 AT05

最終得到更換名稱后的注釋文件

Contig1 EVM gene    6013    6784    .   -   .   ID=AT01G000001
Contig1 EVM mRNA    6013    6784    .   -   .   ID=AT01G000001.mRNA1;Parent=AT01G000001
Contig1 EVM exon    6013    6447    .   -   .   ID=AT01G000001.exon1;Parent=AT01G000001.mRNA1
Contig1 EVM CDS 6013    6447    .   -   0   ID=AT01G000001.cds1;Parent=AT01G000001.mRNA1
Contig1 EVM exon    6530    6784    .   -   .   ID=AT01G000001.exon2;Parent=AT01G000001.mRNA1
Contig1 EVM CDS 6530    6784    .   -   0   ID=AT01G000001.cds2;Parent=AT01G000001.mRNA1
Contig1 EVM gene    7765    9021    .   +   .   ID=AT01G000002
Contig1 EVM mRNA    7765    9021    .   +   .   ID=AT01G000002.mRNA1;Parent=AT01G000002
Contig1 EVM exon    7765    8006    .   +   .   ID=AT01G000002.exon1;Parent=AT01G000002.mRNA1
Contig1 EVM CDS 7765    8006    .   +   0   ID=AT01G000002.cds1;Parent=AT01G000002.mRNA1

參考資料
http://www.chenlianfu.com/?p=1323
http://daler.github.io/gffutils/