SAMPLES=["lodi", "branching", "e1-101", "ranch9", "cstrain"]

rule all:
    input: 
       expand("{sample}_R1.fastq.gz", sample=SAMPLES),
       expand("{sample}_R2.fastq.gz", sample=SAMPLES),
       expand("{sample}_peaks_w_mt.txt", sample=SAMPLES),
       expand("{sample}_peaks_wo_mt.txt", sample=SAMPLES)

##download_genome: download the reference mitochondrial genome using efetch
rule download_genome:
   conda: "env.yml"    # conda environment to use
   output:
      "MT880588.fasta" # output fasta file
   shell: """
      efetch             \
         -db nuccore     \
         -id  MT880588.1 \
         -format fasta > {output}
   """


##download_reads: downloads fastq files for different strains. URLs were obtained with https://sra-explorer.info/.
# for a better version for this rule, see http://ivory.idyll.org/blog/2020-snakemake-hacks-collections-files.html
rule download_reads:
   output: 
      lodir1="lodi_R1.fastq.gz", 
      lodir2="lodi_R2.fastq.gz",
      ranch9r1="ranch9_R1.fastq.gz", 
      ranch9r2="ranch9_R2.fastq.gz",
      e1_101r1="e1-101_R1.fastq.gz",
      e1_101r2="e1-101_R2.fastq.gz",
      branchingr1="branching_R1.fastq.gz", 
      branchingr2="branching_R2.fastq.gz",
      cstrainr1="cstrain_R1.fastq.gz",
      cstrainr2="cstrain_R2.fastq.gz",
   shell: """
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/000/SRR1448470/SRR1448470_1.fastq.gz -o {output.lodir1}
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/000/SRR1448470/SRR1448470_2.fastq.gz -o {output.lodir2}

      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/008/SRR1448468/SRR1448468_1.fastq.gz -o {output.e1_101r1}
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/008/SRR1448468/SRR1448468_2.fastq.gz -o {output.e1_101r2}

      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/004/SRR1448454/SRR1448454_1.fastq.gz -o {output.ranch9r1}
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/004/SRR1448454/SRR1448454_2.fastq.gz -o {output.ranch9r2}

      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/003/SRR1448453/SRR1448453_1.fastq.gz -o {output.branchingr1}
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/003/SRR1448453/SRR1448453_2.fastq.gz -o {output.branchingr2}

      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/000/SRR1448450/SRR1448450_1.fastq.gz -o {output.cstrainr1}
      curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR144/000/SRR1448450/SRR1448450_2.fastq.gz -o {output.cstrainr2}
    """


##filter_mtDNA: remove reads with k-mer matching the mitochondrial genome
rule filter_mtDNA:
  conda: "env.yml"    # conda environment to use
  input:
     r1="{sample}_R1.fastq.gz",    # left reads of the pairs
     r2="{sample}_R2.fastq.gz",    # right reads of the pairs
     mt="MT880588.fasta"           # mitochondrial genome downloaded with rule download_genome
  output:
     r1="{sample}_mtless_R1.fastq.gz",    # left reads of the pairs that do not have a k-mer matching MT880588.fasta
     r2="{sample}_mtless_R2.fastq.gz"     # right reads of the pairs that do not have a k-mer matching MT880588.fasta
  params:
     kmer=29    # k-mer size to use
     cpus=4     # how many cpus to use              
  shell: """
     bbduk.sh            \
        in1={input.r1}   \
        in2={input.r2}   \
        out1={output.r1} \
        out2={output.r2} \
        ref={input.mt}   \
        k={params.kmer}  \
        threads={params.cpus}
  """


##estimate_all: estimate nuclear genome size using all reads
rule estimate_all:
  conda: "env.yml"    # conda environment to use
  input: 
     r1="{sample}_R1.fastq.gz",    # left reads of the pairs
     r2="{sample}_R2.fastq.gz"     # right reads of the pairs
  output:
     "{sample}_peaks_w_mt.txt"     # output log with estimated genome size
  params:
     kmer=29    # k-mer size to use
     cpus=4     # how many cpus to use 
  shell: """
     kmercountexact.sh  \
        in1={input.r1}  \
        in2={input.r2}  \
        k=29            \
        k={params.kmer} \
        threads={params.cpus}
  """


##estimate_filt: estimate genome size using filtered reads, after removing those with a k-mer matching the mitochondrial genome
rule estimate_filt:
  conda: "env.yml"    # conda environment to use
  input:
     r1="{sample}_mtless_R1.fastq.gz",    # left reads of the pairs that do not have a k-mer matching MT880588.fasta
     r2="{sample}_mtless_R2.fastq.gz"     # right reads of the pairs that do not have a k-mer matching MT880588.fasta
  output:
     "{sample}_peaks_wo_mt.txt"     # output log with estimated genome size
  params:
     kmer=29    # k-mer size to use
     cpus=4     # how many cpus to use 
  shell: """
     kmercountexact.sh  \
        in1={input.r1}  \
        in2={input.r2}  \
        peaks={output}  \
        k={params.kmer} \
        threads={params.cpus}
  """


