#!/usr/bin/python


'''
SnakeMake pipeline to run ArrayScan2CNV
'''

import glob
import sys
import os
import shutil
import fnmatch
from snakemake.utils import R



configfile: "config.yaml"

scriptDir = config['repo_scripts']
gtc_dir = config['gtc_dir']
gtcFiles = []
SAMPLE_IDS = []
sampIdToGtcDict = {}
for root, dirnames, filenames in os.walk(gtc_dir):
    for filename in fnmatch.filter(filenames, '*.gtc'):
        gtc = os.path.join(root, filename)
        sampId = filename[:-4]
        SAMPLE_IDS.append(sampId)
        sampIdToGtcDict[sampId] = gtc
        gtcFiles.append(gtc)



def getGtc(wildcards):
    sampId = wildcards.sample_id
    return sampIdToGtcDict[sampId]


rule all:
    input:
        expand('gsrc/{sample_id}.gsrc.bed', sample_id = SAMPLE_IDS),
        expand('PennCnvOut/{sample_id}.PennCnvOut.txt', sample_id = SAMPLE_IDS)



rule gtc2PennCNV:
    input:
        gtc = getGtc,
        bpm = config['bpm']
    output:
        'PennCnvIn/{sample_id}.PennCNV_in.txt'
    shell:
        'module load python/2.7.5;python ' + scriptDir + '/gtc2PennCNV.py {input.gtc} {input.bpm} {output}'



rule gsrc2Cnv:
    input:
        'PennCnvIn/{sample_id}.PennCNV_in.txt'
    output:
        txt = 'gsrc/{sample_id}.gsrc.txt',
        bed = 'gsrc/{sample_id}.gsrc.bed'
    run:
        shell('module load R/3.3.0')
        R("""
            library(illuminaio)
            library(gsrc)
            filewpath <- "{input}"
            out_txt <- "{output.txt}"
            out_bed <- "{output.bed}"
            sample.name<-unlist(strsplit(filewpath,"\\."))[1]
            file.df<-read.delim(filewpath,sep="\t",header=T)
            dat<-list("snps" = c(),"chr" = c(), "pos" = c(), "baf" = c(), "rratio" = c(), "samples" = c(), "geno" = c())
            dat$snps<-file.df$Name
            dat$chr<-file.df$Chr
            dat$pos<-file.df$Position
            dat$baf<-file.df$B.Allele.Freq
            dat$rratio<-file.df$Log.R.Ratio
            dat$samples<-sample.name
            dat$geno<-file.df$Gtype

            #call segments-makes cna
            seg<-segm(dat)

            #gsrc cnv calling-makes cnv
            del<- -0.1
            dup<-0.1
            cnv.call<-cnv(seg,del=del,dup=dup)
            cna.bed<-data.frame("chr"=cnv.call$cna$chrom,"s"=cnv.call$cna$loc.start,"e"=cnv.call$cna$loc.end,
               "num.probes"=cnv.call$cna$num.mark,"avg.probe.val"=cnv.call$cna$seg.mean)
            cna.bed$state<-"NULL"
            del<- -0.1
            dup<-0.1
            for(i in 1:nrow(cna.bed)){
                valtmp<-cna.bed[i,"avg.probe.val"]
                val<-ifelse(valtmp<del,-1,ifelse(valtmp>dup,1,0))
                cna.bed[i,"state"]<-val
            }
            write.table(cna.bed, file = out_txt, quote=F, row.names=F, col.names=T)

            write.table(cbind(cna.bed[,1:3], cna.bed[,6]),file = out_bed, quote=F, row.names=F,col.names=F)

        """)



rule pennCnv:
    input:
        'PennCnvIn/{sample_id}.PennCNV_in.txt'
    output:
        'PennCnvOut/{sample_id}.PennCnvOut.txt'
    params:
        hmm = config['hmm'],
        pfb = config['pfb']
    log:
        'logs/{sample_id}.log'
    shell:
        'module load PennCNV/2015-v1.0.3;detect_cnv.pl -test -hmm {params.hmm} -pfb {params.pfb} -log {log} -out {output}'
