Snakefile

################## Import libraries ##################

import pandas as pd
import os
import sys
from subprocess import call
import itertools
from snakemake.utils import R


################## Configuration file and PATHS##################

configfile: "config.yaml"

WORKING_DIR         = config["working_dir"]
RESULT_DIR          = config["results_dir"]
DATA_DIR            = config["data_dir"]
GENOME_FASTA_URL    = config["genome_fasta_url"]
GENE_GTF_URL        = config["gene_gtf_url"]

#units = pd.read_table(config["units"], dtype=str).set_index(["bed"], drop=False)

#BED = units.index.get_level_values('bed').unique().tolist()
units = pd.read_table(config["units"], dtype=str).set_index(["sample"], drop=False)

SAMPLES = units.index.get_level_values('sample').unique().tolist()

###############
# Helper Functions
###############
def get_fastq(wildcards):
    return units.loc[(wildcards.samples), ["fq1", "fq2"]].dropna()


##############
# Wildcards
##############
wildcard_constraints:
    sample = "[A-Za-z0-9]+"

wildcard_constraints:
    unit = "[A-Za-z0-9]+


################## DESIRED OUTPUT ##################
# Here we define the outputs of rules we want the pipeline to produce.
# The varialble defined here is then used in the `rule all`.
# Fastqc rule need to be modified to accept any sample
# HMMRATAC will not work with the test dataset.

FASTQC              =       expand(RESULT_DIR     + "fastqc/sample_{numbers}_{R}_fastqc.html",numbers = ['8','12','4_3','4_1'], R=['R1', 'R2']),
FORWARD_READS       =       expand(WORKING_DIR    + "trimmed/{samples}_forward.fastq.gz", samples = SAMPLES),
REVERSE_READS       =       expand(WORKING_DIR    + "trimmed/{samples}_reverse.fastq.gz", samples = SAMPLES),
TRIMMED_FASTQC      =       expand(RESULT_DIR     + "trimmed_fastqc/{samples}_{direction}_fastqc.html", samples = SAMPLES, direction=['forward', 'reverse']),
MAPPED              =       expand(WORKING_DIR    + "mapped/{samples}.bam", samples = SAMPLES),
UNMAPPED            =       expand([WORKING_DIR   + "unmapped/{samples}.fq." + str(i) +".gz" for i in range(1,2)], samples = SAMPLES),
MAP_SORTED          =       expand(WORKING_DIR    + "sort/{samples}.sorted.bam", samples = SAMPLES),
DEDUP               =       expand(WORKING_DIR    + "dedup/{samples}.dedup.bam", samples = SAMPLES),
STATS               =       expand(WORKING_DIR    + "dedup/{samples}.dedup.stats", samples = SAMPLES),
SORTED_INDEXED      =       expand(WORKING_DIR    + "sort/{samples}.sorted.bam.bai", samples = SAMPLES),
GENOME_INFO         =       expand(WORKING_DIR    + "genome_info/{samples}.genome.info", samples  = SAMPLES),
HMMRATAC            =       expand(RESULT_DIR     + "hmmratac/{samples}_peaks.gappedPeak", samples = SAMPLES),
NAME_LOG            =       expand(RESULT_DIR     + "peaks/{samples}.log", samples = SAMPLES),
NARROWPEAK          =       expand(RESULT_DIR     + "macs2/{samples}_peaks.narrowPeak", samples = SAMPLES),
COVERAGE_TRACK      =       expand(RESULT_DIR     + "bamCoverage/{samples}.bw", samples = SAMPLES),
BIGWIGSUMMARY       =        RESULT_DIR     + "bigwigsummary/multiBigwigSummary.npz",
PCAPLOT             =        RESULT_DIR     + "PCA/PCA_PLOT.pdf"
HEATMAP             =       RESULT_DIR + "heatmap/heatmap_reference_point_genes.pdf"


container: "docker://continuumio/miniconda3:4.4.10"

rule all:
    input:
        FASTQC,
        FORWARD_READS,
        REVERSE_READS,
        TRIMMED_FASTQC,
        MAPPED,
        UNMAPPED,
        MAP_SORTED,
        #DEDUP,
        #STATS,
        #COVERAGE_TRACK,
        #BIGWIGSUMMARY,
        #PCAPLOT,
        #SORTED_INDEXED,
        #NARROWPEAK,
        #HMMRATAC,
        HEATMAP

    message : "Analysis is complete!"
    log:
        "log/rule_all.log"
    shell:""

################## INCLUDE RULES ##################
include: "rules/external_data.smk"
include: "rules/pre-processing.smk"
include: "rules/macs2.smk"
include: "rules/deeptools.smk"
include: "rules/hmmratac.smk"