script:
    "../scripts/do_UMAP.py"

script:
    "../scripts/do_PCA.py"

script:
    "../scripts/plot_dimred.py"

script:
    "../scripts/split_data.py"

script:
    "../scripts/filter_regions.py"

script:
    "../scripts/normalize_TMM.py"

script:
    "../scripts/normalize_CQN.py"

script:
    "../scripts/select_HVR.py"

script:
    "../scripts/plot_mean_var.py"

script:
    "../scripts/get_consensus_regions.py"

script:
    "../scripts/quantify_support_sample.py"

script:
    "../scripts/quantify_support_aggregate.py"

script:
    "../scripts/quantify_counts_sample.py"

script:
    "../scripts/quantify_counts_aggregate.py"

shell:
    """
    uropa -p {params.results_dir}/gencode -i {input.gencode_config} -t {threads} -l {params.results_dir}/uropa.gencode.log
    uropa -p {params.results_dir}/reg -i {input.reg_config} -t {threads} -l {params.results_dir}/uropa.reg.log
    """

shell:
    """
    export PATH="{params.homer_bin}:$PATH";

    resources/homer/bin/annotatePeaks.pl {input.consensus_regions} {config[genome]} \
        > {output.homer_annotations} \
        2> {output.homer_annotations_log};
    """ 

run:
    # load and format uropa gencode results
    gencode_characterization=pd.read_csv(input.gencode_results,sep='\t')
    gencode_characterization = gencode_characterization.set_index("peak_id")
    gencode_characterization.loc[gencode_characterization['feature']=='transcript','feat_type']='transcript:'+gencode_characterization.loc[gencode_characterization['feature']=='transcript','transcript_type']
    gencode_characterization.loc[gencode_characterization['feature']=='gene','feat_type']='gene:'+gencode_characterization.loc[gencode_characterization['feature']=='gene','gene_type']
    gencode_characterization['length']=gencode_characterization['peak_end']-gencode_characterization['peak_start']
    gencode_characterization=gencode_characterization[['peak_chr','peak_start','peak_end','length','feat_anchor','distance','relative_location','feat_type','gene_id','gene_name','name']]
    gencode_characterization.columns=['chr','start','end','length','feat_anchor','distance','location','feat_type','gene_id','gene_name','characterization']
    gencode_characterization.loc[gencode_characterization['characterization'].isna(),'characterization']='NONE'
    gencode_characterization=gencode_characterization.add_prefix('gencode_')

    # load and format uropa regulatory build results
    reg_characterization=pd.read_csv(input.reg_results,sep='\t')
    reg_characterization = reg_characterization.set_index('peak_id')[['feature','ID']]
    reg_characterization.columns=['reg_feature','reg_feature_id']
    reg_characterization.loc[reg_characterization['reg_feature'].isna(),'reg_feature']='reg_NONE'
    reg_characterization=reg_characterization.add_prefix('regulatoryBuild_')

import os
import pandas as pd
import pybedtools as bedtools

# configuration
sloop_extension=snakemake.params["sloop_extension"]
output = snakemake.output[0]
blacklist_file= snakemake.params["blacklist_file"]
chrom_file = snakemake.params["chrom_file"]
results_dir = snakemake.params["results_dir"]
# load annotation file and get sample names
annotations = pd.read_csv(snakemake.input[0], index_col=0)
annotations=annotations[(annotations['pass_qc']>0)]
# save filtered annotation file
annotations.to_csv(snakemake.output[1])

peakfiles = [os.path.join(results_dir,"{}".format(sample),"peaks", "{}_summits.bed".format(sample)) for sample in annotations.index]

output_bed = None

if not os.path.exists(os.path.split(output)[0]):
    os.mkdir(os.path.split(output)[0])

for peakfile in peakfiles:
    peak_bed = bedtools.BedTool(peakfile)
    if (blacklist_file is not None):
        peak_bed=peak_bed.intersect(blacklist_file,v=True, wa=True)

    peak_bed = peak_bed.slop(g=chrom_file, b=sloop_extension)

    if (output_bed is None):
        output_bed = peak_bed
    else:
        output_bed = output_bed.cat(peak_bed,force_truncate=True)

output_bed.saveas(output)

peaks = bedtools.BedTool(output).sort(faidx=chrom_file).to_dataframe(names=['CHR','START','END'],dtype={'START':int,'END':int})
peaks['ID'] = peaks.index.format(formatter=(lambda x: "CONS{:011d}".format(x)))
bedtools.BedTool().from_dataframe(peaks).saveas(output)

import pybedtools
import os
import time
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# utility functions
def calc_norm_factors(counts_df, method):
    assert method in ['TMM', 'TMMwsp', 'RLE', 'upperquartile']
    from rpy2.robjects import numpy2ri, pandas2ri, r
    numpy2ri.activate()
    pandas2ri.activate()
    r.source(os.path.join('workflow','scripts','utils_calcNormFactors.R'))
    return r.calcNormFactors(counts_df, method=method)

def cpm(counts_matrix, feature_len=1, norm_factors=1, log=False, pseudocount=0.5, libsize_pseudocount=1, million=1e6):
    """
    Counts per million normalization
    If feature_len=1 then you get counts per million
    Setting pseudocount=0.5, libsize_pseudocount=1 ensures results equivalent to LIMMA-voom
    """
    rpk = counts_matrix.astype(float) / feature_len
    cpm = million * (rpk + pseudocount) / (rpk.sum(axis=0) * norm_factors + libsize_pseudocount)

    return np.log2(cpm) if log else cpm

#### configurations

# input
filtered_counts=pd.read_csv(snakemake.input[0], index_col=0)

# output
output_data=snakemake.output[0]

# normalize filtered data by selected standard method
norm_method = 'TMM'

# determine norm factors
norm_factors = calc_norm_factors(filtered_counts, method=norm_method)
norm_factors = pd.Series(norm_factors, index=filtered_counts.columns, name='norm_factor')

# calculate normalized counts
normalized_lcpm = cpm(filtered_counts, norm_factors=norm_factors, log=True, pseudocount=0.5, libsize_pseudocount=1)

# save normalized counts
normalized_lcpm.to_csv(os.path.join(output_data))

import os
import pandas as pd
import numpy as np
from itertools import compress

# visualization
import seaborn as sns
import matplotlib.pyplot as plt

# # utils for enumerating string lists
# from collections import defaultdict
# from itertools import count
# from functools import partial


#### configurations

# annot=observations x annotations
annot=pd.read_csv(snakemake.input[0], index_col=0)

# type of dimensionality reduction
dimred=snakemake.params["dimred"]

# variables=columns of annot to plot
variables=snakemake.params['variables']

#  label=plot title and file name
label = snakemake.params["label"]

# results folder
results_dir = snakemake.params["results_dir"]

# if 3 samples or less UMAP could not be performed
if (dimred=="UMAP" and annot.shape[0]<4):
    from pathlib import Path
    for variable in variables:
        Path(os.path.join(results_dir, dimred+"_"+label+"_"+variable+".svg")).touch()
    import sys
    sys.exit()

# if 2 samples or less PCA only consists of one PC
if (dimred=="PCA" and annot.shape[0]<3):
    from pathlib import Path
    for variable in variables:
        Path(os.path.join(results_dir, dimred+"_"+label+"_"+variable+".svg")).touch()
    import sys
    sys.exit()

# data=observations x features
data_df=pd.read_csv(snakemake.input[1], index_col=0)

#plot parameters
color_dict=None
alpha=1


# plot data in 2D
# sns.set(color_codes=True)
for variable in variables:
    fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))
    unique_vals = annot[variable].unique()

#     if len(unique_vals)==1 or (sum(pd.isna(unique_vals))>0):
#         print("{} only one value or contains NaNs".format(variable))
#         continue

    unique_vals = unique_vals[pd.notna(unique_vals)]

    if all([isinstance(i, (str, bool, np.bool_)) for i in unique_vals]):
        print('discrete variable ', variable)

        if color_dict==None:
#             colors = plt.cm.get_cmap("tab20").colors
            cm = plt.get_cmap('gist_ncar')
            colors=[cm(1.*i/len(unique_vals)) for i in range(len(unique_vals))]
        else:
            colors = [color_dict[x] for x in unique_vals]

        # plot data by unique value
        for g, c in zip(unique_vals, colors):
            c = [c]  # to remove the warnings
            axes.scatter(
                data_df.loc[annot.index[annot[variable] == g].tolist(), '0',],
                data_df.loc[annot.index[annot[variable] == g].tolist(), '1',],
                label=g,
                marker=".",
                c=c,
                alpha=alpha,
            )

            # centroids by mean
            axes.scatter(
                data_df.loc[annot.index[annot[variable] == g].tolist(), '0',].mean(),
                data_df.loc[annot.index[annot[variable] == g].tolist(), '1',].mean(),
                marker="s",
                alpha=0.5,
                label=str(g) + " centroid",
                c=c,
            )
            axes.text(
                data_df.loc[annot.index[annot[variable] == g].tolist(), '0',].mean(),
                data_df.loc[annot.index[annot[variable] == g].tolist(), '1',].mean(),
                s=g,
                horizontalalignment="center",
                verticalalignment="bottom",
                alpha=0.5,
            )

        # edit and position legend
        handles, labels = axes.get_legend_handles_labels()
        legend_idx = ["centroid" in label for label in labels]
        handles = list(compress(handles, legend_idx))
        labels = list(compress(labels, legend_idx))
        axes.legend(handles, labels, loc="center left", bbox_to_anchor=(1.05, 0.5))

    elif all([isinstance(i, (int, float, np.int, np.float, np.int64)) for i in unique_vals]):
        print('continous variable ', variable)


        # check if enumerated ie numerical -> probably not needed anymore
#         if annot[variable].dtype != 'float64' and annot[variable].dtype != 'int64':
#             label_to_number = defaultdict(partial(next, count(1)))
#             annot[variable]=[label_to_number[label] for label in annot[variable]]

        cmap = sns.cubehelix_palette(as_cmap=True)
        points = axes.scatter(
            data_df.loc[:, '0',],
            data_df.loc[:, '1',],
            marker=".",
            c=annot[variable],
            s=50, 
            cmap=cmap,
            alpha=alpha,
        )
        fig.colorbar(points)

    else:
        print("variable type not-detected for {}".format(variable))
        continue

    # show and save figure
    x_postfix=''
    y_postfix=''

    if dimred=='PCA':
        explained_variance=list(pd.read_csv(os.path.join(results_dir,"PCA_{}_{}_explained_variance.csv".format(snakemake.wildcards["split"],snakemake.wildcards["step"])), index_col=0)['0'])
        x_postfix=' ({:.1f}%)'.format(100*explained_variance[0])
        y_postfix=' ({:.1f}%)'.format(100*explained_variance[1])

    axes.set_xlabel(dimred+" 1"+x_postfix)
    axes.set_ylabel(dimred+" 2 "+y_postfix)
    plt.title(dimred+" of " + label + " "+variable)
    plt.show()

    # save plot
    fig.savefig(
        fname=os.path.join(results_dir, dimred+"_"+label+"_"+variable+".svg"),
        format="svg",
        dpi=300,
        bbox_inches="tight",
    )

import os
import pandas as pd

# configuration
output = snakemake.output[0]
results_dir = snakemake.params["results_dir"]

# load annotation file and get sample names
annotations = pd.read_csv(snakemake.input[0], index_col=0)

results=[]

for sample in annotations.index:
    print(sample)
    result=pd.read_csv(os.path.join(results_dir,"{}".format(sample),"mapped", "{}_quantification_counts.csv".format(sample)))
    results.append(result)

# results = [item for item in results if item is not None]
results = pd.concat(results)
results.set_index(results.columns[0]).T.to_csv(output)

import os
import pandas as pd
import pybedtools as bedtools

# configuration
output = snakemake.output[0]
chrom_file = snakemake.params["chrom_file"]
results_dir = snakemake.params["results_dir"]
sample=snakemake.wildcards["sample"]

elements_to_quantify = bedtools.BedTool(snakemake.input[0])

bamfile=snakemake.input[1]

print("Processing "+sample)
try:
    result= elements_to_quantify.coverage(b=bamfile,sorted=True,g=chrom_file).to_dataframe(
                names=["CHR", "START", "END", "ID", sample, "NA1", "NA2", "NA3"],
                dtype={sample: int},
                usecols=['ID', sample],
                index_col='ID').T
    result.to_csv(output)
except Exception as e:
    print("Error occured while processing sample "+sample)
    pd.DataFrame(0,index=elements_to_quantify.index,columns=[sample]).T.to_csv(output)

import os
import pandas as pd

# configuration
output = snakemake.output[0]
results_dir = snakemake.params["results_dir"]

# load annotation file and get sample names
annotations = pd.read_csv(snakemake.input[0], index_col=0)

results=[]

for sample in annotations.index:
    print(sample)
    result=pd.read_csv(os.path.join(results_dir,"{}".format(sample),"peaks", "{}_quantification_support.csv".format(sample)))
    results.append(result)

# results = [item for item in results if item is not None]
results = pd.concat(results)
results.set_index(results.columns[0]).T.to_csv(output)

import os
import pandas as pd
import pybedtools as bedtools

# configuration
output = snakemake.output[0]
chrom_file = snakemake.params["chrom_file"]
sample=snakemake.wildcards["sample"]

consensus_peaks = bedtools.BedTool(snakemake.input[0])
consensus_peaks_df = bedtools.BedTool(snakemake.input[0]).to_dataframe().set_index('name')

peakfile=snakemake.input[1]
result = pd.DataFrame(0,index=consensus_peaks_df.index,columns=[sample])

try:
    if (peakfile is not None):
        sample_peaks = bedtools.BedTool(peakfile)
        result = consensus_peaks.intersect(
            sample_peaks,
            g=chrom_file, 
            wa=True,
            c=True
        ).to_dataframe(
            index_col='name',
            usecols=[3,4],
            names=['name',sample]
        )
except Exception as e:
    print("Error occured while processing sample "+sample)
finally:
    result.T.to_csv(output)

import os
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# for robust dispersion norm calculation
from statsmodels import robust
import warnings
warnings.filterwarnings("ignore")

#### configurations

# input
data=pd.read_csv(snakemake.input[0], index_col=0)

# parameters
HVR_percentage=snakemake.params['HVR_percentage']
top_n = int((data.shape[0]/100)*HVR_percentage)

# output
output_data=snakemake.output[0]
output_plot=snakemake.output[1]

#####  determine highly variable regions (HRV) based on normalized dispersion
## inspired by “highly_variable_genes” function from scanpy 1.5.1 (Python 3.6.10) with the “flavor” argument set to “cell_ranger”
## but instead of dispersion=var/mean we use dispersion=std to keep it meaningful for negative mean values
## overlap of region selecrtion in test data was 96.6%

mean = data.mean(axis=1)
# var = data.var(axis=1)

# now actually compute the dispersion
mean[mean == 0] = 1e-12  # set entries equal to zero to small value
dispersion = data.std(axis=1) #var / mean

# all of the following quantities are "per-gene" here
df = pd.DataFrame()
df['means'] = mean
df['dispersions'] = dispersion

# bin regions by their mean in 20 bins
df['mean_bin'] = pd.cut(df['means'], np.r_[
    -np.inf,
    np.percentile(df['means'], np.arange(10, 105, 5)),
    np.inf
])
# group regions by mean_bin
disp_grouped = df.groupby('mean_bin')['dispersions']
# determine median (robust) dispersion by group
disp_median_bin = disp_grouped.median()

# the next line raises the warning: "Mean of empty slice"
with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    # calculate median absolute deviation (mad; robust measure of variance) per group
    disp_mad_bin = disp_grouped.apply(robust.mad)
    # determine "normalized" dispersion via (dispersion-dispersion_group_median)/dispersion_mad
    df['dispersions_norm'] = (df['dispersions'].values
        - disp_median_bin[df['mean_bin'].values].values
        ) / disp_mad_bin[df['mean_bin'].values].values

dispersion_norm = df['dispersions_norm'].values.astype('float32')



HVR_idx=df.index[df['dispersions_norm'].rank(ascending=False)<=top_n]
# HVR_idx=std.index[std.rank(ascending=False)<=top_n]
data_HVR = data.loc[HVR_idx, :]
data_HVR.to_csv(output_data)

# make plots describing the HVR selection
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
plt.subplots_adjust(wspace=0.5, hspace=0.5)

# plot region variability
ax = sns.scatterplot(df['dispersions_norm'].rank(ascending=False),
                     df['dispersions_norm'],
                     ax=axs[0], 
                     rasterized=True,
                    )
ax = sns.scatterplot(df.loc[HVR_idx,'dispersions_norm'].rank(ascending=False),
                     df.loc[HVR_idx,'dispersions_norm'],
                     ax=axs[0], 
                     rasterized=True,
                     color='r'
                    )
ax.set_xlabel('rank')
ax.set_ylabel('normalized dispersion')
sns.despine()
ax.set_title("Regions ranked by normalized dispersion")

# mean-dispersion plot (highlighting selected regions)
ax = sns.scatterplot(df.loc[:,'means'],
                     df.loc[:,'dispersions_norm'],
                     ax=axs[1], 
                     rasterized=True,
                    )
ax = sns.scatterplot(df.loc[HVR_idx,'means'], 
                     df.loc[HVR_idx,'dispersions_norm'],
                     ax=axs[1], 
                     rasterized=True,
                     color='r'
                    )

ax.set_xlabel('mean')
ax.set_ylabel('normalized dispersion')
sns.despine()
ax.set_title("Mean to normalized dispersion relationship")

fig.suptitle('Top {} HVRs by binned normalized dispersion'.format(top_n), fontsize=14)

plt.savefig(
            fname=output_plot,
            format="svg",
            dpi=300,
            bbox_inches="tight",
            )

plt.show()
plt.close()

import os
import pandas as pd
import numpy as np

#### configurations

# input
data_counts=pd.read_csv(snakemake.input[0], index_col=0)
annotation=pd.read_csv(snakemake.input[1], index_col=0)

# parameters
project_path = snakemake.params["project_path"]
split_by = snakemake.params["split_by"]

#####  make output folders, split and save data
for split in list(annotation[split_by].unique()):
    # make split folder
    if not os.path.exists(os.path.join(project_path,split)):
        os.mkdir(os.path.join(project_path,split))

    # split counts
    split_counts = data_counts.loc[:, annotation[(annotation[split_by]==split)].index]
    # save split counts
    split_counts.to_csv(os.path.join(project_path,split,"{}_counts.csv".format(split)))

    # split annotatios
    split_annot = annotation.loc[annotation[split_by]==split, :]
    # save split annotatios
    split_annot.to_csv(os.path.join(project_path,split,"{}_annotation.csv".format(split)))