Hi all,
I am interested in pulling out "CDS" fasta sequences from a GFF3 file. My GFF3 file lists one or more fasta sequences at the end.
I looked into using the gffutils python package but according to this documentation, it does not deal with fasta sequences listed at the end of the GFF3 file.
Am I stuck writing my own parser, or is there a better solution?
Thanks for any suggestions!
EDITED to add a snippet of my GFF file:
##gff-version 3
##sequence-region Consensus_10_consensus_sequence_1 1 645439
Consensus_10_consensus_sequence_1 feature gene 551 2104 . + . name=gene
Consensus_10_consensus_sequence_1 feature CDS 551 2104 . + 0 name=YALI0C06512p CDS
...
Consensus_9_consensus_sequence_4 feature mRNA 9625 9891 . + . name=mRNA
Consensus_9_consensus_sequence_4 feature CDS 9625 9891 . + 0 name=YALI0A21351p CDS
##FASTA
>Consensus_10_consensus_sequence_1 <unknown description>
TGCCACCTCCAAATTAACTCTCGCTTATTTCTTGTACCTGTCATATCACGTGATGTAGCT
TCCCAATCAAGAGCGGATCCTGCCTGTTTGGCTGCGTGGGTTTGCGTCTTCTTTCCGTTT
GAAGCAGTGGTATTATTCCCCCATTGTGCCAAAAGTAATGCTGAAAAGATGCCCACGAAT
>Consensus_10_consensus_sequence_2 <unknown description>
CCGAAACCACAGCCATGATCAGAGTCACTCCTATTCCAACCCCGCCAACTGCCGTGGGGT
ACACGTTATATTGGGTAGGTGTGTATCCTTGAGACTTGAGCCAGGAGATCATGGAAGGTT
GGGTGCTGGCAATACAGGTGTTATTGTAGCAAAGAAAGATGAGAGAGAAGAAATAAATGT
GCCAGGTTTTCAAAGAGCGTTTCAAAACTTGCAAGAACGGCTCTTCTTTGTTACCGGTAT
CGCCGATTTGTGCTCTTCTCTCCTTGGCTATAGCAATGTCTTCCTTGGTGAAGTACCAGC
TGGTAGTTGTCTCCGGTGTGTTTGGGTTGACGAACATGGTGTAAAGTGCCACTGGAAATG
>Consensus_10_consensus_sequence_2 <unknown description>
CACTCTCAAGCATTTAGGAACTTGTCAAGAGGTTCAAAGGTTGGAACTTGCAACTGAACT
GATCGCAACATAATCACCGCTATTAAACCCTGATTACAAGTGCTTCTGATTGCATCCACA
GTTCATTTCCATGGGCTAGGCTATACGAAAATACAAGGATTAGAAACTATATACAATTGA
CTCTGCAATCTTTCCCGCTAAACGGTGGTGTGGTTATGACCTGGCTCGTGTTCATGGCCG
...
•
written
by
Lina F
3 answers
You should just split your file into two separate files, a GFF and FASTA. Then you can do:
As you have the coordinates of your sequence of interest, you can use samtools faidx.
You first have to make an index of your reference genome using samtools faidx itself after which you can then use the same command to extract your sequence of interest..the usage example available from the command line of samtools faidx is helpful..
If the FASTA sequence is already in the attributes field, it may be a simple parsing step to excise that info.
bedtools getfasta can do that
use -s to force strandness
http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html
bedtools getfasta -fi yourgenome.fa -fo out.fa -bed yourgff.bed -s
Log in to answer this question.
What does that mean?
Please post a snippet of your GFF3 file.