Does WASP allele-specific mapping work with long-reads as input?

[dmanahan1]

Hello,

Question as written in the subject line. I have long-read RNA-Seq data (aligned with minimap2 and flair) and I need to do allele-specific mapping. From what I understand, WASP is designed to work with bowtie2 as an aligner.

[sparthib]

Hi @dmanahan1,

Have you been able to figure this out?

Thanks,
Sowmya

[dmanahan1]

Hello Sowmya,

Some good news and bad news. Bad news, no, I did not figure this out. Or at least, not with WASP.

Talking to colleagues who used WASP for allele-specific mapping with short-reads, it seemed that they didn't use WASP directly, but rather it was a feature of STAR that implements WASP. I am trying to recall off the top of my head, but the feature was basically tagging reads as coming from one allele or another and then sorting them into allele-specific BAMs. As far as I know, there are no long-read mappers that implement WASP functions.

Good news is that I did end up successfully parsing my BAMs corresponding to their respective alleles using whatshap haplotag followed by whatshap split. In my case, I was sorting reads from F1 hybrids into their respective parental alleles and had a phased VCF as reference. I had to phase my F1 VCFs using this reference with whatshap phase.

I understand that my analysis' goals and starting data may be different enough, but replying with this information in case it is useful to you.

-David

[sparthib]

Hi @dmanahan1, thank you, this is helpful to know!

Could I ask what type of data you are creating your VCFs from? short-read-WGS? I am assuming both the WGS and long-read RNA-Seq was generated from the same sample? Please correct me if I am wrong.

I am also curious what percentage of your reads were haplotagged using whatshap?

Thanks,
Sowmya

[dmanahan1]

The genomic VCF containing fixed differences that were used for tagging was previously generated by a colleague in the lab. It was probably generated using GATK variant caller.

The VCF of my long-read expression data was generated with longshot. Longshot by design tries to do haplotype-aware variant calling and phasing, so you'd expect it to handle the allele-specific reads on its own. However, it tries to haplotype based on all variants, and for my data it was only able to phase 6% of my reads. Longshot was not able to take in that reference VCF with fixed differences (that I described above) as input. Hence, I just called variants in my long-read expression data and generated a VCF with longshot, and then pushed that VCF through whatshap to tag the reads in my BAM and then parse the reads into allele-specific BAMs.

I was able to split about 10% of my F1 hybrid reads with this method, though that number will change depending on how diverged your haplotypes are. For what it's worth, my colleagues had the same rate of allele-specific mapping using short-read data from the same animals (processed by WASP implemented via STAR).

[sparthib]

Thanks so much @dmanahan1 , this is very helpful!

Just so I am getting this right, you ended up using longshot on your ONT RNASeq data to variant call, and then passed the VCF and bam files to whatshap.

When you say splitting 10% of your reads, do you mean that you ended up phasing 10% of your reads?

Sowmya

[dmanahan1]

Yes, that's right. I was working with hybrids between two populations of the same species, so I guess I am not too surprised that that many reads could be assigned to one haplotype vs another.

[sparthib]

Thank you! Is this work published anywhere? Curious what your next steps were.

[dmanahan1]

Sorry, not yet! But soon. I am writing up the manuscript now.

[sparthib]

No worries, good luck!!

[gmcvicker]

Hi, In principal WASP should be able to work with long reads. It can work with any aligner and is not tied to bowtie2. As you note above, STAR has an implementation of the WASP read filtering algorithm and should be fine. A limitation is that if the reads are very long and there are a large number of unphased SNPs overlapping reads, then WASP can run into problems because it has to remap all possible allelic combinations of the read (2*N, where N is the number of SNPs). Currently, reads that overlap > 10 SNPs that are polymorphic within your sample are discarded.

[gmcvicker]

You might also consider trying WASP2, which should be easier to use and is faster.
https://github.com/mcvickerlab/WASP2

[sparthib]

Thanks for pointing me to WASP2 @gmcvicker!