Hello again,
in the pipeline for which I use SpecHLA, I have recently added testing using parallel execution on downsampled alignments. During that addition, I noticed that for one file, a WES tumor sample, the HLA calling was non-determinitstic when run with more than one core. I compared outputs from runs started with identical parameters and found, that the main difference was in the FermiKit de-novo assembly of reads aligning to the respective HLA region.
Relevant section in the realignment script:
|
fermi2.pl unitig -s1k -t $thread -T 10 -2 -l $rlen -p $outdir/prefix2 $outdir/extract.fa > $outdir/prefix2.mak |
This is called in the main genotyping script here:
|
bash $dir/../run.assembly.realign.sh $sample $outdir/$sample.merge.bam $outdir 70 $assemble_region ${num_threads:-5} |
In their Limitations section, the FermiKit developer states that it is known that parallel exectution of FermiKit may lead to different assemblies.
Specifically modifying the realignment script to always run fermikit single-threaded did lead to me no longer finding the alternative genotypes. I did check larger real world runs and did not find similar inconsistencies, however, for me it would be important to rule out as many sources of non-determinism in my pipeline as possible.
Proposed fix
I would propose to add an option to the genotyping script to specifically run the assembly step single-threaded while letting everything else still run in parallel. This would allow users to rule out this source of non-determinism while still allow for parallel execution of all other steps.
Hello again,
in the pipeline for which I use SpecHLA, I have recently added testing using parallel execution on downsampled alignments. During that addition, I noticed that for one file, a WES tumor sample, the HLA calling was non-determinitstic when run with more than one core. I compared outputs from runs started with identical parameters and found, that the main difference was in the FermiKit de-novo assembly of reads aligning to the respective HLA region.
Relevant section in the realignment script:
SpecHLA/script/run.assembly.realign.sh
Line 39 in 9c480a3
This is called in the main genotyping script here:
SpecHLA/script/whole/SpecHLA.sh
Line 244 in 9c480a3
In their Limitations section, the FermiKit developer states that it is known that parallel exectution of FermiKit may lead to different assemblies.
Specifically modifying the realignment script to always run fermikit single-threaded did lead to me no longer finding the alternative genotypes. I did check larger real world runs and did not find similar inconsistencies, however, for me it would be important to rule out as many sources of non-determinism in my pipeline as possible.
Proposed fix
I would propose to add an option to the genotyping script to specifically run the assembly step single-threaded while letting everything else still run in parallel. This would allow users to rule out this source of non-determinism while still allow for parallel execution of all other steps.