Skip to content

lrtk ALIGN fails in bam_stats step #7

Description

@Leelaa30

I'm using LRTK to align stLFR reads to a reference genome using the following command:

lrtk ALIGN -FQ1 /scratch/project/adna/Papaya/Murdoch_sequencing/Papaya_stlfr/PAPHI001.R1.fq.gz -FQ2 /scratch/project/adna/Papaya/Murdoch_sequencing/Papaya_stlfr/PAPHI001.R2.fq.gz -R /scratch/project/adna/Papaya/Murdoch_sequencing/papaya_ref_grade/PAPHI003/hap2/PAPHI003_hap2.fasta -O /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/PAPHI001/PAPHI001.rg.bam -RG "@RG\tID:PAPHI001.E150016243\tSM:PAPHI001\tLB:PAPHI001.E150016243\tPU:E150016243\tPL:MGI\tPM:DNBseq_T7\tCN:MU_SABC" -T 10 -P stLFR 

However, the alignemnt fails with this error:

The output directory is: /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/PAPHI001
@RG\tID:PAPHI001.E150016243\tSM:PAPHI001\tLB:PAPHI001.E150016243\tPU:E150016243\tPL:MGI\tPM:DNBseq_T7\tCN:MU_SABC
[2026-03-04 15:32:12.339734] sort_barcoded_FQ1 starts
[2026-03-04 15:32:45.481001] sort_barcoded_FQ1 ends
[2026-03-04 15:32:45.481069] sort_unbarcoded_FQ1 starts
[2026-03-04 15:33:42.542408] sort_unbarcoded_FQ1 ends
[2026-03-04 15:33:42.542474] sort_barcoded_FQ2 starts
[2026-03-04 15:34:15.231935] sort_barcoded_FQ2 ends
[2026-03-04 15:34:15.232026] sort_unbarcoded_FQ2 starts
[2026-03-04 15:35:12.082798] sort_unbarcoded_FQ2 ends
[2026-03-04 15:35:12.086547] EMA_alignment starts
BWA initialization...
[M::bwa_idx_load_from_disk] read 0 ALT contigs
Processing reads...
[2026-03-04 15:35:12.554866] EMA_alignment ends
[2026-03-04 15:35:12.554966] BWA_alignment starts
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[main] Version: 0.7.18-r1243-dirty
[main] CMD: bwa mem -1 -R @RG\tID:PAPHI001.E150016243\tSM:PAPHI001\tLB:PAPHI001.E150016243\tPU:E150016243\tPL:MGI\tPM:DNBseq_T7\tCN:MU_SABC -t 8 /scratch/project/adna/Papaya/Murdoch_sequencing/papaya_ref_gr
ade/PAPHI003/hap2/PAPHI003_hap2.fasta /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/PAPHI001/tmp.unbarcoded.1.fq /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/PAP
HI001/tmp.unbarcoded.2.fq
[main] Real time: 0.098 sec; CPU: 0.099 sec
[2026-03-04 15:35:12.659517] BWA_alignment ends
[2026-03-04 15:35:12.659581] merge_bam starts
[2026-03-04 15:35:12.667025] merge_bam ends
[2026-03-04 15:35:12.667086] sort_bam starts
[2026-03-04 15:35:12.675469] sort_bam ends
[2026-03-04 15:35:12.675539] Mark_bam starts
[2026-03-04 15:35:15.320830] Mark_bam ends
[2026-03-04 15:35:15.320900] index_bam starts
[2026-03-04 15:35:15.328093] index_bam ends
[2026-03-04 15:35:15.328277] bam_stat starts
Illegal division by zero at /home/leelamanoharan/miniforge3/envs/genomics/lib/python3.9/site-packages/script/AlignStat_WGS.pl line 84.
[2026-03-04 15:35:15.366629] bam_stat fails: perl /home/leelamanoharan/miniforge3/envs/genomics/lib/python3.9/site-packages/script/AlignStat_WGS.pl -in /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/
aligned_reads/PAPHI001/PAPHI001.rg.bam -sam PAPHI001.rg.bam -samtools /home/leelamanoharan/miniforge3/envs/genomics/bin/samtools -outdir /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/P
APHI001 -move /scratch/project/adna/leela/FB_SNP_calling_10_10_2025/aligned_reads/PAPHI001

Also, note that the documentation didn't say anything about indexing the genome with bwa, so I had to figure it out myself.

I would appreciate any help or suggestions.

Thanks, Leela

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions