# gffread: GFaSeqGet errors on coordinate overhang

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
-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? ## 1 Answer 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; myusage = q/Usage:
cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
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.