HaplotypeCaller on single-end data

[SimonaSecomandi]

Hello,

I am trying to use HaplotypeCaller on "single-end" reads. Specifically, I have ancient DNA (aDNA) data that were sequence as paired-end, but I merge the overlapping reads before aligning the data to the reference genome to increase read length (it's a common step in aDNA data preprocessing).

Would elprep haplotype caller work on these kind of reads? If I add the flag "--single-end" to the command the error is:
2025/02/04 15:33:16 Error: Haplotypecaller (sfm) not supported on single end reads.

Many thanks,

Simona

[caherzee]

Hi Simona,

You can try running with "elprep filter" instead of "elprep sfm". Please be aware that we did not test variant calling for single-end reads. There is some logic in the haplotype caller algorithm that performs certain checks assuming there are paired-end reads, but I think it could still work - but I do not have data to verify this.

From my understanding, the original haplotype caller algorithm can work with single-end reads, but the results are less accurate compared to variant calling for paired-end reads.

We do not currently support the combination of "elprep sfm" (elprep split/filter/merge tools) with single-end reads + haplotype caller. It would require some work on the command-line interface.

The "elprep filter" mode loads all data into memory for running, so, depending on your data set size, this may require more RAM.

If you decide to try with "elprep filter", can you let me know if that works?

Another option may be to use the GATK haplotypeCaller and use elPrep for all other steps.

I am not familiar with ancient DNA. The elprep pipeline expects some standard snp databases as inputs. I do not know if these would be the same for ancient DNA. Are there any other tools for ancient DNA you usually use?

Thank you,

Best,
Charlotte

[SimonaSecomandi]

Hi Charlotte,

thank you for your detailed answer.

I am pretty new to the field of aDNA too and I wanted to test the GATK pipeline. Usually aDNA studies use samtools mpileup or bcftools mpileup for variant calling, especially because the samples are low coverage. However, I have 20-30X of coverage for each of my aDNA samples after merging of overlapping reads, mapping, duplication removal, mapQ filter, and post-mortem damage BAM base recalibration. Therefore, I wanted to use GATK to genotype my high-coverage aDNA samples. Since I'm also working with a pangenome reference and I saw that in the 2023 Chicken pangenome paper they used elprep, I wanted to test it to improve the run time of the analysis.

From my understanding, elprep is calling GATK 4 modules such as HaplotypeCaller, but GATK4 can handle single-end reads (maybe with less accuracy than paired-end) while elprep sfm can't. Is that correct?
And why using elprep filter would be beneficial? I though the difference with elprep sfm was only the splitting by chromosome.

Moreover, for the standard GATK pipeline I am not using a standard SNP dataset as input, while for elprep is it mandatory even when not running BQSR?

Many thanks!

Best,
Simona

[caherzee]

Hi Simona,

elPrep is not calling GATK4. It is its own implementation which produces identical results to GATK4, but executes significantly faster. GATK4 can indeed handle single-end reads, but it warns about lower accuracy of its algorithms compared to paired-end reads.

elPrep sfm currently cannot perform haplotype calling with single-end reads, because we would have to implement support for this. It is currently conflicting with the way we split the input, which is not as simple as splitting by chromosome because ... of reasons explained in our paper :) Our other single-end read users do not use haplotype caller, so it is not something I have looked at.

elPrep filter just runs the pipeline in memory and may function correctly with single-end reads, but I have not tested this. I would be happy to hear your results with this. Normally the advantage of "elprep filter" is faster execution, but it loads all data in memory.

It is not necessary to use any SNP dataset as input. I was just curious about your application.

Thanks!

Best,
Charlotte