title: Snakemake: How do I use a function that takes in a wildcard and returns a value?
body: <p>I have cram(bam) files that I want to split by read group. This requires reading the header and extracting the read group ids.</p>
  
  <p>I have this function which does that in my Snakemake file:</p>
  
  <pre><code>def identify_read_groups(cram_file):
      import subprocess
      command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
      read_groups = subprocess.check_output(command, shell=True)
      read_groups = read_groups.split('\n')[:-1]
      return(read_groups)
  </code></pre>
  
  <p>I have this rule all:</p>
  
  <pre><code>rule all:
  input:
      expand('cram/RG_bams/{sample}.RG{read_groups}.bam', read_groups=identify_read_groups('cram/{sample}.bam.cram'))
  </code></pre>
  
  <p>And this rule to actually do the split:</p>
  
  <pre><code>rule split_cram_by_rg:
  input:
      cram_file='cram/{sample}.bam.cram',
      read_groups=identify_read_groups('cram/{sample}.bam.cram')
  output:
      'cram/RG_bams/{sample}.RG{read_groups}.bam'
  run:
      import subprocess
      read_groups = open(input.readGroupIDs).readlines()
      read_groups = [str(rg.replace('\n','')) for rg in read_groups]
      for rg in read_groups:
          command = 'samtools view -b -r ' + str(rg) + ' ' + str(input.cram_file) + ' &gt; ' + str(output)
          subprocess.check_output(command, shell=True)
  </code></pre>
  
  <p>I get this error when doing a dry run:</p>
  
  <pre><code>[E::hts_open_format] fail to open file 'cram/{sample}.bam.cram'
  samtools view: failed to open "cram/{sample}.bam.cram" for reading: No such file or directory
  TypeError in line 19 of /gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile:
  a bytes-like object is required, not 'str'
  File "/gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile", line 37, in &lt;module&gt;
    File "/gpfs/gsfs5/users/mcgaugheyd/projects/nei/mcgaughey/EGA_EGAD00001002656/Snakefile", line 19, in identify_read_groups
  </code></pre>
  
  <p>{sample} isn't being passed to the function.</p>
  
  <p>How do I solve this problem? I'm open to other approaches if I'm not doing this in a 'snakemake-ic' way.</p>
  
  <p>==============</p>
  
  <h1><strong>EDIT 1</strong></h1>
  
  <p>Ok, the first set of examples I gave had many many issues.</p>
  
  <p>Here's a better (?) set of code, which I hope demonstrates my issue.</p>
  
  <pre><code>import sys
  from os.path import join
  
  shell.prefix("set -eo pipefail; ")
  
  def identify_read_groups(wildcards):
      import subprocess
      cram_file = 'cram/' + wildcards + '.bam.cram'
      command = 'samtools view -H ' + cram_file + ' | grep ^@RG | cut -f2 | cut -f2 -d":" '
      read_groups = subprocess.check_output(command, shell=True)
      read_groups = read_groups.decode().split('\n')[:-1]
      return(read_groups)
  
  SAMPLES, = glob_wildcards(join('cram/', '{sample}.bam.cram'))
  RG_dict = {}
  for i in SAMPLES:
      RG_dict[i] = identify_read_groups(i)
  
  rule all:
      input:
          expand('{sample}.boo.txt', sample=list(RG_dict.keys()))
  
  rule split_cram_by_rg:
      input:
          file='cram/{sample}.bam.cram',
          RG = lambda wildcards: RG_dict[wildcards.sample]
      output:
          expand('cram/RG_bams/{{sample}}.RG{input_RG}.bam') # I have a problem HERE. How can I get my read groups values applied here? I need to go from one cram to multiple bam files split by RG (see -r in samtools view below). It can't pull the RG from the input.
      shell:
          'samtools view -b -r {input.RG} {input.file} &gt; {output}'
  
  
  rule merge_RG_bams_into_one_bam:
      input:
          rules.split_cram_by_rg.output
      output:
          '{sample}.boo.txt'
      message:
          'echo {input}'
      shell:
          'samtools merge {input} &gt; {output}' #not working
          """
  </code></pre>
  
  <p>==============</p>
  
  <h1><strong>EDIT 2</strong></h1>
  
  <p>Getting MUCH closer, but currently struggling with expand properly building the lane bam files and keeping the wildcards</p>
  
  <p>I'm using this loop to create the intermediate file names:</p>
  
  <pre><code>for sample in SAMPLES:
      for rg_id in list(return_ID(sample)):
          out_rg_bam.append("temp/lane_bam/{}.ID{}.bam".format(sample, rg_id))
  </code></pre>
  
  <p><code>return_ID</code> is a function which takes the sample wildcard and returns a list of the read groups the sample contains </p>
  
  <p>If I use <code>out_rg_bam</code> as an input for a merge rule, then ALL of the files get combined into a merged bam, instead of being split by <code>sample</code>. </p>
  
  <p>If I use <code>expand('temp/realigned/{{sample}}.ID{rg_id}.realigned.bam', sample=SAMPLES, rg_id = return_ID(sample))</code> then rg_id gets applied to each sample. So if I have two samples (a,b) , with read groups (0,1) and (0,1,2), I end up with a0, a1, a0, a1, a2 and b0, b1, b0, b1, b2. </p>