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:
- Genome assembled from nanopore long reads, corrected using Illumina cDNA reads
- 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