gffread: GFaSeqGet errors on coordinate overhang

Question

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?

Answer 1

Basically, the source of the problem is issues with hard/soft clipped bases in short-reads RNAseq reads alignments to reference genome in post-transcriptome assembly steps.

Truncating the coordinates in the gtf file to match the reference is a solution. For this, I will share a perl solution (non-official) that will do the job.

perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf > transcripts_fixed.gtf

Where:

• genome.fa.fai is the fasta index of the genome file
• transcripts_to_fix.gtf is the transcript set out of the assembler with issues
• transcripts_fixed.gtf is the new transcript file with the repaired end coordinates

Content of cufftrim.pl:

#!/usr/bin/perl
use strict;
use Getopt::Std;

my $usage = q/Usage:
 cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
umask 0002;
getopts('o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
if ($outfile) {
  open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
  select(OUTF);
  }
# --
my $fai=shift(@ARGV) || die ("$usage\nSamtools fasta index required!\n");
my $gtf=shift(@ARGV) || die ("$usage\nCufflinks' transcripts file required!\n");
my %slen;

open(FAI, $fai) || die("Error opening $fai\n");
open(GTF, $gtf) || die("Error opening $gtf\n");
while (<FAI>) {
 chomp;
 next if m/^#/;
 my @t=split(/\t/);
 die("Parsing error: $t[1] length for sequence $t[0] ?\n") 
    unless $t[1]>10;
 $slen{$t[0]}=int($t[1]);
}
close(FAI);

# now process the GTF file
my ($tcount, $ecount); #trim counts
while (<GTF>) {
 my @t=split(/\t/);
 my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
 if ($t[4]>$l) {
   $t[4]=$l;
   $t[2] eq 'exon' ? ++$ecount : ++$tcount;
 }
 print join("\t", @t);
}
close(GTF);
print STDERR "$tcount transcripts, $ecount exons adjusted.\n";
# --
if ($outfile) {
 select(STDOUT);
 close(OUTF);
 }

Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution.