Execute snakemake for iterative mapping

user_11425932 0
Viewed 85

I would like to use Snakemake for iterative mapping of my samples. I don't know beforehand how many times the samples have to be remapped. For certain samples it will probably be 2-3 times, for others 10 times. If I understand correctly, Snakemake cannot use while loops, but maybe some kind of checkpoint is possible?

Basically what I want to do in this loop is calling the fasta sequence for my Illumina reads each time until it doesn't change anymore. (this is done using bowtie2 > samtools view > samtools mpileup > bcftools call > bcftools view > bcftools index > bcftools consensus)

I wrote previously a bash script to do this, but a Snakemake could really speed up this process. In the bash script, I used a R-script I wrote that counts the number of differences between the new and old fasta file. If this number is = 0, then it can stop the loop, if it is not the same then it has to rerun the steps above. Ideally, it should have a minimum of 3 loops and a max of 15 loops.

If someone could help, that would be fantastic! Thanks in advance

snakemake bioinformatics mapping iteration loops
1 Answers

Since your sequence of steps bowtie2 > samtools view > samtools mpileup ... has to be done serially, i.e. it cannot be parallelized for a given input fastq, I think I would do what you describe in your bash script. Something like, largely pseudocode:

rule all:
    input:
        ['a.consensus.fastq', 'b.consensus.fastq', ...]

rule iterative_mapping:
    input:
        fq= '{sample}.fastq',
    output:
        fq= '{sample}.consensus.fastq',
    run:
        i = 0
        while i < 15 and fastq has changed:
            shell("bowtie2 > samtools view ...")
            shell("check_fastq_has_changed.R ...")

This will do the iterative mapping in parallel for the different {sample}s. Does it make sense?

user_11144530
Related