shell.executable("/bin/bash")
shell.prefix("source ~/.bash_profile; ")
import os
import datetime
""" Snakemake pipeline for RNA editing detection from RNA-Seq data """

configfile: "config.yaml"
  
DATA = config["DATA"]

LIB = "../src"
DBASES = "../dbases"
RESULTS = "../results"
DOCS = "../docs"

PICARD = config["PICARD"]
GATK = config["GATK"]
SNPEFF = config["SNPEFF"]
PANDOC = config["PANDOC"]
TRANSCRIPTS = config["TRANSCRIPTS"]
GENOME = config["GENOME"]
GENOME_DIR = config["GENOME_DIR"]
ANNOTATION_DIR = config["ANNOTATION_DIR"]
FASTQ = config["FASTQ"]
CHROM_SIZES = config["CHROM_SIZES"]
POS_REF = config["pos_ref"]
POS_ALT = config["pos_alt"]
RESULT_DATE = ["2018-08-06"]
VARIANT_TYPES = ["variant_allele_counts_by_strand"]
FQ_BASENAME = [os.path.basename(x) for x in FASTQ]

def comp(x):

  bp_dict = {"A":"T",
             "T":"A",
             "G":"C",
             "C":"G"}

  return bp_dict[x]

C_POS = [comp(x) for x in POS_REF]
C_ALT = [comp(x) for x in POS_ALT]
FOREBRAIN_FQ = [x for x in FASTQ if "brainrest" in x]

rule all:
  input:
    #### trim brainrest samples ####
    expand("{data}/raw_data/{fq_fastq}_{read}_001.fastq.gz",
      data = DATA, 
      fq_fastq = FOREBRAIN_FQ, 
      read = ["R1", "R2"]),

    #### index and align data with two pass star mapping ####
    GENOME_DIR + "/Genome",
    expand("{data}/star/{fastq}_2pass_Aligned.sortedByCoord.out.bam",
    data=DATA, fastq=FASTQ),
    expand("{data}/star/{fastq}_2pass_Aligned.sortedByCoord.out.bam.bai",
    data=DATA, fastq=FASTQ),
    expand("{data}/star/{fastq}_split.bam",
    data=DATA, fastq=FASTQ),

   #### run base recalibration 2x and variant calling 3x
    expand("{data}/vcf/uncalibrated/{fastq}_firstbootstrap.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/uncalibrated/filtered/{fastq}_hardfiltered.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/uncalibrated/merged_uncalibrated_filtered.vcf.gz",
    data=DATA),
    expand("{data}/vcf/recalibrated/{fastq}_recal_1x.table",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/{fastq}_secondbootstrap.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/filtered/{fastq}_hardfiltered_1x.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/merged_recalibrated_hardfiltered_1x.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/{fastq}_recal_2x.table",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/{fastq}_thirdbootstrap.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/filtered/{fastq}_hardfiltered_2x.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/recalibrated/merged_recalibrated_hardfiltered_2x_snps.vcf.gz",
    data=DATA, fastq=FASTQ),

   # #### get stats on vcfs pre and post filtering
    expand("{data}/vcf/recalibrated/merged_recalibrated_2xpass.vcf.gz",
    data=DATA),
    expand("{data}/vcf/recalibrated/filtered/dp10/{fastq}_hardfiltered_2x_snps.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/uncalibrated/filtered/dp10/{fastq}_filtered.vcf.gz",
    data=DATA, fastq=FASTQ),
    expand("{data}/vcf/{calibration}/stats/filtered/{fastq}_stats.txt",
    data=DATA, fastq=FASTQ, calibration=['recalibrated', 'uncalibrated']),

   # #### annotate variants ####
    ANNOTATION_DIR + "gtf_exons.bed.gz.tbi",
    ANNOTATION_DIR + "gtf_transcripts.bed.gz",

   # #### get counts, kallisto, and bigwigs #####
    expand("{data}/featurecounts/count_summary.tsv",
    data=DATA),
    expand('{data}/bigwigs/{fastq}.pos.bw', 
    data=DATA, fastq=FASTQ ),
    expand('{data}/bigwigs/{fastq}.neg.bw', 
    data=DATA, fastq=FASTQ ),
    
   # DBASES + "/kallisto_idx/Ictidomys_tridecemlineatus.spetri2.denovo.kallisto.idx",
   # expand("{data}/kallisto_aln/{fastq}/abundance.tsv", 
   # data=DATA, fastq=FASTQ),

   # #### get editing counts and parse out variants into all possible snvs ####

    expand("{data}/vcf/variant_allele_counts_by_strand/strand_counts/{fastq}_strand_counts.txt.gz",
    data=DATA, fastq=FASTQ ),
    expand("{data}/vcf/variant_allele_counts_by_strand/filtered.bed",
    data=DATA),
    
   # # candidate for ugliest snakemake code ever written below 
   # # generate counts for all possible combinations of variants, 
   # # but not A to A for example
   # 
    expand(expand("{data}/vcf/{{variant_types}}/{pos_ref}_{pos_alt}_alleles/counts/{{fastq}}_counts.txt.gz",
    zip, data=[DATA] * len(POS_REF), pos_ref = POS_REF, pos_alt = POS_ALT), 
    fastq = FASTQ, variant_types = VARIANT_TYPES),
  #  
  # # expand(expand("{data}/vcf/{{variant_types}}/{pos_ref}_{pos_alt}_alleles_trimmed/counts/{{fastq}}_counts.txt.gz",
  # # zip, data=[DATA] * len(POS_REF), pos_ref = POS_REF, pos_alt = POS_ALT), 
  # # fastq = FASTQ, variant_types = VARIANT_TYPES),
  #  
    expand(expand("{data}/vcf/{{variant_types}}/{pos_ref}_{pos_alt}_alleles/filtered_select_variants.bed.gz",
    zip, data=[DATA] * len(POS_REF), pos_ref = POS_REF, pos_alt = POS_ALT), 
    variant_types = VARIANT_TYPES),
    
    expand("{data}/vcf/variant_allele_counts_by_strand/annotated_variants.bed.gz",
    data=DATA),
  #  
  # # #### run Rmarkdown for Diff editing ####
    expand(expand("{data}/../results/{{result_date}}/{{variant_types}}/{pos_ref}_{pos_alt}_alleles/{pos_ref}_{pos_alt}_alleles.html", 
    zip , data = [DATA] * len(POS_REF), pos_ref = POS_REF, pos_alt =
    POS_ALT), result_date = RESULT_DATE, 
    variant_types = VARIANT_TYPES ),
    
   # expand(expand("{data}/../results/{{result_date}}/{{variant_types}}/{pos_ref}_{pos_alt}_alleles_trimmed/{pos_ref}_{pos_alt}_alleles_trimmed.html", 
   # zip , data = [DATA] * len(POS_REF), pos_ref = POS_REF, pos_alt =
   # POS_ALT), result_date = RESULT_DATE, 
   # variant_types = VARIANT_TYPES ),


   # #### find hyperedited regions ####
   # expand("{data}/star/{fastq}_2pass_Unmapped.out.mate{read}.fastq.gz",data = DATA,
   # fastq = FASTQ, read = ["1", "2"]),
   # expand("{data}/hyperedits/{fastq}_2pass_Unmapped.out.mate{read}_filtered.fastq.gz",
   # data = DATA, fastq=FASTQ, read = ["1", "2"]),
   # 
   # #expand(expand("{{data}}/hyperedits/{{fastq}}_2pass_Unmapped.out.mate{{read}}_filtered_{ref}to{alt}.fastq.gz",
   # #zip, ref = POS_REF, alt = POS_ALT),
   # #data = DATA, fastq=FASTQ, read = ["1", "2"]),

   # expand(DBASES + "/hyperediting/ensemb85_genome_{ref}to{alt}.fasta.bwt", zip,
   # ref = POS_REF, alt = POS_ALT), 
   # 
   # #expand(expand("{{data}}/hyperedits/{{fastq}}_{{read}}_{ref}to{alt}_gnome{c_ref}to{c_alt}.bam",
   # #zip, ref = POS_REF*2, alt = POS_ALT*2, c_ref = C_POS + POS_REF, c_alt= C_ALT + POS_ALT), 
   # #data = DATA, fastq=FASTQ, read = ["1", "2"]), 
   # 
   # #expand(expand("{{data}}/hyperedits/bed/{{type}}/{ref}_to_{alt}/{{fastq}}_{{read}}_{f_ref}to{f_alt}_gnome{c_ref}to{c_alt}.bed",
   # #zip, ref = POS_REF*4, alt = POS_ALT*4, f_ref = POS_REF*2 + C_POS*2, f_alt = POS_ALT*2 + C_ALT*2, c_ref = (C_POS + POS_REF)*2, c_alt= (C_ALT + POS_ALT)*2), 
   # #data = DATA, fastq=FASTQ, type = ["edits", "reads"], read = ["1", "2"]),
   # 
   # #expand(expand("{{data}}/hyperedits/bed/{{type}}/filtered/{ref}_to_{alt}/{{fastq}}_{{read}}_{f_ref}to{f_alt}_gnome{c_ref}to{c_alt}.bed",
   # #zip, ref = POS_REF*4, alt = POS_ALT*4, f_ref = POS_REF*2 + C_POS*2, f_alt = POS_ALT*2 + C_ALT*2, c_ref = (C_POS + POS_REF)*2, c_alt= (C_ALT + POS_ALT)*2), 
   # #data = DATA, fastq=FASTQ, type = ["edits", "reads"], read = ["1", "2"]),
   #  
   # expand("{data}/hyperedits/site_counts/A_G_counts/{fastq}_counts.txt.gz",
   # data = DATA, fastq=FASTQ),

    #### check editing site positions """"
    #expand("{data}/edit_sites/edit_positions/{fastq}_editpos.txt",
    #data = DATA, fastq=FASTQ),
    #expand("{data}/edit_sites/nonsig_edit_positions/{fastq}_editpos.txt",
    #data = DATA, fastq=FASTQ),

include: "rules/star_idx.snake"
include: "rules/star_two_pass.snake"
include: "rules/mark_duplicates.snake"
include: "rules/base_recalibration.snake"
include: "rules/featurecounts.snake"
include: "rules/make_bigwigs.snake"
include: "rules/vcf_stats.snake"
include: "rules/parse_strandedness.snake"
include: "rules/run_reports.snake"
include: "rules/bam_stats.snake"
#include: "rules/hyperediting.snake"
include: "rules/qc_check.snake"
include: "rules/trim_forebrain.snake"
