Answers.Body: <p>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. </p>
  
  <p>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. </p>
  
  <pre><code>perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf &gt; transcripts_fixed.gtf
  </code></pre>
  
  <p>Where:</p>
  
  <ul>
  <li><code>genome.fa.fai</code> is the fasta index of the genome file</li>
  <li><code>transcripts_to_fix.gtf</code> is the transcript set out of the assembler with issues</li>
  <li><code>transcripts_fixed.gtf</code> is the new transcript file with the repaired end coordinates</li>
  </ul>
  
  <p>Content of <code>cufftrim.pl</code>:</p>
  
  <pre><code>#!/usr/bin/perl
  use strict;
  use Getopt::Std;
  
  my $usage = q/Usage:
   cufftrim.pl genome.fa.fai transcripts.gtf &gt; transcripts.trim.gtf
  /;
  umask 0002;
  getopts('o:') || die($usage."\n");
  my $outfile=$Getopt::Std::opt_o;
  if ($outfile) {
    open(OUTF, '&gt;'.$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 (&lt;FAI&gt;) {
   chomp;
   next if m/^#/;
   my @t=split(/\t/);
   die("Parsing error: $t[1] length for sequence $t[0] ?\n") 
      unless $t[1]&gt;10;
   $slen{$t[0]}=int($t[1]);
  }
  close(FAI);
  
  # now process the GTF file
  my ($tcount, $ecount); #trim counts
  while (&lt;GTF&gt;) {
   my @t=split(/\t/);
   my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
   if ($t[4]&gt;$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);
   }
  </code></pre>
  
  <p>Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution. </p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T14:55:38.870Z
Answers.LastActivityDate: 2018-02-27T14:55:38.870Z
Answers.UserDisplayName: aechchiki
Body: <p><strong><em>Disclaimer</strong>: I had this issue posted on <a href="https://groups.google.com/forum/#!topic/tuxedo-tools-users/i79gfggEBAQ" rel="nofollow noreferrer">Tuxedo Tools users group</a> and shared it on <a href="https://twitter.com/aechchiki/status/953632647074516992" rel="nofollow noreferrer">Twitter</a>, 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.</em></p>
  
  <p>Using <a href="http://cole-trapnell-lab.github.io/cufflinks/file_formats/" rel="nofollow noreferrer">gffread</a> utility to convert transcripts from GTF as output of <a href="http://cole-trapnell-lab.github.io/cufflinks/cufflinks/index.html" rel="nofollow noreferrer">Cufflinks</a> to FASTA, I get some <code>GFaSeqGet</code> errors, as follows (here is the full log file):</p>
  
  <pre><code>$ 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)!
  </code></pre>
  
  <p>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.</p>
  
  <p>The alignment file (SAM) was generated using <a href="https://ccb.jhu.edu/software/hisat2/index.shtml" rel="nofollow noreferrer">Hisat2</a>:</p>
  
  <pre><code>$ 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
  </code></pre>
  
  <p>Then sorted:</p>
  
  <pre><code>$ cat out.sam | sort -k 3,3 -k 4,4n -T . &gt; out_sorted.sam
  </code></pre>
  
  <p>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. </p>
  
  <p>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:</p>
  
  <pre><code>$ cat transcripts.gtf | awk '$3=="transcript"{print $0}' | wc -l
  231139
  $ cat transcripts.fasta | grep '^&gt;' | wc -l
  218075
  </code></pre>
  
  <p>Do you know what would be a cause of this problem?</p>
  
Tags: <rna-seq><phylogenetics><software-recommendation><base-clipping>
Title: gffread: GFaSeqGet errors on coordinate overhang
ThreadText: # gffread: GFaSeqGet errors on coordinate overhang
  
  (Asked by: username_0 on 2018-02-27)
  
  ***Disclaimer***: I had this issue posted on [Tuxedo Tools users group](https://groups.google.com/forum/#!topic/tuxedo-tools-users/i79gfggEBAQ) and shared it on [Twitter](https://twitter.com/username_0/status/953632647074516992), 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](http://cole-trapnell-lab.github.io/cufflinks/file_formats/) utility to convert transcripts from GTF as output of [Cufflinks](http://cole-trapnell-lab.github.io/cufflinks/cufflinks/index.html) 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](https://ccb.jhu.edu/software/hisat2/index.shtml):
  
  ```
  $ 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?
  
  **Tags:** rna-seq, phylogenetics, software-recommendation, base-clipping
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2018-02-27. Score: 2 | Accepted 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;
  
  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.
  
  
  ---
  
  