I got a genome (
data/genome/genome.fasta) and braker-based genome annotation (
data/genome/annotation.gff3), now I would like to get sequences of one transcript per gene (let's say the longest one).
Cufflinks has a
gffread program to get the transcriptome out of the genome and the annotation. Something like
gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa
Will do the job. However, it is pulling multiple transcripts per gene. I read the help page twice, found parameters that could potentially do what I would like to (collapsing alternative transcripts into single one), but even when I specify my candidate parameters
-M -K I get still
.t2 sequences for some genes.
The task is surprisingly simple, and I got already a custom script solution, but it just feels wrong. There must be a better way!
What is a reasonable/elegant way to get nucleotide sequences of genes out of a genome and annotation? Am I missing something with cufflinks?
Related, not exactly identical:
- How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms (dealing more with terminology of Trinity)
- When I used Kallisto/Salmon for differential expression analysis, I collapsed transcripts into genes while importing to R via