This workflow stems from the discussion with @sandergoossens , the PPT was shared on the email and further details were shared here #32

Hey Tim,

I'm updating the modules for the tool paths on the cluster - as soon as I'm done, I'll ping you here so that you can start testing.

@abhi18av @LennertVerboven
please see https://github.com/TORCH-Consortium/xbs-nf/blob/master/conda_envs/setup_conda_envs.sh line 22
rm -rf resistance_db_who ./
note that as this point in the script you are located in xbs-nf/conda_envs/resistance_db_who, this command specifies the removal of two things: 1. the current directory where you are located ./ which is odd as this is not possible (you would have to go a dir up) and 2. what is presumably a directory resistance_db_who, which would mean a dir exists in this location with exactly the same name as the parent dir, guess this is not what was meant.

I am just running some other dataset atm and encountering quite a few issues with corrupted fastq files.
It would be great if we could check the file integrity of each fastq.gz with gzip -t [file] somewhere before fastqc.

Since conda is often a cause of problems on a cluster - we should update the pipeline to have

• a script to download and setup the necessary conda environments within the specified folder

• a conda_local profile which relies on this conda environment

On this thread we can keep track of the enhancements which would be nice to have.

Relevant discussions

#50
#64

@LennertVerboven , could you please confirm (in isolation) that the latest tb-profiler (v4) is working as expected?

3 annotations (DP/...) work with 4 gaussians

GATK_VARIANT_RECALIBRATOR
--max-gaussians 

Always reduce the value from 4

MERGE_WF:RESISTANCE_ANALYSIS:TBPROFILER_VCF_PROFILE__LOFREQ crashes because, when run in parallel, multiple processes are reading and writing the same files in /conda/xbs-nf-env-1-d99876e5fea88a1c4bd18887d111ae27/share/tbprofiler/

temporary solution is to keep on restarting the run, each time completing more samples. Decreasing queueSize and increasing errorStrategy retries is likely to be a good temporary solution as well.

I've created this issue to track the progress and feedback for the MERGE_WF. Will ping with further/updates on this thread.

• Change the relevant references in code
• Change the name of the repo

This is primarily a reminder for myself to include code so that GATK will correctly handle MIXED and MNP variants, we've encountered this rare issue elsewhere.

This task builds upon the data shared in #20

This task about identifying all the individual tools which have been used and to add the NF stubs to them.

So far, I've identified these many.

(base) bash-5.0$ tree -L 2
.
├── bcftools
│   └── view
├── bwa
│   └── mem
├── clusterpicker
├── delly
│   └── call
├── fastqc
├── gatk
│   ├── apply_vqsr
│   ├── base_recalibrator
│   ├── collect_wgs_metrics
│   ├── combine_gvcfs
│   ├── flag_stat
│   ├── genotype_gvcfs
│   ├── haplotype_caller
│   ├── index_feature_file
│   ├── mark_duplicates
│   ├── merge_vcfs
│   ├── select_variants
│   ├── variant_recalibrator
│   └── variants_to_table
├── iqtree
├── lofreq
│   ├── call
│   └── indelqual
├── samtools
│   ├── index
│   ├── merge
│   └── stats
├── snpdists
├── snpeff
├── snpsites
├── tbprofiler
│   ├── collate
│   └── vcf_profile
└── utils
    └── gunzip

The inspiration for this task is the way nf-core uses the MultiQC reports to provide a single-file summary stats about various steps.

Following up from #80 (comment)

the refname entry is not transferred from the quanttb output to this file hence the columns/ data are unaligned

Hi guys,

This task is related to defining the default biologically-informed publish directories as opposed to the current tool-specific directories.

This is an independent task and can be accomplished by customizing the results_dir section for any specific process within your pbs.config file.

I've added an example for this here https://github.com/abhi18av/xbs-nf/blob/master/conf/pbs.config#L60

This would be a useful information for retreacability of analysis and results.

process TEST {
    output:
    ...
    path(".*command.{sh,out,err,log}")

    ...
}

And in the nextflow (>= v21.10.0 ) we can specify this in the config file as shown below

params {
    outdir = "${projectDir}/results/"
}




process {
    withName: "TEST" {
        publishDir = [[path: "${params.outdir}/test/", enabled: 'true', pattern: "temp.txt"],
                      [path: "${params.outdir}/test/logs/", enabled: 'true', pattern: "*command.{sh,out,err,log}"]]
    }
}

testing testing.

One major issue we have to tackle before release is the unstable VQSR, which can crash or works sub-optimal given variant annotations wit insufficient variation.
Solution:

1. Run VQSR with all annotations.
2. Extract informativeness of each annotation.
3. Rerun several times each time excluding the weakest annotation.
4. Pick the best model based on 99.9 sensitivity tranche VQSLOD score.

Once #2 is done, it's time to think of the right sub-workflows and the overall workflow.

Hello @everyone (@abhi18av @vrennie @LennertVerboven ).

At the moment we identify multiple infections with QuantTB so that such cases can be excluded from subsequent analyses. Unfortunately we see that this program is not working as well as we would like it to, resulting in false positive multiple infection calls and therefore perfectly good samples excluded from the rest of the analyses.

After some discussion here at Antwerp it has become clear that QuantTB needs to be replaced by a better process, rather than fixing QuantTB which seems a can of worms. The good thing is that we already produce a better output to infer multiple infection, TBprofiler when run on the LoFreq output produces a json with frequencies of the identified (sub)lineages.

Lennert has already written a script to extract the necessary information from this json: https://github.com/TORCH-Consortium/xbs-nf/blob/master/bin/multiple_infection_filter.py

Note that in this particular scenario the multiple infection information is reported by TBprofiler, which in turn is run on the LoFreq output, which in turn is run on the mapped (bam) files.

Additional advantages of this proposal:

• QuantTB is not compatible with many OS's, excluding it would make our pipeline easier to set-up on different systems.
• By the time the pipeline has inferred the multiple infection with TBprofiler all the other stats for selection criteria have also been inferred, it is then possible to report all selection criteria in one summary file (rather than the 4 we have now)
• This in turn makes the whole workflow of the pipeline much cleaner as all samples go through the same flow up to the selection criteria assessment, there is no need for a rejected workflow.
• The output from TBprofiler is much cleaner and easier to interpret: there is clear lineage identification (rather than cryptic QuantTB codes) and associated percentages.
• If this proposal combined with the TBprofiler fix (#104) than both flows can be incorporated at once and solve all major issues for release.

If all stats for selection criteria can be collated and assessed at one point then the best would be to do this before haplotype caller.

@abhi18av could you please have a look at this proposal so that we can discuss?

Tim

Labels feature in Nextflow can be used to simplify the overall management of hardware directives for a profile.

Some good examples are

• nf-core pipelines (which have dynamic computation of hardware resources)
• https://github.com/hoelzer-lab/rnaflow

At the moment we use Coll2018 for excluding DR markers from phylogenetic inference, which is good but not perfect.
https://github.com/TORCH-Consortium/xbs-nf/blob/4db84e6bc4256484ab713206ffd60207a1e4cdb9/workflows/merge_wf.nf#L37
https://github.com/TORCH-Consortium/xbs-nf/blob/4db84e6bc4256484ab713206ffd60207a1e4cdb9/workflows/merge_wf.nf#L54
Ideally we would use the WHO list for this purpose as well, this means WHO would have to be translated to VCF.

We should accomodate the use-case where people only have SRA IDs.

As per the discussion with Lennert, we can use some default values in this case.

As discussed in the TORCH meeting sometimes files are empty or corrupted.

We should add a tool that checks file size and corruption

As discussed in the call earlier today, it's high time we replace the content of https://github.com/abhi18av/xbs-nf/blob/master/modules/utils/sample_stats.nf with a Python script.

Hi @TimHHH and @LennertVerboven ,

Continuing the discussion from #80 (comment) about using custom database in tb-profiler, I am currently testing on Azure/AWS with the custom containers for XBS-nf (https://github.com/abhi18av/xbs-nf/tree/master/containers).

Have you seen this problem before?

Using gff file: /opt/conda/share/tbprofiler/resistance_db_who.gff
Using ref file: /opt/conda/share/tbprofiler/resistance_db_who.fasta
Using barcode file: /opt/conda/share/tbprofiler/resistance_db_who.barcode.bed
Using bed file: /opt/conda/share/tbprofiler/resistance_db_who.bed
Using json_db file: /opt/conda/share/tbprofiler/resistance_db_who.dr.json
Using version file: /opt/conda/share/tbprofiler/resistance_db_who.version.json
Using variables file: /opt/conda/share/tbprofiler/resistance_db_who.variables.json


Can't find /opt/conda/share/tbprofiler/resistance_db_who.variables.json // <--- CORE PROBLEM

The resistance_db_who.variables.json in question is present here

We need to load the TBProfiler database in TBProfiler when creating its virtual environment using the following commands

To create the database from a list of variants
tb-profiler create_db --prefix <new_library_name>

Load the newly created database into TBProfiler
tb-profiler load_library --prefix <new_library_name>

Hi @TimHHH and @LennertVerboven ,

As a result of #29, the current state of pipeline has two stages

• before the sample bams are merged, samples with different libraries

[b5/73822f] process > TEST:QUANTTB_QUANT (EXIT-RIF.FS031-2.L2.A464.1.1.1)        [100%] 13 of 13, cached: 13 ✔
[29/8dd134] process > TEST:MAP_WF:FASTQC (EXIT-RIF.FS031-2.L2.A464.1.1.1)        [100%] 13 of 13, cached: 13 ✔
[d7/425f90] process > TEST:MAP_WF:BWA_MEM (EXIT-RIF.FS031-2.L3.A464.1.1.1)       [100%] 13 of 13, cached: 13 ✔

• after the bams are merged via SAMTOOLE_MERGE

[81/803656] process > TEST:CALL_WF:SAMTOOLS_MERGE (EXIT-RIF.FS031-2)             [100%] 10 of 10, cached: 10 ✔
[c0/b3165c] process > TEST:CALL_WF:GATK_MARK_DUPLICATES (EXIT-RIF.FS007-2)       [100%] 10 of 10, cached: 10 ✔
[be/4095fc] process > TEST:CALL_WF:SAMTOOLS_INDEX (EXIT-RIF.EC090-209)           [100%] 10 of 10, cached: 10 ✔
[4e/775254] process > TEST:CALL_WF:GATK_HAPLOTYPE_CALLER (EXIT-RIF.FS031-2)      [100%] 10 of 10, cached: 10 ✔
[27/e96f03] process > TEST:CALL_WF:LOFREQ_CALL__NTM (EXIT-RIF.EC550A-330)        [100%] 10 of 10, cached: 10 ✔
[19/c59a69] process > TEST:CALL_WF:LOFREQ_INDELQUAL (EXIT-RIF.EC083-202)         [100%] 10 of 10, cached: 10 ✔
[5e/7d8b8d] process > TEST:CALL_WF:SAMTOOLS_INDEX__LOFREQ (EXIT-RIF.EC083-202)   [100%] 10 of 10, cached: 10 ✔
[c9/b90c8c] process > TEST:CALL_WF:LOFREQ_CALL (EXIT-RIF.EC083-202)              [100%] 10 of 10, cached: 10 ✔
[50/db759f] process > TEST:CALL_WF:LOFREQ_FILTER (EXIT-RIF.EC088-207)            [100%] 10 of 10, cached: 10 ✔
[32/edb65b] process > TEST:CALL_WF:DELLY_CALL (EXIT-RIF.FS007-2)                 [100%] 10 of 10, cached: 10 ✔
[d0/e68561] process > TEST:CALL_WF:BCFTOOLS_VIEW (EXIT-RIF.FS031-2)              [100%] 10 of 10, cached: 10 ✔
[3a/26bbd3] process > TEST:CALL_WF:GATK_INDEX_FEATURE_FILE (EXIT-RIF.FS007-2)    [100%] 10 of 10, cached: 10 ✔
[6a/d9b220] process > TEST:CALL_WF:SAMTOOLS_STATS (EXIT-RIF.FS007-2)             [100%] 10 of 10, cached: 10 ✔
[14/76d561] process > TEST:CALL_WF:GATK_COLLECT_WGS_METRICS (EXIT-RIF.EC090-209) [100%] 10 of 10, cached: 10 ✔
[d0/b21681] process > TEST:CALL_WF:GATK_FLAG_STAT (EXIT-RIF.EC550A-330)          [100%] 10 of 10, cached: 10 ✔

The problem that arises is that now the output of QUANTTB_QUANT has a different shape than that of any process after the SAMTOOLS_MERGE - as shown below

# QuantTB
[EXIT-RIF.EC083-202.L1.A254.1.1.1, /home/labbactfiocruz/projects/xbs-nf/work/8b/fc36dbf6a1cd9682e84b1f25d3d787/output/EXIT-RIF.EC083-202.L1.A254.1.1.1.quant.txt]

# WgsMetrics
[EXIT-RIF.EC090-209, /home/labbactfiocruz/projects/xbs-nf/work/14/76d561b1a276bb56ffa272e0af8942/EXIT-RIF.EC090-209.WgsMetrics.txt]

# FlagStat
[EXIT-RIF.EC083-202, /home/labbactfiocruz/projects/xbs-nf/work/0d/2cf006fe1def80c1006f13cf4a5bea/EXIT-RIF.EC083-202.FlagStat.txt]

# SamtoolStats
[EXIT-RIF.EC090-209, /home/labbactfiocruz/projects/xbs-nf/work/72/9508ce671b9a145552de3cf64c33a6/EXIT-RIF.EC090-209.SamtoolStats.txt]

Which means that the code used to merge all of the stats outputs based on sampleID is not working since the

https://github.com/abhi18av/xbs-nf/blob/master/workflows/call_wf.nf#L137-L141

[-        ] process > TEST:CALL_WF:UTILS_SAMPLE_STATS                            -
[-        ] process > TEST:CALL_WF:UTILS_COHORT_STATS                            -

What is the course of action you suggest?

This release relates to #11 and includes the following workflows

1. MAP_WF
2. QUANTTB_QUANT
3. CALL_WF

Once we have understood the behavior of the platform and permissions (on-the-fly conda env) on the cluster, I will give the finishing touches and push the MERGE_WF as well.

With #4 , the basic stubs have been added, but they need to be customized as per each module.

The optimal way to complete the next iteration is to make sure all the process and file stubs are correct and working as per the original design of XBS.

The samples which do not pass the relative abundance filter are eliminated from analysis too soon. They are currently removed from analysis before mapping, they should instead be run through the entire map and call workflows and only be eliminated before merging (same as samples that are eliminated due to insufficient coverage etc.)

I have regenerated the EXIT-RIF expansion set for those working with few or clonal samples.
On our system it can be found here: /home/shared_data_ghi/CWGSMTB/xbs-nf_output/EXIT-RIF_expansion_set/
We'll have to find a place on line or on github where we can deposit this 191MB file, any ideas?

Notes after the meeeting on 12-Jul-2022

• We want to capture the stats of all failing samples (from QUANTTB step) into the stats file (cohort stats).
• Not to let these samples reach MERGE_WF

I finally had a chance to have a closer look at the EXIT-RIF run.
See Tower.nf

Good to see that the the whole thing now runs on pbs, only 10 failed jobs, and I don't think that had anything to do with XBS-nf, but rather a system issue (storage was likely full).

As for the run/results, I have to say that, despite the smallish issues remaining things looks excellent now, the run was surprisingly fast (<2days) and there were no major issues in the general data flow of the pipeline. 👍

The nature of the remaining issues also mean that these can all be tested on smaller datasets, no need to run this big data set again, which is great.

Issues remaining:

1. 🤷 The run ended with a tower communication error, hence tower thinks the run is still ongoing:

Duration    : 1d 20h 22m 26s
CPU hours   : 4'822.2 (3.6% cached, 0% failed)
Succeeded   : 8'393
Cached      : 93
Failed      : 10

WARN: Unexpected HTTP response.
Failed to send message to https://api.tower.nf -- received
- status code : 400
- response msg: Oops... Unable to process request - Error ID: 4Jrmw1EpnkVWYOkEBSx4hx
- error cause : Oops... Unable to process request - Error ID: 4Jrmw1EpnkVWYOkEBSx4hx

2. VQSR, two issues here:

   • ✅ Only 1 Gaussian was used, this should be 4 unless VQSR crashes.
   • 🕵️ VQSR was performed for several annotation combinations but the best one was not selected based on the VQSLOD score (-156.6738 for -an DP -an AS_MQ -an AS_QD -an AS_FS)
   • ✅ *.R.pdf and *.tranches.pdf are not published to the results directory

3. joint.quanttb_cohort_stats.tsv two issues:

   • 🚧 the refname entry is not transferred from the quanttb output to this file hence the columns/ data are unaligned
   • 🚧 like in joint.cohort_stats.tsv there should be a column MULTIPLE_INFECTION_THRESHOLD_MET
   • 🚧 for a sample like quanttb/quant/output/EXIT-RIF.EC104-256.L1.A271.1.1.1.quant.txt there are actually two output lines, both should be reported, rather than one, do make sure that this does not mess up the selection criteria.

4. ✅ TBprofiler is identifying lineage but not Drug Resistance, likely some database issue.

So, these should be very fixable. Let us know who will tackle which issues.
Thanks.

GATK has seen some important updates particularly in light of log4shell, would be good to update our version to the latest 4.2.4.1

• Choose a doc framework (public website)
• Make no assumptions on the background of user
• Provide a script for generating the initial samplesheet

• HPC (PBS @TimHHH @vrennie /SLURM @abhi18av )
• Server (CLIMB <16 GB RAM default> @vrenni ; KHAOS <16 CPUs> @abhi18av )
• Docker (Cloud Azure) @abhi18av
• Workstation (find a machine) @LennertVerboven

@abhi18av to share the instructions for using pipeline URL and updating configs.

Given two samples with the same name before a point can result in an unwanted sample merge, e.g.:
Votintseva2017.614406.m
Votintseva2017.614406.1
In this case both different samples are interpreted as one, namely Votintseva2017.614406
Hence we should make a note in the manual that points are not allowed in the sample sheet.

Lennert is correct, this is for merging multiple libraries from the same sample. This may happen when one library is not too good and a better one also available, which happens some times, but this can also happen when we decide we simply want more/deeper data for a sample, this also happens some times. And there are more scenarios. It is important enough that this step should be implemented now please, we will encounter this issue in the dataset for the paper. Please see my post in the testing section on the impact of library/sample names.
cheers.

Originally posted by @TimHHH in #26 (comment)

Following up from the discussion on #80 (comment)

like in joint.cohort_stats.tsv there should be a column MULTIPLE_INFECTION_THRESHOLD_MET

// GATK_SELECT_VARIANTS__INTERVALS(GATK_INDEX_FEATURE_FILE.out, params.drgenes_list) in the CALL workflow is not yet enabled, is the -L ${intervalsFile} \\ available yet? Or should I create one?

Keep track of which strain has the highest relative abundance and list that strain in the qc summary

https://github.com/TORCH-Consortium/xbs-nf/blob/1ff0abfb4b5f641ba75ecd280ce398df749bfa22/bin/quanttb_qc.py#L114

Since the analysis would be used with small differences, it'd be quite handy to be able to see the resolved configs printed in the console.

I noticed that the quanttb tool (https://github.com/AbeelLab/quanttb) doesn't have a corresponding bioconda package - however it seems to be available in a separate conda channel here https://anaconda.org/jemunro/quanttb

The benefit of bioconda is that the corresponding docker and singularity containers are automatically built and are downloadable from trustable sources.