args = commandArgs(trailingOnly=TRUE)
# 1: counts file in
# 2: prefix for deseq out
# 3: name of contrast


library("DGCA")
library("dplyr")
library("yaml")
library("readr")
library("tidyr")
library("stringr")
library("magrittr")

xpr_in <- args[1]
nombre <- args[2]
cor_pre <- args[3]




#	build the sample DF
trammel <- read_yaml("config.yaml")
data_sets.df <- plyr::ldply(trammel$data_sets, data.frame)
data_sets.df$name <- as.factor(data_sets.df$name)
data_sets.df$day<- as.factor(data_sets.df$day)
data_sets.df$subgroups<- as.factor(data_sets.df$subgroups)
data_sets.df$rep<- as.factor(data_sets.df$rep)
data_sets.df$housing<- as.factor(data_sets.df$housing)
data_sets.df$genotype<- as.factor(data_sets.df$genotype)
data_sets.df$tissue<- as.factor(data_sets.df$tissue)

#	build the correlation experiment
corr.df <- plyr::ldply(trammel$correlations, data.frame) %>% filter(name == nombre )

# prep the sample metadata
sbgrp <- (corr.df %>% select(filt))[1,] %>% as.character
#coldata <- data_sets.df %>% as.data.frame() %>% ungroup()%>% filter(subgroups == sbgrp)
#coldata <- data_sets.df %>% ungroup() %>% group_by(name) %>%  mutate(genotype = paste0(as.character(genotype), collapse = "," ))  %>% filter(subgroups == sbgrp)
coldata <- data_sets.df%>% group_by(name) %>%  mutate(genotype = paste0(unique(as.character(genotype)), collapse = "." )) %>% ungroup()  %>% filter(subgroups == sbgrp)  %>% as.data.frame() 

rownames(coldata) <- coldata$name
dsgn.df <- coldata %>% select(name, corr.df$vars %>% as.character() )%>% mutate(present =1) %>% spread(key=genotype, value=present, fill=0) 
rownames(dsgn.df) <- dsgn.df$name
dsgn.df %<>% select(-c("name"))
dsgn.mat <- dsgn.df %>% as.matrix()




#xpr_in <- "expression/rotund.vs_dm6main.dm6_genes.mapspliceMulti.MpBC.rpkm"

#xpr_in <- "expression/rotund.vs_dm6main.dm6_genes.mapspliceMulti.MpBC.rpkm"


#load and format the count data
xpr.df <- read_delim(xpr_in, "\t", escape_double = FALSE, trim_ws = TRUE)
xpr.df %<>% as.data.frame()
rownames(xpr.df) <- xpr.df$Geneid
xpr.df %<>% select(rownames(dsgn.df)) 

ddcor_res <- ddcorAll(inputMat = xpr.df, design = dsgn.mat, compare = c("wt", "rn"),adjust = "none", heatmapPlot = FALSE, nPerm = 0)

#Write to TSV
write_delim(ddcor_res, paste(cor_pre, "corr","tsv", sep="."), delim = "\t")

args = commandArgs(trailingOnly=TRUE)
# 1: counts file in
# 2: prefix for deseq out
# 3: name of contrast


library("DESeq2")
library("dplyr")
library("yaml")
library("readr")
library("tidyr")
library("stringr")

counts_in <- args[1]
dsq_pre <- args[2]
nombre <- args[3]
config <- args[4]

#	build the sample DF
#trammel <- read_yaml("config.yaml")
trammel <- read_yaml(config)
data_sets.df <- plyr::ldply(trammel$data_sets, data.frame)
data_sets.df$name <- as.factor(data_sets.df$name)
data_sets.df$day<- as.factor(data_sets.df$day)
data_sets.df$subgroups<- as.factor(data_sets.df$subgroups)
data_sets.df$rep<- as.factor(data_sets.df$rep)
data_sets.df$housing<- as.factor(data_sets.df$housing)
data_sets.df$genotype<- as.factor(data_sets.df$genotype)
data_sets.df$tissue<- as.factor(data_sets.df$tissue)

#	build the contrast experiment
contrasts.df <- plyr::ldply(trammel$contrasts, data.frame) %>% filter(name == nombre )

# prep the sample metadata
sbgrp <- (contrasts.df %>% select(filt))[1,] %>% as.character
#coldata <- data_sets.df %>% as.data.frame() %>% ungroup()%>% filter(subgroups == sbgrp)
#coldata <- data_sets.df %>% ungroup() %>% group_by(name) %>%  mutate(genotype = paste0(as.character(genotype), collapse = "," ))  %>% filter(subgroups == sbgrp)
coldata <- data_sets.df%>% group_by(name) %>%  mutate(genotype = paste0(unique(as.character(genotype)), collapse = "." )) %>% ungroup()  %>% filter(subgroups == sbgrp) %>% unique() %>% as.data.frame() 

rownames(coldata) <- coldata$name
coldata <- coldata %>% select(contrasts.df$vars %>% as.character() )

#load and format the count data
counts.df <- read_delim(counts_in, "\t", escape_double = FALSE, trim_ws = TRUE)
counts.df.gath <- counts.df %>% select(-c(Chr, Start, End, Strand, Length)) %>% gather(key = "sample", value = "count", -one_of("Geneid"))
counts.df.sprud <- counts.df.gath %>% spread(key = sample, value = count) %>% as.data.frame() 
rownames(counts.df.sprud) <- counts.df.sprud$Geneid
counts.df.sprud <- counts.df.sprud %>%  select(rownames(coldata))

# convert counts, metadata, design into DESeq data object
#counts.dds <- DESeqDataSetFromMatrix(countData = counts.df.sprud, colData = coldata, design = ~ housing) 
desgn <- as.character(contrasts.df$design[1])
counts.dds <- DESeqDataSetFromMatrix(countData = counts.df.sprud, colData = coldata, design = eval(parse(text=desgn)))

#prefilter to remove low-count rows
keep <- rowSums(counts(counts.dds)) >= 10
counts.dds <- counts.dds[keep,]

#relevel 
#counts.dds$condition <- relevel(counts.dds$condition, ref = "untreated")
for (i in seq(nrow(contrasts.df))){
	var <- contrasts.df[i,]$vars %>% as.character()
	ruf <-contrasts.df[i,]$ref %>% as.character()
	counts.dds[[var]] <- relevel(counts.dds[[var]], ref = ruf)
}

#Run DESeq
counts.dds.dsq <- DESeq(counts.dds)

# run the vanilla DE
counts.dds.res <- results(counts.dds.dsq) %>% as.data.frame()
counts.dds.res$geneid <- rownames(counts.dds.res)

# Effect-size shrinkage
coff <- resultsNames(counts.dds.dsq)[2]
counts.dds.res.shrunk <- lfcShrink(counts.dds.dsq, coef=coff, type="apeglm") %>% as.data.frame()
#explore shrinkage types??
counts.dds.res.shrunk$geneid <- rownames(counts.dds.res.shrunk)

# scale by gene length
counts.dds.res<-inner_join(counts.dds.res, counts.df%>% select(Geneid, Length), by=c("geneid"="Geneid")) %>% mutate(expression = baseMean/Length) %>% select(-c(Length))
counts.dds.res.shrunk<-inner_join(counts.dds.res.shrunk, counts.df%>% select(Geneid, Length), by=c("geneid"="Geneid")) %>% mutate(expression = baseMean/Length) %>% select(-c(Length))

#Write to TSV
counts.dds.res.full <- inner_join(counts.dds.res %>%  as.data.frame(), counts.dds.res.shrunk %>%  as.data.frame(), by =c("geneid"="geneid"), suffix=c(".fresh",".shrunk")) %>% select(geneid, everything())
write_delim(counts.dds.res.full, paste(dsq_pre,"de",sep="."), delim = "\t")


###		Now run the itemized DE

#get the factors in the design
desgn.splt <- (desgn %>% str_replace_all("~","")%>% str_replace_all(" ","") %>% strsplit("+", fixed=TRUE))[[1]]

#first, collect the geneids and the basemeans: just run the first factor and its first nonref level
fact <- desgn.splt[1]
ref_lev <- (contrasts.df %>% filter(vars == fact))$ref
all_levs<-levels(counts.dds.dsq[[fact]])
alt_levs<- all_levs[all_levs != ref_lev]

basic.df <- results(counts.dds.dsq, contrast=c(fact, as.character(ref_lev), alt_levs[1])) %>% as.data.frame() %>% select(c("baseMean"))
basic.df$geneid <- rownames(basic.df)
basic.df <- select(basic.df, c("geneid", "baseMean"))

#for each factor in factors....

dds.dsq.byFactor.df <- select(basic.df, c("geneid") )


for (fact in desgn.splt){

	ref_lev <- as.character((contrasts.df %>% filter(vars == fact))$ref)
	all_levs<-levels(counts.dds.dsq[[fact]])
	alt_levs<- all_levs[all_levs != ref_lev]


#running basemean by level:
#https://support.bioconductor.org/p/63567/

	#	# gather baseMeans for each level in factor.levels:
	factor.baseMeanPerLvl <- sapply( levels(counts.dds.dsq[[fact]]), function(lvl) rowMeans( counts(counts.dds.dsq,normalized=TRUE)[,counts.dds.dsq[[fact]] == lvl, drop=F] ) ) %>% as.data.frame()

	names(factor.baseMeanPerLvl ) <- paste0(paste("baseMean.",fact,".", sep=""), names(factor.baseMeanPerLvl))
	factor.baseMeanPerLvl$geneid <- rownames(factor.baseMeanPerLvl)
	factor.baseMeanPerLvl <-inner_join(factor.baseMeanPerLvl, counts.df%>% select(Geneid, Length), by=c("geneid"="Geneid"))

	for (levi in all_levs){

		factor.baseMeanPerLvl$x <- factor.baseMeanPerLvl[[paste("baseMean",fact,levi,sep=".")]]
		factor.baseMeanPerLvl[[paste("expression",fact,levi,sep=".")]] <- factor.baseMeanPerLvl$x/factor.baseMeanPerLvl$Length

	}

	factor.baseMeanPerLvl <- factor.baseMeanPerLvl %>% select(-c("x", "Length"))

	dds.dsq.byFactor.df <- full_join(dds.dsq.byFactor.df, factor.baseMeanPerLvl, by = c("geneid"="geneid"))

#	#for each alt_level in factor.levels:
	for (aLev in alt_levs){

		#dds.res_factor_level <- results(counts.dds.dsq, contrast=c(fact, ref_lev, aLev))
		dds.res_factor_level.shrunk <- lfcShrink(counts.dds.dsq, coef=paste(fact, aLev, "vs", ref_lev, sep = "_"), type="apeglm") %>% as.data.frame()

		dds.res_factor_level.shrunk$geneid <- rownames(dds.res_factor_level.shrunk)
		dds.res_factor_level.shrunk<-inner_join(dds.res_factor_level.shrunk, counts.df%>% select(Geneid, Length), by=c("geneid"="Geneid")) %>% mutate(expression = baseMean/Length) %>% select(-c(Length))
		rownames(dds.res_factor_level.shrunk) <- dds.res_factor_level.shrunk$geneid
		dds.res_factor_level.shrunk <- dds.res_factor_level.shrunk %>% select(-c(geneid))
		names(dds.res_factor_level.shrunk ) <- paste0(names(dds.res_factor_level.shrunk),paste(".",fact,".vs_",aLev,".apeglm", sep=""))
		dds.res_factor_level.shrunk$geneid <- rownames(dds.res_factor_level.shrunk)

		dds.dsq.byFactor.df <- full_join(dds.dsq.byFactor.df, dds.res_factor_level.shrunk, by = c("geneid"="geneid"))
#		counts.dds.res_factor_level <- results(counts.dds.dsq, contrast=c(factor1, ref_level, alt_level1)) 
#		counts.dds.res_factor_level$geneid <- rownames(counts.dds.res_factor_level)

	}

}

#Write to TSV
dds.dsq.byFactor.df <- dds.dsq.byFactor.df %>% select(geneid, everything())
write_delim(dds.dsq.byFactor.df, paste(dsq_pre, "itemized.de", sep="."), delim = "\t")

args = commandArgs(trailingOnly=TRUE)
# 1: counts file in
# 2: prefix for deseq out
# 3: name of contrast


library("dplyr")
#library("yaml")
library("readr")
library("tidyr")
library("magrittr")

#library("stringr")

counts_in <- args[1] # counts from featureCount
#counts_in <- "counts/all.vs_dm6main.dm6_genes.mapspliceMulti.MpBC.counts"
aln_meta <- args[2] # metadata for alignments (to get million mapped reads)
#aln_meta <- "summaries/alignments.vs_dm6main.mapspliceMulti.summary"
exp_out <- args[3] # output file prefix to write expression
#exp_out <- "expression.potato"

counts.df <- read_delim(counts_in, "\t", escape_double = FALSE, trim_ws = TRUE)

RPK.df <- counts.df %>% select(-c(Chr, Start, End, Strand)) %>% mutate(Length = Length/1000) # convert bases to kb
#turn counts into per-kb rates
#https://stackoverflow.com/a/48693978
RPK.df[,-(1:2)] %<>% sapply(`/`, RPK.df[,2])
RPK.df.gath <- RPK.df  %>% gather(key = "sample", value = "count", -one_of("Geneid", "Length"))


#collect the million reads mapped (proper pairs, given the current featureCounts settings)
aln_meta.df <- read_delim(aln_meta, "\t", escape_double = FALSE, trim_ws = TRUE,col_names = F)
names(aln_meta.df)<- c("ref_genome","aligner","sample","measure","value")

aln_meta.df %<>% filter(measure == "properly_paired_count") %<>% select(-c("measure", "aligner", "ref_genome"))


RPKM.df.gath <- inner_join(RPK.df.gath, aln_meta.df, by =c("sample"="sample"))

RPKM.df.gath %<>% mutate(MMR = value/1000000, expression = count/MMR )  %>% select(c("Geneid", "sample", "expression"))

RPKM.df.sprud <- RPKM.df.gath %>% spread(key = "sample", value = "expression")

write_delim(RPKM.df.sprud, paste(exp_out,".rpkm", sep=""), delim = "\t", col_names = T)





RPK.df <- counts.df %>% select(-c(Chr, Start, End, Strand)) %>% mutate(Length = Length/1000) # convert bases to kb
# TPM definition https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/
#				https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/
RPK.df[,-(1:2)] %<>% sapply(`/`, RPK.df[,2])

RPK.gath.df <- RPK.df %>% select(-c("Length")) %>% gather(key = "sample", value = "rpk", -one_of("Geneid"))
RPK.gath.grupt.sum.df <- RPK.gath.df  %>% group_by(sample) %>% summarise(total=sum(rpk))

TPM.df <- inner_join(RPK.gath.df, RPK.gath.grupt.sum.df, by=c("sample"="sample")) %>% mutate(TPM = rpk/(total/1000000)) %>% select(-c("rpk","total")) %>% spread(key=sample, value=TPM) 

write_delim(TPM.df, paste(exp_out,".tpm", sep=""), delim = "\t", col_names = T)

import json
#	https://stackoverflow.com/questions/19483351/converting-json-string-to-dictionary-not-list
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("json_in", help="JSON file to be condensed")
parser.add_argument("flat_out", help="flatfile summary")
parser.add_argument("-t", "--tag", help="line-name for the processed JSON", default=None)
args = parser.parse_args()

def prune_jason(jsn):
	sparse_dict = {}
	sparse_dict = dict([(k,{}) for k in ['prefiltered', 'postfiltered', 'filtration', 'duplication']])

	for i in ['gc_content','q20_rate', 'q30_rate', 'read1_mean_length','total_reads']:
		try:
			sparse_dict['prefiltered'][i] =  jsn['summary']['before_filtering'][i] 
		except KeyError:
			sparse_dict['prefiltered'][i] = "NA"

		try:
			sparse_dict['postfiltered'][i] =  jsn['summary']['after_filtering'][i] 
		except KeyError:
			sparse_dict['postfiltered'][i] = "NA"

	for i in ['adaptor_trimmed_reads']:
		try: 
			sparse_dict['filtration'][i] = jsn['adapter_cutting']['adapter_trimmed_reads']
		except KeyError:
			sparse_dict['filtration'][i] = "NA"

	for i in ['corrected_reads','low_quality_reads','passed_filter_reads']:
		try:
			sparse_dict['filtration'][i] = jsn['filtering_result'][i]
		except KeyError:
			sparse_dict['filtration'][i] = "NA"

	for i in ['rate']:
		sparse_dict['duplication'][i] = jsn['duplication'][i]

	return sparse_dict

def print_pruned(prn):
	lines2write = [ [x, y, prn[x][y]]  for x in prn.keys() for y in prn[x]]
	if args.tag:
		[ ell.insert(0, args.tag) for ell in lines2write ]
	phial_out = open(args.flat_out, 'w')
	for preline in lines2write:
		field_count = len(preline)
		line = ("%s" + "\t%s"*(field_count-1) + "\n") % tuple(preline)
		phial_out.write(line)
	phial_out.close()

jason_file=open(args.json_in)
jason_str = jason_file.read()
jason = json.loads(jason_str)

jason_condensed = prune_jason(jason)
print_pruned(jason_condensed)

args = commandArgs(trailingOnly=TRUE)

library("readr")
library("biomaRt")
library("org.Dm.eg.db")
library("topGO")
library("dplyr")
library("yaml")
library("stringr")
library("magrittr")
data("geneList")

deSeq_file <- args[1]
nombre <- args[2]
go_out <- args[3]
config <- args[4]

#nombre <- "hausWtVsMut"


sig_thresh <- 0.01

#hausWtVsMut_vs_dm6main_dm6_genes_mapspliceMulti_MpBC_itemized <- read_delim("diff_expr/hausWtVsMut/hausWtVsMut.vs_dm6main.dm6_genes.mapspliceMulti.MpBC.itemized.de", "\t", escape_double = FALSE, trim_ws = TRUE)
DESeq.itemized<- read_delim(deSeq_file, "\t", escape_double = FALSE, trim_ws = TRUE)

trammel <- read_yaml(config)
#trammel <- read_yaml("config.yaml")
#	build the contrast experiment
contrasts.df <- plyr::ldply(trammel$contrasts, data.frame) %>% filter(name == nombre )
#get the factors in the design
desgn <- as.character(contrasts.df$design[1])
desgn.splt <- (desgn %>% str_replace_all("~","")%>% str_replace_all(" ","") %>% strsplit("+", fixed=TRUE))[[1]]

query_genes.list <- DESeq.itemized$geneid %>% unique() %>% as.list()

#ensembl_dm = useMart("ensembl",dataset="dmelanogaster_gene_ensembl")
marty <- useDataset("dmelanogaster_gene_ensembl",  useMart("ensembl",  host = "useast.ensembl.org") )
G_list <- getBM(attributes= c("flybase_gene_id", "ensembl_gene_id", "external_gene_name", "go_id", "name_1006"), mart= marty, useCache= FALSE) %>% mutate(go_name = name_1006) %>% select(-c("name_1006"))
#G_list <- getBM(attributes= c("flybase_gene_id", "ensembl_gene_id", "external_gene_name", "go_id", "name_1006"), mart= marty) %>% mutate(go_name = name_1006) %>% select(-c("name_1006"))


### ### ### ### ### ### ### ### 
### plz fix this
go_bad <- c("GO:0110109","GO:0120176","GO:0120177","GO:0120170","GO:0062023")
bad_boiz <- G_list %>% filter(go_id %in% go_bad) %>% dplyr::select("flybase_gene_id")
bad_boiz <- rbind(bad_boiz, "FBgn0050046", "FBgn0033710",  "FBgn0050203", "FBgn0026721")
#mask the baddies
G_list.unq <- G_list %>% filter(!flybase_gene_id %in% bad_boiz$flybase_gene_id ) %>% dplyr::select("flybase_gene_id") %>% unique()
DESeq.itemized %<>% filter(!(geneid %in% bad_boiz$flybase_gene_id))
### ### ### ### ### ### ### ### 

#present_gene_list <- G_list.unq$flybase_gene_id %in% DESeq.itemized$geneid %>% as.numeric()
#names(present_gene_list) <- G_list.unq$flybase_gene_id

jeanOnt.mega.df <- as.data.frame(c())

#for each factor in factors....

for (fact in desgn.splt){

	ref_lev <- (contrasts.df %>% filter(vars == fact))$ref
	#all_levs<-levels(counts.dds.dsq[[fact]])

	#subset the DE data to the factor
#	DESeq.itemized.factor <- DESeq.itemized[,!is.na(str_locate(names(DESeq.itemized), paste(fact,".vs_" ,sep=""))[,1])]
	DESeq.itemized.factor <- DESeq.itemized[,!is.na(str_match(names(DESeq.itemized), paste(fact,".vs_(.*?).apeglm" ,sep=""))[,1])]

	#pull out the alt levels from the data
#	alt_levs<- (str_split_fixed(names(DESeq.itemized.factor), paste(fact,".vs_", sep=""), n=2)[,2] %>% str_split_fixed("\\.", n=2))[,1] %>% unique()

	noms <- str_match(names(DESeq.itemized), paste(fact,".vs_(.*?).apeglm" ,sep=""))[,2]
	alt_levs <- unique(noms[!is.na(noms)])
	DESeq.itemized.factor$geneid <- DESeq.itemized$geneid

	#	#for each alt_level in factor.levels:

	for (alt in alt_levs) {
		print(alt)
		#subset the DE data to that level
		DESeq.itemized.level <- DESeq.itemized.factor[,!is.na(str_locate(names(DESeq.itemized.factor), paste(fact,"\\.vs_",alt,"\\." ,sep=""))[,1])]
		DESeq.itemized.level$geneid <- DESeq.itemized.factor$geneid
		#pull out the genes & significance values
		sig_gene_list <- DESeq.itemized.level[paste("padj.",fact,".vs_",alt,".apeglm",sep="")][[1]]# < sig_thresh)[,1] #%>% as.numeric()
		sig_gene_list[is.na(sig_gene_list)] <- 1
		names(sig_gene_list) <- DESeq.itemized.level$geneid

		jeanOnt.df <- as.data.frame(c())

		for (ont in c("MF","BP","CC") ){

			#topDiffGenes uses p cutoff of 0.01
			GOdata <- new("topGOdata", ontology = ont, allGenes = sig_gene_list, geneSel = topDiffGenes, annot = annFUN.org, mapping = "org.Dm.eg.db", ID="Ensembl")
			fishy <- runTest(GOdata, algorithm = "classic", statistic = "fisher")
			kissy <- runTest(GOdata, algorithm = "classic", statistic = "ks")

			goRes.tmp <- GenTable(GOdata, classicFisher =fishy, classicKS = kissy,  topNodes = 50) %>% mutate(ontology = as.factor(ont), classicFisher = as.numeric(classicFisher), classicKS = as.numeric(classicKS)) %>% select(-c(Term), Term) %>% as_tibble() 
			goRes.tmp %<>% left_join(G_list%>% select(c("go_id", "go_name")) %>% unique(), by =c("GO.ID" = "go_id")) %>% select(-c("Term"))
			jeanOnt.df <- rbind(jeanOnt.df, goRes.tmp)
			#store & bind the topGO results. 

		}

		jeanOnt.df$variable <- as.factor(fact)
		jeanOnt.df$alt_level <- as.factor(alt)

		jeanOnt.mega.df <- rbind(jeanOnt.mega.df,jeanOnt.df)

	}

}

write.table(jeanOnt.mega.df, file=args[3], row.names=FALSE, col.names = TRUE, sep = "\t", quote=F)

import pysam
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("-i", "--samfile_in", help="file in SAM format to ")
# parser.add_argument("-i", "--idxstat_in", help="samtools idxstat report")
# parser.add_argument("-g", "--split_coverage", help="bedtools genomecov report, with split reads")
# parser.add_argument("-G", "--spanning_coverage", help="bedtools genomecov report, with spanning reads")
# parser.add_argument("-d", "--depthstats_in", help="samtools depth report")
# parser.add_argument("-D", "--dupestats_in", help="samtools markdup stats report")
# parser.add_argument( "-m", "--mapping_multiplicity", help="histogram of mapping multiplicity scraped from the IH:: tags")
# parser.add_argument( "-c", "--mapped_count", help="count of mapped reads")
parser.add_argument("-o", "--samfile_out", help="flatfile summary")
# parser.add_argument("-t", "--tag", help="line-name for the flatfile", default=None)
args = parser.parse_args()

sam_in = args.samfile_in
sammy = pysam.AlignmentFile(sam_in,"rb")  


sam_out = args.samfile_out
samuel = pysam.AlignmentFile(sam_out, "w", template=sammy)


# cigartuples encoding:
# M = 0
# I = 1
# D = 2
# N = 3
# S = 4


for read in sammy.fetch():
	cig = read.cigar
	nu_cig = []
	suffix = []
	if cig[0][0] == 4:
		cig = cig[1:]
	if cig[-1][0] == 4:	#discard the softclips
		cig = cig[:-1]
	front_pad_tally = 0
	while cig[0][0] in [0,1,2]:
		front_pad_tally += cig[0][1]
		cig.pop(0)
	back_pad_tally = 0
	while cig[-1][0] in [0,1,2]:
		back_pad_tally += cig[-1][1]
		cig.pop(-1)

	nu_cig.extend([(3,front_pad_tally-1),(0,1)]) # new cigar starts with 5' end clipped down to the actual splice

	# if cig[0][0] == 0: #take the terminal match blocks
	# 	num_match = cig[0][-1]
	# 	prefix = [(3,num_match-1),(0,1)]
	# 	nu_cig.extend(prefix) #add a bunch of Ns to pad ; ignore 5' and 3' read ends. 
	# 	cig = cig[1:]
	# if cig[-1][0] == 0:
	# 	num_match = cig[-1][1]
	# 	suffix = 
	# 	cig = cig[:-1]

	rill_tally = 0
	mids = []
	for rillo in cig:
		if rillo[0] in [0,1,2]:
			rill_tally += rillo[1] #consolidate exons
		else:
			if rill_tally > 0:
				mids.extend([(0,1),(3,rill_tally-2),(0,1)]) #write hollowed exon
			mids.extend([rillo]) # write intron
			rill_tally=0 # reset tally
	nu_cig.extend(mids)
	nu_cig.extend([(0,1),(3,back_pad_tally-1),]) # end cigar with splice and Ns to 3' end

	read.cigar = nu_cig
	read.seq = "*"
	samuel.write(read)

samuel.close()

run:
	print("%s%s" % tuple(["FASTQ/", wildcards.path]))
	try:
		sra = sra_by_fqpath["%s%s/" % tuple(["FASTQ/", wildcards.path])]
		shell(""" mkdir -p FASTQ/{wildcards.path}/ """)
		shell("""
			fasterq-dump  --split-3 --outdir FASTQ/{wildcards.path}/ {sra}
		""")
	except KeyError:
		raise KeyError("Sample is listed as empirical but no reads found and no SRA to download!" )