Snakemake with custom script to split bedgraph by contig

user_15613766 0
Viewed 88

I'm quite new to bioinformatics and Snakemake both, but I'm trying to put together an automated pipeline for Tn-seq data analysis.

I've written a script that reads in a .bedgraph file and outputs separate files for each contig, as I'd like to analyze each contig separately. I've written the code to output a file with a basename of the input file + the contig name:

input_handle = FILE
path = PATH

import csv
import re

contigs = {}

with open(input_handle) as data:
    data_reader = csv.reader(data, delimiter='\t')
    contigs = {row[0] for row in data_reader}

for c in contigs:
    with open(input_handle) as data:
        data_reader = csv.reader(data, delimiter='\t')
        out_file = path + re.search(r".+\/(.+)(?=.bedgraph)", input_handle).group(1) + "-" + c + ".bedgraph"
        f_out = open(out_file, 'w')
        for row in data_reader:
            if row[0] == c:
                f_out.write("\t".join(row)+"\n")

I'm struggling to figure out how to incorporate this into Snakemake appropriately. I'm also pretty unclear on how to then incorporate that output within my script.

EDIT: I had previously said I thought I should be using dynamic but after looking around more it appears that is deprecated.

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
---->   "split_beds/{sample}/?????"
    script:
        "scripts/split_bed.py"

calls:

input_bedgraph = snakemake.input[0]

import csv
import re

contigs = {}

with open(input_bedgraph) as data:
    data_reader = csv.reader(data, delimiter='\t')
    contigs = {row[0] for row in data_reader}

for c in contigs:
    with open(input_bedgraph) as data:
        data_reader = csv.reader(data, delimiter='\t')
---->   out_file = snakemake.output[0]
        f_out = open(out_file, 'w')
        for row in data_reader:
            if row[0] == c:
                f_out.write("\t".join(row)+"\n")

If anyone could point my in the right direction it would be much appreciated! The tutorial was great for getting started and I've got plenty of rules that work well before this, but I'm a bit lost at this point.

EDIT2:

I've confirmed that simply using the following with an "all" rule expanding the split_beds folder works to make a single file with a single contig from each, so the script works fine... Just need to be able to do multiple output into different folders...

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        "split_beds/{sample}.bedgraph"
    script:
        "scripts/split_bed.py"
snakemake bioinformatics python
2 Answers

After getting some help from someone who knows what they're doing the problem became a bit easier to solve.

I did end up using checkpoints (thanks @Luigi). I believe the problem was probably simplified by using an awk one liner instead of a script but using the code below along with a modified all rule ended up getting me all of the files that I needed.

checkpoint split_bed:
    input:
        "bedgraphs/{sample}.bedgraph"
    output:
        directory("split_beds/{sample}/")
    shell:
        "mkdir split_beds/{wildcards.sample}; awk \'{{print $0 > \"split_beds/{wildcards.sample}/{wildcards.sample}_\"$1\".bedgraph\"}}\' {input}"

def aggregate_input(wildcards):
    checkpoint_output = checkpoints.split_bed.get(**wildcards).output[0]
    return expand("split_beds/{sample}/{contig}.bedgraph",
           sample=wildcards.sample,
           contig=glob_wildcards(os.zpath.join(checkpoint_output, "{contig}.bedgraph")).contig)

rule aggregate:
    input:
        aggregate_input
    output:
        "split_beds/{sample}/{contig}"

And the "all" rule is:

rule all:
    input:
        expand("split_beds/{sample}/", sample=config["samples"])
user_156137660

Do you need to use the contigs downstream in your snakemake pipeline? If so, dynamic() is definitely the way to go. This related question is a very accessible intro to dynamic(). It'd potentially look something like this:

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        dynamic("split_beds/{sample}/{contig}.fasta")
    script:
        "scripts/split_bed.py"

rule process_contigs:
    input:
        dynamic("split_beds/{sample}/{contig}.fasta")
    ...

If my understanding is correct, snakemake understands this as "process_contigs needs some number of files with names matching the pattern 'split_beds/{sample}/{contig}.fasta', which will be produced by split_bed".

Alternately, if you just want to splort out your split contigs and don't need those files processed individually within snakemake, a less elegant way is to just touch a file once you're done processing the contigs, then use that as the output for the snakemake rule (using pathlib).

Your sample processing code:

input_handle = FILE
path = PATH

...

from pathlib import Path

...

for c in contigs:

    ...

Path('split_beds/status/{sample}_completed.txt').touch()

Your rule:

rule split_bed:
    input:
        "bam_coverage/{sample}.bedgraph"
    output:
        "split_beds/status/{sample}_completed.txt"
    script:
        "scripts/split_bed.py"

I would only really do the second one if the split contigs are a terminal point of your pipeline and you don't want to bother with making another rule to take the dynamic contig files as input.

Edit:

I've just remembered that checkpoints exist. I think these replaced dynamic()--checkpoints basically allow snakemake to re-evaluate the DAG at certain points of the workflow, which I think is what you need, since you'll be making a variable number of individual contig files for each sample. Thus, the DAG can't be finalized at the beginning of the run (I think this is what checkpoints are designed to solve). I haven't used them myself so I'm not going to mislead you by trying to psuedocode an example, but hopefully this points you in the right direction!

user_34403020
Related