tools to reconcile experimental transcripts with reference annotation

Question

Looking for tools to reconcile alignment file of experimental transcripts mapped to genome (SAM/BAM) with the reference transcriptome annotation (GTF) from Ensembl (organism: D. melanogaster).

The aim would be to check which transcripts reported in the alignment file are also reported in the reference annotation, generating a new annotation (e.g., GTF) including the new transcripts as seen in the alignment file (experimental transcripts to genome).

Any idea?

-- Edit:

Up to now, I tried an in-house script that does the following:

1. generate a BED file from the BAM file with the exon coordinates (start/end) per transcript - using BEDOPS v2p4p27
2. convert the reference GTF to BED to retrieve the coordinates per annotated feature
3. compare the experimental vs the reference BED intervals
4. infer the suite of annotated exons per experimental transcript, thus, ideally, the isoform it corresponds to

However, with this simple approach, I only get around 20% of the experimental transcripts that can be identified in all their length using reference annotation. This might be due either to a too simple approach (e.g. if the start/end of the exon in the experimental transcript is shifted by 1 base as compared to the reference, it won't be detected), or to a real biological signal (since I am also using data from long-read sequencing technologies).

Answer 1

I've never tried this myself, so I don't know how easy this is...

One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You should then be able to get the information you want from the transfrag codes.

Answer 2

As per my answer to @_julien_roux on twitter:

Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA: http://pasapipeline.github.io/#A_ComprehensiveTranscriptome

We did something similar in much simpler organisms in our recent paper: http://dx.doi.org/10.1186/s12864-017-3505-0

Your biggest problems are going to be alternative spliced gene models, the large amount of fragmentation inherent in de novo transcriptome assembly and largely underreported non-coding RNAs in most genomes. It was easier for us as the dictyostelids have few and small introns - fly is going to be a different ball-game. Be prepared to spend a lot of time optimising parameters before finding a result you're happy with (maybe) ;)

You may think the fly genome pretty good and complete, but be prepared to find lots of errors, corrections and novelty. We were very surprised at the amount of new annotation we were able to find in D. discoideum - an organism where every single gene had been manually curated!

Good luck!

Answer 3

Thanks to @Devon_Ryan's suggestion, I had a go at using GMAP with our own de-novo assembly. Here's what I started with:

1. Genome assembled from nanopore long reads, corrected using Illumina cDNA reads
2. genome-guided Trinity transcriptome assembly, with reads grouped based on mapping to the assembled genome, then assembled by Trinity

The transcriptome assembly was filtered to only include transcripts that had reasonable coverage from cDNA reads, and further filtered to cluster similar transcripts and only use the longest isoform. This process is detailed in our paper. If I were to do that again, I'd probably create a supertranscript containing all exons, rather than the longest isoform.

I first created an index for the genome using gmap_build:

gmap_build -k 15 -j 8 -d 'MINBRA' \
  Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta -D MINBRA_GMAP

I then used gmap to map the transcriptome to the genome, generating gene models as a GFF3 file:

gmap -t 10 -d MINBRA -D ./MINBRA_GMAP \
  -A locusTagged_CDLI_CD98LMOHC50_TBNOCFED.tran.fasta \
  -f 2 > Nb_TrinTran_vs_Canu.gff3

Here's an example from the annotation file that was produced by GMAP:

...
tig00022147 MINBRA  gene    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path1;Name=MINBRA_TGG_483_c1_g1_i1
tig00022147 MINBRA  mRNA    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path1;coverage=100.0;identity=99.9;matches=1452;mismatches=0;indels=2;unknowns=0
tig00022147 MINBRA  exon    413327  414778  99  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00022147 MINBRA  CDS 413417  413557  100 +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00001947 MINBRA  gene    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path2;Name=MINBRA_TGG_483_c1_g1_i1
tig00001947 MINBRA  mRNA    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path2;coverage=100.0;identity=98.7;matches=1435;mismatches=1;indels=18;unknowns=0
tig00001947 MINBRA  exon    37526   38961   98  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00001947 MINBRA  CDS 37614   37750   97  +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00022147 MINBRA  gene    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.path1;Name=MINBRA_TGG_489_c13_g1_i2
tig00022147 MINBRA  mRNA    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.path1;coverage=100.0;identity=99.3;matches=907;mismatches=0;indels=6;unknowns=0
tig00022147 MINBRA  exon    464550  465456  99  +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.exon1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 913 1 .
tig00022147 MINBRA  CDS 464909  465053  98  +   0   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.cds1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 507 361 .
###
...

GMAP told me on standard output what transcripts were not mappable to the genome (136 of 36830, or about 0.4%).

I've now used STAR to map these cDNA reads to the transcripts (using the GFF3 annotation file for initial evidence) in order to confirm that the mapping carried out by GMAP is reasonable. I ended up throwing out another 7.4% of transcripts. First, I generated a genome index using the generated gff3 file:

~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
   --runMode genomeGenerate --genomeDir ./star_genome \
   --genomeFastaFiles Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta \ 
   --sjdbGTFfile Nb_TrinTran_vs_Canu.gff3 \
   --sjdbGTFtagExonParentTranscript Parent

Followed by mapping my RNASeq experiments using STAR. STAR can sort the genome mapping, but not the transcriptome mapping, so an additional sort/index was done on that:

for x in $(seq 1 9); do
  ~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
     --readFilesCommand zcat --genomeDir ./star_genome \
     --readFilesIn *-1563-0${x}-*.fastq.gz \
     --outFileNamePrefix mapped/1563-0${x}_vs_BNOCFED \
     --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM
  samtools sort mapped/1563-0${x}.toTranscriptome.*.bam \
    > mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
  samtools index mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
done

BAM hits were collated into a single file:

(for x in $(seq 1 9); 
    do samtools idxstats 1563-0${x}.toTranscriptome.sortedByCoord.out.bam | \
      perl -pe "s/^/1563-0${x}\t/"; 
 done) | gzip  > readCounts_all_vs_transcriptome.tsv.gz

Which was then processed using R to identify unmapped transcript annotations:

#!/usr/bin/Rscript
data.df <-
    read.delim("readCounts_all_vs_transcriptome.tsv.gz", header=FALSE,
               col.names=c("Sample","Transcript",
                           "Length","Mapped","Unmapped"));
gene.counts <- tapply(data.df$Mapped + data.df$Unmapped,
                      list(data.df$Transcript, data.df$Sample), sum);
sum(rowSums(gene.counts)==0) / nrow(gene.counts);
TranscriptsToRemove <- rownames(gene.counts)[rowSums(gene.counts) == 0];
cat(TranscriptsToRemove, file="transcripts_to_remove.txt", sep="\n");

The unmapped annotations were then excluded from the gff3 file:

grep -v -F -f mapped/transcripts_to_remove.txt Nb_TrinTran_vs_Canu.gff3 \
  > mappable_Nb_TrinTran_vs_Canu.gff3