Hello all! I have a technical question related to RNA-seq variant calling using STAR and mutect2 or Haplotypecaller. Considering that mutect2 filter by default reads with a MAPQ < 20, and STAR multimapped reads have always a MAPQ < 20 (due to the used formula

int(-10 * log10(1 - 1/NH))

). Should not be better to set the

--outFilterMultimapNmax 1

to avoid saving multimapped reads that will never be used by mutect2? It will be useful to save memory and reducing storage cost. Is there any contraindication that I am not considering?

That does sound smart. what's your take @Nicolas Vannieuwkerke

Sure, sounds good, as long as it doesn't break haplotypecaller

Thanks a lot for the quick reply!

Be careful with STAR multimappers -- the default value of

--outFilterMultimapScoreRange

is

1

(and as far as I can see this is not being changed in the pipeline). This means that a read may map perfectly to one locus and with 1 mismatch to one or more other loci, and it will be counted as multimapping to all these loci. Maybe that makes some sense for RNAseq, but for variant calling I consider that suboptimal, because you have the choice of not counting multimappers at all (and throwing out perfectly valid reads that way) or of counting them and thus detecting "variants" that are just mapping errors. In my opinion this STAR parameter should be set to

0

in this pipeline.

Thanks @Thomas Danhorn! It sounds a very good point.