How to parse GFF3 file and get start and end coordinates in Python

Question

[rimjhim.roy.ch]

I am new to Python. I want to create a new GFF3 file from an existing file by filtering out features that are smaller than 1 kb.

Is there a way to parse a GFF3 file containing feature information of semicolon separated tags as well as start and end of the feature, and create a new file altogether?

My gff file looks like below:

seq_chr1    S-MART    match    158337    160567    .    -    .    Superfamily=LINE;Target=RIL-Map20 356 2619;ID=ms1_seq_chr1_RIL-Map20;Order=TE;Class=Unknown;Identity=93.9881;Name=ms1_seq_chr1_RIL-Map20

Answer 1

[Ryan Dale]

As your analysis gets more complex, you might find gffutils useful for accessing GFF attributes. It has a fairly robust GFF/GTF parser. Here's an example that only keeps >1kb features that are in the LINE superfamily (as you suggest is your goal in a comment on another answer).

	from gffutils.iterators import DataIterator
	input_filename = 'example.gff'
	output_filename = 'output.gff'
	with open(output_filename, 'w') as fout:
	for feature in DataIterator(input_filename):
	# len() works to get the length of a feature in bp
	if len(feature) < 1000:
	continue
	# Attribute access with dictionary-like syntax
	if 'LINE' not in feature['Superfamily']:
	continue
	fout.write(str(line) + '\n')

view raw gffparse_example.py hosted with ❤ by GitHub

(edit: disclaimer, I'm the primary author)

Answer 2

[Daniel]

This is a simple perl script to filter by size which I had lying around which could do that, but python would be just as simple to learn for the long run. I recommend www.pythonforbiologists.com

#!/usr/bin/perl
use strict;
use warnings;

my $usage = "## USAGE:  filter_csv.pl infile.txt [min value]\n";

my $infile = $ARGV[0] or die $usage;
my $filter = $ARGV[1] or die $usage;

open INFILE, "<", "$infile" or die "$usage | check $infile\n";

my $n=0;

while (my $line = <INFILE>){
    $n++;
    if ($n == 1){
        ## Don't process on header line (Maybe you don't need this bit)
        print $line;
    }else{
        my @l = split(/\t/, $line);
        if ($l[4] - $l[3] >= $filter){
            print $line;
        }
    }
}

You would run it like so:

./filter_csv.pl my_file.gff 1000

You could add in a filter for the right term by using a similar if clause and adding it in the size loop.

EDIT: For example (but not tested):

while (my $line = <INFILE>){
    $n++;
    if ($n == 1){
        ## Don't process on header line (Maybe you don't need this bit)
        print $line;
    }else{
       my @l = split(/\t/, $line);
       if ($l[4] - $l[3] >= $filter) && ($line ~= /$term/`){
            print $line;
        }
    }
}

Answer 3

[Daniel Standage]

Since you asked for Python, this simple script should work.

import sys
for line in sys.stdin:
    line = line.rstrip()
    values = line.split("\t")
    start, end = int(values[3]), int(values[4])
    length = end - start + 1
    if length >= 1000:
        print line

It was not tested on a file containing comments or pragmas, but that should be easy to account for.

In your question, you mention semicolon-separated attributes. If you're trying to further process the data (rather than just filtering out lines < 1 kb), you'll need to be more specific about what you need.

[rimjhim.roy.ch]

Thanks for your reply. I implemented the method and it worked. Yes, I also need to keep everything with attribute SuperFamily=LINES and filter out everything else. I hoped there was an implementation with Biopython, but I guess I can work with re and search.

[Daniel Standage]

No need for regular expressions or the search function there: a simple change will suffice. Change if length >= 1000: to if length >= 1000 and "SuperFamily=LINES" in values[8]:.