Bioinformatics Asked by Felipe Almeida on August 21, 2021
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 know 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 in the gene feature.
Any ideas?
I think this simple script that uses pyranges solves it. Half of the below is just setup to make a minimal reproducible example.
# pip install pyranges
# or
# conda install -c bioconda pyranges
import pandas as pd
from io import StringIO
import pyranges as pr
#########
# setup #
#########
gff_name = "gtfo.gtf"
contents = StringIO("""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 gene 22901 45904 . + . ID=Gene2.g;Name=Gene2.g
Scaffold_1 JGI mRNA 22901 45904 . + . ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
Scaffold_1 JGI five_prime_UTR 22901 23284 . + . ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
Scaffold_1 JGI CDS 23285 23423 . + 0 ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ""")
pd.read_table(contents, sep="s+", header=None).to_csv(gff_name, sep="t", index=False, header=False)
genes_to_keep = set(["Genecv11000001m.g"])
############
# solution #
############
df = pr.read_gff3(gff_name, as_df=True)
def all_descendants(df, genes_to_keep):
old_len, new_len = -1, len(genes_to_keep)
while old_len != new_len:
genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
old_len, new_len = new_len, len(genes_to_keep)
return df[df.ID.isin(genes_to_keep)]
desired_gff = all_descendants(df, genes_to_keep)
print(desired_gff)
# Chromosome Source Feature Start End Score Strand Frame ID Name longest Parent
# 0 Scaffold_1 JGI gene 22900 45904 . + . Genecv11000001m.g Genecv11000001m.g NaN NaN
# 1 Scaffold_1 JGI mRNA 22900 45904 . + . PAC4GC:50510902 Genecv11000001m 1 Genecv11000001m.g
# 2 Scaffold_1 JGI five_prime_UTR 22900 23284 . + . PAC4GC:50510902.five_prime_UTR.1 NaN NaN PAC4GC:50510902
# 3 Scaffold_1 JGI CDS 23284 23423 . + 0 PAC4GC:50510902.CDS.1 NaN NaN PAC4GC:50510902
# 4 Scaffold_1 JGI CDS 24030 24062 . + 2 PAC4GC:50510902.CDS.2 NaN NaN PAC4GC:50510902
Correct answer by The Unfun Cat on August 21, 2021
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.
Answered by James Hawley on August 21, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP