Skip to content

Experiment/enable ntm variant calling #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
bobturneruk wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Experiment/enable ntm variant calling #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bf66caa
feat: support reference input in channels for ntm variant calling
bobturneruk Nov 12, 2025
eee063c
test: update for ntm variant calling channels
bobturneruk Nov 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,10 @@ This is included in the `data` folder. Just run `get_clair3_models.sh` to downlo
Can run with bcftools or with clair3
```bash
nextflow run . --workflow bcftools --input_dir test_data/assemblers --publish_dir results \
--ref-fasta data/h37rv_20231215.fa.gz
--ref-fasta data/NC_000962.3.fa.gz

nextflow run . --workflow clair3 --input_dir test_data/assemblers --publish_dir results \
--ref-fasta data/h37rv_20231215.fa.gz --clair3_models_dir data/clair3_models \
--ref-fasta data/NC_000962.3.fa.gz --clair3_models_dir data/clair3_models \
--basecalling_model [email protected]
```

Expand Down
0 data/h37rv_20231215.fa.gz → data/NC_000962.3.fa.gz
View file
Open in desktop
File renamed without changes.
0 data/h37rv_20231215.fa.gz.fai → data/NC_000962.3.fa.gz.fai
View file
Open in desktop
File renamed without changes.
69 changes: 44 additions & 25 deletions main.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,22 +26,32 @@ workflow {
exit(0)
}

input_files = Channel.fromPath("${params.input_dir}/${params.input_single_suffix}", checkIfExists: true)
input_files = Channel
.fromPath("${params.input_dir}/${params.input_single_suffix}", checkIfExists: true)
.ifEmpty { error("cannot find any reads matching ${params.input_single_suffix} in ${params.input_dir}") }
.map { it -> tuple(it.getName().replaceFirst(/(?i)\.(fastq|fq)\.gz$/, ""), it) }

input_files.take(3).view()

ref = Channel.fromPath("${params.ref_fasta}").first()
read_ch = input_files.map { it ->
[
it[0],
it[1],
file(params.ref_fasta),
params.reference,
params.accession,
]
}

clair3_models_dir = params.clair3_models_dir.startsWith("/")
? params.clair3_models_dir
: "${projectDir}/${params.clair3_models_dir}"

if (params.workflow == "clair3") {
rundial(input_files, ref, clair3_models_dir, params.basecalling_model)
rundial(read_ch, clair3_models_dir, params.basecalling_model)
}
else if (params.workflow == "bcftools") {
rundial_with_bcftools(input_files, ref)
rundial_with_bcftools(read_ch)
}
else {
error("Please provide a valid workflow: bcftools or clair3")
Expand All @@ -52,74 +62,83 @@ workflow {
workflow rundial {
take:
fastq_files
ref
clair3_models_dir
dorado_model

main:

ref_name_ch = fastq_files.map { it -> [it[0], it[2]] }

ref_path_ch = fastq_files.map { it -> [it[0], it[4]] }

fastq_files = fastq_files.map { it -> tuple(it[0], it[1], it[4]) }

clair3_model = get_clair3_model(dorado_model)

valid_model_ch = fastq_files
.combine(clair3_model)
.filter { it -> it[2] != "failed" }
.map { it -> tuple(it[0], it[1]) }
.filter { it -> it[3] != "failed" }
.map { it -> tuple(it[0], it[1], it[2]) }
invalid_model_ch = fastq_files
.combine(clair3_model)
.filter { it -> it[2] == "failed" }
.map { it -> tuple(it[0], it[1]) }
.filter { it -> it[3] == "failed" }
.map { it -> tuple(it[0], it[1], it[2]) }

valid_model_ch.take(1).view { "Appropriate Clair3 Model Found" }
invalid_model_ch.take(1).view { "No appropriate Clair3 Model Found" }

// If no appropriate clair3 model is found, use bcftools
rundial_with_bcftools(invalid_model_ch, ref)
rundial_with_bcftools(invalid_model_ch)

// If appropriate clair3 model is found, run clair3
model_path = clair3_model.map { it -> clair3_models_dir + "/${it}.tar.gz" }
model_path.view { "Using Clair3 Model: ${it}" }

minimap2(valid_model_ch, ref)
clair3(minimap2.out.sorted_alignment, ref, model_path)
minimap2(valid_model_ch)
clair3(minimap2.out.sorted_alignment, model_path)

clair3_filter_params = Channel.fromPath("${moduleDir}/process/clair3_filter_params.yml").first()
apply_filters_clair3(clair3.out.vcf, clair3_filter_params, "clair3_")

// Still use normal params with bcftools support vcf
filter_params = Channel.fromPath("${moduleDir}/process/filter_params.yml").first()
call_snps(minimap2.out.sorted_alignment, ref)
call_snps(minimap2.out.sorted_alignment)
apply_filters(call_snps.out.gvcf, filter_params, "bcftools_")


consensus_params = Channel.fromPath("${moduleDir}/process/clair3_consensus_params.yml").first()
calls = apply_filters.out.filtered_gvcf.join(apply_filters_clair3.out.filtered_gvcf)
make_clair3_consensus(calls, ref, consensus_params)
.join(ref_path_ch)
calls.view { "Calls channel: ${it}" }
make_clair3_consensus(calls, consensus_params)

emit:
alignment = rundial_with_bcftools.out.alignment.concat(minimap2.out.sorted_alignment)
gvcf = rundial_with_bcftools.out.gvcf.concat(clair3.out.vcf)
final_fasta = rundial_with_bcftools.out.final_fasta.concat(make_clair3_consensus.out.final_fasta)
variable_length_fasta = rundial_with_bcftools.out.variable_length_fasta.concat(make_clair3_consensus.out.variable_length_fasta)
final_vcf = rundial_with_bcftools.out.final_vcf.concat(make_clair3_consensus.out.final_vcf)
full_vcf = rundial_with_bcftools.out.full_vcf.concat(make_clair3_consensus.out.full_vcf)
creation_report_json = rundial_with_bcftools.out.creation_report_json.concat(make_clair3_consensus.out.report_json)
alignment = rundial_with_bcftools.out.alignment.concat(minimap2.out.sorted_alignment).join(ref_name_ch).map { it -> tuple(it[0], it[1], it[3], null, it[2]) }
gvcf = rundial_with_bcftools.out.gvcf.concat(clair3.out.vcf).join(ref_name_ch)
final_fasta = rundial_with_bcftools.out.final_fasta.concat(make_clair3_consensus.out.final_fasta).join(ref_name_ch)
variable_length_fasta = rundial_with_bcftools.out.variable_length_fasta.concat(make_clair3_consensus.out.variable_length_fasta).join(ref_name_ch)
final_vcf = rundial_with_bcftools.out.final_vcf.concat(make_clair3_consensus.out.final_vcf).join(ref_name_ch)
full_vcf = rundial_with_bcftools.out.full_vcf.concat(make_clair3_consensus.out.full_vcf).join(ref_name_ch)
creation_report_json = rundial_with_bcftools.out.creation_report_json.concat(make_clair3_consensus.out.report_json).join(ref_name_ch)
}


workflow rundial_with_bcftools {
take:
fastq_files
ref

main:
filter_params = Channel.fromPath("${moduleDir}/process/filter_params.yml").first()
consensus_params = Channel.fromPath("${moduleDir}/process/consensus_params.yml").first()

minimap2(fastq_files, ref)
call_all(minimap2.out.sorted_alignment, ref)
minimap2(fastq_files)
call_all(minimap2.out.sorted_alignment)
apply_filters(call_all.out.gvcf, filter_params, "")

calls = apply_filters.out.filtered_gvcf
make_consensus(calls, ref, consensus_params)
calls = calls.join(fastq_files).map { it -> tuple(it[0], it[1], it[3]) }
calls.view { "Calls channel: ${it}" }
make_consensus(calls, consensus_params)

emit:
alignment = minimap2.out.sorted_alignment
Expand Down
6 changes: 5 additions & 1 deletion nextflow.config
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,13 @@ params {
help = ''
input_dir = ''
workflow = 'clair3'
ref_fasta = 'data/h37rv_20231215.fa.gz'
ref_fasta = 'data/NC_000962.3.fa.gz'
clair3_models_dir = 'data/clair3_models'
basecalling_model = '' // e.g. "dna_r9.4.1_450bps_sup_prom"
reference_genome_suffix = '.fa.gz'
reference_genomes_dir = 'test_data/ref_data/'
accession = 'NC_000962.3'
reference = 'Mycobacterium tuberculosis'
}

profiles {
Expand Down
6 changes: 2 additions & 4 deletions process/consensus.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,7 @@ process make_consensus {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path(filtered_gvcf)
path reference
tuple val(sample_name), path(filtered_gvcf), path(reference)
path consensus_params

output:
Expand Down Expand Up @@ -76,8 +75,7 @@ process make_clair3_consensus {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path("filtered.gvcf.gz"), path("clair3.vcf.gz")
path reference
tuple val(sample_name), path("filtered.gvcf.gz"), path("clair3.vcf.gz"), path(reference)
path consensus_params

output:
Expand Down
22 changes: 10 additions & 12 deletions process/variant_calling.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,10 @@ process minimap2 {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path(fq)
path reference
tuple val(sample_name), path(fq), path(reference)

output:
tuple val(sample_name), path("final.bam"), emit: sorted_alignment
tuple val(sample_name), path("final.bam"), path(reference), emit: sorted_alignment

script:
"""
Expand All @@ -43,8 +42,7 @@ process call_snps {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path(sorted_alignment)
path reference
tuple val(sample_name), path(sorted_alignment), path(reference)

output:
tuple val(sample_name), path("calls.gvcf.gz"), emit: gvcf
Expand Down Expand Up @@ -89,8 +87,7 @@ process call_all {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path(sorted_alignment)
path reference
tuple val(sample_name), path(sorted_alignment), path(reference)

output:
tuple val(sample_name), path("calls.gvcf.gz"), emit: gvcf
Expand All @@ -99,8 +96,10 @@ process call_all {
"""
echo "Running bcftools mpileup in parallel"

gunzip -c ${reference} > ref.fasta # Handle gzipped (but not bgzipped) reference input

# script also makes indexes
multithreaded_bcftools -a ${sorted_alignment} -r ${reference} \
multithreaded_bcftools -a ${sorted_alignment} -r ref.fasta \
-o pileup.bcf -t ${task.cpus} \
--settings "-x -Q 10 -a INFO/SCR,INFO/ADR,INFO/ADF,FORMAT/SP,FORMAT/AD -h100 -M10000"

Expand All @@ -109,10 +108,10 @@ process call_all {
# Report all alleles during calling
bcftools call --threads ${task.cpus} --ploidy 1 -m -A -V indels pileup.bcf \
| bcftools filter -e 'DP==0' \
| bcftools norm -d exact -f ${reference} -Oz -o snps.gvcf.gz
| bcftools norm -d exact -f ref.fasta -Oz -o snps.gvcf.gz
bcftools call --threads ${task.cpus} --ploidy 1 -m -A -v -V snps pileup.bcf \
| bcftools filter -e 'DP==0' \
| bcftools norm -d exact -f ${reference} -Oz -o indels.gvcf.gz
| bcftools norm -d exact -f ref.fasta -Oz -o indels.gvcf.gz
bcftools index snps.gvcf.gz
bcftools index indels.gvcf.gz
bcftools concat -a snps.gvcf.gz indels.gvcf.gz -o calls.gvcf
Expand Down Expand Up @@ -166,8 +165,7 @@ process clair3 {
pod label: "run_id", value: "${params.run_id}"

input:
tuple val(sample_name), path(bam)
path reference
tuple val(sample_name), path(bam), path(reference)
path model

output:
Expand Down
4 changes: 2 additions & 2 deletions test_data/clair3_consensus/clair3.vcf
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
##fileformat=VCFv4.2
##source=Clair3
##clair3_version=1.0.4
##cmdline=/home/ubuntu/miniconda3/envs/clair3/bin/run_clair3.sh --bam_fn=/home/ubuntu/pipelines/sundial/MHall_comparison/assembly/ont/4.sorted.bam --ref_fn=/home/ubuntu/pipelines/sundial/data/h37rv_20231215.fa.gz --threads=8 --platform=ont --model_path=/home/ubuntu/miniconda3/envs/clair3/bin/models/r941_prom_sup_g5014 --include_all_ctgs --print_ref_calls --gvcf --no_phasing_for_fa --output=clair_out/4
##reference=/home/ubuntu/pipelines/sundial/data/h37rv_20231215.fa.gz
##cmdline=/home/ubuntu/miniconda3/envs/clair3/bin/run_clair3.sh --bam_fn=/home/ubuntu/pipelines/sundial/MHall_comparison/assembly/ont/4.sorted.bam --ref_fn=/home/ubuntu/pipelines/sundial/data/NC_000962.3.fa.gz --threads=8 --platform=ont --model_path=/home/ubuntu/miniconda3/envs/clair3/bin/models/r941_prom_sup_g5014 --include_all_ctgs --print_ref_calls --gvcf --no_phasing_for_fa --output=clair_out/4
##reference=/home/ubuntu/pipelines/sundial/data/NC_000962.3.fa.gz
##FILTER=<ID=PASS,Description="All filters passed">
##FILTER=<ID=LowQual,Description="Low quality variant">
##FILTER=<ID=RefCall,Description="Reference call">
Expand Down
9 changes: 3 additions & 6 deletions tests/nextflow/process_call_all.nf.test
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,7 @@ nextflow_process {
}
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand All @@ -30,8 +29,7 @@ nextflow_process {
}
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand All @@ -48,8 +46,7 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/not_exist.fastq.gz")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/not_exist.fastq.gz"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand Down
9 changes: 3 additions & 6 deletions tests/nextflow/process_call_snps.nf.test
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,7 @@ nextflow_process {
}
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand All @@ -30,8 +29,7 @@ nextflow_process {
}
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand All @@ -48,8 +46,7 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/not_exist.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[0] = ['sample', file("$projectDir/test_data/assemblers/not_exist.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
"""
}
}
Expand Down
10 changes: 4 additions & 6 deletions tests/nextflow/process_clair3.nf.test
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,8 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[2] = "$projectDir/data/clair3_models/r1041_e82_400bps_sup_v430.tar.gz"
input[0] = ['sample', file("$projectDir/test_data/assemblers/final.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
input[1] = "$projectDir/data/clair3_models/r1041_e82_400bps_sup_v430.tar.gz"
"""
}
}
Expand All @@ -25,9 +24,8 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/assemblers/non_exist.bam")]
input[1] = file("$projectDir/data/h37rv_20231215.fa.gz")
input[2] = "$projectDir/data/clair3_models/r1041_e82_400bps_sup_v430.tar.gz"
input[0] = ['sample', file("$projectDir/test_data/assemblers/non_exist.bam"), file("$projectDir/data/NC_000962.3.fa.gz")]
input[1] = "$projectDir/data/clair3_models/r1041_e82_400bps_sup_v430.tar.gz"
"""
}
}
Expand Down
10 changes: 4 additions & 6 deletions tests/nextflow/process_clair3_consensus.nf.test
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,8 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/clair3_consensus/bcftools.vcf"), file("$projectDir/test_data/clair3_consensus/clair3.vcf")]
input[1] = file("$projectDir/test_data/clair3_consensus/ref.fasta.gz")
input[2] = file("$projectDir/test_data/clair3_consensus/consensus_params.yml")
input[0] = ['sample', file("$projectDir/test_data/clair3_consensus/bcftools.vcf"), file("$projectDir/test_data/clair3_consensus/clair3.vcf"), file("$projectDir/test_data/clair3_consensus/ref.fasta.gz")]
input[1] = file("$projectDir/test_data/clair3_consensus/consensus_params.yml")
"""
}
}
Expand All @@ -31,9 +30,8 @@ nextflow_process {
when {
process{
"""
input[0] = ['sample', file("$projectDir/test_data/non_existent.vcf")]
input[1] = file("$projectDir/test_data/test_ref.fasta.gz")
input[2] = file("$projectDir/test_data/clair3_consensus/consensus_params.yml")
input[0] = ['sample', file("$projectDir/test_data/non_existent.vcf"), file("$projectDir/test_data/test_ref.fasta.gz")]
input[1] = file("$projectDir/test_data/clair3_consensus/consensus_params.yml")
"""
}
}
Expand Down
Loading
Loading