1
$\begingroup$

I would like to subset a GFF file (gene and nested features) from a gene list.

The GFF file looks like this

##gff-version 3
Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903

And a have the target genes in a list such as

Genecv11033552m
Genecv11003131m
Genecv11036683m
Genecv11012576m
Genecv11003654m
Genecv11012587m

I now that is possible to subset the gff using grep -f gene_list.txt <gff_file>. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five_prime_UTR, CDS, three_prime_UTR).

This happens because the ID in CDS and UTR features are the same of the mRNA and not as the gene feature.

Any ideas?

$\endgroup$
1
$\begingroup$

Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython. Dealing with nested data in bash or with some fancy awk script is probably going to trickier to write than it's worth.

You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with mRNA in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.

If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.

Repeat until you reach the end of the GFF file or your list of target genes is empty.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I have actually done it in python using the "flow" you suggested. Thank you! $\endgroup$ – Felipe Almeida Jul 21 at 0:33
  • $\begingroup$ Answer as a python script for this task specifically $\endgroup$ – Felipe Almeida Jul 21 at 12:34
  • $\begingroup$ That script looks good! My only suggestion is that, depending on the GFF file you're loading, you may get slowed down when trying to load the entire file as a pandas DataFrame. If you are trying to parse a multi-gigabyte GFF file, going line-by-line keeps the memory usage low. Otherwise, what you've done looks good $\endgroup$ – James Hawley Jul 21 at 15:04
  • $\begingroup$ Precisely, I have already felt that pain :). I am actually working on a better solution for this. Thanks for the help! $\endgroup$ – Felipe Almeida Jul 21 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.