Disclaimer: I had this issue posted on Tuxedo Tools users group and shared it on Twitter, but could not get an answer to this from the developers, nor find documentation of this issue online. So, I'll share my solution here below and see if some of you have a better solution to add.
Using gffread utility to convert transcripts from GTF as output of Cufflinks to FASTA, I get some GFaSeqGet
errors, as follows (here is the full log file):
$ gffread transcripts.gtf -g ../genome.fa -w transcripts.fasta
Error (GFaSeqGet): end coordinate (32519) cannot be larger than sequence length 32512
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (6327) cannot be larger than sequence length 6326
Error (GFaSeqGet): end coordinate (4651) cannot be larger than sequence length 4648
Error (GFaSeqGet): subsequence cannot be larger than 3551
Error getting subseq for CUFF.29386.2 (1..3568)!
That's a bit weird since I used the same genome file as the one used to generate the alignment file required by Cufflinks, and the coordinates overhang does not make a lot of sense to me.
The alignment file (SAM) was generated using Hisat2:
$ hisat2
--rna-strandness RF
--dta-cufflinks
-x ./hisat_idx
-q
-1 ../1A_1.fq.gz,../1B_1.fq.gz,../2A_1.fq.gz,../2B_1.fq.gz,../3A_1.fq.gz,../3B_1.fq.gz
-2 ../1A_2.fq.gz,../1B_2.fq.gz,../2A_2.fq.gz,../2B_2.fq.gz,../3A_2.fq.gz,../3B_2.fq.gz
-S out.sam
Then sorted:
$ cat out.sam | sort -k 3,3 -k 4,4n -T . > out_sorted.sam
My question is: are only the few sequences reported in the log affected? Or is it a more general problem? The answer seem to be the second one.
Out of gffread, alongside the log file, I actually got a FASTA file with no additional error messages during its generation (537M out of 404M GTF file) - but there are around 13k transcripts which are present in the GTF but seem to be missing from the FASTA:
$ cat transcripts.gtf | awk '$3=="transcript"{print $0}' | wc -l
231139
$ cat transcripts.fasta | grep '^>' | wc -l
218075
Do you know what would be a cause of this problem?