snakemake script access to stdin / stdout for stream processing

user_2325562 0
Viewed 492

For a Snakemake workflow, I need to manipulate tags in many BAM files, and would like to process these by piping them through a script (using the Snakemake script: directive). The specific way I'm doing this is with pysam stream processing.

infile = pysam.AlignmentFile("-", "rb")
outfile = pysam.AlignmentFile("-", "wb", template=infile)

for s in infile:
    (flowcell, lane) = s.query_name.split(':')[0:2]
    rg_id = ".".join([flowcell, lane])
    s.set_tag('RG',rg_id,'Z')
    outfile.write(s)

This script works well standalone, but I haven't been able to figure out how to integrate it via the snakemake script directive. I prefer this way to minimize IO and RAM usage.

Edit: Resorted to direct loading to fix the RG tag.

# parameters passed from snakemake
bam_file = snakemake.input[0]
fixed_bam_file = snakemake.output[0]

bamfile  = pysam.AlignmentFile(bam_file, "rb")
fixed_bamfile = pysam.AlignmentFile(fixed_bam_file, "wb", template = bamfile)

for i, read in enumerate(bamfile.fetch()):
    (flowcell, lane) = read.query_name.split(':')[0:2]
    rg_id = ".".join([flowcell, lane])
    read.set_tag('RG',rg_id,'Z')
    fixed_bamfile.write(read)
    if not (i % 100000):
        print("Updated the read group for {} reads to {}".format(i, rg_id))

bamfile.close()
fixed_bamfile.close()

EDIT: Snakemakes run: and shell: directives set the workdir: directory, while the script: directive operates relative to the directory where the Snakefile was executed (keeping everything nice and tidy). Hence the problem of putting a stream processor under script:.

pysam snakemake
2 Answers

Using shell instead of script directive:

rule all:
    input:
        expand('{sample}_edited.bam'), sample=['a', 'b', 'c']


rule somename:
    input:
        '{sample}.bam'
    output:
        '{sample}_edited.bam'
    shell:
        '''
        cat {input} > python edit_bam.py > {output}
        '''
user_39982520

@Krischan it seems you found a solution already and if so maybe good to post it as an answer.

Alternatively, you can use the object {workflow} to get the directory of the Snakefile and from there construct the path to your python script. If your directory structure is:

./
├── Snakefile
├── data
│   └── sample.bam
└── scripts
    └── edit_bam.py

The Snakefile may look like:

rule all:
    input:
        'test.tmp',

rule one:
    input:
        'sample.bam',
    output:
        'test.tmp',
    shell:
        r"""
        cat {input} \
        | {workflow.basedir}/scripts/edit_bam.py >  {output}
        """

Executed with snakemake -d data ...

It seems the workflow object is not documented but check this thread Any way to get the full path of the Snakefile within the Snakefile?

user_11144530
Related