Has anyone ever dived deep into how deterministic

bwa

,

deacon

and

picard

are? Or is bored and interested in a puzzle? 1. I just updated a custom

deacon

module in my pipeline to the latest release in nf-core modules. The underlying container is the same, but there are some slight differences in prefix checks, a tab got changed into four spaces, etc. The main command itself is unchanged. 2. When I run the pipeline on my test data, the

fastq

files filtered by

deacon

files no longer match with the previous ones. I verified with

seqkit sum

that the actual reads are the same, so it must just be the order in which they appear that differs. I saw some speed differences in the

deacon

logs too. 3. The bam files produced by

bwa

also differ, but

gatk CompareSAMs

and

samtools flagstat

do not entirely agree on how they differ, and the latter differences disappear when sorting the file (even though it should already have been sorted in my pipeline now that I think about it). Flagstat sees the largest amount of difference after the bam files have been merged per-lane by

picard

. EDIT: moving logs to thread to save space.

I was a picard developer a decade ago, and we worked hard to ensure that identical input produced consistent output. I currently use a lot of picard tools, but not MarkDuplicates, and I find that still to be the case. I'm a nextflow novice, and one of the things I stumbled over is that a process that combines multiple inputs doesn't get the inputs in a consistent order. Maybe this is an already-solved problem, but this is how I addressed it:

// Collect a channel of tuple(meta, value) in which meta.collectIndex is an integer on which to sort.
// Return the collected values in order by collectIndex.  Return value is a channel containing a single list item.
def collectInOrder(inChannel) {
    def ret = inChannel
        .map { meta, value ->
            tuple(meta.collectIndex, value)
        }
        .toSortedList { a, b -> a[0] <=> b[0] }
        .map { sorted ->
            sorted.collect { pair ->
                pair[1]
            }
        }

    return ret
}

Honestly, my old brain glazes over when I look at the above code, but it appears to do the right thing.

I don't know anything about

deacon

, but if you run it multi-threaded and the output is not explicitly sorted, the relative speed of the threads would determine the order of the reads in the output.

Thanks for the input and advice, @Alec Wysoker! Adding explicit ordering to channels sounds like it's worth exploring! After further consideration, I think the non-deterministic results indeed do not stem from any of the GATK tools, but from the upstream

bwa

step. @Thomas Danhorn: ooh, that sounds plausible! --- I've done some more sleuthing and believe that the differences I'm seeing are due to a documented behaviour of

bwa

. 1. Deacon's output is deterministic in terms of the number of filtered reads, but their order in the fastq files is random (due to multi-threading). Of course, R1 / R2 reads are kept in the same order in their respective files. 2.

bwa

can produce slight differences in output depending on the order of the reads. See https://github.com/malariagen/pipelines/issues/41, https://www.biostars.org/p/238628/, https://github.com/CCDG/Pipeline-Standardization/issues/2 and https://github.com/nf-core/sarek/issues/101. Supposedly the hidden option

-K 1000000

(which I'm using, just like #CGFUX04HZ), should make

bwa

more deterministic, but apparently it's not enough in this case. 3. The downstream outputs of e.g.

picard

consequently differ. Additionally, even if the

bam

outputs were identical, the random order of per-lane bam files in the nextflow channel could lead to different combined bam files (same content, but different order and header). Lastly, I've made sure that the

bam

files are coordinate-sorted in my pipeline before feeding them to

picard

. I've added the output of some smaller test files below:

└─▶ zcat hostremoval/deacon/filter/fastq/106264-002-079-LANE_L008-FC_22LJC2LT4_1.fq.gz | head -n 10000 > sample_L8_R1.fq.gz
└─▶ zcat hostremoval/deacon/filter/fastq/106264-002-079-LANE_L008-FC_22LJC2LT4_2.fq.gz | head -n 10000 > sample_L8_R2.fq.gz

└─▶ for i in *.fq.gz; do seqkit sort ${i} > "$(basename ${i} .fq.gz)".sort.fq.gz; done

└─▶ bwa mem -Y -K 100000000 \
  -R '@RG\tID:106264-002-079.22NY35LT3.L007\tSM:106264-002-079\tLB:106264-002-079\tPL:ILLUMINA\tPU:22NY35LT3.L007' \
  PlasmoDB-68_Pfalciparum3D7_Genome \
  sample_L8_R1.fq.gz sample_L8_R2.fq.gz \
  > sample_L8.K.bam

└─▶ bwa mem -Y \
  -R '@RG\tID:106264-002-079.22NY35LT3.L007\tSM:106264-002-079\tLB:106264-002-079\tPL:ILLUMINA\tPU:22NY35LT3.L007' \
  PlasmoDB-68_Pfalciparum3D7_Genome \
  sample_L8_R1.fq.gz sample_L8_R2.fq.gz \
  > sample_L8.bam

└─▶ bwa mem -Y -K 100000000 \
  -R '@RG\tID:106264-002-079.22NY35LT3.L007\tSM:106264-002-079\tLB:106264-002-079\tPL:ILLUMINA\tPU:22NY35LT3.L007' \
  PlasmoDB-68_Pfalciparum3D7_Genome \
  sample_L8_R1.sort.fq.gz sample_L8_R2.sort.fq.gz \
  > sample_L8.sort.K.bam

└─▶ bwa mem -Y \
  -R '@RG\tID:106264-002-079.22NY35LT3.L007\tSM:106264-002-079\tLB:106264-002-079\tPL:ILLUMINA\tPU:22NY35LT3.L007' \
  PlasmoDB-68_Pfalciparum3D7_Genome \
  sample_L8_R1.sort.fq.gz sample_L8_R2.sort.fq.gz \
  > sample_L8.sort.bam

└─▶ for i in *.bam; do printf "\n$i\n"; samtools flagstat $i; done

sample_L8.K.bam
5090 + 0 in total (QC-passed reads + QC-failed reads)
5000 + 0 primary
0 + 0 secondary
90 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5068 + 0 mapped (99.57% : N/A)
4978 + 0 primary mapped (99.56% : N/A)
5000 + 0 paired in sequencing
2500 + 0 read1
2500 + 0 read2
4762 + 0 properly paired (95.24% : N/A)
4974 + 0 with itself and mate mapped
4 + 0 singletons (0.08% : N/A)
52 + 0 with mate mapped to a different chr
45 + 0 with mate mapped to a different chr (mapQ>=5)

sample_L8.bam
5090 + 0 in total (QC-passed reads + QC-failed reads)
5000 + 0 primary
0 + 0 secondary
90 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5068 + 0 mapped (99.57% : N/A)
4978 + 0 primary mapped (99.56% : N/A)
5000 + 0 paired in sequencing
2500 + 0 read1
2500 + 0 read2
4762 + 0 properly paired (95.24% : N/A)
4974 + 0 with itself and mate mapped
4 + 0 singletons (0.08% : N/A)
52 + 0 with mate mapped to a different chr
45 + 0 with mate mapped to a different chr (mapQ>=5)

sample_L8.sort.K.bam
5090 + 0 in total (QC-passed reads + QC-failed reads)
5000 + 0 primary
0 + 0 secondary
90 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5068 + 0 mapped (99.57% : N/A)
4978 + 0 primary mapped (99.56% : N/A)
5000 + 0 paired in sequencing
2500 + 0 read1
2500 + 0 read2
4762 + 0 properly paired (95.24% : N/A)
4974 + 0 with itself and mate mapped
4 + 0 singletons (0.08% : N/A)
50 + 0 with mate mapped to a different chr
44 + 0 with mate mapped to a different chr (mapQ>=5)

sample_L8.sort.bam
5090 + 0 in total (QC-passed reads + QC-failed reads)
5000 + 0 primary
0 + 0 secondary
90 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5068 + 0 mapped (99.57% : N/A)
4978 + 0 primary mapped (99.56% : N/A)
5000 + 0 paired in sequencing
2500 + 0 read1
2500 + 0 read2
4762 + 0 properly paired (95.24% : N/A)
4974 + 0 with itself and mate mapped
4 + 0 singletons (0.08% : N/A)
50 + 0 with mate mapped to a different chr
44 + 0 with mate mapped to a different chr (mapQ>=5)

The numbers for "mate mapped to a different chr" are ever so slightly off depending on the order of the reads in the fastq files, regardless of the

-K

option.