Fastq input strategy to BALSAMIC

[mathiasbio]

Defining paths used in discussion:

• casefolder ../production/cancer/cases/[caseid]/
• analysisfolder ../production/cancer/cases/[caseid]/analysis
• fastqinputfolder ../production/cancer/cases/[caseid]/fastqs

background

Currently (v11.2.0) fastqs from different lanes are concatenated by CG into the fastqinputfolder such as:
HNN57DRX2_ACCXXXXX_S12_L001_R1_001.fastq.gz
HNN57DRX2_ACCXXXXX_S12_L001_R2_001.fastq.gz
HNN57DRX2_ACCXXXXX_S12_L002_R1_001.fastq.gz
HNN57DRX2_ACCXXXXX_S12_L002_R2_001.fastq.gz
---------->>>
concatenated_ACCXXXXX_XXXXXX_R_1.fastq.gz
concatenated_ACCXXXXX_XXXXXX_R_2.fastq.gz

With the pull request here:
#1090
This concatenation will be removed, and instead the start of the pipeline will be parallelised to run per lane. This means primarily trimming and mapping.

When the concatenation step is removed from CG the intention is to let it instead produce links in the fastqinputfolder to fastq-files in housekeeper. These links are given a new name in CG which differs from the original fastq-name, example:

original: HYMT5DSX3_[LIMS-ID]S297_L001_R1_001.fastq.gz
link: 1_171015_HYMT5DSX3[LIMS-ID]_XXXXXX_R_1.fastq.gz

With the update in this PR: #1069
The path to the fastqfiles in housekeeper will be added to the singularity bind path so that rules in snakemake will have access to them.

the question

What we are discussing now is what the best strategy is for how to input and use this fastq info for the "map per lane" PR #1090

As an example below is a list of fastqlinks for a T / N case in the fastqinputfolder:
1_171015_HYMT5DSX3_[LIMS-ID-TUMOR]XXXXXX_R_1.fastq.gz
2_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_1.fastq.gz
3_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_1.fastq.gz
4_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_1.fastq.gz
1_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_2.fastq.gz
2_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_2.fastq.gz
3_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_2.fastq.gz
4_171015_HYMT5DSX3[LIMS-ID-TUMOR]XXXXXX_R_2.fastq.gz
1_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_1.fastq.gz
2_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_1.fastq.gz
3_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_1.fastq.gz
4_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_1.fastq.gz
1_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_2.fastq.gz
2_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_2.fastq.gz
3_171015_HYMYKDSX3[LIMS-ID-NORMAL]XXXXXX_R_2.fastq.gz
4_171015_HYMYKDSX3[LIMS-ID-NORMAL]_XXXXXX_R_2.fastq.gz

option 1

As it is implemented right now, the path to the fastq-files (fastqinputfolder) is input into BALSAMIC when you configure the case.
This is then included in the [case].json file which the analysis in BALSAMIC is using to configure the workflow.

In the balsamic.smk file a function is imported:

def get_fastq_info(samplename, fastq_dir):
    """
    input:
    output:
    """
    # Fastq suffixes to look for
    fastq_suffixes = {'1': {'fwd': '_1.fastq.gz', 'rev': '_2.fastq.gz'},
                      '2': {'fwd': '_R1_001.fastq.gz', 'rev': '_R2_001.fastq.gz'}}

    fastq_dict = {}
    fastq_dict["fastqpair_patterns"] = {}

    # Prepare Fastq Variables
    for suffix in fastq_suffixes:
        fwd_suffix = fastq_suffixes[suffix]["fwd"]
        rev_suffix = fastq_suffixes[suffix]["rev"]
        fwd_fastqs = glob.glob(f"{fastq_dir}/*{samplename}*{fwd_suffix}")
        if fwd_fastqs:
            for fwd_fastq in fwd_fastqs:
                fastqpair_pattern = os.path.basename(fwd_fastq).replace(fwd_suffix, "")
                if fastqpair_pattern in fastq_dict["fastqpair_patterns"]:
                    LOG.error(f"Fastq name conflict. Fastq pair pattern {fastqpair_pattern} already assigned to "
                              f"dictionary for sample: {samplename}.")
                fastq_dict["fastqpair_patterns"][fastqpair_pattern] = {}
                fastq_dict["fastqpair_patterns"][fastqpair_pattern]["fwd"] = fwd_fastq
                fastq_dict["fastqpair_patterns"][fastqpair_pattern]["rev"] = fwd_fastq.replace(fwd_suffix, rev_suffix)

    return fastq_dict

It takes as input the path to the fastqinputdir and the samplename and returns a dictionary populated with "fastqpatterns" which is the name of the fastq-files with the fastq-suffix and read-number removed, and which holds the paths to the fwd and rev reads (for that sample). For example:
["fastqpair_patterns"][4_171015_HYMYKDSX3_{LIMS-ID-NORMAL}_XXXXXX_R]
["fwd"][fastqpath]
["rev"][fastqpath]

This information is returned to a sample_dict in the balsamic.smk

# Prepare sample_dict
sample_dict = {}
for sample in config["samples"]:
    sample_type = config["samples"][sample]["type"]
    if sample_type == "tumor":
        tumor_sample = sample
        sample_dict[tumor_sample] = get_fastq_info(tumor_sample, fastqinput_dir)
        sample_dict[tumor_sample]["sample_type"] = "TUMOR"
    else:
        normal_sample = sample
        sample_dict[normal_sample] = get_fastq_info(normal_sample, fastqinput_dir)
        sample_dict[normal_sample]["sample_type"] = "NORMAL"

And used in the rules as a fastqpattern wildcard.

Example trimming (start):

rule fastp_umi:
    input:
        fastq_r1 = lambda wildcards: fastq_dict[wildcards.fastqpattern]["fwd"],
        fastq_r2 = lambda wildcards: fastq_dict[wildcards.fastqpattern]["rev"]
    output:
        fastq_r1 = temp(fastq_dir + "{fastqpattern}_1.umi_optimized.fastq.gz"),
        fastq_r2 = temp(fastq_dir + "{fastqpattern}_2.umi_optimized.fastq.gz"),
        json = qc_dir + "fastp/{fastqpattern}_fastp_umi.json",
        html = qc_dir + "fastp/{fastqpattern}_fastp_umi.html",

Example aligning per lane:

rule sentieon_align_sort:
    input:
        ref = config["reference"]["reference_genome"],
        fastq_r1 = Path(fastq_dir, "{fastqpattern}_1.fp.fastq.gz").as_posix(),
        fastq_r2 = Path(fastq_dir, "{fastqpattern}_2.fp.fastq.gz").as_posix(),
        refidx = expand(config["reference"]["reference_genome"] + ".{prefix}", prefix=["amb","ann","bwt","pac","sa"])
    output:
        bamout = Path(bam_dir, "{sample}_align_sort_{fastqpattern}.bam").as_posix()

Example unifying into sample wildcard in dedup step:

rule sentieon_dedup:
    input:
        bamfiles = lambda wildcards: sample_dict[wildcards.sample]["bam"]["align_sort_bamlist"]
    output:
        bam = Path(bam_dir, "{sample}.dedup.bam").as_posix(),

As a side-note, this preparation of the fastqpattern dict can be moved to the previous CONFIG step

option 2

Instead of relying on the fastqinputfolder and sample-name:
CG or BALSAMIC config, will provide in the [case].json configfile to the balsamic.smk a complemented dictionary with information about flowcells and lanes which can be parsed in the snakemake workflow with some hard-coded components with an expected fastq-name.

The [case].json configfile would be modified something like this:

Now it looks like this:

    "samples": {
        "[LIMS-ID]": {
            "type": "tumor"
        },
        "[LIMS-ID]"": {
            "type": "normal"
        }
    },

Into:

    "samples": {
        "[LIMS-ID]": {
            "type": "tumor"
            "flowcells": {"HYMT5DSX3": lanes: [1,2], "HYMT5XXXX": lanes: [1,2,3,4]}
        },
        "[LIMS-ID]"": {
            "type": "normal"
            "flowcells": {"ABCAC5DSX3": lanes: [1,2,3]}
        }
    },

I am at the moment not sure how to implement this in snakemake but I think a strategy could be to create a function like this:

def get_mapping_names(samplename, sample_dict, bam_dir):
    """
    input:
    output:
    """
    alignments = {}
    alignments["align_sort_bamlist"] = []
    for sample in sample_dict[samplename]:
        alignments["align_sort_bamlist"].append(
            bam_dir + "{sample}_align_sort_{flowcell}.{lane}.bam".format(sample=samplename, flowcell=flowcell, lane=lane))
    return alignments

Keep this rule:

rule sentieon_dedup:
    input:
        bamfiles = lambda wildcards: sample_dict[wildcards.sample]["bam"]["align_sort_bamlist"]
    output:
        bam = Path(bam_dir, "{sample}.dedup.bam").as_posix(),

Modify this:

rule sentieon_align_sort:
    input:
        ref = config["reference"]["reference_genome"],
        fastq_r1 = Path(fastq_dir, "{sample}_fc_{flowcell}_l_{lane}_1.fp.fastq.gz").as_posix(),
        fastq_r2 = Path(fastq_dir, "{sample}_fc_{flowcell}_l_{lane}_1.fp.fastq.gz").as_posix(),
        refidx = expand(config["reference"]["reference_genome"] + ".{prefix}", prefix=["amb","ann","bwt","pac","sa"])
    output:
        bamout = Path(bam_dir, "{sample}_align_sort_{flowcell}.{lane}.bam").as_posix()

And this to something like:

rule fastp_umi:
    input:
        fastq_r1 = fastq_input_dir +  {lane}_171015_{flowcell}_{sample}_XXXXXX_R_1.fastq.gz, 
        fastq_r2 = fastq_input_dir +  {lane}_171015_{flowcell}_{sample}_XXXXXX_R_2.fastq.gz
    output:
        fastq_r1 = temp(fastq_dir + "{sample}_fc_{flowcell}_l_{lane}.umi_optimized.fastq.gz"),
        fastq_r2 = temp(fastq_dir + "{sample}_fc_{flowcell}_l_{lane}.umi_optimized.fastq.gz"),
        json = qc_dir + "fastp/{fastqpattern}_fastp_umi.json",
        html = qc_dir + "fastp/{fastqpattern}_fastp_umi.html",

comment

Bare in mind these are examples and either option could be changed code-wise. I think the major decision is if we want to hard-code fastq-names in the rules, or rely on a larger more expansive wildcard (fastqpatterns).

[ivadym]

✅ #1176