shell:
  """
  mkdir -p {params.dir}
  ln -s {input} {output}
  """

shell: """
  # Get encoding of Phred Quality Scores
  module load python
  encoding=$(python {params.script_dir}/phred_encoding.py {input.R1})
  echo "Detected Phred+${{encoding}} ASCII encoding"

  module load bbtools/38.87
  bbtools bbmerge-auto in1={input.R1} in2={input.R2} qin=${{encoding}} \
  ihist={output} k=62 extend2=200 rem ecct -Xmx900G
      """

shell:
  """
  module load bismark/0.23.0
  module load samtools/1.15
  cd {params.dir}
  samtools view -hb {input.F1} | samtools sort -n -@ {threads} -O BAM -o {output.T1}
  deduplicate_bismark --paired --bam --outfile {params.prefix} {output.T1}
  samtools view -T {params.genome} -@ 16 -h -C {output.B1} > {output.C1}
  samtools flagstat -@ {threads} {output.B1} > {output.B2}
  """

shell:
  """
  mkdir -p {params.outdir}
  module load bismark samtools bowtie
  bismark_methylation_extractor --paired-end --no_overlap --multicore 8 --gzip --report --bedGraph --counts --buffer_size 100G --no_header --cytosine_report --output {params.outdir} --scaffolds --genome_folder {params.bismark_index} {input.bam}
  """

shell:
  """
  module load bismark/0.23.0
  mkdir -p {params.dir}
  cd {params.dir}
  cp reference_genome {params.dir}/genome.fa
  bismark_genome_preparation --verbose --parallel {threads} {params.dir} #--single_fasta
  """

shell:
  """
  module load bismark/0.23.0
  mkdir -p {params.dir}
  bismark --multicore {threads} --temp_dir /lscratch/$SLURM_JOBID/ {params.command} --output_dir {params.dir} --genome {params.genome_dir} -1 {input.F1} -2 {input.F2}
  """

shell:
  """
  module load rseqc/4.0.0
  mkdir -p {params.dir}
  infer_experiment.py -r {params.bedref} -i {input.file1} -s 1000000 > {output.out1}
  read_distribution.py -i {input.file1} -r {params.bedref} > {output.out2}
  """

shell:
  """
  module load rseqc/4.0.0
  mkdir -p {params.dir}
  inner_distance.py -i {input.bam} -r {params.genemodel} -k 10000000 -o {params.prefix}
  """

shell:
  """
  mkdir -p {params.dir1}
  mkdir -p {params.dir2}
  module load R
  Rscript {params.script_dir}/get_methy.R {input} {wildcards.samples} {params.cutoff} {output}
  """

run:
  df_ref=pd.read_csv(params.map_table,sep='\t',header=None)
  df_ref.columns=['chromosome','start','end','cgid']
  df_ref=df_ref.loc[(df_ref['chromosome'].isin(CHRS)),]
  dfm=pd.read_csv(input[0])
  dfm=pd.merge(df_ref,dfm,on=['chromosome','start','end'],how='inner')
  dfm=dfm.drop(labels=['chromosome','start','end'],axis=1)
  dfm=dfm.set_index('cgid')
  dfm.to_csv(output[0])

shell:
  """
  module load python
  cd {params.dir}
  python {params.script_dir}/deconvolve2.py --atlas_path {params.ref} --residuals {input}  > {output}  2>&1
  """

run:
  dfm=pd.read_csv(input[0])
  for f in input[1:]:
  	df=pd.read_csv(f)
  	dfm=pd.merge(dfm,df,on='cgid',how='outer')
  dfm.to_csv(output[0],index=False)

shell:
  """
  module load python
  cd {params.dir}
  python {params.script_dir}/deconvolve2.py --atlas_path {params.ref} --residuals {input}
  """

shell:
  """
  module load bedtools crossmap
  crossmap bed {params.lift_file} {input.bed} {output.graph}
  """

shell:
  """
  module load bedtools
  bedtools sort -i {input.graph} | bedtools intersect -wo -a {params.markers} -b stdin -sorted | awk '$6-$5==1 {print $0}' | awk 'NF{NF-=1};1' > {output.sort}
  """

shell:
  """
  module load R
  Rscript {params.script_dir}/aggregate_over_regions.R {input.sort} {output.tsv}
  """

shell:
  """
  module load R
  Rscript ${params.script_dir}/tissues_of_origin_v2.R \
  {input.bed}  \
  {params.reference_markers} \
  {output.tsv} \
  {params.reference_IDs} \
  {params.sampleName} \
  TRUE \
  FALSE \
  TRUE \
  colon_5/hsc_5/hsc_2/spleen_1/spleen_2/spleen_3/spleen_4/
  """

shell:
  """
  module load R
  mkdir -p {params.dir}
  Rscript {params.script_dir}/bsseq_lm.R {params.chr} {input.bizfile} {output.bed1} {params.sample_prop} {params.cov}
  """

shell:
  """
  mkdir -p {params.dir}
  module load R/4.1
  Rscript {params.script_dir}/mvp_plots.R {input.dir} {params.fdr_cutoff} {params.highlight} {output.miami} {output.violin}
  ## {output.pie_cpg} {output.pie_genic}
  """

shell:
  """
  module load R
  mkdir -p {params.dir}
  zcat {params.prefix}*.regions-p.bed.gz | grep -v "^#" > {output.p1}
  Rscript {params.script_dir}/combp_Manhatten.R {output.p1} {output.f3}
  Rscript {params.script_dir}/combp_GREAT.R {output.p1} {output.f1} {output.o1} {output.o2}
  """

shell:
  """
  mkdir -p {params.dir}
  ls {params.dir}/{wildcards.group}_*_betas_pval.bed > {output.LIST}
  module load R
  Rscript {params.script_dir}/bsseq_Manhatten.R {output.LIST} {output.MAN}
  """