# TODO: Add comment
# 
# Author: E.Korsching  2026 version
###############################################################################



de.wrap <- function(x, groups, coll, out, wpath, test="nbW", cond1=c("condition","A","N"), stati="de_overview"){
	# DE wrapper function (2025-01)
	#  is only working with selcted experiments in DE analysis
	# x: one or more data name(s), vector with character names - even composed like "a.sam.mt$mir$sts"
	#    all will be processed with one to many DE analyses given in coll (needs that they all have the same column names available)
	#    for all parameters see input R script  e.g. 0data_structure_wrap.R
	# groups: column names building a group in vector joint together to a list
	# coll: group status and selected groups, list of vectors with two status names and two group names 
	# out: result names of comparisons per input element and comparison, vector
	# wpath: relative output path,   test: one of  DESeq > DESeqMean > nbW > LRT,  cond1: see de.fn()
	# stati: name of statistic output file
	xlen <- length(x)
	clen <- length(coll)
	olen <- length(out)
	if((clen*xlen)!=olen){ stop("length of parameter x,coll,out are not appropriate") }
	dir.create(file.path(wpath, test))
	wpath <- paste(wpath, test, sep="/")		# paste(format(Sys.time(), "0_%Y_%m_%d_%M%S_")
	fConn <- file(paste(wpath,"/",stati,".tsv",sep=""), open="w")
	writeLines("START", fConn)
	close(fConn)
	k <- 1	# results
	for(i in 1:xlen){	# data elements
		bl <- eval(parse(text=x[i]))
		if(is.null(bl)){
			stop(paste("problem with data name",x[i],sep=" "))
		}
		for(j in 1:clen){	# calculation(s)
			blx <- coll[[j]]
			bl1 <- groups[[blx[3]]]
			bl1len <- length(bl1)
			bl2 <- groups[[blx[4]]]
			bl2len <- length(bl2)
			dat <- de.pre(bl[,c(bl1,bl2)], wpath=wpath, sfile=stati) # preprocessing of subset
			
			len <- bl1len+bl2len
			a1 <- vector("character",len)
			a1[1:bl1len] <- blx[1]
			a1[(bl1len+1):len] <- blx[2]
			
			blout <- paste(out[k], test, sep=".")
			#if(j==5){ return(list(dat=dat,blx=blx,bl1=bl1,bl1len=bl1len,bl2=bl2,bl2len=bl2len,len=len,a1=a1,blout=blout)) }
			res <- de.fn(dat, a1=a1, wpath=wpath, wfile=blout, sfile=stati, anno=NULL, cond1=cond1, sign=F, test=test, ret=T)		# DE
			assign(blout, res, pos=1)
			k <- k+1
		}
	}
	return()
}


de.pre <- function(x, wpath, sfile){
	# customize and remove empty rows
	# x: data frame / matrix
	x <- as.matrix(x)
	dnl <- length(dimnames(x)[[2]])
	nrx <- nrow(x)
	fConn <- file(paste(wpath,"/",sfile,".tsv",sep=""), open="a")
	writeLines("## pre-customization", fConn)
	writeLines(paste("\nnumber of columns in ",dnl," rows ",nrx,sep=""), fConn)
	# check
	isn <- sum(is.na(x))
	writeLines(paste("\nis NA number         ",isn,sep=""), fConn)
	iszero <- sum(x==0)
	writeLines(paste("\nis ZERO number       ",iszero," / ",dnl*nrx,sep=""), fConn)
	# remove all '0' rows
	a <- apply(x,1,sum)
	x <- x[!(a==0),]
	nrx <- nrow(x)
	writeLines(paste("\nremaining rows       ",nrx,sep=""), fConn)
	# adjust
	x[x<1 & x>=0] <- as.integer(1)		# consider also DESeq in the 'poscounts' mode (in estimateSizeFactors) as an alternative (prod(x)^1/n of n positive numbers, gene wise)
	cadj <- sum(x<1 & x>=0)
	writeLines(paste("\nremaining <1 ?       ",cadj,sep=""), fConn)
	mode(x) <- "integer"
	writeLines(paste("\nstorage.mode of x ?  ",storage.mode(x),sep=""), fConn)
	rxx <- range(x)
	writeLines(paste("\nrange of x           ",rxx[1]," - ",rxx[2],sep=""), fConn)
	writeLines("", fConn)
	close(fConn)
	return(x)
}


# DE
de.fn <- function(countdata, a1, wpath, wfile, sfile, anno=NULL, cond1=NULL, sign=F, test="DESeq", ret=F){
	# DESeq function
	# countdata: data matrix, samples in columns
	# a1: column data - conditions, vector e.g. c("A","N","N","N","A","N","A","C", ...)
	# wpath+wfile: path & file name results, sfile: report file w/o extension
	# anno: annotation: e.g. miR names - needs the row names of the countdata file
	# cond1: the pairwise condition which should become a result: contrast: e.g. c("condition","A","N")  fc=A/N
	# test: standard: DESeq, other: DESeqMean, nbW, LRT see below
	# ret: return results or not
	coldata <- data.frame(		# numbers with as.factor()
			row.names=dimnames(countdata)[[2]],
			condition=as.factor(a1),
			type=as.factor(a1)
	)
	cat("\n")
	#print(coldata)
	cat("\n nrow coldata ",nrow(coldata))
	
	## QC on filtered data
	cat("\n countdata<0  ",sum(countdata[countdata<0]))
	cat("\n countdata==0 ",sum(countdata[countdata==0]))
#	return(list(countdata=countdata,coldata=coldata))
	
	## create deseq data structure for differential analysis
	# countdata:  rows: transcript IDs  cols: sample names
	#          DESeqDataSetFromMatrix
	dds.mat <- DESeqDataSetFromMatrix(
			countData = countdata,
			colData = coldata,
			design = ~ condition)
	#		design = ~ batch + condition)		# err -> Please read the vignette section 'Model matrix not full rank':  seems to be no work around
	
	dds.mat <- estimateSizeFactors(dds.mat)		# dds <- estimateSizeFactors(dds, type = "poscounts") -or- dds <- DESeq(dds, sfType="poscounts")
	
	# check properties
	cat("\n nrow dds.mat ",nrow(dds.mat),"\n")
	
	## differential
	if(test=="DESeq"){
		dds.mat <- DESeq(dds.mat)
	}else if(test=="DESeqMean"){
		dds.mat <- DESeq(dds.mat, fitType="mean")		# fitType mean for sparse miR data -- sometimes even this is crashing
	}else if(test=="nbW"){
		dds.mat <- estimateDispersionsGeneEst(dds.mat)
		dispersions(dds.mat) <- mcols(dds.mat)$dispGeneEst
		dds.mat <- nbinomWaldTest(dds.mat)
	}else if(test=="LRT"){
		dds.mat <- estimateDispersionsGeneEst(dds.mat)
		dispersions(dds.mat) <- mcols(dds.mat)$dispGeneEst
		dds.mat <- nbinomLRT(dds.mat, reduced= ~1)
	}
	
	## results 1
	# list names
	if(is.null(cond1)){
		cat("\n dds.mat obj ",resultsNames(dds.mat),"\n")
		cat("\n >>> choose a contrast argument and restart\n")
		return()
	}else{
		cat("\n chosen ",cond1,"\n")
	}
	
	# A-B
	res.A.N <- results(dds.mat, contrast=cond1)
	# order by adjusted p-value
	res.A.N <- as.data.frame(res.A.N)
	res.A.N[is.na(res.A.N$pvalue), "pvalue"] <- 1	# technically substitute NA by 1 for intersection
	res.A.N[is.na(res.A.N$padj), "padj"] <- 1	# technically substitute NA by 1 for intersection
	res.A.N <- res.A.N[order(res.A.N$pvalue), ]
	
	# number of significants
	yall <- nrow(res.A.N)
	y1true <- sum(res.A.N$pvalue<0.05)
	y2true <- sum(res.A.N$padj<0.05)
	cat("\n <0.05 :  FALSE  TRUE")
	cat("\n sign. p      ",yall-y1true," ",y1true)
	cat("\n sign. padj   ",yall-y2true," ",y2true)
	
	# write number of significants in a file
	fConn <- file(paste(wpath,"/",sfile,".tsv",sep=""), open="a")
	writeLines(paste(wfile, "pF", "pT", "padjF", "padjT", sep="\t"), fConn)
	writeLines(paste("", (yall-y1true), y1true, (yall-y2true), y2true, sep="\t"), fConn)
	writeLines("", fConn)
	close(fConn)
	
	# add additional annotations - needs row.names
	if(!is.null(anno)){
		res.A.N.df <- merge(as.data.frame(res.A.N), as.data.frame(anno), by="row.names", sort=F)
		row.names(res.A.N.df) <- res.A.N.df[,1]
		res.A.N.df <- res.A.N.df[,-1]
		res.A.N <- res.A.N.df
	}
	
	# create a  data+results  dataframe
	dcm <- as.data.frame(counts(dds.mat, normalized=T)-1)		# subtract the added 1 from de.pre()
	dcm <- round(dcm, 0)
	sti <- as.data.frame(res.A.N)
	#sti <- round(sti, 2)
	res.A.N.df <- merge(sti, dcm, by="row.names", sort=F)
	row.names(res.A.N.df) <- res.A.N.df[,1]
	res.A.N.df <- res.A.N.df[,-1]
	res.A.N.df <- res.A.N.df[order(res.A.N.df$padj), ]	# it is - but to assure
	# keep only sign.
	if(sign){
		res.A.N.df <- res.A.N.df[res.A.N.df$padj<0.05, ]
	}
	
	# write data DE used/adjusted counts
	write.table(x=res.A.N.df, file=paste(wpath,"/",wfile,".tsv",sep=""), sep="\t", dec=".", col.names=NA)		# NA means colnames with first blank
	cat("\n\n")
	if(ret){ return(res.A.N.df) }else{ return() }
}



# Extract the lines with number of significant factors from the statistic file
de.overview <- function(p, f, ret=F){
	# extract from wrap.de() report file the lines with number of significants
	# p: path name, f: file name with extension, ret: additionally return table to workspace
	con <- file(paste(p,"/",f,sep=""), "r")
	out <- "\tpF\tpT\tpadjF\tpadjT"
	while( TRUE ){
		line <- readLines(con, n=1)
		if(length(line)==0){ break }
		logi <- grepl(pattern="^.+\\tpF\\tpT\\tpadjF\\tpadjT", x=line, perl=T)
		if( logi ){
			txt <- regmatches(x=line, regexpr(pattern="(?<=^).+?(?=\t)", text=line, perl=T))
			line.2 <- readLines(con, n=1)
			out <- rbind(out, paste(txt, line.2, sep=""))
		}
	}
	close(con)
	write.table(out, paste(p,"/0res_",f,sep=""), append=F, quote=F, sep="", row.names=F, col.names=F)
	out <- read.table(paste(p,"/0res_",f,sep=""), header=T, sep="\t", quote="", stringsAsFactors=F)
	out <- out[,-1]
	if(ret){ return(out) }
}
# de.overview(p="res20250113/DESeq", f="de_overview_DESeq.tsv", ret=T)