Hi all, question implying

picard

and

samtools

I am trying to upgrade a lab-inherited, half-constucted workflow. It originally runs:

bwa mem -C -M -t {params.threads} -H {input.sam_header} {params.tmp_ref} {input.fastq1} {input.fastq2} > {params.tmp_sam}
samtools view -@ {params.threads} -S -b {params.tmp_sam} > {params.tmp_bam}
samtools view -@ {params.threads} -bh -q 30 -F 4 -F 8 -F 0X0100 {params.tmp_bam} > {output}
picard FixMateInformation -Xmx4g -Djava.io.tmpdir=tmp I={output} O={output2} SORT_ORDER=coordinate

would

samtools

fixmate

+

sort

create the same output? AFAIK,

picard FixMateInformation

adds MC, recomputes UQ and MQ, but only the first is added by

samtools fixmate

-m. Would that be problematic in following steps of a DNA-seq pipeline (e.g. GATK)?

bwa mem -C -M -t {params.threads} -H {input.sam_header} {params.ref} {input.fastq1} {input.fastq2} \
  | samtools view -@ {params.threads2} -bh -q 30 -F 4 -F 8 -F 0X0100 \
  | samtools fixmate -@ {params.threads2} -m - - \
  | samtools sort -@ {params.threads2} -m 8G -o {output}

Thanks a lot!!

I checked a CRAM from the #CGFUX04HZ pipeline (which of course used GATK), it does not have

UQ

, so I'am pretty sure that is not required.

MC

and

MQ

are there, but since they are just "helpful" but redundant (you can look at the mate to get that information), I would hope that GATK does not actually require them. (AFAIK, GATK does not require a specific mapper, so it should not expect any but the most basic SAM tags.) The definitive proof is of course to just try out your proposed tool chain and see if the downstream processes work with it. Don't know if your workflow provides that option, but with #CGFUX04HZ you can use the

--step

parameter to provide BAMs/CRAMs at any of the processing steps, so you could make a test file and feed it to the pipeline to see if it would work. 🤞 But unless your workflow can do something that #CGFUX04HZ can't do, I would just use the curated, community-supported pipeline ...

Thanks a lot for the reply, Thomas! I’ll definitely try PS: Indeed I’d use sarek, but the pipeline is for scDNA-seq data, in combination with other modalities. That’s why 😅

So you are doing variant calling on single-cell data? How is that going? I tried that once with an scATACseq data set, and the coverage was just too sparse to get much of anything (I was also looking for somatic mutations, which are much harder than germline with sparse data). I didn't bother with a single-cell specific method and just treated the data as bulk and ran sarek on it, mainly to see if anything was salvageable.

Sorry for the delay in replying, Thomas. During the last years I’ve done several variant calling on multiple sc omics data, but I’m cautious to say until which extend one can use the results. Here are some of my thoughts: • One can do variant calling on scRNA-seq (GEX library, no PCR enrichment of specific cDNAs), but only for sample demultiplexing. I used cellsnp-lite + vireo on the bam files after alignment and cell calling in that sense. • One can do variant calling on scRNA-seq for enriched libraries of concrete transcript-derived cDNA molecules. Good examples of that is NPM1 variant calling in AML or MAESTER with mitochondrial mutations • Also, mtDNA genomes are variant called on scATAC-seq (mtscATAC-seq and REDEEM), both by enriching with concrete primers • New technique GIFT allows genotyping of transcripts with high accuracy, but needs specific probes to be designed • In any way, I’d only believe variant calling (with same aims of variant calling on bulk DNA) on cells by profiling with specific scDNA-seq methods (SMART-PTA, targeted scDNA-seq, scWGA, scWGS) Hope it helps 😄

OK, thank you! The data set I looked at was regular 10x "multiomics" (scRNAseq + scATACseq). I find it interesting that you can call mitochondrial variants while that "chromosome" is normally filtered out/blacklisted for other analyses, but it makes sense.