Hi all! We recently ran a small experiment comparing variant calls across different nf-core pipelines and mappers using the same sample. Specifically: • nf-core/sarek (germline workflow, BWAMEM2 mapper) • nf-core/rnavar (STAR mapper) Both were run with HaplotypeCaller followed by SnpEff annotation. Surprisingly, we observed substantial differences in the number of variants identified between the two approaches. nf-sarek identified 871553 variants, whereas nf-rnavar identified 22407 variants. Has anyone seen similar behavior before, or do you have ideas on why such discrepancies might arise? We’d really appreciate insights on potential sources of divergence. Thanks a lot in advance!

Was this a WGS sample? In that case this makes a lot of sense since RNA calling only happens on a subset of the genome

Apologies, the question was missing this critical context; it was indeed a DNA-seq WGS sample

yes I'd say this makes sense in this case 😛

Thank you! We'll try next with a WES sample, we should expect similar magnitudes in terms ammount of variants between the two runs right?

I'd say so, it might be a little different depending on your used WES regions but it should be fairly similar I think

Or you could potentially give a WES bed file to sarek?

Interesting idea

Thanks for the suggestion! I'll get my hands on a WES sample and re-run both pipelines, creating the WES BED file for sarek

Check out https://www.biorxiv.org/content/10.1101/2022.09.18.508404v1

image.png

Do you have stats about the variant calling afterwards?

Totally! We know STAR and BWA were designed for very different purposes, but we still wanted to run them side by side to see what discordances show up when using otherwise similar germline workflows on the same sample.

Nope I didn't test that

Hello again! I tested what was suggested here (supplying the same intervals BED file to both sarek and rnavar pipelines) and now the results from both pipelines are much more similar. I'm attaching an "overlap barplot" to illustrate the similarity of the results. There are still variants called only by one of the two pipelines and not by the other, but before supplying the same intervals file, the overlap was minimal. To me, these new results look reassuring, but wanted to consult them here with you as well because there might be some obvious thing I'm missing that makes the overlap not total or higher. What do you think? Thanks in advance!

What reference did you align the data to?

For both pipelines I used the GATK.GRCh38. For sarek I supplied it with the

genome

param but for rnavar I could not since it does not support this reference (I also wanted to ask why?) so I supplied the individual reference files (fasta, GTF, etc) that I took from sarek

Does this reference include alt, HLA, decoys etc?

It does, it actually includes tons of decoys

I’m mostly worried about alt contigs. Depending on your aligner and settings, reads can (primarily) align to a different contig, resulting in a different coverage profile and different variant calls downstream

Got it, it also has alt contigs and HLA loci, to account for natural variation, so this makes a lot of sense 🤔

since you’re benchmarking for concordance, and not accuracy, I would suggest you use a no-alt reference

Thanks for the idea, I'll do that and hopefully results match even more!

Nice

I might not be following quite well but one can compare two VCFs in variantbenchmarking. one VCF can be considered as truth the other can be query. Matching/mismatched variants will be output then too (I have upset plots and indel distribution plots)

Hello all! I'm writing in this thread because the question relates to this topic! I noticed that nf-core/sarek and nf-core/rnavar use two different approaches for variant filtering. 1. nf-core/rnavar --> VariantFiltration (Hard-filtering variants based) 2. nf-core/sarek --> CNNScoreVariants (Annotating variants with a score derived from a Convolutional Neural Network (CNN) model) + FilterVariantTranches (Applying filtering based on a variant's quality score (the one from CNNScoreVariants)) Is there a reason behind the use of this two different workflows? Thanks a lot in advance!

That's a very good question, and I feel like I should know the answer, and I am ashamed to say that I do not, so I'll ping @Nicolas Vannieuwkerke first

I don't know either, it was like that when I made my first contribution to the pipeline 😅