Input function for technical / biological replicates in snakemake

Viewed 27
I´m currently trying to write a Snakemake workflow that can check automatically via a sample.tsv file if a given sample is a biological or technical replicate. And then use in this case at some point of my workflow a rule to merge technical / biological replicates.

My tsv file looks like this:

sample unit_bio unit_tech fq1 fq2
bCalAnn1 1 1 /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_1_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_1_R2.fastq.gz
bCalAnn1 1 2 /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_2_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_2_R2.fastq.gz
bCalAnn2 1 1 /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_1_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_1_R2.fastq.gz
bCalAnn2 1 2 /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_2_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_2_R2.fastq.gz
bCalAnn2 2 1 /home/assembly_downstream/data/arima_HiC/bCalAnn2_2_1_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn2_2_1_R2.fastq.gz
bCalAnn2 3 1 /home/assembly_downstream/data/arima_HiC/bCalAnn2_3_1_R1.fastq.gz /home/assembly_downstream/data/arima_HiC/bCalAnn2_3_1_R2.fastq.gz

My Pipeline looks like this:

import pandas as pd
import os
import yaml

configfile: "config.yaml"


samples = pd.read_table(config["samples"], dtype=str)

rule all:
    input:
        expand(config["arima_mapping"] + "final/{sample}_{unit_bio}_{unit_tech}.bam", zip, 
            sample=samples["sample"], unit_bio=samples["unit_bio"], unit_tech=samples["unit_tech"])



..
some rules
..

rule add_read_groups:
    input:
        config["arima_mapping"] + "paired/{sample}_{unit_bio}_{unit_tech}.bam"
    output:
        config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam"
    params:
        platform = "ILLUMINA",
        sampleName = "{sample}",
        library = "{sample}",
        platform_unit ="None"
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/paired_read_groups/{sample}_{unit_bio}_{unit_tech}.log"
    shell:
        "picard AddOrReplaceReadGroups I={input} O={output} SM={params.sampleName} LB={params.library} PU={params.platform_unit} PL={params.platform} 2> {log}"

rule merge_tech_repl:
    input:
        config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam" 
    output:
        config["arima_mapping"] + "merge_tech_repl/{sample}_{unit_bio}_{unit_tech}.bam"
    params:
        val_string = "SILENT" 
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/merged_tech_repl/{sample}_{unit_bio}_{unit_tech}.log"
    threads:
        2 #verwendet nur maximal 2
    shell:
        "picard MergeSamFiles -I {input} -O {output} --ASSUME_SORTED true --USE_THREADING true --VALIDATION_STRINGENCY {params.val_string} 2> {log}"


rule mark_duplicates:
    input:
        config["arima_mapping"] + "merge_tech_repl/{sample}_{unit_bio}_{unit_tech}.bam" if config["tech_repl"] else config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam"
    output:
        bam = config["arima_mapping"] + "final/{sample}_{unit_bio}_{unit_tech}.bam",
        metric = config["arima_mapping"] + "final/metric_{sample}_{unit_bio}_{unit_tech}.txt"
    #params:
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/mark_duplicates/{sample}_{unit_bio}_{unit_tech}.log"
    shell:
        "picard MarkDuplicates I={input} O={output.bam} M={output.metric} 2> {log}"

At the moment I have set a boolean in a config file that tells the mark_duplicates rule whether to take its input from the add_read_group or the merge_technical_replicates rule. This is of course not optimal since it could be that some samples may have duplicates (of any numbers) while others don´t. Therefore I want to create a syntax that checks the tsv table if a given sample name and unit_bio number are identical while the unit_tech number is different (and later analog to this for biological replicates), thus merging these specific samples while nonduplicate samples skip the merging rule.

Reading up on the snakemake documentation I think the proper way to do this would be via an input function. According to the docs the base structure would look like this:

def check_tech_repl(wildcards):
    return[... a list of input files depending on given wildcards ...]

After some attempts my main problem I struggle with is the table checking. As far as I know snakemake takes templates as input and checks for all samples that match this. But I would need to check the table for every sample that shares (e.g. for technical replicate) the sample name and the unit_bio number take all these samples and give them as input for the first rule run. Then I would have to take the next sample which was not already part of a previous run to prevent merging the same samples multiple times.

I guess I will just have to write a script for that but I was wondering if there is a snakemake way of solving this?

1 Answers

Therefore I want to create a syntax that checks the tsv table if a given sample name and unit_bio number are identical while the unit_tech number is different (and later analog to this for biological replicates), thus merging these specific samples while nonduplicate samples skip the merging rule.

rule gather_techdups_of_a_biodup:
    output: "{sample}/{unit_bio}"
    input: gather_techdups_of_a_biodup_input_fn
    shell: "true"  # Fill this in

rule gather_biodips_of_a_techdup:
    output: "{sample}/{unit_tech}"
    input: gather_biodips_of_a_techdup_input_fn
    shell: "true"  # Fill this in

After some attempts my main problem I struggle with is the table checking. As far as I know snakemake takes templates as input and checks for all samples that match this. But I would need to check the table for every sample that shares (e.g. for technical replicate) the sample name and the unit_bio number take all these samples and give them as input for the first rule run. Then I would have to take the next sample which was not already part of a previous run to prevent merging the same samples multiple times.

The logic you describe here can be implemented in the gather_techdups_of_a_biodup_input_fn and gather_biodips_of_a_techdup_input_fn functions above. For example, read your sample TSV file with pandas, filter for wildcards.sample and wildcards.unit_bio (or wildcards.unit_tech), then extract columns fq1 and fq2 from the filtered data frame.

Related