Minimap2: Mapping cDNA reads

[lh3]

Update: the "cdna" branch has been merged into "master" and subsequently deleted.

---

Original post:

This is a feature request from twitter users. A very preliminary implementation is available at the cdna branch. Some notes:

• perform global alignment (not local hits to each exon)
• introns reported as long deletions (which made possible by the 2-piece gap penalty)
• allow up to 100kb intron by default (longer threshold reduces performance)
• deletions longer than 68bp are considered as introns
• find small exons if nearby exons are long
• splicing boundaries are not accurate
• the cDNA mode is slower due to more DP-based alignment

To use it,

git clone https://github.com/lh3/minimap2
git checkout -b cdna remotes/origin/cdna
make
./minimap2 -a -x cdna ref.fa reads.fa
./minimap2 -a -x cdna -k14 -G 50000 ref.fa reads.fa # parameter tuning

cdna is an experimental branch. If it is not totally crap, I will refine it and merge it to the master. Your feedbacks are warmly welcomed and greatly appreciated. I don't have a good sense about what is a good cDNA alignment for now.

[gringer]

I keep getting index errors in Tablet when working with SAM files produced by minimap2. These only seem to be a problem with Tablet, and I can't see anything obvious in the SAM records to indicate what is wrong.

Sorted BAM file is attached (mapped to Dec. 2011 (GRCm38/mm10) assembly of the mouse genome).

mm2_1.2M_vs_Mmus_mtDNA.bam.gz

[nickschurch]

Just tried this on some oxford nanopore data. was >160x faster than gmap (wow!). The command I used was:

$home/cluster_installs/minimap2/minimap2 -t 12 -ax map-ont /ensembl/release_35/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa sample01_tissue2243.dna.fasta

On out cluster with 12 threads this took just 19s to align 150k long reads. By and large the alignments look very similar to both BWA and gmap with the exception that it seems to struggle with short (i.e.O100bp) exons, which it seems to miss completely. An excellent and fast first step though!

Another problem for us RNA guys (and especially for sharing the data with wet-lab colleagues via a genome browser) is the lack of connectivity between read fragments. Yes they have the same read ID, but because they are fragmented into different reads the genome browsers don show the intron connectivity between the block. The is true of BWA alignments too. It'd be nice to concatenate the fragmented read alignments into one entry in the sam/bam file after the fact.

@gringer I get the same problem with gmap bam output on this data in IGB and the really makes things chug. I suspect these tools are struggling to deal with the very-long cigar strings for the long read data...

[lh3]

@gringer I can open your BAM in IGV. If you can convert SAM to BAM, it is usually the downstream tools to blame for errors like this. SAM to BAM converters do sufficient integrity check to make sure SAM is correctly formatted.

@nickschurch could you try the following cdna mode?

git clone https://github.com/lh3/minimap2
git checkout -b cdna remotes/origin/cdna
make
./minimap2 -a -x cdna -k14 -G 50000 ref.fa reads.fa # parameter tuning

It should partly addresses the two issues you have raised. On the first one, the "cdna" mode will try to align short exons between two long exons, though short exons at the ends of mRNAs still have problems. For these, reducing -k and -w will help. These two options trade performance for sensitivity to short hits. On the second problem, the "cdna" mode will produce much more contiguous alignment in comparison to "map-ont".

If you see the "cdna" mode often miss short exons or fragment alignment, could you send me a few reads/SAM lines? That will greatly help me to tune the parameters. I am cautiously optimistic that minimap2 can be adapted to work with RNA-seq.

[lh3]

By the way, the "cdna" mode will be several times slower than "map-ont" due to the extra work to retain long introns.

[nickschurch]

A great improvement! The -cdna -k14 -G 50000 run looks (based on eyeballing a few examples) much better for shorter exons. Runtime approximately doubled to ~40s (so only ~80x gmap! - same cluster and threads).

It's still struggling with shorter exons (O10s). so it might need even more tweaking. I'm happy with giving you access to the reads to optimise with but we're not sharing them around publicly yet. Wheres the best place to drop you and link to them privately?

I actually really like the idea of representing the splices as long deletions in the sequence by the way - its really quite a simple and accurate description and is pretty intuitive when viewing the alignments in a genome browser.

[nickschurch]

Also, can't see any entry in the man page or the help string for the -G flag - I'm guessing it's different from the -g INT flag?

[lh3]

@nickschurch Thanks! It would be great if you may hand-pick some problem-causing reads, gzip them and send me via email. I don't have other better ways to transfer data privately. Sorry.

Exons of tens of bp will be hard. One thing you could try is to reduce the intron-open cost with -q:

minimap2 -ax cdna -k14 -q 2,40 -G 50000

where the first number is the gap-open cost and the second number is the intron-open cost in effect. In the context of rna-seq mapping, -G is the max intron length. It is set to 100kb by default. I guess A. thaliana does not have many large introns, so suggested a smaller number. Documentations will be updated when this cdna branch gets more mature.

I actually really like the idea of representing the splices as long deletions in the sequence by the way

The latest "cdna" branch will represent introns with the "N" CIGAR operator. This, I believe, is the convention among RNA-seq mappers. Nonetheless, it is also easy to add a command-line option to switch between "N" and "D" for introns.

[nickschurch]

OK, I'll have a look at doing that. I'm on leave until Next Tuesday though so it'll probably be wed/thurs next week to sort them out...

[lh3]

@nickschurch No hurry. I have found some real data to tune minimap2. Thank you anyway.

[gringer]

It seems that the issue may have been in Tablet not knowing what to do when it sees unexpected strings in the fields of a mapped read. Here's the response I got from Gordon Stephen (Tablet developer):

I've been looking into this issue and I believe I have found the cause in code which creates features from CIGAR strings to add to the feature tracks in Tablet. It looks like that code hadn't been updated to handle the sequence of a SAM/BAM read being set to * when that was added to the SAM/BAM spec, unlike the main block of CIGAR parsing code. We're currently testing a solution and will look to push out a bug fix release of Tablet next week.

[lh3]

This cdna branch is concluded: it has been merged into master. The next release of minimap2 will support spliced alignment. As of now, minimap2 still doesn't find precise exon boundaries, but if you allow ~10bp differences, it is slightly more sensitive and specific than gmap on 2D nanopore reads from Byrne et al (2017), and is ~30X faster against the mouse genome. Probably gmap can be tuned for improved accuracy. I haven't tried as I am not familiar with its internal. I will play with how to find precise boundaries, but the feature won't come soon.

The following screenshot shows minimap2 (track r289.bam) and gmap alignment, as well as the truth isoforms. They are very close.

[lh3]

The master branch now tries to find canonical GT-AG.

[nickschurch]

Thanks Heng, this looks great!

On Mon, 14 Aug 2017 at 02:59 Heng Li [e-mail] wrote:

The master branch now tries to find canonical GT-AG.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/lh3/minimap2/issues/12#issuecomment-322084050, or mute
the thread
https://github.com/notifications/unsubscribe-auth/ACHCLhxoeYadcqB_jsdsKrt5HeUnYotfks5sX6ocgaJpZM4Ow92H
.

[gringer]

If you can get it to be codon-sensitive as well, choosing intron lengths (and/or making insertions/deletions) as to reduce the number of stop codons within coding regions of transcripts, I can't see why I'd ever need another cDNA mapper.

... who am I kidding? There aren't any existing mappers that do that anyway (that I'm aware of), so it's not like there's any better competition in that area.

[lh3]

Intron length is unrelated to stop codons.

[lh3]

Ah, I see what you mean, but that would need a genewise-like model which is vey slow. The gain is only marginal because 1) it doesn't work with utrs or non coding rna; 2) the gene finding part may introduce errors as well; 3) transcript alignment can be made fairly accurate - it is just that devs haven't paid enough efforts to optimize for noisy reads.

[nickschurch]

One thing that would be useful for rapid viewing of the alignments in a genome browser like IGB/IGV would be the option to write out the alignments as a gff/gtf file (also an option in gmap). At the moment the browsers really chug with the long read alignments because they are having to parse the long and error-ridden cigar strings. Being able to output the alignments as gff without the insertion/deletion information, alongside the sam file, makes the downstream viewing much more rapid and means we don't have to drill into the alignment details for a read unless we figure its important...

Actually, I'm a bit surprised there isn't a good intron-away bam-to-gff/gtf converter about somewhere, but I can't seem to find one. Maybe writing one independently would be more useful that adding an option to minimap2. I might look into that....

[gringer]

You can do bedtools bamtobed, followed by bedtools bedtobam to remove the CIGAR and sequence information from reads.

[nickschurch]

Yes, but I don't think that preserves the intron information though, right?
Bamtobed definitely ends up with just Exon blocks...

On 16 Aug 2017 7:22 pm, "David Eccles (gringer)" [e-mail]
wrote:

You can do bedtools bamtobed, followed by bedtools bedtobam to remove the
CIGAR and sequence information from reads.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/lh3/minimap2/issues/12#issuecomment-322857485, or mute
the thread
https://github.com/notifications/unsubscribe-auth/ACHCLnUT-iMEFSlv8faaryJ6wt-JNIqSks5sYzNvgaJpZM4Ow92H
.

[lh3]

Minimap2-2.1 is released with the support of spliced alignment. The technical report evaluating the performance will appear on arXiv on the next Monday (or you can compile the latex source in the tex/ directory).

I am closing this issue. You can open new ones if you see bugs or request features. Thank you for all the invaluable input! I really appreciate.