More adaptable mapping wrapper.

PRO-seq/PRO-seq.snakefile

@@ -55,6 +55,8 @@ from libworkflows import read_htseq_counts, sum_htseq_counts
 from libworkflows import read_feature_counts, sum_feature_counts
 from smincludes import rules as irules
 
+alignment_settings = {"bowtie2": ""},
+
 #TRIMMERS = ["cutadapt", "fastx_clipper"]
 TRIMMERS = ["cutadapt"]
 COUNTERS = ["feature_count"]
@@ -312,6 +314,7 @@ rule map_on_genome:
         sam = temp(OPJ(output_dir, "{trimmer}", aligner, "mapped_C_elegans", "{lib}_{rep}_{type}_on_C_elegans.sam")),
         nomap_fastq = OPJ(output_dir, "{trimmer}", aligner, "not_mapped_C_elegans", "{lib}_{rep}_{type}_unmapped_on_C_elegans.fastq.gz"),
     params:
+        aligner = aligner,
         index = index,
         settings = "",
     message:

small_RNA-seq/small_RNA-seq.snakefile

@@ -129,6 +129,7 @@ from libworkflows import filter_combinator, sum_feature_counts, sum_htseq_counts
 from smincludes import rules as irules
 strip = str.strip
 
+alignment_settings = {"bowtie2": "-L 6 -i S,1,0.8 -N 0"},
 
 # Positions in small RNA sequences for which to analyse nucleotide distribution
 #POSITIONS = ["first", "last"]
@@ -733,8 +734,9 @@ rule map_on_genome:
         sam = temp(OPJ(output_dir, aligner, "mapped_C_elegans", "{lib}_{rep}", "%s_on_C_elegans.sam" % size_selected)),
         nomap_fastq = OPJ(output_dir, aligner, "not_mapped_C_elegans", "{lib}_{rep}_%s_unmapped_on_C_elegans.fastq.gz" % size_selected),
     params:
+        aligner = aligner,
         index = index,
-        settings = "-L 6 -i S,1,0.8 -N 0",
+        settings = alignment_settings[aligner],
     message:
         "Mapping {wildcards.lib}_{wildcards.rep}_%s on C. elegans genome." % size_selected
     benchmark:

snakemake_wrappers/map_on_genome/wrapper.py

 from snakemake.shell import shell
 
-cmd = """
-genome_dir="${{HOME}}/Genomes"
-genome="C_elegans"
-cmd="bowtie2 --seed 123 -t {snakemake.params.settings} --mm -x {snakemake.params.index} -U {snakemake.input.fastq} --no-unal --un-gz {snakemake.output.nomap_fastq} -S {snakemake.output.sam}"
-        echo ${{cmd}} > {snakemake.log.log}
-        eval ${{cmd}} 1>> {snakemake.log.log} 2> {snakemake.log.err}
+def mapping_command(aligner):
+    """This function returns the shell commands to run given the *aligner*."""
+    if aligner == "hisat2":
+        shell_commands = """
+cmd="hisat2 -p {snakemake.threads} --dta --seed 123 -t {snakemake.params.settings} --mm -x {snakemake.params.index} -U {snakemake.input.fastq} --no-unal --un-gz {snakemake.output.nomap_fastq} -S {snakemake.output.sam}"
+echo ${{cmd}} 1> {log.log}
+eval ${{cmd}} 1>> {log.log} 2> {log.err}
 """
+    elif aligner == "bowtie2":
+        shell_commands = """
+cmd="bowtie2 -p {snakemake.threads} --seed 123 -t {snakemake.params.settings} --mm -x {snakemake.params.index} -U {snakemake.input.fastq} --no-unal --un-gz {snakemake.output.nomap_fastq} -S {snakemake.output.sam}"
+echo ${{cmd}} > {snakemake.log.log}
+eval ${{cmd}} 1>> {snakemake.log.log} 2> {snakemake.log.err}
+"""
+    elif aligner == "crac":
+        shell_commands = """
+cmd="crac --nb-threads {snakemake.threads} --summary %s --all %s {snakemake.params.settings} -i {snakemake.params.index} -r {snakemake.input.fastq} --sam {snakemake.output.sam}"
+echo ${{cmd}} > {snakemake.log.log}
+eval ${{cmd}} 1>> {snakemake.log.log} 2> {snakemake.log.err}
+# TODO: extract non mappers from the sam output
+> {snakemake.output.nomap_fastq}
+""" % (snakemake.output.nomap_fastq.split("_unmapped")[0] + "_summary.txt", snakemake.output.nomap_fastq.split("_unmapped")[0] + "_output.txt")
+    else:
+        raise NotImplementedError("%s is not yet handled." % aligner)
+    return shell_commands
 
-shell(cmd)
+shell(mapping_command(snakemake.params.aligner))