inputdir="/fast3/ranaseq/results/1/FASTQ"
analysisname="/fast3/ranaseq/results/1/analysis_340"
reportsdir="/fast3/ranaseq/scripts/reports"
methodsdir="/fast3/ranaseq/scripts/analysis"
inputquant="quant.genes.sf"
alignparams=rbind(c("1","A","","","0","0"))
files=data.frame(id=c(100007,100005,100003,100001),name=c("SRR5615379_1","SRR5615378_1","SRR5615377_1","SRR5615376_1"),samplename=c("Nrf2_KO_2","Nrf2_KO_1","WT_2","WT_1"))
samples=c("Nrf2_KO_1","Nrf2_KO_2","WT_1","WT_2")
groups=c(0,0,1,1)
groupa="Samples A"
groupb="Samples B"
paired=FALSE
method="DESeq2"
test="Wald"
fitType="parametric"
cutoff=0.05
tabgen_file="/fast3/ranaseq/genomes/Mus_musculus.GRCm38.cdna.all.fa.tabgen"
fdbdir="/fast3/ranaseq/genomes/functional/"
taxon_id="10090"

source(paste0(methodsdir,'/functions.R'))
load(file=paste0(analysisname,"/tpms.RData"))
load(file=paste0(analysisname,"/counts.RData"))

if(nrow(counts)==1)
  stop("Only one gene is expressed. Please use the raw counts data for your analyses. RaNA-Seq is conceived for the analysis of whole transcriptome.")

setwd(analysisname)

suppressMessages(library(parallel))
no_cores <- detectCores() / 2

## SERE coefficent functions
SERE <- function(observed){
  #calculate lambda and expected values
  laneTotals <- colSums(observed)
  total <- sum(laneTotals)
  fullObserved <- observed[rowSums(observed)>0,]
  fullLambda <- rowSums(fullObserved)/total
  fullLhat <- fullLambda > 0
  fullExpected<- outer(fullLambda, laneTotals)

  #keep values
  fullKeep <- which(fullExpected > 0)

  #calculate degrees of freedom (nrow*(ncol -1) >> number of parameters - calculated (just lamda is calculated >> thats why minus 1)
  #calculate pearson and deviance for all values
  oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/ fullExpected[fullKeep] # pearson chisq test
  dfFull <- length(fullKeep) - sum(fullLhat!=0)

  sqrt(sum(oeFull)/dfFull)
}


## Represent tpms
rma <- tpms


## PCA
rmaPCA <- prcomp(t(rma))

eigs <- rmaPCA$sdev^2
#pctExplained = (eigs[1] / sum(eigs)) + (eigs[2] / sum(eigs))
pctExplained1 = paste0("PC1 (",format(eigs[1] / sum(eigs)*100,digits=4),"% of variance explained)")
pctExplained2 = paste0("PC2 (",format(eigs[2] / sum(eigs)*100,digits=4),"% of variance explained)")


## Correlation between samples
size <- ncol(counts)
indices <- as.matrix(data.frame(x=rep(1:size, each=size), y=rep(1:size, size)))
tabSERE <- rep(0,nrow(indices))
unrep <- indices[,1]<indices[,2]
indices <- indices[unrep, , drop=FALSE]
cl <- makeCluster(no_cores)
clusterExport(cl=cl, varlist=c("counts","SERE","indices"))
res <- parSapply(cl, 1:nrow(indices),function(i){
  SERE(counts[,indices[i,]])
})
stopCluster(cl)
tabSERE[unrep] <- res
dim(tabSERE) <- c(size, size)
tabSERE[upper.tri(tabSERE)] <- t(tabSERE)[upper.tri(tabSERE)]
colnames(tabSERE) <- rownames(tabSERE) <- colnames(counts)
corsamples <- signif(tabSERE,3)

if(!identical(groups,FALSE)){
  groupnames <- c(groupa,groupb)
  groups <- groupnames[groups+1]
  names(groups) <- samples
}else{
  groups <- NULL
}

suppressMessages(library(RJSplot))

extractJSON <- function(directory){
  html <- scan(file = paste0(directory,"/index.html"), what = character(0), sep = "\n", quiet = TRUE)
  json <- gsub("(<script type=\"application/json\" id=\"data\">)|(</script>)","",html[13])
  write(json,paste0(directory,".json"))
  unlink(directory, recursive=TRUE)
}

directory <- paste0(analysisname,"/QC_PCAplot_rjs")
suppressMessages(scatterplot_rjs(signif(rmaPCA$x[,1],3),signif(rmaPCA$x[,2],3),col=groups,xlab=pctExplained1,ylab=pctExplained2,id=row.names(rmaPCA$x),plot=FALSE,dir=directory))
extractJSON(directory)

directory <- paste0(analysisname,"/QC_boxplot_rjs")
suppressMessages(boxplot_rjs(rma,outline=FALSE,col=groups,xlab="",ylab="Expression Value",plot=FALSE,dir=directory))
extractJSON(directory)

directory <- paste0(analysisname,"/QC_ExpressedGenes_rjs")
suppressMessages(barplot_rjs(colSums(rma>1),xlab="",ylab="Expressed Genes",plot=FALSE,dir=directory))
extractJSON(directory)

if(is.null(groups)){
  metadata <- NULL
}else{
  metadata <- matrix(groups)
}
directory <- paste0(analysisname,"/QC_heatmap_rjs")
suppressMessages(heatmap_rjs(corsamples,metadata,color="Reds",plot=FALSE,dir=directory))
extractJSON(directory)

# PDF Report
create_report('QCreport')