gff3 to CDS fasta

Question: gff3 to CDS fasta

2.7 years ago by

n.caillou • 20 wrote:

Given a genomic fasta file and a GFF3 file like this (this is just one transcript):

groupI  AUGUSTUS        gene    24128   27467   1       +       .       ID=g2
groupI  AUGUSTUS        transcript      24128   27467   .       +       .       ID=g2.t1;Parent=g2
groupI  AUGUSTUS        start_codon     24128   24130   .       +       0       Parent=g2.t1
groupI  AUGUSTUS        CDS     24128   24311   .       +       0       ID=g2.t1.cds;Parent=g2.t1
groupI  AUGUSTUS        CDS     25287   25354   .       +       2       ID=g2.t1.cds;Parent=g2.t1
groupI  AUGUSTUS        CDS     27267   27467   .       +       0       ID=g2.t1.cds;Parent=g2.t1
groupI  AUGUSTUS        stop_codon      27465   27467   .       +       0       Parent=g2.t1

Is there a simple way to extract a fasta of the CDS (or of the translated sequence) for each transcript?

I have already seen:

Extract Cds Fastas From A Gff Annotation + Reference Sequence where some people have posted and amended perl scripts to do it
http://biopython.org/wiki/Gene_predictions_to_protein_sequences a script base that uses a non-standard part of Biopython
also I have tried to use cufflinks's gffread like this: gffread file.gff3 -g genome.fa -y proteins.fa but the proteins sequence are nonsense, although the lengths are consistent.. (Am I doing something wrong?)
edit: as well as bedtools getfasta, but I get each element's sequence separately (see my comments below)

Isn't there something more simple and reliable than the above custom scripts?

sequence gff3 genome • 5.6k views

ADD COMMENT • link •

modified 2.7 years ago by A. Domingues • 2.1k • written 2.7 years ago by n.caillou • 20

Did you try bedtools getfasta ?

ADD REPLY • link modified 2.7 years ago • written 2.7 years ago by Santosh Anand ♦ 4.9k

I did, see my answer to WouterDeCoster below. I cannot get per-transcript CDS sequences, I just get each element (gene, transcript, start/end codons, and each CDS fragment) separately. Is this the expected result?

ADD REPLY • link written 2.7 years ago by n.caillou • 20

I see what you mean. This might help you A: Gencode Gtf To Bed12 Or Psl

ADD REPLY • link written 2.7 years ago by Santosh Anand ♦ 4.9k

2.7 years ago by

A. Domingues • 2.1k

Dresden, Germany

A. Domingues • 2.1k wrote:

Are you happy with R/Bioconductor? There is a solution here which could solve your problem. Look for item #6 of 'Alignments (and genomic annotations)'. I will transcribe the relevant (but incomplete) bits of code here for convenience, but all credits go to Martin Morgan.

## cds for all transcripts
cdsByTx <- cdsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "tx")

## reference genome
library(BSgenome.Hsapiens.UCSC.hg19)

dna <- extractTranscriptSeqs(BSgenome.Hsapiens.UCSC.hg19, cdsByTx)
protein <- translate(dna)

ADD COMMENT • link modified 2.7 years ago • written 2.7 years ago by A. Domingues • 2.1k

If I am going to do a non-trivial amount of scripting I'd probably go with Biopython. Also, would this work with any fasta/gff3 pair of files?

ADD REPLY • link written 2.7 years ago by n.caillou • 20

Assuming the gff3 is well formatted, that is, contains the CDS information, yes. Of course, when it comes to bioinformatics assuming a file is well formatted and that things will go as expected it a big ask so ymmv.

ADD REPLY • link written 2.7 years ago by A. Domingues • 2.1k

2.7 years ago by

WouterDeCoster ♦ 40k

Belgium

WouterDeCoster ♦ 40k wrote:

I would suggest to try if bedtools getfasta fit your needs.

ADD COMMENT • link written 2.7 years ago by WouterDeCoster ♦ 40k

Thanks, but actually, I have also tried getfasta. And I can only get the sequences of each element separately.

Would the -split option do what I want? Apparently it will only work if the input is BED. How do I convert my GFF3 to BED? Which other tool do I have to use besides getfasta?

edit: The converter is bedops's gff2bed, I suppose. It didn't help, though (same result).

ADD REPLY • link modified 2.7 years ago • written 2.7 years ago by n.caillou • 20

And can't you pre filter your gff to only contain 'transcript' features, e.g. using grep "transcript" yourfile.gff > onlytranscript.gff?

ADD REPLY • link written 2.7 years ago by WouterDeCoster ♦ 40k

I can do that to have all the CDS and just them (that's part of the "shell tricks" in my answer below), but then I would have to piece them together. Definitely possible but not simple either.

ADD REPLY • link written 2.7 years ago by n.caillou • 20

But the transcript annotation in a gff is the full transcript, entire CDS of a gene, right? I don't understand what you would need to paste together.

ADD REPLY • link written 2.7 years ago by WouterDeCoster ♦ 40k

As I understand the question of n.caillou, he wants to extract coding sequences, not transcript sequence (i.e. he needs to paste together exons for multi-exon genes)...

ADD REPLY • link written 18 months ago by al-ash • 100

2.7 years ago by

Daniel Standage ♦ 3.9k

Davis, California, USA

Daniel Standage ♦ 3.9k wrote:

The xtractore program from the AEGeAn Toolkit should work for this. Set --type=CDS on the command line. For CDSs and other discontiguous features, xtractore pastes the fragments together correctly, taking into account strand, etc.

ADD COMMENT • link modified 2.7 years ago • written 2.7 years ago by Daniel Standage ♦ 3.9k

This would work, thanks. But isn't there something more standard? It is very impractical for me to install new software.

ADD REPLY • link written 2.7 years ago by n.caillou • 20

I'm not sure I understand your constraints. If you are uncomfortable evaluating custom Perl/Python scripts and if installing new software is impractical for you, you may be in the wrong field!

I'm not trying to be witty or smug, I'm being completely frank. The bioinformatics state of the art unfortunately provides very few clean solutions; almost always you must get your hands a bit dirty to piece together a working solution.

ADD REPLY • link written 2.7 years ago by Daniel Standage ♦ 3.9k

Hello Daniel,

I used xtractore to extract CDS from a gff file. It does extract CDS taking into account strand and so on, but it does not paste CDS fragments belonging to the same gene together. Am I missing sth?

Also, the -w option does not seem to work properly:

xtractore -t CDS -w 0 -o timema_cristinae.fa timema_cristinae.gff Timema_cristinae_scaffolds.fas

Setting -w to zero gives me an empty output fasta file. In the documentation it is mentioned: "-w|--width: INT width of each line of sequence in the Fasta output; default is 80; set to 0 for no formatting"

Anyways, thanx for your program, it would take 2-3 days to write a CDS extractor myself :P

ADD REPLY • link modified 14 months ago • written 14 months ago by panosprov • 0

2.7 years ago by

n.caillou • 20

n.caillou • 20 wrote:

Using bedops and bedtools and shell tricks I can get a TSV file like:

g2.t1.cds   ATGAAGACGACGCTTGC...
g2.t1.cds   AGCGGAGTACCACGTGT...
g2.t1.cds   CACCACCGCGACCCC...TGA

But it is just a pile of hacks, and it is not even finished.

ADD COMMENT • link modified 2.7 years ago • written 2.7 years ago by n.caillou • 20

Please log in to add an answer.

Similar posts • Search »