# Snakemake workflow for creating a matrix of
# coverage values for genomic features and flanking regions 

# Usage (snakemake --cores should reflect available cores):
# conda env create --file environment.yaml --name ChIPseq_mapping
# source activate ChIPseq_mapping
# snakemake -p --cores 48
# source deactivate

import os

# To make the per_base_coverage rule work with a shell script invoked using the "shell" directive,
# we need to determine the base path of Snakefile since we expect the scripts directory to be there as well
SRCDIR = srcdir("")

# Specify config file parameters
configfile: "config.yaml"
sample               = config["SAMPLE"]
superfamFeaturesName = config["FEATURES"]["superfamName"]
reference            = config["MAPPING"]["reference"]
refbase              = os.path.basename(reference)
mode                 = config["MAPPING"]["mode"]
flankName            = config["COVERAGE"]["flankName"]
binName              = config["COVERAGE"]["binName"]

# Specify the desired end target file(s)
rule all:
    input:
        #expand("matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_matrix_bin{binName}_flank{flankName}.tab",
        #       sample = sample,
        #       refbase = refbase,
        #       mode = mode,
        #       binName = binName,
        #       flankName = flankName),
        #expand("matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_ranLoc_matrix_bin{binName}_flank{flankName}.tab",
        #       sample = sample,
        #       refbase = refbase,
        #       mode = mode,
        #       binName = binName,
        #       flankName = flankName),
        expand("matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_matrix_bin{binName}_flank{flankName}.tab",
               sample = sample,
               refbase = refbase,
               mode = mode,
               superfamFeaturesName = superfamFeaturesName,
               binName = binName,
               flankName = flankName),
        expand("matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_ranLoc_matrix_bin{binName}_flank{flankName}.tab",
               sample = sample,
               refbase = refbase,
               mode = mode,
               superfamFeaturesName = superfamFeaturesName,
               binName = binName,
               flankName = flankName)

# Use computeMatrix function from deepTools to create a matrix of
# coverage values for genomic features and flanking regions, and
# for equivalent random loci and flanking regions
rule compute_matrix_superfamFeatures:
    """Create matrices of coverage values for genomic features and random loci"""
    input:
        featuresFile = "/home/ajt200/analysis/wheat/featureProfiles/TEs/superfamilies/iwgsc_refseqv1.0_TransposableElements_2017Mar13_superfamily_{superfamFeaturesName}.bed",
        ranLocFile   = "/home/ajt200/analysis/wheat/featureProfiles/TEs/superfamilies/iwgsc_refseqv1.0_TransposableElements_2017Mar13_superfamily_{superfamFeaturesName}_randomLoci.bed",
        BW           = "../{mode}/bw/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm.bw"
    output:
        featuresGZ  = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_matrix_bin{binName}_flank{flankName}.gz",
        featuresTAB = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_matrix_bin{binName}_flank{flankName}.tab",
        ranLocGZ  = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_ranLoc_matrix_bin{binName}_flank{flankName}.gz",
        ranLocTAB = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_ranLoc_matrix_bin{binName}_flank{flankName}.tab"
    params:
        bodyLength  = config["COVERAGE"]["bodyLength"],
        startLabel  = config["COVERAGE"]["startLabel"],
        endLabel    = config["COVERAGE"]["endLabel"],
        upstream    = config["COVERAGE"]["upstream"],
        downstream  = config["COVERAGE"]["downstream"],
        binSize     = config["COVERAGE"]["binSize"],
        sortRegions = config["COVERAGE"]["sortRegions"]
    log:
        features = "logs/computeMatrix/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_matrix_bin{binName}_flank{flankName}.log",
        ranLoc   = "logs/computeMatrix/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_{superfamFeaturesName}_ranLoc_matrix_bin{binName}_flank{flankName}.log"
    threads:
        config["THREADS"]  
    shell:
        "(computeMatrix scale-regions"
        " --regionsFileName {input.featuresFile}"
        " --scoreFileName {input.BW}"
        " --outFileName {output.featuresGZ}"
        " --outFileNameMatrix {output.featuresTAB}"
        " --regionBodyLength {params.bodyLength}"
        " --startLabel {params.startLabel}"
        " --endLabel {params.endLabel}"
        " --upstream {params.upstream}"
        " --downstream {params.downstream}"
        " --binSize {params.binSize}"
        " --sortRegions {params.sortRegions}"
        " --numberOfProcessors {threads}) 2> {log.features}; "
        "(computeMatrix scale-regions"
        " --regionsFileName {input.ranLocFile}"
        " --scoreFileName {input.BW}"
        " --outFileName {output.ranLocGZ}"
        " --outFileNameMatrix {output.ranLocTAB}"
        " --regionBodyLength {params.bodyLength}"
        " --startLabel {params.startLabel}"
        " --endLabel {params.endLabel}"
        " --upstream {params.upstream}"
        " --downstream {params.downstream}"
        " --binSize {params.binSize}"
        " --sortRegions {params.sortRegions}"
        " --numberOfProcessors {threads}) 2> {log.ranLoc}"

# Use computeMatrix function from deepTools to create a matrix of
# coverage values for genomic features and flanking regions, and
# for equivalent random loci and flanking regions
rule compute_matrix_features:
    """Create matrices of coverage values for genomic features and random loci"""
    input:
        featuresFile = "/home/ajt200/analysis/wheat/featureProfiles/TEs/iwgsc_refseqv1.0_TransposableElements_2017Mar13.bed",
        ranLocFile   = "/home/ajt200/analysis/wheat/featureProfiles/TEs/iwgsc_refseqv1.0_TransposableElements_2017Mar13_randomLoci.bed",
        BW           = "../{mode}/bw/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm.bw"
    output:
        featuresGZ  = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_matrix_bin{binName}_flank{flankName}.gz",
        featuresTAB = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_matrix_bin{binName}_flank{flankName}.tab",
        ranLocGZ   = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_ranLoc_matrix_bin{binName}_flank{flankName}.gz",
        ranLocTAB  = "matrices/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_ranLoc_matrix_bin{binName}_flank{flankName}.tab"
    params:
        bodyLength  = config["COVERAGE"]["bodyLength"],
        startLabel  = config["COVERAGE"]["startLabel"],
        endLabel    = config["COVERAGE"]["endLabel"],
        upstream    = config["COVERAGE"]["upstream"],
        downstream  = config["COVERAGE"]["downstream"],
        binSize     = config["COVERAGE"]["binSize"],
        sortRegions = config["COVERAGE"]["sortRegions"]
    log:
        features = "logs/computeMatrix/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_matrix_bin{binName}_flank{flankName}.log",
        ranLoc   = "logs/computeMatrix/{sample}_MappedOn_{refbase}_lowXM_{mode}_sort_norm_allTEs_ranLoc_matrix_bin{binName}_flank{flankName}.log"
    threads:
        config["THREADS"]  
    shell:
        "(computeMatrix scale-regions"
        " --regionsFileName {input.featuresFile}"
        " --scoreFileName {input.BW}"
        " --outFileName {output.featuresGZ}"
        " --outFileNameMatrix {output.featuresTAB}"
        " --regionBodyLength {params.bodyLength}"
        " --startLabel {params.startLabel}"
        " --endLabel {params.endLabel}"
        " --upstream {params.upstream}"
        " --downstream {params.downstream}"
        " --binSize {params.binSize}"
        " --sortRegions {params.sortRegions}"
        " --numberOfProcessors {threads}) 2> {log.features}; "
        "(computeMatrix scale-regions"
        " --regionsFileName {input.ranLocFile}"
        " --scoreFileName {input.BW}"
        " --outFileName {output.ranLocGZ}"
        " --outFileNameMatrix {output.ranLocTAB}"
        " --regionBodyLength {params.bodyLength}"
        " --startLabel {params.startLabel}"
        " --endLabel {params.endLabel}"
        " --upstream {params.upstream}"
        " --downstream {params.downstream}"
        " --binSize {params.binSize}"
        " --sortRegions {params.sortRegions}"
        " --numberOfProcessors {threads}) 2> {log.ranLoc}"
