diff --git a/assets/magma_samplesheet_test.csv b/assets/magma_samplesheet_test.csv
new file mode 100644
index 0000000..7266c86
--- /dev/null
+++ b/assets/magma_samplesheet_test.csv
@@ -0,0 +1,4 @@
+Study,Sample,Library,Attempt,R1,R2,Flowcell,Lane,Index Sequence
+MAGMA,SRR26331590,1,1,https://zenodo.org/records/20671479/files/SRR26331590_R1.fastq.gz,https://zenodo.org/records/20671479/files/SRR26331590_R2.fastq.gz,UNNAMED,1,ATGCATGC
+MAGMA,SRR26331595,1,1,https://zenodo.org/records/20671479/files/SRR26331595_R1.fastq.gz,https://zenodo.org/records/20671479/files/SRR26331595_R2.fastq.gz,UNNAMED,1,ATGCATGC
+MAGMA,SRR26331599,1,1,https://zenodo.org/records/20671479/files/SRR26331599_R1.fastq.gz,https://zenodo.org/records/20671479/files/SRR26331599_R2.fastq.gz,UNNAMED,1,ATGCATGC
diff --git a/assets/schema_input.json b/assets/schema_input.json
index 1d6befa..3d4db8d 100644
--- a/assets/schema_input.json
+++ b/assets/schema_input.json
@@ -34,6 +34,9 @@
"errorMessage": "FastQ file for reads 2 cannot contain spaces and must have extension '.fq.gz' or '.fastq.gz'"
}
},
- "required": ["sample", "fastq_1"]
+ "anyOf": [
+ { "required": ["sample", "fastq_1"] },
+ { "required": ["Sample", "R1"] }
+ ]
}
}
diff --git a/bin/convert_ismapper_to_vcf.py b/bin/convert_ismapper_to_vcf.py
new file mode 100755
index 0000000..7883246
--- /dev/null
+++ b/bin/convert_ismapper_to_vcf.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import argparse
+import csv
+import glob
+from Bio import SeqIO
+
+# Define the VCF header
+vcf_header = """##fileformat=VCFv4.2
+##source=ISMapper
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##INFO=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMPLE
+"""
+
+# Function to read the reference genome
+def read_reference_genome(reference_file):
+ reference_sequences = {}
+ reference_sequence_length = 0
+ with open(reference_file, 'r') as ref_file:
+ # Read the reference genome using Biopython
+ for record in SeqIO.parse(ref_file, "fasta"):
+ reference_sequences[record.id] = record.seq
+ reference_sequence_length = len(record.seq)
+ return reference_sequences, reference_sequence_length
+
+# Function to read the transposable element information
+def read_transposable_elements(query_file):
+ te_info = {}
+ with open(query_file, 'r') as query_file:
+ # Read the multifasta using Biopython
+ for record in SeqIO.parse(query_file, "fasta"):
+ te_name = record.id
+ te_info[te_name] = len(record.seq)
+
+ # NOTE: Renable in case we need to dump this on disk.
+ # with open("temp_transposable_elements.csv", mode='w', newline='') as file:
+ # writer = csv.writer(file)
+ # writer.writerow(["name", "length"])
+ # for te_name, te_length in te_info.items():
+ # writer.writerow([te_name, te_length])
+ return te_info
+
+# Function to convert ISMapper data to VCF format
+def convert_is_mapper_to_vcf(is_mapper_dir, vcf_file, reference_sequences, te_info):
+
+ is_mapper_txt_file = glob.glob(is_mapper_dir + "*.txt")[0]
+
+ with open(is_mapper_txt_file, 'r') as infile_is_mapper, open(vcf_file, 'w') as outfile:
+ # Write the VCF header
+ outfile.write(vcf_header)
+
+ # Read the ISMapper TSV file
+ reader_is_mapper = csv.DictReader(infile_is_mapper, delimiter='\t')
+
+ for row in reader_is_mapper:
+ chrom = "NC-000962-3-H37Rv" # Using the specified chromosome
+ pos = int(row['x']) # Convert position to integer
+ region_id = row['region']
+ ref = reference_sequences[chrom][pos - 1] # Extract the reference allele from the reference sequence
+ orientation = row['orientation']
+ # FIXME Hard-code the name of this specific element
+ te_name = 'IS6110'
+ te_length = te_info.get(te_name, 'NA')
+ alt = f"{ref}[<{te_name},{orientation}>:{te_length}[" # Use transposable element, orientation, and its length
+ qual = '.'
+ filter_status = 'PASS'
+ info = (
+ f"SVTYPE=BND;"
+ f"Orientation={orientation};"
+ f"Gap={row['gap']};"
+ f"Call={row['call']};"
+ f"Percent_ID={row['percent_ID']};"
+ f"Percent_cov={row['percent_cov']};"
+ f"Left_gene={row['left_gene']};"
+ f"Left_description={row['left_description']};"
+ f"Left_strand={row['left_strand']};"
+ f"Left_distance={row['left_distance']};"
+ f"Right_gene={row['right_gene']};"
+ f"Right_description={row['right_description']};"
+ f"Right_strand={row['right_strand']};"
+ f"Right_distance={row['right_distance']};"
+ f"Gene_interruption={row['gene_interruption']}"
+ )
+ format_field = "GT:AD:DP:GQ:PL"
+ sample_field = "1:10,10:10:99:1800"
+
+ # Write the VCF entry
+ vcf_entry = f"{chrom}\t{pos}\t{region_id}\t{ref}\t{alt}\t{qual}\t{filter_status}\t{info}\t{format_field}\t{sample_field}\n"
+ outfile.write(vcf_entry)
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Convert ISMapper output to VCF format.')
+ parser.add_argument('--is_mapper_dir', required=True, help='Path to the ISMapper output directory.')
+ parser.add_argument('--reference_file', required=True, help='Path to the reference genome file.')
+ parser.add_argument('--query_file', required=True, help='Path to the transposable elements multifasta file.')
+ parser.add_argument('--output_vcf_file', required=True, help='Path to the output VCF file.')
+
+ # Step 3: Parse the command-line arguments
+ args = parser.parse_args()
+
+ # Step 4: Replace hardcoded file paths with variables
+ is_mapper_dir = args.is_mapper_dir
+ vcf_file = args.output_vcf_file
+ reference_file = args.reference_file
+ query_file = args.query_file
+
+ # Read the reference genome sequence
+ reference_sequences, _ = read_reference_genome(reference_file)
+
+ # Read the transposable element information
+ query_info = read_transposable_elements(query_file)
+
+ # Run the conversion function
+ convert_is_mapper_to_vcf(is_mapper_dir, vcf_file, reference_sequences, query_info)
diff --git a/bin/fastq_cohort_validation.py b/bin/fastq_cohort_validation.py
new file mode 100755
index 0000000..37eeccb
--- /dev/null
+++ b/bin/fastq_cohort_validation.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import glob
+import argparse
+import csv
+import json
+
+import pandas as pd
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Summarize the input FASTQ validation report')
+
+ parser.add_argument('magma_samplesheet', metavar='magma_samplesheet', type=str, help='')
+
+ parser.add_argument('merged_fastq_reports', metavar='merged_fastq_reports', type=str, help='')
+
+ parser.add_argument('magma_analysis', metavar='magma_analysis', type=str, help='')
+
+
+ args = vars(parser.parse_args())
+
+ # ============================================
+ # Parse the validation reports for exact sample names which passed/failed
+ # ============================================
+
+ # Load the JSON file into a dictionary
+ with open(args['merged_fastq_reports'], 'r') as f:
+ fastq_report_dict = json.load(f)
+
+ with open(args['magma_samplesheet'], 'r') as f:
+ magma_samplesheet_json = json.load(f)
+
+ fastq_report_keys_list = list(fastq_report_dict.keys())
+
+ magma_analysis_json = {}
+
+ for elem in magma_samplesheet_json:
+ elem["fastq_report"] = {}
+
+ if elem['R1'] is not None:
+ fastq_1_name = elem['R1'].split("/")[-1]
+ elem["fastqs_approved"] = True
+ if fastq_1_name in fastq_report_keys_list:
+ elem["fastq_report"][fastq_1_name] = {"file": fastq_report_dict[fastq_1_name]}
+ else:
+ elem["fastq_report"][fastq_1_name] = {"fastq_utils_check": "failed"}
+ elem["fastqs_approved"] = False
+
+ if elem['R2'] != "" :
+ fastq_2_name = elem['R2'].split("/")[-1]
+ if fastq_2_name in fastq_report_keys_list:
+ elem["fastq_report"][fastq_2_name] = {"file": fastq_report_dict[fastq_2_name]}
+ else:
+ elem["fastq_report"][fastq_2_name] = {"fastq_utils_check": "failed"}
+ elem["fastqs_approved"] = False
+
+ magma_analysis_json[elem["magma_sample_name"]] = elem
+
+ with open(args['magma_analysis'], 'w') as f:
+ json.dump(magma_analysis_json, f, indent=4, ensure_ascii= False)
diff --git a/bin/generate_merged_cohort_stats.py b/bin/generate_merged_cohort_stats.py
new file mode 100755
index 0000000..83e3c44
--- /dev/null
+++ b/bin/generate_merged_cohort_stats.py
@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import argparse
+import pandas as pd
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Summarize the CALL_WF and MINOR_VARIANTS_ANALYSIS_WF analysis results')
+ parser.add_argument('--relabundance_approved_tsv', default="approved_samples.relabundance.tsv", metavar='relabundance_approved_tsv', type=str, help='File enlisting the approved samples from MINOR_VARIANTS_ANALYSIS_WF')
+ parser.add_argument('--relabundance_rejected_tsv', default="rejected_samples.relabundance.tsv", metavar='relabundance_rejected_tsv', type=str, help='File enlisting the rejected samples from MINOR_VARIANTS_ANALYSIS_WF')
+ parser.add_argument('--call_wf_cohort_stats_tsv', default="joint.cohort_stats.tsv", metavar='call_wf_cohort_stats_tsv', type=str, help='File enlisting the cohort results of CALL_WF')
+ parser.add_argument('--output_file', default="joint.merged_cohort_stats.tsv", metavar='output_file', type=str, help='Name of the output file merged cohort statistics')
+ args = vars(parser.parse_args())
+
+ # Read the TSV files into dataframes
+ df_cohort_stats = pd.read_csv(args['call_wf_cohort_stats_tsv'], sep="\t")
+ df_cohort_stats.columns = df_cohort_stats.columns.str.strip()
+ df_cohort_stats['SAMPLE'] = df_cohort_stats['SAMPLE'].str.strip()
+ df_cohort_stats = df_cohort_stats.set_index('SAMPLE')
+
+ df_approved_relabundance_stats = pd.read_csv(args['relabundance_approved_tsv'], sep="\t")
+ df_approved_relabundance_stats.columns = df_approved_relabundance_stats.columns.str.strip()
+ df_approved_relabundance_stats['SAMPLE'] = df_approved_relabundance_stats['SAMPLE'].str.strip()
+ df_approved_relabundance_stats = df_approved_relabundance_stats.set_index('SAMPLE')
+
+ df_rejected_relabundance_stats = pd.read_csv(args['relabundance_rejected_tsv'], sep="\t")
+ df_rejected_relabundance_stats.columns = df_rejected_relabundance_stats.columns.str.strip()
+ df_rejected_relabundance_stats['SAMPLE'] = df_rejected_relabundance_stats['SAMPLE'].str.strip()
+ df_rejected_relabundance_stats = df_rejected_relabundance_stats.set_index('SAMPLE')
+
+ # Join the datasets
+ df_relabundance_stats_concat = pd.concat([df_approved_relabundance_stats, df_rejected_relabundance_stats])
+ df_joint_cohort_stats = df_cohort_stats.join(df_relabundance_stats_concat, how="outer")
+
+ # Reorder the columns
+ df_joint_cohort_stats.columns = df_joint_cohort_stats.columns.str.strip()
+ new_cols = ['AVG_INSERT_SIZE', 'MAPPED_PERCENTAGE', 'RAW_TOTAL_SEQS', 'AVERAGE_BASE_QUALITY', 'MEAN_COVERAGE', 'SD_COVERAGE', 'MEDIAN_COVERAGE', 'MAD_COVERAGE', 'PCT_EXC_ADAPTER', 'PCT_EXC_MAPQ', 'PCT_EXC_DUPE', 'PCT_EXC_UNPAIRED', 'PCT_EXC_BASEQ', 'PCT_EXC_OVERLAP', 'PCT_EXC_CAPPED', 'PCT_EXC_TOTAL', 'PCT_1X', 'PCT_5X', 'PCT_10X', 'PCT_30X', 'PCT_50X', 'PCT_100X', 'LINEAGES', 'FREQUENCIES', 'MAPPED_NTM_FRACTION_16S', 'MAPPED_NTM_FRACTION_16S_THRESHOLD_MET', 'COVERAGE_THRESHOLD_MET', 'BREADTH_OF_COVERAGE_THRESHOLD_MET', 'RELABUNDANCE_THRESHOLD_MET', 'ALL_THRESHOLDS_MET']
+ df_final_cohort_stats = df_joint_cohort_stats[new_cols]
+
+ # Impute the NaN value after join
+ df_final_cohort_stats['RELABUNDANCE_THRESHOLD_MET'] = df_final_cohort_stats['RELABUNDANCE_THRESHOLD_MET'].fillna(0)
+
+ # Prepare for boolean operation
+ df_final_cohort_stats['MAPPED_NTM_FRACTION_16S_THRESHOLD_MET'] = df_final_cohort_stats['MAPPED_NTM_FRACTION_16S_THRESHOLD_MET'].fillna(0).astype('Int64')
+ df_final_cohort_stats['COVERAGE_THRESHOLD_MET'] = df_final_cohort_stats['COVERAGE_THRESHOLD_MET'].fillna(0).astype('Int64')
+ df_final_cohort_stats['BREADTH_OF_COVERAGE_THRESHOLD_MET'] = df_final_cohort_stats['BREADTH_OF_COVERAGE_THRESHOLD_MET'].fillna(0).astype('Int64')
+ df_final_cohort_stats['RELABUNDANCE_THRESHOLD_MET'] = df_final_cohort_stats['RELABUNDANCE_THRESHOLD_MET'].fillna(0).astype('Int64')
+
+ # Derive the final threshold using Boolean operations
+ df_final_cohort_stats['ALL_THRESHOLDS_MET'] = (
+ df_final_cohort_stats['MAPPED_NTM_FRACTION_16S_THRESHOLD_MET'].apply(lambda x: bool(x) if pd.notna(x) else False) &
+ df_final_cohort_stats['COVERAGE_THRESHOLD_MET'].astype('bool') &
+ df_final_cohort_stats['BREADTH_OF_COVERAGE_THRESHOLD_MET'].astype('bool') &
+ df_final_cohort_stats['RELABUNDANCE_THRESHOLD_MET'].astype('bool')
+ )
+ df_final_cohort_stats['ALL_THRESHOLDS_MET'] = df_final_cohort_stats['ALL_THRESHOLDS_MET'].replace({True: 1, False: 0})
+
+ # Write the final dataframe to file
+ df_final_cohort_stats.to_csv(args['output_file'], sep="\t")
\ No newline at end of file
diff --git a/bin/multiple_infection_filter.py b/bin/multiple_infection_filter.py
new file mode 100755
index 0000000..ea3b08b
--- /dev/null
+++ b/bin/multiple_infection_filter.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import os
+import json
+import csv
+import glob
+
+import argparse
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Analyse resistance output from MAGMA pipeline')
+ parser.add_argument('indir', metavar='indir', type=str, help='The directory containing the LoFreq TBProfiler output')
+ parser.add_argument('relative_abundance_threshold', metavar='relative_abundance_threshold', type=float, help='Minimum relative abundance of the majority strain required to process the sample')
+
+ args = vars(parser.parse_args())
+ samples = []
+ json_files = glob.glob(args['indir'] + "/" + "*.json")
+ for file_name in json_files:
+ with open(file_name) as json_file:
+ samples.append(json.load(json_file))
+ accepted = [['SAMPLE', 'LINEAGES', 'FREQUENCIES', 'RELABUNDANCE_THRESHOLD_MET']]
+ rejected = [['SAMPLE', 'LINEAGES', 'FREQUENCIES', 'RELABUNDANCE_THRESHOLD_MET']]
+ for sample in samples:
+ sublins = [i for i in sample['lineage'] if i['lineage'] in sample['sub_lineage'].split(';')]
+ lins = [i['lineage'] for i in sublins]
+ fracs = [i['fraction'] for i in sublins]
+ if not lins:
+ rejected.append([sample['id'], 'None', 'None', 0])
+ elif max(fracs) < args['relative_abundance_threshold']:
+ rejected.append([sample['id'], ';'.join(lins), ';'.join(['{:.0%}'.format(i) for i in fracs]), 0])
+ else:
+ accepted.append([sample['id'], ';'.join(lins), ';'.join(['{:.0%}'.format(i) for i in fracs]), 1])
+
+ with open('approved_samples.relabundance.tsv', 'w') as csvfile:
+ csvwriter = csv.writer(csvfile, delimiter='\t', quotechar='"', quoting=csv.QUOTE_MINIMAL)
+ csvwriter.writerows(accepted)
+
+ with open('rejected_samples.relabundance.tsv', 'w') as csvfile:
+ csvwriter = csv.writer(csvfile, delimiter='\t', quotechar='"', quoting=csv.QUOTE_MINIMAL)
+ csvwriter.writerows(rejected)
diff --git a/bin/preprocess_multiqc_input.py b/bin/preprocess_multiqc_input.py
new file mode 100755
index 0000000..0bbc4d4
--- /dev/null
+++ b/bin/preprocess_multiqc_input.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+
+def prepare_cohort_stats(merged_cohort_stats_file,output_name):
+
+ cohort_stats_output = []
+ with open(merged_cohort_stats_file) as tsv_file:
+ tsv_reader = csv.reader(tsv_file, delimiter='\t')
+ cohort_stats_output.append("\t".join(next(tsv_reader, None)))
+ for row in tsv_reader:
+ # replace colums 25 - 31
+ for column in range(26,31):
+ # replace 1 by APPROVED in "*_MET" and 0 by REJECTED
+ row[column] = row[column].replace("1","APPROVED").replace("0","REJECTED")
+ cohort_stats_output.append("\t".join(row))
+ with open(output_name, 'w') as f:
+ f.write("\n".join(cohort_stats_output))
+
+
+def prepare_distance_matrix(matrix_file, output_file):
+ samples_to_remove =["MTb_L1.SAMN10185847",
+ "MTb_L2.SAMEA1877219",
+ "MTb_L3.SAMEA1877181",
+ "MTb_L410.ERR216945",
+ "MTb_L43.ERR1193883",
+ "MTb_L5.SAMEA1877169",
+ "MTb_L6.SAMEA1877150",
+ "MTb_L7.ERR1971849",
+ "MTb_L8.SRR10828835",
+ "MTb_L9.ERR4162024",
+ "MTB_L10.ERR2707158",
+ "Mcanettii.ERR5104570"]
+
+ with open(matrix_file, "r") as f:
+ lines = [line.strip().split("\t") for line in f.readlines()]
+ header = lines[0]
+ filtered_indices = [i for i, sample in enumerate(header) if sample not in samples_to_remove]
+ filtered_matrix = [[""]+[header[i] for i in filtered_indices]]
+ for line in lines[1::]:
+ sample_name = line[0]
+ if sample_name not in samples_to_remove:
+ filtered_matrix.append([line[0]]+[line[i+1] for i in filtered_indices])
+
+ with open(output_file, "w", newline="") as f:
+ writer = csv.writer(f, delimiter="\t")
+ writer.writerows(filtered_matrix)
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Prepare TSV files for multiqc report generation.')
+ parser.add_argument('--merged_cohort_stats', required=True, help='Path to the merged cohort status TSV file.')
+
+ parser.add_argument('--output_prefix',required=False,default="prepared", type=str)
+ parser.add_argument('--skip_merge_analysis', action="store_false", help='If used, specify that merge analysis skipped.')
+ parser.add_argument('--distance_matrix', required=False, help='Path to the distance matrix TSV file')
+
+
+
+ args = parser.parse_args()
+ prepare_cohort_output_name = args.output_prefix + "_" + args.merged_cohort_stats
+ prepare_cohort_stats(args.merged_cohort_stats,prepare_cohort_output_name)
+
+ if args.skip_merge_analysis:
+ print("merge was not skipped")
+ prepare_distance_matrix_output_name = args.output_prefix + "_" + args.distance_matrix
+ prepare_distance_matrix(args.distance_matrix, prepare_distance_matrix_output_name)
+
+
diff --git a/bin/reduce_annotations.py b/bin/reduce_annotations.py
new file mode 100755
index 0000000..2d89978
--- /dev/null
+++ b/bin/reduce_annotations.py
@@ -0,0 +1,209 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import os
+import sys
+import argparse
+import re
+import json
+
+
+#####################################################
+# Arg processing
+#####################################################
+
+
+def validate_file(f):
+ if not os.path.exists(f):
+ raise argparse.ArgumentTypeError(
+ "\n\nERROR: The file {0} does not exist.".format(f)
+ )
+ return f
+
+
+def parse_args(args=None):
+ parser = argparse.ArgumentParser(
+ description="Extract and reshape the annotations for variant recalibration",
+ epilog="Example usage: reduce_annotations.py -i annotations.log -o new.annotations.txt",
+ )
+
+ parser.add_argument(
+ "-i",
+ "--input",
+ dest="input_annotations_file",
+ required=True,
+ type=validate_file,
+ help="input file",
+ metavar="FILE",
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output",
+ dest="output_annotations_file",
+ required=True,
+ help="output file",
+ metavar="FILE",
+ )
+
+ parser.add_argument(
+ "-t",
+ "--tranches",
+ dest="input_tranches_file",
+ required=True,
+ type=validate_file,
+ help="input tranches file",
+ metavar="FILE",
+ )
+
+ parser.add_argument(
+ "-d",
+ "--output_tranches_and_annotations",
+ dest="output_annotations_and_tranches_data",
+ required=True,
+ help="output annotations and tranches data",
+ metavar="FILE",
+ )
+
+ args = parser.parse_args()
+ return args
+
+
+#####################################################
+# Extract annotations
+#####################################################
+
+
+def extract_reduce_write_annotations(input_filename, output_filename):
+
+ annotations_order_marker = "VariantDataManager - Annotation order is"
+
+ with open(input_filename) as input_annotations_file:
+ lines = input_annotations_file.readlines()
+ annotations_order_string = list(
+ filter(lambda a_line: a_line.find(annotations_order_marker) != -1, lines)
+ )[0]
+ annotations_order_list = (
+ annotations_order_string.strip()
+ .split(":")[-1]
+ .replace("[", "")
+ .replace("]", "")
+ .split(",")
+ )
+ print(
+ "ANNOTATIONS_ANALYSIS: Initial ordered annotations list: ",
+ annotations_order_list,
+ )
+
+ eliminated_annotation_string = (
+ "-an " + " -an ".join(annotations_order_list[:-1]).strip()
+ )
+ print("ANNOTATIONS_ANALYSIS: Dropping : " + annotations_order_list[-1])
+ print(
+ "ANNOTATIONS_ANALYSIS: Final annotation string : ",
+ eliminated_annotation_string,
+ )
+
+ with open(output_filename, "w") as output:
+ output.write(eliminated_annotation_string)
+
+ return eliminated_annotation_string
+
+
+#####################################################
+# Tranches and annotations file
+#####################################################
+
+
+def tranches_and_annotations(
+ input_annotations_file, input_tranches_filename, output_tranches_file
+):
+
+ annotations_order_marker = "VariantDataManager - Annotation order is"
+
+ with open(input_annotations_file) as input_annotations_file:
+ lines = input_annotations_file.readlines()
+ annotations_order_string = list(
+ filter(lambda a_line: a_line.find(annotations_order_marker) != -1, lines)
+ )[0]
+ annotations_order_list = (
+ annotations_order_string.strip()
+ .split(":")[-1]
+ .replace("[", "")
+ .replace("]", "")
+ .split(",")
+ )
+
+ with open(input_tranches_filename) as input_tranches_file:
+ lines = input_tranches_file.readlines()
+ minVQSLod_score = lines[-2].split(",")[5]
+ print(
+ "TRANCHES_ANALYSIS: minVQSLod score (for initial ordered annotations list) corresponding to targetTruthSensitivity: 99.90 is => ",
+ minVQSLod_score,
+ )
+
+ with open(output_tranches_file, "w") as output:
+ annotations_order_joined_string = ",".join(annotations_order_list)
+
+ annotations_count = re.search(r"\w*ANN\d\w*", input_tranches_filename).group()
+
+ annotations_dict = {
+ "minVQSLod": minVQSLod_score,
+ "annotationsOrder": annotations_order_joined_string,
+ "annotationsCount": annotations_count,
+ }
+
+ json.dump(annotations_dict, output)
+
+
+#####################################################
+# Main
+#####################################################
+
+
+def main(args=None):
+ args = parse_args(args)
+
+ extract_reduce_write_annotations(
+ args.input_annotations_file, args.output_annotations_file
+ )
+
+ tranches_and_annotations(
+ args.input_annotations_file,
+ args.input_tranches_file,
+ args.output_annotations_and_tranches_data,
+ )
+
+
+#####################################################
+#####################################################
+
+if __name__ == "__main__":
+ sys.exit(main())
+
+#####################################################
+#####################################################
diff --git a/bin/reduce_strand_bias.py b/bin/reduce_strand_bias.py
new file mode 100755
index 0000000..c87f1c3
--- /dev/null
+++ b/bin/reduce_strand_bias.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+
+import argparse
+import vcf
+from scipy.stats import binom_test, binomtest
+
+import warnings
+warnings.filterwarnings("ignore", category=DeprecationWarning)
+
+def filter_vcf_file(vcf_input, vcf_output, pval):
+ # Open the VCF file for reading
+
+ vcf_reader = vcf.Reader(open(vcf_input, 'r'))
+
+ # Create a list to store the filtered records
+
+ filtered_records = []
+
+ # Iterate over each record (variant) in the VCF file
+
+ for record in vcf_reader:
+
+ # Extract the DP4 field from the INFO dictionary
+
+ dp4_values = record.INFO.get('DP4')
+
+ # Perform binomial test on the last two values
+
+ if dp4_values:
+
+ # The binom_test tests the null hypothesis that ALT variants are fairly distributed across the read directions, the hypothetical success probability of which is 0.5, p=0.5 reflects this hypothetical probability not the significance level.
+ p_value = binom_test([dp4_values[-2], dp4_values[-1]], n=sum(dp4_values[-2:]), p=0.5)
+
+ # Check if p-value is above or equal to 0.05
+ if p_value >= pval:
+ filtered_records.append(record)
+
+ # Write the filtered records to a new VCF file
+ vcf_writer = vcf.Writer(open(vcf_output, 'w'), vcf_reader)
+
+ for record in filtered_records:
+ vcf_writer.write_record(record)
+
+ # Close the VCF writer
+ vcf_writer.close()
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Reduce strand bias in VCF file')
+ parser.add_argument('pval', metavar='pval', type=str, help='The pval to be used for filtering')
+ parser.add_argument('input', metavar='input_vcf', type=str, help='The initial VCF file')
+ parser.add_argument('output', metavar='output_vcf', type=str, help='The output VCF filed')
+ args = vars(parser.parse_args())
+
+ input_pval = float(args['pval'])
+ input_vcf_file = args['input']
+ output_vcf_file = args['output']
+
+ filter_vcf_file(input_vcf_file, output_vcf_file, input_pval)
diff --git a/bin/reformat_lofreq.py b/bin/reformat_lofreq.py
new file mode 100755
index 0000000..475ec14
--- /dev/null
+++ b/bin/reformat_lofreq.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import ast
+import argparse
+
+import pandas as pd
+
+header_formats = """##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+##FORMAT=
+"""
+
+def read_vcf(filename):
+ with open(filename, 'r') as vcf:
+ header = ''
+ for line in vcf:
+ if line[:2] != '##':
+ break
+ else:
+ header += line
+ try:
+ df = pd.read_csv(vcf, header=None, sep='\t')
+ df.columns = line[:-1].split('\t')
+ not_empty = True
+ except pd.errors.EmptyDataError as e:
+ df = pd.DataFrame(columns=line[:-1].split('\t'))
+ not_empty = False
+ return df, header, not_empty
+
+def write_vcf(filename, df, header):
+ with open(filename, 'w') as vcf:
+ vcf.write(header)
+ vcf.write(header_formats)
+ df.to_csv(vcf, sep='\t', index=False)
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Analyse resistance output from the MAGMA pipeline')
+ parser.add_argument('lofreq_vcf_file', metavar='lofreq_vcf_file', type=str, help='The input lofreq vcf file')
+ parser.add_argument('lofreq_sample_name', metavar='lofreq_sample_name', type=str, help='The sample name')
+ parser.add_argument('outfile', metavar='outfile', type=str, help='The name of the output VCF file')
+
+ args = vars(parser.parse_args())
+
+ vcf, header, not_empty = read_vcf(args['lofreq_vcf_file'])
+ header = '\n'.join([i for i in header.split('\n') if 'lofreq' not in i])
+ if not_empty:
+ vcf['FORMAT'] = 'GT:AD:DP:GQ:PL'
+
+ for idx, row in vcf.iterrows():
+ info = [ast.literal_eval(i.split('=')[1]) for i in row['INFO'].split(';')[:4]]
+ ref_dp = sum(info[3][:2])
+ alt_dp = sum(info[3][2:])
+ GT = 1
+ AD = '{},{}'.format(ref_dp, alt_dp)
+ DP = sum(info[3])
+ GQ = 99
+ PL = '1800,0'
+ vcf.loc[idx, args['lofreq_sample_name']] = '{}:{}:{}:{}:{}'.format(GT,AD,DP,GQ,PL)
+
+ write_vcf(args['outfile'], vcf, header)
diff --git a/bin/rename_vcf_chrom.py b/bin/rename_vcf_chrom.py
new file mode 100755
index 0000000..c5c8e81
--- /dev/null
+++ b/bin/rename_vcf_chrom.py
@@ -0,0 +1,74 @@
+#! /usr/bin/env python
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+'''Original author Jody Phelan at https://github.com/jodyphelan/pathogen-profiler/blob/master/scripts/rename_vcf_chrom.py'''
+import sys
+import argparse
+import subprocess
+
+def errlog(x,ext=False):
+ sys.stderr.write('\033[91m' + str(x) + '\033[0m' + '\n')
+ if ext==True:
+ quit(1)
+
+def cmd_out(cmd,verbose=1):
+ cmd = "set -u pipefail; " + cmd
+ if verbose==2:
+ sys.stderr.write("\nRunning command:\n%s\n" % cmd)
+ stderr = open("/dev/stderr","w")
+ elif verbose==1:
+ sys.stderr.write("\nRunning command:\n%s\n" % cmd)
+ stderr = open("/dev/null","w")
+ else:
+ stderr = open("/dev/null","w")
+ try:
+ res = subprocess.Popen(cmd,shell=True,stderr = stderr,stdout=subprocess.PIPE)
+ for l in res.stdout:
+ yield l.decode().rstrip()
+ except:
+ errlog("Command Failed! Please Check!")
+ raise Exception
+ stderr.close()
+
+def main(args):
+ generator = cmd_out("bcftools view " + args.vcf) if args.vcf else sys.stdin
+ convert = dict(zip(args.source,args.target))
+ for l in generator:
+ if l[0]=="#":
+ sys.stdout.write(l.strip()+"\n")
+ else:
+ row = l.strip().split()
+ row[0] = convert[row[0]]
+ sys.stdout.write("\t".join(row)+"\n")
+
+parser = argparse.ArgumentParser(description='tbprofiler script',formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument('--vcf',type=str,help='')
+parser.add_argument('--source',nargs="+",type=str,help='')
+parser.add_argument('--target',nargs="+",type=str,help='')
+parser.set_defaults(func=main)
+args = parser.parse_args()
+args.func(args)
diff --git a/bin/sample_stats.py b/bin/sample_stats.py
new file mode 100755
index 0000000..9593621
--- /dev/null
+++ b/bin/sample_stats.py
@@ -0,0 +1,96 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import ast
+import argparse
+import re
+
+import pandas as pd
+
+re_mapped_p = re.compile(r'\d* mapped \((.*)%\)')
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Process the sample stats')
+ parser.add_argument('--sample_name', dest='sample_name', required=True, metavar='sample_name', type=str, help='The sample name')
+ parser.add_argument('--flagstat_file', dest='flagstat_file', required=True, metavar='flagstat_file', type=str, help='The flag stats file')
+ parser.add_argument('--samtoolsstats_file', dest='samtoolsstats_file', required=True, metavar='samtoolsstats_file', type=str, help='The samtools stats file')
+ parser.add_argument('--wgsmetrics_file', dest='wgsmetrics_file', required=True, metavar='wgsmetrics_file', type=str, help='The WGS metrics file')
+ parser.add_argument('--ntmfraction_file', dest='ntmfraction_file', required=True, metavar='ntmfraction_file', type=str, help='The NTM fraction file')
+
+ parser.add_argument('--cutoff_median_coverage', metavar='cutoff_median_coverage', default=10, type=float, help='The median coverage cutoff threshold')
+ parser.add_argument('--cutoff_breadth_of_coverage', metavar='cutoff_breadth_of_coverage', default=0.9, type=float, help='The breadth of coverage cutoff threshold')
+ parser.add_argument('--cutoff_ntm_fraction', metavar='cutoff_ntm_fraction', default=0.2, type=float, help='The NTM fraction cutoff threshold')
+
+## NOTE: This is computed by the multiple_infection_filter script
+# parser.add_argument('--cutoff_rel_abundance', metavar='cutoff_rel_abundance', default=0.8, type=float, help='The relative abundance cutoff threshold')
+
+ args = vars(parser.parse_args())
+
+ with open(args['wgsmetrics_file']) as f:
+ for line in f:
+ if '## METRICS CLASS' in line:
+ rows = [f.readline().strip(), f.readline().strip()]
+ wgsmetrics = pd.DataFrame([rows[1].split('\t')], columns=rows[0].split('\t'))
+ with open(args['ntmfraction_file']) as f:
+ ntm_fraction = float(f.read().strip())
+ with open(args['samtoolsstats_file']) as f:
+ for line in f:
+ if 'insert size average' in line:
+ ins_size = float(line.strip().split('\t')[2])
+ if 'raw total sequences' in line:
+ total_seqs = int(line.strip().split('\t')[2])
+ if 'average quality' in line:
+ avg_qual = float(line.strip().split('\t')[2])
+ with open(args['flagstat_file']) as f:
+ for line in f:
+ m = re_mapped_p.match(line)
+ if m:
+ mapped_p = float(m[1])
+
+ if int(wgsmetrics.loc[0, 'MEDIAN_COVERAGE']) >= args['cutoff_median_coverage']:
+ coverage_threshold_met = 1
+ else:
+ coverage_threshold_met = 0
+
+ if float(wgsmetrics.loc[0, 'PCT_1X']) >= args['cutoff_breadth_of_coverage']:
+ breadth_of_coverage_threshold_met = 1
+ else:
+ breadth_of_coverage_threshold_met = 0
+
+ if ntm_fraction <= args['cutoff_ntm_fraction']:
+ ntm_fraction_threshold_met = 1
+ else:
+ ntm_fraction_threshold_met = 0
+
+ if coverage_threshold_met and breadth_of_coverage_threshold_met and ntm_fraction_threshold_met:
+ all_thresholds_met = 1
+ else:
+ all_thresholds_met = 0
+
+ with open('{}.stats.tsv'.format(args['sample_name']), 'w') as f:
+ f.write('\t'.join([str(i) for i in [args['sample_name'], ins_size, mapped_p, total_seqs, avg_qual] + list(wgsmetrics.loc[0, ['MEAN_COVERAGE', 'SD_COVERAGE', 'MEDIAN_COVERAGE', 'MAD_COVERAGE', 'PCT_EXC_ADAPTER', 'PCT_EXC_MAPQ', 'PCT_EXC_DUPE', 'PCT_EXC_UNPAIRED', 'PCT_EXC_BASEQ', 'PCT_EXC_OVERLAP', 'PCT_EXC_CAPPED', 'PCT_EXC_TOTAL', 'PCT_1X', 'PCT_5X', 'PCT_10X', 'PCT_30X', 'PCT_50X', 'PCT_100X']]) + [ntm_fraction, ntm_fraction_threshold_met, coverage_threshold_met, breadth_of_coverage_threshold_met, all_thresholds_met]]))
+ f.write('\n')
diff --git a/bin/samplesheet_validation.py b/bin/samplesheet_validation.py
new file mode 100755
index 0000000..716f605
--- /dev/null
+++ b/bin/samplesheet_validation.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+
+import re
+import argparse
+import csv
+import json
+
+from sys import exit
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Run the MAGMA pipeline samplesheet validation')
+ parser.add_argument('input_file', metavar='input_file', type=str, help='The input sample file')
+ parser.add_argument('output_file', metavar='output_file', type=str, help='The validate output sample file')
+ args = vars(parser.parse_args())
+
+ name_re = re.compile('^[a-zA-Z0-9\-_]*$')
+
+ standard_magma_fields = ['Study', 'Sample', 'Library', 'Attempt', 'R1', 'R2', 'Flowcell', 'Lane', 'Index Sequence']
+
+ # Read the CSV file using the csv library
+ with open(args['input_file'], 'r') as infile, open(args['output_file'], 'w', newline='') as outfile:
+ reader = csv.DictReader(infile)
+
+ writer = csv.DictWriter(outfile,
+ fieldnames=standard_magma_fields + ['magma_sample_name', 'magma_bam_rg_string'],
+ quoting=csv.QUOTE_NONE, # Disable quoting
+ escapechar='\\') # Use backslash as escape character
+ writer.writeheader()
+
+ samplesheet_data = []
+ fail = False
+ for row in reader:
+
+ if 'Library' not in reader.fieldnames or row['Library'] == "" or row['Library'] == None:
+ row['Library'] = "1"
+
+ if 'Attempt' not in reader.fieldnames or row['Attempt'] == "" or row['Attempt'] == None:
+ row['Attempt'] = "1"
+
+ if 'Flowcell' not in reader.fieldnames or row['Flowcell'] == "" or row['Flowcell'] == None:
+ row['Flowcell'] = "1"
+
+ if 'Lane' not in reader.fieldnames or row['Lane'] == "" or row['Lane'] == None:
+ row['Lane'] = "1"
+
+ # if 'IndexSequence' not in reader.fieldnames or row['IndexSequence'] == "":
+ # row['IndexSequence'] = "1"
+
+ if 'Index Sequence' not in reader.fieldnames or row['Index Sequence'] == "" or row['Index Sequence'] == None:
+ row['Index Sequence'] = "1"
+
+ if 'Study' not in reader.fieldnames or row['Study'] == "":
+ row['Study'] = "MAGMA"
+
+ # Perform validation checks
+ if 'Study' in reader.fieldnames and not name_re.match(row['Study']):
+ print(f'Row {reader.line_num}: {row["Study"]} Illegal character in STUDY id')
+ fail = True
+ if not name_re.match(row['Sample']):
+ print(f'Row {reader.line_num}: {row["Sample"]} - Illegal character in SAMPLE id')
+ fail = True
+ if row['R1'] == row['R2']:
+ print(f'Row {reader.line_num}: {row["R1"]}, {row["R2"]} - DUPLICATED fastq file specified')
+ fail = True
+
+ # Create the magma_sample_name column
+ row[
+ 'magma_sample_name'] = f"{row['Study']}.{row['Sample']}.L{row['Library']}.A{row['Attempt']}.{row['Flowcell']}.{row['Lane']}.{row['Index Sequence']}"
+
+ # Create the magma_bwa_rg_string column
+ row[
+ 'magma_bam_rg_string'] = f"@RG\\tID:{row['Flowcell']}.{row['Lane']}\\tSM:{row['Study']}.{row['Sample']}\\tPL:illumina\\tLB:lib{row['Library']}\\tPU:{row['Flowcell']}.{row['Lane']}.{row['Index Sequence']}"
+
+ # Write the validated row to the output file
+ writer.writerow(row)
+
+ # Append the row to the samplesheet_data list
+ samplesheet_data.append(row)
+
+ # Convert the samplesheet_data list to JSON and write to a file
+ with open('magma_samplesheet.json', 'w') as f:
+ json.dump(samplesheet_data, f, indent=4, ensure_ascii=False)
+
+ if not fail:
+ print('Samplesheet format validation checks PASSED')
+ exit(0)
+ else:
+ print('Samplesheet format validation checks FAILED')
+ exit(1)
diff --git a/bin/select_best_annotations.py b/bin/select_best_annotations.py
new file mode 100755
index 0000000..3ce9a63
--- /dev/null
+++ b/bin/select_best_annotations.py
@@ -0,0 +1,193 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import os
+import sys
+import argparse
+import json
+from pathlib import Path
+
+
+#####################################################
+# Arg processing
+#####################################################
+
+
+def validate_path(f):
+ if not os.path.exists(f):
+ raise argparse.ArgumentTypeError(
+ "\n\nERROR: The path {0} does not exist.".format(f)
+ )
+ return f
+
+
+def parse_args(args=None):
+ parser = argparse.ArgumentParser(
+ description="Select the best set of annotations based on the minimal minVQSLod score",
+ epilog="Example usage: select_best_annotations.py --input_directory annotations_and_tranches_json_files --output_directory best_annotation_files --output_json_file_name all_annotations_data.json",
+ )
+
+ parser.add_argument(
+ "-i",
+ "--input_directory",
+ dest="input_folder_with_json_files",
+ required=True,
+ type=validate_path,
+ help="input directory",
+ metavar="FOLDER",
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output_directory",
+ dest="output_folder_with_best_annotations",
+ required=True,
+ type=validate_path,
+ help="output directory",
+ metavar="FOLDER",
+ )
+
+ parser.add_argument(
+ "-j",
+ "--output_json_file_name",
+ dest="output_json_file_name",
+ required=True,
+ help="output directory",
+ metavar="FOLDER",
+ )
+
+ args = parser.parse_args()
+ return args
+
+
+#####################################################
+# Find and move files related to the best annotations
+#####################################################
+
+
+def move_files(annotations_dict, output_directory):
+
+ best_annotations_count = annotations_dict["annotationsCount"]
+
+ p = Path(".")
+ grep_pattern = "*" + best_annotations_count + "*"
+
+ best_ann_tranches_files = list(p.glob("./tranches_files/" + grep_pattern))
+ print("BEST ANNOTATIONS TRANCHES FILE: ", best_ann_tranches_files)
+
+ for tranches_file in best_ann_tranches_files:
+ tranches_file.rename("./" + output_directory + "/" + tranches_file.name)
+
+ best_ann_recalibrated_vcf_files = list(p.glob("./recal_vcf_files/" + grep_pattern))
+ print("BEST ANNOTATIONS RECALIBRATED VCF FILES: ", best_ann_recalibrated_vcf_files)
+
+ for recal_vcf_file in best_ann_recalibrated_vcf_files:
+ recal_vcf_file.rename("./" + output_directory + "/" + recal_vcf_file.name)
+
+
+def read_annotations_json_file(annotation_data_json_file_path):
+
+ with open(annotation_data_json_file_path, "r") as annotation_json_data_file:
+ annotation_data_dict = json.load(annotation_json_data_file)
+
+ return annotation_data_dict
+
+
+def collect_annotations_data(base_folder, list_of_files):
+
+ annotations_dict_list = []
+
+ for annotation_json_file in list_of_files:
+ annotation_data_json_file_path = base_folder + "/" + annotation_json_file
+
+ annotations_dict_list.append(
+ read_annotations_json_file(annotation_data_json_file_path)
+ )
+
+ return annotations_dict_list
+
+
+def find_best_annotations(input_folder, json_file_name):
+
+ # NOTE: Explore the use of path.iterdir
+ list_of_annotations_files = os.listdir(input_folder)
+
+ all_annotations_dict_list = collect_annotations_data(
+ input_folder, list_of_annotations_files
+ )
+
+ print("ALL ANNOTATIONS DATA: ", all_annotations_dict_list)
+
+ with open(json_file_name, "w") as fp:
+ json.dump(all_annotations_dict_list, fp, indent=4)
+
+ tentative_best_annotations = all_annotations_dict_list[0]
+ tentative_max_minvqslod_score = float(tentative_best_annotations["minVQSLod"])
+
+ for annotations_dict in all_annotations_dict_list[1:]:
+
+ candidate_minvqslod = float(annotations_dict["minVQSLod"])
+
+ tentative_max_minvqslod_score = max(
+ tentative_max_minvqslod_score, candidate_minvqslod
+ )
+
+ best_annotations_dict = list(
+ filter(
+ lambda a_dict: float(a_dict["minVQSLod"]) == tentative_max_minvqslod_score,
+ all_annotations_dict_list,
+ )
+ )[0]
+
+ return best_annotations_dict
+
+
+#####################################################
+# Main
+#####################################################
+
+
+def main(args=None):
+ args = parse_args(args)
+
+ best_annotations_dict = find_best_annotations(
+ args.input_folder_with_json_files, args.output_json_file_name
+ )
+
+ print("BEST ANNOTATIONS: ", best_annotations_dict)
+
+ move_files(best_annotations_dict, args.output_folder_with_best_annotations)
+
+
+#####################################################
+#####################################################
+
+if __name__ == "__main__":
+ sys.exit(main())
+
+#####################################################
+#####################################################
diff --git a/bin/summarize_resistance.py b/bin/summarize_resistance.py
new file mode 100755
index 0000000..49d2810
--- /dev/null
+++ b/bin/summarize_resistance.py
@@ -0,0 +1,225 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+import glob
+import json
+import os
+import re
+import argparse
+
+import pandas as pd
+import numpy as np
+
+from tqdm import tqdm
+
+alphabet = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
+# who_v2_source = 'https://www.who.int/publications/i/item/9789240082410'
+
+drugs = {'amikacin', 'bedaquiline', 'capreomycin', 'clofazimine', 'cycloserine', 'delamanid', 'ethambutol',
+ 'ethionamide', 'imipenem', 'isoniazid', 'kanamycin', 'levofloxacin', 'linezolid', 'meropenem', 'moxifloxacin',
+ 'para_aminosalicylic_acid', 'pretomanid', 'prothionamide', 'pyrazinamide', 'rifabutin', 'rifampicin',
+ 'streptomycin', 'terizidone'}
+
+unknown_positions = 2.5
+
+map_resistance_class = {
+ 0: 'R',
+ 1: 'Very high probability',
+ 2: 'High probability',
+ 3: 'Moderate probability',
+ 4: 'Low probability',
+ 5: 'Very low probability',
+ 6: 'S',
+ -1: 'Unknown'
+}
+
+map_confidence_who = {
+ 'Assoc w R': 0,
+ 'Assoc w R - Interim': 1,
+ 'Uncertain significance': -1,
+ 'Not assoc w R - Interim': 5,
+ 'Not assoc w R': 6,
+}
+
+
+def create_resistance_df(sample_res, method):
+ pt_df = pd.DataFrame(columns=['Drug', 'Variant', 'Resistance interpretation','Source', 'Source notation', 'Type', 'Frequency', 'Literature']).set_index(['Drug', 'Variant'])
+
+ # DR variants from the MAGMA analysis. Added after LoFreq DR variants to give priority to MAGMA frequencies.
+ for var in sample_res['dr_variants']:
+ gene = var['gene_name']
+ var_repr = '{}_{}'.format(gene, var['change'])
+
+ # Add variants with an unlisted classification and overwrite later if necessary
+ for drug in var['gene_associated_drugs']:
+ drug_name = drug.lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [unknown_positions, 'non-Catalogue', '', var['type'], '{:.0%}'.format(var['freq']), 'Manually curated']
+
+ # Overwrite the variant classification for drugs which have a WHO classification last as to overrule unlisted classifications
+ for drug in var['drugs']:
+ drug_name = drug['drug'].lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [map_confidence_who[drug['confidence']], 'Catalogue', drug['original_mutation'], var['type'], '{:.0%}'.format(var['freq']), drug['source']]
+
+ # Other variants from MAGMA analysis.
+ for var in sample_res['other_variants']:
+ gene = var['gene_name']
+ var_repr = '{}_{}'.format(gene, var['change'])
+
+ # Add variants with an unlisted classification and overwrite later if necessary
+ for drug in var['gene_associated_drugs']:
+ drug_name = drug.lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [unknown_positions, 'non-Catalogue', '', var['type'], '{:.0%}'.format(var['freq']), 'Manually curated']
+
+ # Overwrite the variant classification for drugs which have a WHO classification last as to overrule unlisted classifications
+ if 'annotation' in var:
+ for annotation in var['annotation']:
+ pt_df.loc[(annotation['drug'].lower().replace(' ', '_'), var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [map_confidence_who[annotation['confidence']], 'Catalogue', annotation['original_mutation'], var['type'], '{:.0%}'.format(var['freq']), annotation['source']]
+
+ pt_df['Method'] = method
+ return pt_df
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Analyse resistance output from magma Pipeline')
+ parser.add_argument('merged_cohort_stats_file', metavar='merged_cohort_stats_file', type=str, help='The file containing the merged cohort statistics')
+ parser.add_argument('major_res_var_dir', metavar='major_res_var_dir', type=str, help='The directory containing the major variants TBProfiler output files')
+ parser.add_argument('minor_res_var_dir', metavar='minor_res_var_dir', type=str, help='The directory containing the minor variants TBProfiler output files')
+ parser.add_argument('struc_res_var_dir', metavar='struc_res_var_dir', type=str, help='The directory containing the structural variants TBProfiler output files')
+ parser.add_argument('summary_output_dir', metavar='summary_output_dir', type=str, help='The directory where the resulting excel sheets should be placed')
+ args = vars(parser.parse_args())
+
+ dir_summary = args['summary_output_dir']
+ if not os.path.exists(dir_summary):
+ os.makedirs(dir_summary)
+
+ dir_summary_json = args['summary_output_dir'] + "/json_format"
+ if not os.path.exists(dir_summary_json):
+ os.makedirs(dir_summary_json)
+
+
+ samples = {}
+
+ dir_major_vars = args['major_res_var_dir'] + "/" + 'results'
+ json_major_vars = glob.glob(dir_major_vars + "/" + "*.json")
+ for file_name in json_major_vars:
+ keys = '.'.join(file_name.split('.')[:-2]).split('/')[-1]
+ samples[keys] = {}
+ with open(file_name) as json_file:
+ samples[keys]['magma'] = json.load(json_file)
+
+ dir_minor_vars = args['minor_res_var_dir'] + "/" + 'results'
+ if os.path.exists(dir_minor_vars):
+ json_minor_vars = glob.glob(dir_minor_vars + "/" + "*.json")
+ for file_name in json_minor_vars:
+ keys = '.'.join(file_name.split('.')[:-2]).split('/')[-1]
+ if keys not in samples:
+ continue
+ with open(file_name) as json_file:
+ samples[keys]['lofreq'] = json.load(json_file)
+
+ dir_struc_vars = args['struc_res_var_dir'] + "/" + 'results'
+ if os.path.exists(dir_struc_vars):
+ json_struc_vars = glob.glob(dir_struc_vars + "/" + "*.json")
+ for file_name in json_struc_vars:
+ keys = '.'.join(file_name.split('.')[:-2]).split('/')[-1]
+ if keys not in samples:
+ continue
+ with open(file_name) as json_file:
+ samples[keys]['delly'] = json.load(json_file)
+
+#===============
+# ADD FILTER FOR SAMPLES PASSING << ALL_THRESHOLDS >>
+#===============
+ stats_df = pd.read_csv(args["merged_cohort_stats_file"], sep="\t")
+ filtered_stats_df = stats_df[stats_df["ALL_THRESHOLDS_MET"] == 1]
+
+ samples_df = pd.DataFrame(list(samples), columns=['full_sample'])
+ filtered_samples_df = samples_df[samples_df["full_sample"].isin(filtered_stats_df["SAMPLE"].to_list())]
+ samples_df = filtered_samples_df.set_index('full_sample').sort_index()
+
+#===============
+# POPULATE THE DATAFRAME
+#===============
+ for patient, sample in tqdm(samples_df.iterrows(), total=samples_df.shape[0]):
+ sample_res = samples[patient]
+
+ pt_df_magma = create_resistance_df(sample_res['magma'], 'XBS')
+
+ if 'lofreq' in sample_res:
+ pt_df_lof = create_resistance_df(sample_res['lofreq'], 'LoFreq')
+ else:
+ pt_df_lof = pd.DataFrame(columns=['Drug', 'Variant']).set_index(['Drug', 'Variant'])
+
+ if 'delly' in sample_res:
+ pt_df_delly = create_resistance_df(sample_res['delly'], 'Delly')
+ else:
+ pt_df_delly = pd.DataFrame(columns=['Drug', 'Variant']).set_index(['Drug', 'Variant'])
+
+ pt_df = pd.concat([pt_df_magma, pt_df_lof, pt_df_delly]).reset_index().drop_duplicates(subset=['Drug', 'Variant']).set_index(['Drug', 'Variant']).sort_index()
+
+ for drug in list(drugs - set([i[0] for i in pt_df.index.values])):
+ pt_df.loc[(drug, 'No variants found'), 'Resistance interpretation'] = 6
+
+ pt_df['Resistance interpretation'].replace(-1, unknown_positions, inplace=True)
+ for drug in pt_df.index.levels[0]:
+ conc = min(pt_df.loc[drug, 'Resistance interpretation'].value_counts().keys())
+ pt_df.loc[drug, 'Conclusion'] = conc
+
+ pt_df = pt_df.reset_index().sort_values(['Conclusion', 'Drug', 'Resistance interpretation', 'Variant'])
+ for column in ['Conclusion', 'Resistance interpretation']:
+ pt_df[column] = pt_df[column].apply(lambda c: map_resistance_class[-1] if c == unknown_positions else None if pd.isna(c) else map_resistance_class[c])
+
+ pt_df = pt_df[['Drug', 'Conclusion', 'Variant', 'Resistance interpretation', 'Type', 'Frequency', 'Method', 'Literature', 'Source notation']]
+
+ # Write a json output
+ json_data = pt_df.to_json()
+ with open(os.path.join(dir_summary_json, '{}.json'.format(patient)), 'w') as f:
+ f.write(json_data)
+
+ # Write the sheet to excel with formatting
+ with pd.ExcelWriter(os.path.join(dir_summary, '{}.xlsx'.format(patient)), engine='xlsxwriter') as writer:
+ pt_df.set_index(['Drug', 'Conclusion', 'Variant']).to_excel(writer, sheet_name='resistance_variants')
+
+ format_sens = writer.book.add_format({'bold': False, 'font_color': 'green'})
+ format_mod = writer.book.add_format({'bold': True, 'font_color': 'orange'})
+ format_res = writer.book.add_format({'bold': True, 'font_color': 'red'})
+ cond_res_0 = {'type': 'cell', 'criteria': '==', 'value': '"R"', 'format': format_res}
+ cond_res_1 = {'type': 'cell', 'criteria': '==', 'value': '"Very high probability"', 'format': format_res}
+ cond_res_2 = {'type': 'cell', 'criteria': '==', 'value': '"High probability"', 'format': format_res}
+ cond_mod_3 = {'type': 'cell', 'criteria': '==', 'value': '"Moderate probability"', 'format': format_mod}
+ cond_sens_4 = {'type': 'cell', 'criteria': '==', 'value': '"Low probability"', 'format': format_sens}
+ cond_sens_5 = {'type': 'cell', 'criteria': '==', 'value': '"Very low probability"', 'format': format_sens}
+ cond_sens_6 = {'type': 'cell', 'criteria': '==', 'value': '"S"', 'format': format_sens}
+
+ # Add formatting to Variants sheet
+ for cond_format in [cond_res_0, cond_res_1, cond_res_2, cond_mod_3, cond_sens_4, cond_sens_5, cond_sens_6]:
+ writer.sheets['resistance_variants'].conditional_format('{}1:{}{}'.format(alphabet[1], alphabet[1], pt_df.shape[0] + 1), cond_format)
+ writer.sheets['resistance_variants'].conditional_format('{}1:{}{}'.format(alphabet[3], alphabet[3], pt_df.shape[0] + 1), cond_format)
+
+ # Autofit the worksheet and hide columns
+ # writer.sheets['resistance_variants'].autofit() # Not available in the version we are using
+ writer.sheets['resistance_variants'].set_column(2, 2, 25, None, None)
+ writer.sheets['resistance_variants'].set_column(5, 5, 25, None, None)
diff --git a/bin/summarize_resistance_mixed_infection.py b/bin/summarize_resistance_mixed_infection.py
new file mode 100755
index 0000000..6a7c5e5
--- /dev/null
+++ b/bin/summarize_resistance_mixed_infection.py
@@ -0,0 +1,236 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+import glob
+import json
+import os
+import re
+import argparse
+
+import pandas as pd
+import numpy as np
+
+from tqdm import tqdm
+
+alphabet = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
+# who_v2_source = 'https://www.who.int/publications/i/item/9789240082410'
+
+drugs = {'amikacin', 'bedaquiline', 'capreomycin', 'clofazimine', 'cycloserine', 'delamanid', 'ethambutol',
+ 'ethionamide', 'imipenem', 'isoniazid', 'kanamycin','levofloxacin', 'linezolid', 'meropenem', 'moxifloxacin',
+ 'para_aminosalicylic_acid', 'pretomanid', 'prothionamide', 'pyrazinamide', 'rifabutin', 'rifampicin',
+ 'streptomycin', 'terizidone'}
+
+unknown_positions = 2.5
+
+map_resistance_class = {
+ 0: 'R',
+ 1: 'Very high probability',
+ 2: 'High probability',
+ 3: 'Moderate probability',
+ 4: 'Low probability',
+ 5: 'Very low probability',
+ 6: 'S',
+ -1: 'Unknown'
+}
+
+map_confidence_who = {
+ 'Assoc w R': 0,
+ 'Assoc w R - Interim': 1,
+ 'Uncertain significance': -1,
+ 'Not assoc w R - Interim': 5,
+ 'Not assoc w R': 6,
+}
+
+# drug_mapping = {
+# 'AMI': 'amikacin',
+# 'BDQ': 'bedaquiline',
+# 'CAP': 'capreomycin',
+# 'CFZ': 'clofazimine',
+# 'DLM': 'delamanid',
+# 'EMB': 'ethambutol',
+# 'ETH': 'ethionamide',
+# 'INH': 'isoniazid',
+# 'KAN': 'kanamycin',
+# 'LEV': 'levofloxacin',
+# 'LZD': 'linezolid',
+# 'MXF': 'moxifloxacin',
+# 'PZA': 'pyrazinamide',
+# 'RIF': 'rifampicin',
+# 'STM': 'streptomycin',
+# 'PTH': 'prothionamide'
+# }
+#
+# drug_mapping_inv = {
+# item[1]: item[0] for item in drug_mapping.items()
+# }
+
+
+def create_resistance_df(sample_res, method):
+ pt_df = pd.DataFrame(columns=['Drug', 'Variant', 'Resistance interpretation','Source', 'Source notation', 'Type', 'Frequency', 'Literature']).set_index(['Drug', 'Variant'])
+
+ # DR variants from the MAGMA analysis. Add after LoFreq DR variants to give priority to MAGMA frequencies.
+ for var in sample_res['dr_variants']:
+ gene = var['gene_name']
+ var_repr = '{}_{}'.format(gene, var['change'])
+
+ # Add all the DR variants with classification and overwrite later if necessary
+ for drug in var['gene_associated_drugs']:
+ drug_name = drug.lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [unknown_positions, 'non-Catalogue', '', var['type'], '{:.0%}'.format(var['freq']), 'Manually curated']
+
+ # Overwrite the variant classification for drugs which have a WHO classification last as to overrule all other classifications
+ for drug in var['drugs']:
+ drug_name = drug['drug'].lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [map_confidence_who[drug['confidence']], 'Catalogue', drug['original_mutation'], var['type'], '{:.0%}'.format(var['freq']), drug['source']]
+
+ # Other variants from MAGMA analysis.
+ for var in sample_res['other_variants']:
+ gene = var['gene_name']
+ var_repr = '{}_{}'.format(gene, var['change'])
+
+ # Add all the other variants wi classification and overwrite later if necessary
+ for drug in var['gene_associated_drugs']:
+ drug_name = drug.lower().replace(' ', '_')
+ pt_df.loc[(drug_name, var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [unknown_positions, 'non-Catalogue', '', var['type'], '{:.0%}'.format(var['freq']), 'Manually curated']
+
+ # Overwrite the variant classification for drugs which have a WHO classification last as to overrule all other classifications
+ if 'annotation' in var:
+ for annotation in var['annotation']:
+ pt_df.loc[(annotation['drug'].lower().replace(' ', '_'), var_repr), ('Resistance interpretation', 'Source', 'Source notation', 'Type', 'Frequency', 'Literature')] = [map_confidence_who[annotation['confidence']], 'Catalogue', annotation['original_mutation'], var['type'], '{:.0%}'.format(var['freq']), annotation['source']]
+
+ pt_df['Method'] = method
+ return pt_df
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description='Analyse resistance output from magma Pipeline')
+ parser.add_argument('merged_cohort_stats_file', metavar='merged_cohort_stats_file', type=str, help='The file containing the merged cohort statistics')
+ parser.add_argument('minor_res_var_dir', metavar='minor_res_var_dir', type=str, help='The directory containing the minor variants TBProfiler output files')
+ parser.add_argument('struc_res_var_dir', metavar='struc_res_var_dir', type=str, help='The directory containing the structural variants TBProfiler output files')
+ parser.add_argument('summary_output_dir', metavar='summary_output_dir', type=str, help='The directory where the resulting CSV files should be placed')
+ args = vars(parser.parse_args())
+
+ dir_summary = args['summary_output_dir']
+ if not os.path.exists(dir_summary):
+ os.makedirs(dir_summary)
+
+ dir_summary_json = args['summary_output_dir'] + "/json_format"
+ if not os.path.exists(dir_summary_json):
+ os.makedirs(dir_summary_json)
+
+ samples = {}
+
+ dir_minor_vars = args['minor_res_var_dir'] + "/" + 'results'
+ json_minor_vars = glob.glob(dir_minor_vars + "/" + "*.json")
+ for file_name in json_minor_vars:
+ keys = '.'.join(file_name.split('.')[:-2]).split('/')[-1]
+ samples[keys] = {}
+ with open(file_name) as json_file:
+ samples[keys]['lofreq'] = json.load(json_file)
+
+
+ dir_struc_vars = args['struc_res_var_dir'] + "/" + 'results'
+ if os.path.exists(dir_struc_vars):
+ json_struc_vars = glob.glob(dir_struc_vars + "/" + "*.json")
+ for file_name in json_struc_vars:
+ keys = '.'.join(file_name.split('.')[:-2]).split('/')[-1]
+ if keys not in samples:
+ continue
+ with open(file_name) as json_file:
+ samples[keys]['delly'] = json.load(json_file)
+
+#===============
+# ADD FILTER FOR SAMPLES FAILING ONLY << RELABUNDANCE THRESHOLD_MET >>
+#===============
+ stats_df = pd.read_csv(args["merged_cohort_stats_file"], sep="\t")
+ filtered_stats_df = stats_df.loc[ (stats_df["RELABUNDANCE_THRESHOLD_MET"]==0) & (stats_df["MAPPED_NTM_FRACTION_16S_THRESHOLD_MET"]==1) & (stats_df["COVERAGE_THRESHOLD_MET"]==1) & (stats_df["BREADTH_OF_COVERAGE_THRESHOLD_MET"]==1)]
+
+ samples_df = pd.DataFrame(list(samples), columns=['full_sample'])
+ filtered_samples_df = samples_df[samples_df["full_sample"].isin(filtered_stats_df["SAMPLE"].to_list())]
+ samples_df = filtered_samples_df.set_index('full_sample').sort_index()
+
+
+#===============
+# POPULATE THE DATAFRAME
+#===============
+ for patient, sample in tqdm(samples_df.iterrows(), total=samples_df.shape[0]):
+ sample_res = samples[patient]
+
+ if 'lofreq' in sample_res:
+ pt_df_lof = create_resistance_df(sample_res['lofreq'], 'LoFreq')
+ else:
+ pt_df_lof = pd.DataFrame(columns=['Drug', 'Variant']).set_index(['Drug', 'Variant'])
+
+ if 'delly' in sample_res:
+ pt_df_delly = create_resistance_df(sample_res['delly'], 'Delly')
+ else:
+ pt_df_delly = pd.DataFrame(columns=['Drug', 'Variant']).set_index(['Drug', 'Variant'])
+
+ pt_df = pd.concat([pt_df_lof, pt_df_delly]).reset_index().drop_duplicates(subset=['Drug', 'Variant']).set_index(['Drug', 'Variant']).sort_index()
+
+ for drug in list(drugs - set([i[0] for i in pt_df.index.values])):
+ pt_df.loc[(drug, 'No variants found'), 'Resistance interpretation'] = 6
+
+ pt_df['Resistance interpretation'].replace(-1, unknown_positions, inplace=True)
+ for drug in pt_df.index.levels[0]:
+ conc = min(pt_df.loc[drug, 'Resistance interpretation'].value_counts().keys())
+ pt_df.loc[drug, 'Conclusion'] = conc
+
+ pt_df = pt_df.reset_index().sort_values(['Conclusion', 'Drug', 'Resistance interpretation', 'Variant'])
+ for column in ['Conclusion', 'Resistance interpretation']:
+ pt_df[column] = pt_df[column].apply(lambda c: map_resistance_class[-1] if c == unknown_positions else None if pd.isna(c) else map_resistance_class[c])
+
+ pt_df = pt_df[['Drug', 'Conclusion', 'Variant', 'Resistance interpretation', 'Type', 'Frequency', 'Method', 'Literature', 'Source notation']]
+
+ # Write a json output
+ json_data = pt_df.to_json()
+ with open(os.path.join(dir_summary_json, '{}.json'.format(patient)), 'w') as f:
+ f.write(json_data)
+
+ # Write the sheet to excel with formatting
+ with pd.ExcelWriter(os.path.join(dir_summary, '{}.xlsx'.format(patient)), engine='xlsxwriter') as writer:
+ pt_df.set_index(['Drug', 'Conclusion', 'Variant']).to_excel(writer, sheet_name='resistance_variants') # Updated sheet name
+
+ format_sens = writer.book.add_format({'bold': False, 'font_color': 'green'})
+ format_mod = writer.book.add_format({'bold': True, 'font_color': 'orange'})
+ format_res = writer.book.add_format({'bold': True, 'font_color': 'red'})
+ cond_res_0 = {'type': 'cell', 'criteria': '==', 'value': '"R"', 'format': format_res}
+ cond_res_1 = {'type': 'cell', 'criteria': '==', 'value': '"Very high probability"', 'format': format_res}
+ cond_res_2 = {'type': 'cell', 'criteria': '==', 'value': '"High probability"', 'format': format_res}
+ cond_mod_3 = {'type': 'cell', 'criteria': '==', 'value': '"Moderate probability"', 'format': format_mod}
+ cond_sens_4 = {'type': 'cell', 'criteria': '==', 'value': '"Low probability"', 'format': format_sens}
+ cond_sens_5 = {'type': 'cell', 'criteria': '==', 'value': '"Very low probability"', 'format': format_sens}
+ cond_sens_6 = {'type': 'cell', 'criteria': '==', 'value': '"S"', 'format': format_sens}
+
+ # Add formatting to Variants sheet
+ for cond_format in [cond_res_0, cond_res_1, cond_res_2, cond_mod_3, cond_sens_4, cond_sens_5, cond_sens_6]:
+ writer.sheets['resistance_variants'].conditional_format('{}1:{}{}'.format(alphabet[1], alphabet[1], pt_df.shape[0]+1), cond_format)
+ writer.sheets['resistance_variants'].conditional_format('{}1:{}{}'.format(alphabet[3], alphabet[3], pt_df.shape[0]+1), cond_format)
+
+ # Autofit the worksheet and hide columns
+ # writer.sheets['resistance_variants'].autofit() # Not available in the version we are using
+ writer.sheets['resistance_variants'].set_column(2, 2, 25, None, None)
+ writer.sheets['resistance_variants'].set_column(5, 5, 25, None, None)
diff --git a/bin/variant_table_to_fasta.py b/bin/variant_table_to_fasta.py
new file mode 100755
index 0000000..bb33818
--- /dev/null
+++ b/bin/variant_table_to_fasta.py
@@ -0,0 +1,55 @@
+#! /usr/bin/env python3
+#
+# Copyright (c) 2021-2024 MAGMA pipeline authors, see https://doi.org/10.1371/journal.pcbi.1011648
+#
+# This file is part of MAGMA pipeline, see https://github.com/TORCH-Consortium/MAGMA
+#
+# For quick overview of GPL-3 license, please refer
+# https://www.tldrlegal.com/license/gnu-general-public-license-v3-gpl-3
+#
+# - You MUST keep this license with original authors in your copy
+# - You MUST acknowledge the original source of this software
+# - You MUST state significant changes made to the original software
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program . If not, see .
+#
+
+import sys
+import argparse
+
+def main(args):
+ table = []
+ with open(args.table, 'r') as table_file:
+ table.append(table_file.readline().strip().split('\t')) # Get the headerline without modifying
+ # Process the actual variants
+ for idx, l in enumerate(table_file):
+ l = l.strip().split('\t')
+ l = [i.replace('*', '-').replace('.', '-') for i in l]
+ if l.count('-')/len(l) < (1-args.cutoff_site_representation):
+ table.append(l)
+ else:
+ pass
+ with open(args.output_fasta, 'w') as fasta_file:
+ for l in list(map(list, zip(*table))):
+ fasta_file.write('>{}\n{}\n'.format(l[0].replace('.GT', ''), ''.join(l[1:])))
+
+
+
+parser = argparse.ArgumentParser(description='tbprofiler script',formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument('table', type=str, help='The input table to convert (stdin if empty)')
+parser.add_argument('output_fasta', type=str, help='The output fasta file')
+parser.add_argument('cutoff_site_representation', type=float, help='Minimum fraction of samples that need to have a call at a site before it is considered')
+parser.set_defaults(func=main)
+args = parser.parse_args()
+args.func(args)
diff --git a/conf/magma_base.config b/conf/magma_base.config
new file mode 100644
index 0000000..55335a6
--- /dev/null
+++ b/conf/magma_base.config
@@ -0,0 +1,72 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ MAGMA mode default parameters (pinned to MAGMA v0.3.0 / f33a913)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ Loaded when params.mode == 'magma'. References bundled under
+ data/magma-references/ (vendored from mycobactopia-org/MAGMA tag v0.3.0).
+----------------------------------------------------------------------------------------
+*/
+
+params {
+ magma_input_samplesheet = null // backward-compat alias; prefer --input
+ magma_vcf_name = 'joint'
+ magma_use_ref_gvcf = true
+ magma_ref_gvcf = "${projectDir}/data/magma-references/ref_gvcfs/LineagesAndOutgroupV2.g.vcf.gz"
+ magma_ref_gvcf_tbi = "${projectDir}/data/magma-references/ref_gvcfs/LineagesAndOutgroupV2.g.vcf.gz.tbi"
+ magma_cutoff_median_coverage = 10
+ magma_cutoff_breadth_of_coverage = 0.90
+ magma_cutoff_rel_abundance = 0.70
+ magma_cutoff_ntm_fraction = 0.20
+ magma_cutoff_site_representation = 0.95
+ magma_cutoff_strand_bias = 0 // 0 = skip strand-bias filter
+ magma_only_validate_fastqs = false
+ magma_skip_merge_analysis = false
+ magma_skip_variant_recalibration = false // VQSR optimisation
+ magma_skip_base_recalibration = true // BQSR — not recommended for Mtb
+ magma_skip_minor_variants_gatk = true // LoFreq preferred
+ magma_skip_phylogeny_and_clustering = false
+ magma_skip_complex_regions = false
+ magma_skip_ntmprofiler = false
+ magma_skip_tbprofiler_fastq = true
+ magma_skip_spotyping = true
+ magma_skip_rdanalyzer = true // experimental, not production-ready
+ magma_skip_ismapper = false
+ magma_iqtree_standard_bootstrap = true
+ magma_iqtree_fast_ml_only = false
+ magma_iqtree_fast_bootstrapped_phylogeny = false
+ magma_iqtree_accurate_ml_only = false
+ magma_ref_fasta_basename = "NC-000962-3-H37Rv"
+ magma_ref_fasta_dir = "${projectDir}/data/magma-references/genome"
+ magma_ref_fasta = "${params.magma_ref_fasta_dir}/${params.magma_ref_fasta_basename}.fa"
+ magma_ref_fasta_gb = "${params.magma_ref_fasta_dir}/${params.magma_ref_fasta_basename}.gb"
+ magma_ref_fasta_dict = "${params.magma_ref_fasta_dir}/${params.magma_ref_fasta_basename}.dict"
+ magma_ref_fasta_amb = "${params.magma_ref_fasta}.amb"
+ magma_ref_fasta_ann = "${params.magma_ref_fasta}.ann"
+ magma_ref_fasta_bwt = "${params.magma_ref_fasta}.bwt"
+ magma_ref_fasta_fai = "${params.magma_ref_fasta}.fai"
+ magma_ref_fasta_pac = "${params.magma_ref_fasta}.pac"
+ magma_ref_fasta_sa = "${params.magma_ref_fasta}.sa"
+ magma_queries_multifasta = "${projectDir}/data/magma-references/regions/IS_Queries_Mtb.multi.fasta"
+ magma_drgenes_list = "${projectDir}/data/magma-references/regions/WHO_Tier1_Tier2_DR.list"
+ magma_rrna_list = "${projectDir}/data/magma-references/regions/rRNA.list"
+ magma_excluded_loci_list = "${projectDir}/data/magma-references/regions/UVP_List_of_Excluded_loci.list"
+ magma_dbsnp_vcf = "${projectDir}/data/magma-references/known/Benavente2015.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_dbsnp_vcf_tbi = "${params.magma_dbsnp_vcf}.tbi"
+ magma_benavente2015_vcf = "${projectDir}/data/magma-references/known/Benavente2015.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_benavente2015_vcf_tbi = "${params.magma_benavente2015_vcf}.tbi"
+ magma_coll2014_vcf = "${projectDir}/data/magma-references/truth/Coll2014.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_coll2014_vcf_tbi = "${params.magma_coll2014_vcf}.tbi"
+ magma_coll2018_vcf = "${projectDir}/data/magma-references/truth/Coll2018.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_coll2018_vcf_tbi = "${params.magma_coll2018_vcf}.tbi"
+ magma_napier2020_vcf = "${projectDir}/data/magma-references/truth/Napier2020.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_napier2020_vcf_tbi = "${params.magma_napier2020_vcf}.tbi"
+ magma_walker2015_vcf = "${projectDir}/data/magma-references/truth/Walker2015.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_walker2015_vcf_tbi = "${params.magma_walker2015_vcf}.tbi"
+ magma_zeng2018_vcf = "${projectDir}/data/magma-references/truth/Zeng2018.UVPapproved.rRNAexcluded.vcf.gz"
+ magma_zeng2018_vcf_tbi = "${params.magma_zeng2018_vcf}.tbi"
+ magma_ref_fasta_rdanalyzer = "${projectDir}/data/magma-references/rdanalyzer/RDs30.fasta"
+ magma_snpeff_genome_id = "Mycobacterium_tuberculosis_h37rv"
+ magma_snpeff_config_file = "${projectDir}/data/magma-references/snpeff/snpEff.config"
+ magma_snpeff_db_dir = "${projectDir}/data/magma-references/snpeff/data/${params.magma_snpeff_genome_id}"
+ magma_should_publish_bam = false // publish BAMs from HaplotypeCaller
+}
diff --git a/conf/magma_containers.config b/conf/magma_containers.config
new file mode 100644
index 0000000..0134259
--- /dev/null
+++ b/conf/magma_containers.config
@@ -0,0 +1,32 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ MAGMA mode container assignments (vendored from MAGMA v0.3.0 conf/abc_cluster.config)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ MAGMA's local modules (python utils + tool wrappers) carry no container
+ directive of their own — they rely on these per-process-name assignments.
+ Required on any containerised executor (Docker / Nomad); without them a
+ process like SAMPLESHEET_VALIDATION fails with "Missing container image".
+ Mirrors MAGMA's validated container layout (magma-container-1/2 + the
+ per-tool biocontainers), including the functional DELLY/BCFTOOLS overrides
+ where magma-container's bundled tool is too old for the nf-core command.
+----------------------------------------------------------------------------------------
+*/
+
+process {
+
+ withName: '.*SPOTYPING.*' { container = "quay.io/biocontainers/spotyping:2.1--hdfd78af_4" }
+ withName: '.*RDANALYZER.*' { container = "quay.io/biocontainers/rd-analyzer:1.01--hdfd78af_0" }
+ withName: '.*TBPROFILER.*' { container = "ghcr.io/torch-consortium/magma/biocontainer-tbprofiler:6.3.0--1" }
+ withName: 'NTMPROFILER.*' { container = "ghcr.io/torch-consortium/magma/biocontainer-ntmprofiler:0.4.0" }
+ withName: 'ISMAPPER.*|GATK.*|LOFREQ.*|DELLY.*|MULTIQC.*|FASTQC.*|UTILS.*|FASTQ.*|SAMPLESHEET.*' {
+ container = "ghcr.io/torch-consortium/magma/magma-container-1:2.0.0"
+ }
+ withName: 'BWA.*|IQTREE.*|SNPDISTS.*|SNPSITES.*|BCFTOOLS.*|BGZIP.*|SAMTOOLS.*|SNPEFF.*|CLUSTERPICKER.*' {
+ container = "ghcr.io/torch-consortium/magma/magma-container-2:1.1.1"
+ }
+
+ // Functional container overrides: the nf-core delly/call + bcftools/view
+ // module commands use flags newer than the tools bundled in magma-container.
+ withName: '.*DELLY_CALL.*' { container = 'quay.io/biocontainers/delly:1.7.3--hd6466ae_0' }
+ withName: '.*BCFTOOLS_VIEW.*' { container = 'quay.io/biocontainers/bcftools:1.21--h3a4d415_1' }
+}
diff --git a/conf/magma_modules.config b/conf/magma_modules.config
new file mode 100644
index 0000000..9c3d664
--- /dev/null
+++ b/conf/magma_modules.config
@@ -0,0 +1,1042 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ Config file for defining DSL2 per-module options and publishing paths
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ Available keys to override module options:
+ ext.args = Additional arguments appended to command in module.
+ ext.args2 = Second set of arguments appended to command in module (multi-tool modules).
+ ext.prefix = File name prefix for output files.
+ ext.when = Condition for running module (null or true = run; false = skip).
+----------------------------------------------------------------------------------------
+*/
+
+process {
+
+ // ================================================================================
+ // Default publishing (overridden per process below)
+ // ================================================================================
+ publishDir = [
+ path: { "${params.outdir}/${task.process.tokenize(':')[-1].tokenize('_')[0].toLowerCase()}" },
+ mode: params.publish_dir_mode,
+ saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+ ]
+
+ // ================================================================================
+ // nf-core modules (FastQC, MultiQC)
+ // ================================================================================
+
+ withName: 'FASTQC' {
+ ext.args = '--quiet'
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/per_sample/fastqc/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'MULTIQC' {
+ ext.args = { params.multiqc_title ? "--title \"$params.multiqc_title\"" : '' }
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/cohort/multiqc/" },
+ mode: params.publish_dir_mode,
+ saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+ ]
+ }
+
+ // ================================================================================
+ // Validation & QC
+ // ================================================================================
+
+ withName: 'SAMPLESHEET_VALIDATION' {
+ publishDir = [
+ path: { "${params.outdir}/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ errorStrategy = 'terminate'
+ }
+
+ withName: 'FASTQ_VALIDATOR' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/per_sample/fastq_validation/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'UTILS_FASTQ_COHORT_VALIDATION' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/cohort/fastq_validation/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'UTILS_MULTIPLE_INFECTION_FILTER' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/cohort/multiple_infection_filter" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ // ================================================================================
+ // Mapping — BWA + SAMTOOLS_SORT
+ // ================================================================================
+
+ withName: 'BWA_MEM' {
+ // -k 100: ~2/3 of typical Mtb read length for optimal contaminant removal
+ ext.args = '-M'
+ ext.args2 = '-k 100'
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/mapped_libraries/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'BWA_MEM_DELLY' {
+ // Default -k 19 for DELLY-compatible alignments
+ ext.args = '-M'
+ ext.args2 = ''
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/mapped_libraries/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ // ================================================================================
+ // SAMTOOLS merge / index / stats
+ // ================================================================================
+
+ withName: 'SAMTOOLS_MERGE.*' {
+ // -f (force-overwrite) was hardcoded in the local module's script;
+ // nf-core uses ext.args. Same MAGMA filename pattern reproduced via
+ // ext.prefix (nf-core appends .${file_type} for the extension).
+ ext.args = '-f'
+ ext.prefix = { "${meta.id}.sorted_reads" }
+ }
+
+ withName: 'SAMTOOLS_MERGE' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/merged_libraries/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'SAMTOOLS_MERGE_DELLY' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/merged_libraries/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'SAMTOOLS_INDEX' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/marked_duplicates/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'SAMTOOLS_INDEX_DELLY' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/marked_duplicates/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'SAMTOOLS_INDEX_LOFREQ' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'SAMTOOLS_STATS' {
+ ext.args = '-F DUP,SUPPLEMENTARY,SECONDARY,UNMAP,QCFAIL'
+ // Reproduce the MAGMA filename pattern so published outputs stay
+ // grep-able as .SamtoolStats.* (nf-core appends `.stats`).
+ ext.prefix = { "${meta.id}.SamtoolStats" }
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/per_sample/mapping/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // GATK per-sample
+ // ================================================================================
+
+ withName: 'GATK_MARK_DUPLICATES' {
+ // nf-core gatk4/markduplicates uses ext.prefix as the OUTPUT BAM filename
+ // (including its .bam extension; default prefix is '${meta.id}.bam'). We
+ // set it to '${meta.id}.dedup_reads.bam' to preserve MAGMA's filename
+ // pattern. The metrics file becomes '${prefix}.metrics' →
+ // '${meta.id}.dedup_reads.bam.metrics' (was '${meta.id}.MarkDupMetrics.txt')
+ // — metrics is publish-only and not consumed downstream, so the rename
+ // is harmless.
+ ext.prefix = { "${meta.id}.dedup_reads.bam" }
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/marked_duplicates/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_MARK_DUPLICATES_DELLY' {
+ ext.prefix = { "${meta.id}.dedup_reads.bam" }
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/marked_duplicates/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_BASE_RECALIBRATOR' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/base_recalibrated/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_BASE_RECALIBRATOR_DELLY' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/base_recalibrated/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_APPLY_BQSR.*' {
+ // nf-core gatk4/applybqsr defaults to CRAM output; force BAM to keep
+ // the MAGMA chain (which always operates on BAM downstream) intact.
+ ext.suffix = 'bam'
+ // Preserve the MAGMA filename pattern: ${meta.id}.baserecalibrated.bam.
+ ext.prefix = { "${meta.id}.baserecalibrated" }
+ }
+
+ withName: 'GATK_APPLY_BQSR' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/base_recalibrated/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_APPLY_BQSR_DELLY' {
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/mapped_sequences/base_recalibrated/delly/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_HAPLOTYPE_CALLER' {
+ // The local module hardcoded -ERC GVCF + -bamout in its script; nf-core
+ // doesn't, so they move into ext.args. -bamout writes the realigned BAM
+ // to a custom filename matching MAGMA's pattern (nf-core would name it
+ // .realigned.bam — we want .haplotype_caller.bam).
+ ext.args = { "-ERC GVCF -ploidy 1 --dont-use-soft-clipped-bases --read-filter MappingQualityNotZeroReadFilter -G StandardAnnotation -G AS_StandardAnnotation --bam-output ${meta.id}.haplotype_caller.bam" }
+ // VCF filename ${meta.id}.g.vcf.gz (nf-core appends .vcf.gz to prefix).
+ ext.prefix = { "${meta.id}.g" }
+ publishDir = [
+ [
+ path: { "${params.outdir}/extra/per_sample/major_variants/haplotype_caller/" },
+ mode: params.publish_dir_mode,
+ pattern: "*vcf*",
+ enabled: false
+ ],
+ [
+ path: { "${params.outdir}/extra/per_sample/major_variants/haplotype_caller/" },
+ mode: params.publish_dir_mode,
+ pattern: "*bam*",
+ enabled: { params.magma_should_publish_bam }
+ ]
+ ]
+ }
+
+ withName: 'GATK_HAPLOTYPE_CALLER_MINOR_VARIANTS' {
+ ext.args = """-ploidy 1 \
+ --minimum-mapping-quality 60 \
+ --min-base-quality-score 20 \
+ --read-filter MappingQualityNotZeroReadFilter \
+ -G StandardAnnotation \
+ --dont-use-soft-clipped-bases \
+ --output-mode EMIT_ALL_ACTIVE_SITES"""
+ ext.when = { !params.magma_skip_minor_variants_gatk }
+ publishDir = [
+ path: { "${params.outdir}/extra/per_sample/minor_variants/haplotype_caller/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_COLLECT_WGS_METRICS' {
+ ext.args = '--READ_LENGTH 0 --COVERAGE_CAP 10000 --COUNT_UNPAIRED'
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/per_sample/coverage/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_FLAG_STAT' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/per_sample/mapping/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // LoFreq minor-variant calling
+ // ================================================================================
+
+ withName: 'LOFREQ_INDELQUAL' {
+ // torch-magma's default_params set arguments="-m 60" here, but its
+ // indelqual module never interpolated them (dead config) — and
+ // `lofreq indelqual` rejects -m ("invalid option -- 'm'"). The actual
+ // torch-magma command was: lofreq indelqual -f … --dindel -o … .
+ // The nf-core module template puts ${args} before the other flags, so
+ // the emitted command becomes
+ // lofreq indelqual --dindel -f … -o ${meta.id}.dindel.bam
+ // which differs from the reference only in flag order (cosmetic).
+ ext.args = '--dindel'
+ ext.prefix = { "${meta.id}.dindel" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'LOFREQ_CALL' {
+ ext.args = '-m 60 --call-indels -Q 30'
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'LOFREQ_CALL_NTM' {
+ ext.args = '-m 60 -Q 20 -a 1'
+ ext.args2 = '1472307-1472307' // region
+ publishDir = [
+ path: { "${params.outdir}/non-tuberculous_mycobacteria/vcf_files/variants/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'LOFREQ_FILTER' {
+ ext.args = '-a 0.20'
+ // Preserve the MAGMA filename pattern. nf-core's module hardcodes
+ // output to ${prefix}.vcf.gz, so the published artifact becomes
+ // ${meta.id}.Filtered_AF.vcf.gz (was ${meta.id}.Filtered_AF.vcf,
+ // uncompressed). LOFREQ_FILTER.out is publish-only — the downstream
+ // UTILS_REFORMAT_LOFREQ chain consumes the RAW LOFREQ_CALL output,
+ // not this filtered one — so the format change is purely cosmetic.
+ ext.prefix = { "${meta.id}.Filtered_AF" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // DELLY structural variants
+ // ================================================================================
+
+ withName: 'DELLY_CALL' {
+ // torch-magma reference: delly call -g [ -u 30 -o .delly.bcf
+ // nf-core after migration (suffix='bcf'):
+ // delly call -u 30 --threads --outfile .delly.bcf --genome ][
+ // Differences are cosmetic — long-form --genome / --outfile vs -g / -o,
+ // plus the added --threads flag — but the underlying delly invocation
+ // (discovery mode, -u 30) produces the same BCF output.
+ ext.args = '-u 30'
+ ext.prefix = { "${meta.id}.delly" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // NTM profiler
+ // ================================================================================
+
+ withName: 'NTMPROFILER_PROFILE' {
+ ext.when = { !params.magma_skip_ntmprofiler }
+ publishDir = [
+ path: { "${params.outdir}/analyses/non-tuberculous_mycobacteria/per_sample/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'NTMPROFILER_COLLATE' {
+ ext.prefix = 'ntmprofiler.collate'
+ ext.when = { !params.magma_skip_ntmprofiler }
+ publishDir = [
+ path: { "${params.outdir}/analyses/non-tuberculous_mycobacteria/cohort/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // BCFtools filtering for SV annotation
+ // ================================================================================
+
+ withName: 'BCFTOOLS_VIEW_ISMAPPER' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/ismapper/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'BCFTOOLS_VIEW_GATK' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/gatk/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'BCFTOOLS_VIEW_DELLY' {
+ // Local module ran 3 commands (view → bgzip → bcftools index). nf-core
+ // does it in one call via `-Ob` (bgzipped BCF output) and
+ // `--write-index=csi` (auto-create the index). Filename pattern
+ // ${meta.id}.filtered.bcf.gz reproduced via ext.prefix + the .bcf.gz
+ // suffix nf-core adds for -Ob.
+ ext.args = """-i 'GT="1/1"' -Ob --write-index=csi"""
+ ext.prefix = { "${meta.id}.filtered" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // GATK index / select variants (per-sample SV)
+ // ================================================================================
+
+ withName: 'GATK_INDEX_FEATURE_FILE_SV' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_SELECT_VARIANTS_INCLUSION' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // ISMapper
+ // ================================================================================
+
+ withName: 'ISMAPPER' {
+ ext.when = { !params.magma_skip_ismapper }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/structural_variants/ismapper/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ // ================================================================================
+ // QC statistics utils
+ // ================================================================================
+
+ withName: 'UTILS_SAMPLE_STATS' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/samples_thresholds/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'UTILS_COHORT_STATS' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/cohort/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'UTILS_MERGE_COHORT_STATS' {
+ publishDir = [
+ path: { "${params.outdir}/QC_statistics/cohort/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // Cohort — GATK joint genotyping
+ // ================================================================================
+
+ withName: 'GATK_COMBINE_GVCFS' {
+ ext.args = '-G StandardAnnotation -G AS_StandardAnnotation'
+ // nf-core hardcodes the .combined.g.vcf.gz suffix to make the file
+ // explicitly identifiable as a multi-sample GVCF. torch-magma used
+ // ${joint}.combined.vcf.gz; the .g insertion is purely cosmetic
+ // (downstream GenotypeGVCFs treats both identically).
+ ext.prefix = { "${meta.id}" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/raw_variant_files/combined/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_GENOTYPE_GVCFS' {
+ ext.args = '-G StandardAnnotation -G AS_StandardAnnotation --sample-ploidy 1'
+ // Preserve MAGMA filename ${joint}.raw_variants.vcf.gz.
+ ext.prefix = { "${meta.id}.raw_variants" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/raw_variant_files/genotyped/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // SnpEff annotation
+ // ================================================================================
+
+ withName: 'SNPEFF' {
+ ext.args = { "-nostats -ud 100 ${params.magma_snpeff_genome_id}" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/raw_variant_files/annotated/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // BGzip
+ // ================================================================================
+
+ withName: 'BGZIP' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/raw_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'BGZIP_LOFREQ' {
+ ext.args = '-f'
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'BGZIP_MINOR_VARIANTS' {
+ ext.args = '-f'
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/minor_variants/combined_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // LoFreq / minor-variant reformat
+ // ================================================================================
+
+ withName: 'UTILS_REFORMAT_LOFREQ' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_INDEX_FEATURE_FILE_LOFREQ' {
+ // Match standard MAGMA: emit explicit `-O .tbi` (minor-variants path).
+ // Input vcf is ${meta.id}.lofreq.reformat.corrected.vcf.gz.
+ ext.args = { "-O ${meta.id}.lofreq.reformat.corrected.vcf.gz.tbi" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/per_sample/minor_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_INDEX_FEATURE_FILE_COHORT' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/raw_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // GATK SelectVariants — cohort
+ // ================================================================================
+
+ withName: 'GATK_SELECT_VARIANTS_SNP' {
+ // The local module took variantType='SNP' and prefix='raw' as positional
+ // val inputs; nf-core uses ext.args / ext.prefix. --select-type-to-include
+ // SNP appended explicitly here. Output filename joint.raw_SNP.vcf.gz.
+ ext.args = '--select-type-to-include SNP --remove-unused-alternates --exclude-non-variants'
+ ext.prefix = { "${meta.id}.raw_SNP" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_SELECT_VARIANTS_INDEL' {
+ // Same migration pattern as SNP. Output filename joint.raw_INDEL.vcf.gz.
+ ext.args = '--select-type-to-include INDEL --remove-unused-alternates --exclude-non-variants --select-type-to-include MNP --select-type-to-include MIXED'
+ ext.prefix = { "${meta.id}.raw_INDEL" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/indel_variant_files/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'GATK_SELECT_VARIANTS_EXCLUSION_SNP' {
+ // ext.intervals_mode = 'exclude' triggers the patched nf-core selectvariants
+ // module to emit --exclude-intervals (instead of the default
+ // --intervals ) for the path passed in the intervals input slot.
+ ext.intervals_mode = 'exclude'
+ ext.prefix = { "${meta.id}.filtered_SNP_exc-rRNA" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_SELECT_VARIANTS_EXCLUSION_INDEL' {
+ ext.intervals_mode = 'exclude'
+ ext.args = '--select-type-to-include MNP --select-type-to-include MIXED'
+ ext.prefix = { "${meta.id}.filtered_INDEL_exc-rRNA" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/indel_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_SELECT_VARIANTS_PHYLOGENY' {
+ ext.args = '--remove-unused-alternates --exclude-non-variants'
+ publishDir = [
+ path: { "${params.outdir}/analyses/phylogeny/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ // ================================================================================
+ // GATK VQSR
+ // ================================================================================
+
+ withName: 'GATK_VARIANT_RECALIBRATOR_SNP' {
+ ext.prefix = 'SNP'
+ ext.args = """--use-allele-specific-annotations \
+ -AS \
+ --target-titv 1.7 \
+ --truth-sensitivity-tranche 100.0 \
+ --truth-sensitivity-tranche 99.9 \
+ --truth-sensitivity-tranche 99.8 \
+ --truth-sensitivity-tranche 99.7 \
+ --truth-sensitivity-tranche 99.6 \
+ --truth-sensitivity-tranche 99.5 \
+ --truth-sensitivity-tranche 99.4 \
+ --truth-sensitivity-tranche 99.3 \
+ --truth-sensitivity-tranche 99.2 \
+ --truth-sensitivity-tranche 99.1 \
+ --truth-sensitivity-tranche 99.0 \
+ --max-gaussians 4 \
+ -mq-cap 60"""
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_VARIANT_RECALIBRATOR_INDEL' {
+ ext.prefix = 'INDEL'
+ ext.args = """-AS \
+ --target-titv 1.8 \
+ --truth-sensitivity-tranche 100.0 \
+ --truth-sensitivity-tranche 99.9 \
+ --truth-sensitivity-tranche 99.8 \
+ --truth-sensitivity-tranche 99.7 \
+ --truth-sensitivity-tranche 99.6 \
+ --truth-sensitivity-tranche 99.5 \
+ --truth-sensitivity-tranche 99.4 \
+ --truth-sensitivity-tranche 99.3 \
+ --truth-sensitivity-tranche 99.2 \
+ --truth-sensitivity-tranche 99.1 \
+ --truth-sensitivity-tranche 99.0 \
+ --max-gaussians 4 \
+ -mq-cap 60"""
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/indel_variant_files/vqsr/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // VQSR annotation-optimization variants (ANN2–ANN7) — used by optimize_variant_recalibration subworkflow.
+ //
+ // Shared SNP-VQSR args: the OPTIMIZE_VARIANT_RECALIBRATION subworkflow drives
+ // GATK_VARIANT_RECALIBRATOR over 6 annotation combinations (ANN7 → ... → ANN2)
+ // in mode=SNP. Without these flags the gatk call dies with
+ // java.lang.ClassCastException: java.util.ArrayList cannot be cast to java.lang.String
+ // at htsjdk.variant.variantcontext.CommonInfo.getAttributeAsDouble
+ // because AS_QD / AS_MQ / AS_FS / AS_SOR / AS_MQRankSum / AS_ReadPosRankSum
+ // are per-allele lists and require --use-allele-specific-annotations / -AS to
+ // be parsed. This matches torch-magma's reference command exactly. The
+ // per-alias blocks below only carry prefix / publishDir; this block carries
+ // the args (multiple matching withName blocks compose in Nextflow).
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN[2-7]' {
+ ext.args = """--use-allele-specific-annotations \
+ -AS \
+ --target-titv 1.7 \
+ --truth-sensitivity-tranche 100.0 \
+ --truth-sensitivity-tranche 99.9 \
+ --truth-sensitivity-tranche 99.8 \
+ --truth-sensitivity-tranche 99.7 \
+ --truth-sensitivity-tranche 99.6 \
+ --truth-sensitivity-tranche 99.5 \
+ --truth-sensitivity-tranche 99.4 \
+ --truth-sensitivity-tranche 99.3 \
+ --truth-sensitivity-tranche 99.2 \
+ --truth-sensitivity-tranche 99.1 \
+ --truth-sensitivity-tranche 99.0 \
+ --max-gaussians 4 \
+ -mq-cap 60"""
+ }
+
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN2' {
+ ext.prefix = 'ANN2'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann2/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN3' {
+ ext.prefix = 'ANN3'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann3/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN4' {
+ ext.prefix = 'ANN4'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann4/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN5' {
+ ext.prefix = 'ANN5'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann5/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN6' {
+ ext.prefix = 'ANN6'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann6/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'GATK_VARIANT_RECALIBRATOR_ANN7' {
+ ext.prefix = 'ANN7'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann7/" }, mode: params.publish_dir_mode ]
+ }
+
+ withName: 'GATK_APPLY_VQSR_SNP' {
+ // -mode SNP is appended to ext.args (was a positional val(analysisMode)
+ // input on the local module — nf-core uses ext.args).
+ // ext.prefix preserves the MAGMA filename pattern.
+ ext.args = '--ts-filter-level 99.90 -AS --exclude-filtered -mode SNP'
+ ext.prefix = { "${meta.id}.filtered_SNP_inc-rRNA" }
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'GATK_APPLY_VQSR_INDEL' {
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/indel_variant_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // VQSR annotation filtering (used by optimize_variant_recalibration subworkflow)
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN2' {
+ ext.prefix = 'ANN2'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann2/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN3' {
+ ext.prefix = 'ANN3'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann3/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN4' {
+ ext.prefix = 'ANN4'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann4/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN5' {
+ ext.prefix = 'ANN5'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann5/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN6' {
+ ext.prefix = 'ANN6'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann6/" }, mode: params.publish_dir_mode ]
+ }
+ withName: 'UTILS_ELIMINATE_ANNOTATION_ANN7' {
+ ext.prefix = 'ANN7'
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr__ann7/" }, mode: params.publish_dir_mode ]
+ }
+
+ withName: 'UTILS_SELECT_BEST_ANNOTATIONS' {
+ ext.when = { !params.magma_skip_variant_recalibration }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/snp_variant_files/vqsr/best_annotations/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // BCFtools merge (cohort)
+ // ================================================================================
+
+ withName: 'BCFTOOLS_MERGE_LOFREQ' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/minor_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'BCFTOOLS_MERGE_DELLY' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/structural_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // GATK MergeVcfs
+ // ================================================================================
+
+ withName: 'GATK_MERGE_VCFS.*' {
+ // Reproduce MAGMA filename: ${joint}.filtered_SNP.RawIndels.vcf.gz.
+ // The nf-core gatk4/mergevcfs module emits ${prefix}.vcf.gz; we set
+ // ext.prefix to include the MAGMA-specific suffix.
+ ext.prefix = { "${meta.id}.filtered_SNP.RawIndels" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/major_variants/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // GATK VariantsToTable
+ // ================================================================================
+
+ withName: 'GATK_VARIANTS_TO_TABLE' {
+ ext.args = '-GF GT'
+ // Filename pattern: joint..tsv (nf-core appends .tsv).
+ ext.prefix = { "${meta.id}.${meta.phylo_prefix}" }
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/multiple_alignment_files/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'UTILS_VARIANT_TABLE_TO_FASTA' {
+ // Filename pattern: joint..fa
+ ext.prefix = { "${meta.id}.${meta.phylo_prefix}" }
+ }
+
+ // ================================================================================
+ // TBprofiler
+ // ================================================================================
+
+ withName: 'TBPROFILER_VCF_PROFILE_COHORT' {
+ ext.args = '--depth 0,0 --af 0,0 --strand 0 --sv_depth 0,0 --sv_af 0,0 --sv_len 100000,50000'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/major_variants_xbs/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_COLLATE_COHORT' {
+ ext.prefix = 'major_variants'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/major_variants_xbs/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_FASTQ_PROFILE' {
+ ext.args = '--csv'
+ ext.when = { !params.magma_skip_tbprofiler_fastq }
+ publishDir = [
+ path: { "${params.outdir}/analyses/others/per_sample/tbprofiler_fastq/" },
+ mode: params.publish_dir_mode,
+ enabled: false
+ ]
+ }
+
+ withName: 'TBPROFILER_FASTQ_COLLATE' {
+ ext.prefix = 'fastq'
+ ext.when = { !params.magma_skip_tbprofiler_fastq }
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/tbprofiler_fastq/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_VCF_PROFILE_LOFREQ' {
+ ext.args = '--depth 0,0 --af 0,0 --strand 0 --sv_depth 0,0 --sv_af 0,0 --sv_len 100000,50000'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/minor_variants_lofreq/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_COLLATE_LOFREQ' {
+ ext.prefix = 'minor_variants'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/minor_variants_lofreq/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_VCF_PROFILE_DELLY' {
+ ext.args = '--depth 0,0 --af 0,0 --strand 0 --sv_depth 0,0 --sv_af 0,0 --sv_len 100000,50000'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/structural_variants_delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'TBPROFILER_COLLATE_DELLY' {
+ ext.prefix = 'structural_variants'
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/structural_variants_delly/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // SpoTyping
+ // ================================================================================
+
+ withName: 'SPOTYPING' {
+ ext.when = { !params.magma_skip_spotyping }
+ publishDir = [
+ path: { "${params.outdir}/analyses/spotyping/results_excel/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'UTILS_CAT_SPOTYPING' {
+ ext.when = { !params.magma_skip_spotyping }
+ publishDir = [
+ path: { "${params.outdir}/analyses/spotyping/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // RDAnalyzer
+ // ================================================================================
+
+ withName: 'RDANALYZER' {
+ ext.when = { !params.magma_skip_rdanalyzer }
+ publishDir = [
+ path: { "${params.outdir}/analyses/others/per_sample/rdanalyzer/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // Phylogeny & clustering
+ // ================================================================================
+
+ withName: 'SNPSITES' {
+ publishDir = [
+ path: { "${params.outdir}/vcf_files/cohort/multiple_alignment_files/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'SNPDISTS' {
+ // Drive the nf-core SNPDISTS module to emit the same scientific output
+ // as torch-magma's local snp-dists wrapper:
+ // torch-magma: `snp-dists joint.variable..fa -b > joint..snp_dists.tsv`
+ // port (nf-core): `snp-dists -b joint.variable..fa > joint..snp_dists.tsv`
+ // Arg order (-b before vs after the alignment) is cosmetic — snp-dists
+ // is order-insensitive and the output file is byte-identical.
+ ext.args = '-b'
+ ext.prefix = { "${meta.id}.${meta.phylo_prefix}.snp_dists" }
+ publishDir = [
+ path: { "${params.outdir}/analyses/snp_distances/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'IQTREE' {
+ // Per-mode iqtree args. Local module had this logic inside the script
+ // (if/else on params); nf-core uses ext.args. Closure picks the same
+ // mode as the local module.
+ ext.args = {
+ params.magma_iqtree_standard_bootstrap ? '-b 1000' :
+ params.magma_iqtree_fast_ml_only ? '-fast' :
+ params.magma_iqtree_fast_bootstrapped_phylogeny ? '-bb 1000 -alrt 1000' :
+ params.magma_iqtree_accurate_ml_only ? '-allnni' :
+ '-allnni'
+ }
+ // Reproduce MAGMA filename pattern: joint..
+ ext.prefix = { "${meta.id}.${meta.phylo_prefix}" }
+ publishDir = [
+ path: { "${params.outdir}/analyses/phylogeny/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'CLUSTERPICKER' {
+ ext.args = [
+ '0', // bootstrap_1
+ '0', // bootstrap_2
+ '0', // max_cluster_size
+ 'gap' // algorithm
+ ].join(' ')
+ publishDir = [
+ path: { "${params.outdir}/analyses/cluster_analysis/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ // ================================================================================
+ // Resistance summary utils
+ // ================================================================================
+
+ withName: 'UTILS_SUMMARIZE_RESISTANCE_RESULTS' {
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+ withName: 'UTILS_SUMMARIZE_RESISTANCE_RESULTS_MIXED_INFECTION' {
+ publishDir = [
+ path: { "${params.outdir}/analyses/drug_resistance/" },
+ mode: params.publish_dir_mode,
+ ]
+ }
+
+}
diff --git a/conf/magma_test.config b/conf/magma_test.config
new file mode 100644
index 0000000..c31720d
--- /dev/null
+++ b/conf/magma_test.config
@@ -0,0 +1,59 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ MAGMA mode test profile
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ Minimal test for --mode magma: 3 downsampled M. tuberculosis samples
+ (EXIT-RIF / PRJNA1026351: SRR26331590/595/599) on Zenodo (doi 10.5281/zenodo.20671479).
+
+ Use as follows:
+ nextflow run nf-core/tbanalyzer -profile test_magma,docker --outdir ]