GFF Parsing

From Biopython
Revision as of 13:08, 26 February 2010 by Chapmanb (Talk | contribs)
Jump to: navigation, search


GFF Parsing






    • Limiting parsing to features of interest.
    • Iterating over portions of the file.


Examining your GFF file

Since GFF is a very general format, it is extremely useful to start by getting a sense of the type of data in the file and how it is structured. GFFExaminer provides an interface to examine and query the file. To examine relationships between features, examine a dictionary mapping parent to child features:

import pprint
from BCBio.GFF import GFFExaminer
in_file = "your_file.gff"
examiner = GFFExaminer()
in_handle = open(in_file)

This file contains a flexible three level description of coding sequences: genes have mRNA trascripts; those mRNA transcripts each contain common features of coding sequence, the CDS itself, exon, intron and 5' and 3' untranslated regions. This is a common GFF structure allowing representation of multiple transcripts:

{('Coding_transcript', 'gene'): [('Coding_transcript', 'mRNA')],
 ('Coding_transcript', 'mRNA'): [('Coding_transcript', 'CDS'),
                                 ('Coding_transcript', 'exon'),
                                 ('Coding_transcript', 'five_prime_UTR'),
                                 ('Coding_transcript', 'intron'),
                                 ('Coding_transcript', 'three_prime_UTR')]}

Another item of interest for designing your parse strategy is understanding the various tags used to label the features. These consist of:

  • gff_id -- The record identifier being described. This will often refer to a chromosome or other scaffold sequence.
  • gff_source -- The source description from the second column of the GFF file, which specifies how a feature was generated.
  • gff_type -- The type of the feature, pulled from the 3rd column of the GFF file.
  • gff_source_type -- All combinations of sources and types in the file.

The available_limits function in the examiner gives you a high level summary of these feature attributes, along with counts for the number of times they appear in the file:

import pprint
from BCBio.GFF import GFFExaminer
in_file = "your_file.gff"
examiner = GFFExaminer()
in_handle = open(in_file)
{'gff_id': {('I',): 159,
            ('II',): 3,
            ('III',): 2,
            ('IV',): 5,
            ('V',): 2,
            ('X',): 6},
 'gff_source': {('Allele',): 1,
                ('Coding_transcript',): 102,
                ('Expr_profile',): 1,
                ('GenePair_STS',): 8,
                ('Oligo_set',): 1,
                ('Orfeome',): 8,
                ('Promoterome',): 5,
                ('SAGE_tag',): 1,
                ('SAGE_tag_most_three_prime',): 1,
                ('SAGE_tag_unambiguously_mapped',): 12,
                ('history',): 30,
                ('mass_spec_genome',): 7},
 'gff_source_type': {('Allele', 'SNP'): 1,
                     ('Coding_transcript', 'CDS'): 27,
                     ('Coding_transcript', 'exon'): 33,
                     ('Coding_transcript', 'five_prime_UTR'): 4,
                     ('Coding_transcript', 'gene'): 2,
                     ('Coding_transcript', 'intron'): 29,
                     ('Coding_transcript', 'mRNA'): 4,
                     ('Coding_transcript', 'three_prime_UTR'): 3,
                     ('Expr_profile', 'experimental_result_region'): 1,
                     ('GenePair_STS', 'PCR_product'): 8,
                     ('Oligo_set', 'reagent'): 1,
                     ('Orfeome', 'PCR_product'): 8,
                     ('Promoterome', 'PCR_product'): 5,
                     ('SAGE_tag', 'SAGE_tag'): 1,
                     ('SAGE_tag_most_three_prime', 'SAGE_tag'): 1,
                     ('SAGE_tag_unambiguously_mapped', 'SAGE_tag'): 12,
                     ('history', 'CDS'): 30,
                     ('mass_spec_genome', 'translated_nucleotide_match'): 7},
 'gff_type': {('CDS',): 57,
              ('PCR_product',): 21,
              ('SAGE_tag',): 14,
              ('SNP',): 1,
              ('exon',): 33,
              ('experimental_result_region',): 1,
              ('five_prime_UTR',): 4,
              ('gene',): 2,
              ('intron',): 29,
              ('mRNA',): 4,
              ('reagent',): 1,
              ('three_prime_UTR',): 3,
              ('translated_nucleotide_match',): 7}}

GFF Parsing

Basic GFF parsing

Generally, the GFF parser works similar to other parsers in Biopython. Calling parse with a handle to a GFF file returns a set of SeqRecord objects corresponding to the various IDs referenced in the file:

from BCBio import GFF
in_file = "your_file.gff"
in_handle = open(in_file)
for rec in GFF.parse(in_handle):
    print rec

Since a GFF file is not broken down into an explicit record structure, this requires reading the entire file, parsing all of the features, and then returning those as records. This will be fine for small files, but for most real life cases you will want to restrict parsing to a set of features of interest or a section of lines at once to conserve memory.

Limiting to features of interest

A GFF file will commonly contain many types of features, and you will be interested in retrieving a subset of these. The limit_info argument to GFF.parse allows exact specification of which features to parse, turn into objects and retrieve. An example is retrieving all coding sequence on chromosome 1:

from BCBio import GFF
in_file = "your_file.gff"
limit_info = dict(
        gff_id = ["chr1"],
        gff_source = ["Coding_transcript"])
in_handle = open(in_file)
for rec in GFF.parse(in_handle, limit_info=limit_info):
    print rec.features[0]

You will get back a single record for chromosome 1 which contains all of the coding features in memory for further manipulation. Depending on your memory requirements and workflow, it may make sense to do analyses over each chromosome or set of features you are interested in.

Iterating over portions of a file

from BCBio import GFF
in_file = "your_file.gff"
in_handle = open(in_file)
for rec in GFF.parse(in_handle, target_lines=1000):
    print rec

Providing initial sequence records

from BCBio import GFF
from Bio import SeqIO
in_seq_file = "seqs.fa"
in_seq_handle = open(in_seq_file)
seq_dict = SeqIO.to_dict(SeqIO.parse(in_seq_handle, "fasta"))
in_file = "your_file.gff"
in_handle = open(in_file)
for rec in GFF.parse(in_handle, base_dict=seq_dict):
    print rec

Writing GFF3

The GFF3Writer takes an iterator of SeqRecord objects, and writes each SeqFeature as a GFF3 line:

  • seqid -- SeqRecord ID
  • source -- Feature qualifier with key "source"
  • type -- Feature type attribute
  • start, end -- The Feature Location
  • score -- Feature qualifier with key "score"
  • strand -- Feature strand attribute
  • phase -- Feature qualifier with key "phase"

The remaining qualifiers are the final key/value pairs of the attribute.

A feature hierarchy is represented as sub_features of the parent feature. This handles any arbitrarily deep nesting of parent and child features.

from BCBio import GFF
from Bio import SeqIO
in_file = ""
out_file = "your_file.gff"
in_handle = open(in_file)
out_handle = open(out_file, "w")
GFF.write(SeqIO.parse(in_handle, "genbank"), out_handle)
Personal tools