# TODO: Add comment
# 
# Author: E.Korsching  2024
###############################################################################



#### customize sam file(s) - several functions

custom.sam01 <- function(sam, cigar.c=6, chr.c=3, pos.c=4, plot=F){
	# customize a sam files
	# remove unwanted CIGAR strings  and  unify same postion/chr by max(length)
	# import by read.table(..., header=F, sep="\t"), analyze the sam table
	# sam: input data frame,  cigar.c,chr.c,pos.c : CIGAR column, chr column, position column
	# plot: F,T: histogram
	nr <- nrow(sam)
	nc <- ncol(sam)
	cigar.uni <- unique(sam[,cigar.c])
	cigar.uni.len <- length(cigar.uni)
	# keep only '0*M' CIGARs
	sam1 <- sam[grepl(pattern="^[0-9]*M$", sam[,cigar.c], perl=T),,drop=F]
	sam1.nr <- nrow(sam1)
	# condense position identical entries
	chr.sp <- split(sam1, sam1[,chr.c])
	chr.sp.len <- length(chr.sp)
	cat("\n rows,columns",nr,nc," unique CIGARs",cigar.uni.len," valid CIGAR rows",sam1.nr," number of chr",chr.sp.len)
	k <- NULL
	for(j in 1:chr.sp.len){
		chr.tmp <- chr.sp[[j]]
		pos.uni <- split(chr.tmp, chr.tmp[,pos.c])
		pos.uni.len <- length(pos.uni)
		for(i in 1:pos.uni.len){
			pos.tmp <- pos.uni[[i]]#int
			tmp.2.nr <- nrow(pos.tmp)
			
			if(tmp.2.nr==1){
				k <- rbind(k, pos.tmp)
			}else{
				tmp.3 <- as.integer( sub(pattern="M", "", pos.tmp[,cigar.c], perl=T) )
				tmp.arr <- which(tmp.3==max(tmp.3),arr.ind=T)
				k <- rbind(k, pos.tmp[tmp.arr[1],])		# tmp.arr could be more than one number...
			}
		}
	}
	# remove additionally overlapping situations
	nr.new <- nrow(k)
	new.c <- nc+1
	k[,new.c] <- k[,cigar.c]
	k[,new.c] <- as.integer( sub(pattern="M", "", k[,new.c], perl=T) )
	k <- k[order(k[,chr.c],k[,pos.c],decreasing=F),]
	l <- NULL
	j <- 1
	m <- 0
	while(T){
		while(T){
			e1 <- k[j, chr.c]==k[(j+1), chr.c]	# pairwise search on overlap
			if(e1){
				e2 <- (k[j, pos.c] + k[j, new.c]) > (k[(j+1), pos.c] +1)	# we allow 2 base overlap
				if(e2){
					l <- c(l,(j+1))
					j <- j+2
				}else{
					j <- j+1
				}
			}else{
				j <- j+1
			}
			if(j>=nr.new | (j+1)>=nr.new){ break }
		}
		if(!is.null(l)){	# reduce k and start further overlap search
			k <- k[-l,]
			nr.new <- nrow(k)
			cat("\n loop stop at - reduction to",j,nr.new)
			l <- NULL
			j <- 1
			m <- m+1
		}else{
			break
		}
	}
	cat("\n overlap search cycle(s)",m)
	pira <- range(k[,new.c])
	cat("\n piRNA seq length range",pira,"\n")
	if(plot){
		source("../0functions/0general/hist.plot.color.R")
		hist.plot.color(k[,new.c], xcolor=F, y=NULL, bin.num=11, x.range=c(23.5,34.5), norm=F, n.factor=1, lty=3, lwd=1,
				xlab="bins of seq length", ylab="", xlim=c(23,35), y.max=NULL, x.at=c(24,28,34), y.at=NULL, ylog=F, bar.width=0.9, offset=0, digits=2,
				col.grad=c("red","blue"), col.f="blue", col.sh=NULL, col.b=NA, cex=1, h.title="piRNA length distribution", add=F, rC=F)
	}
	#k <- k[,-new.c]
	return(k)
}

#date()
#a1 <- custom.sam01(sam=a, cigar.c=6, chr.c=3, pos.c=4, plot=T)
#date()		# 2 min

#which(c(1,1,2,1)==max(c(1,1,2,1)),arr.ind=T)
#sum(sapply(a2,nrow))	# 57620




#### add BED file strand information to existing SAM table

bed.strand.2.sam <- function(bed, sam.tab, b.strand.c, b.name.c, s.name.c, s.chr.c, s.pos.c, sam.strand=T){
	# import BED file and add strand information to a existing SAM table
	# join by NAME
	# bed: path/file to BED file [1	10531	10559	piR-hsa-7566157	6	-]
	# sam.tab: SAM table in R env [piR-hsa-5213441, 16, 1, 514376, 255, 28M, *, 0, 0, CAGTCAGTTTCTCTTCTTGGATTCAGAA, *, NH:i:1, HI:i:1, NM:i:0, MD:Z:28]
	# b.strand.c: strand column in BED, b.name.c,s.name.c: name columns, s.chr,s.pos: chr and position
	# sam.strand: T: return an extended sam format, F: return only name and strand
	bed <- read.table(bed, header=F, sep="\t")
	cat("\n BED rows   in",nrow(bed))
	bed <- bed[!duplicated(bed[,b.name.c]), ]
	cat("\n BED rows !dup",nrow(bed))
	cat("\n SAM rows   in",nrow(sam.tab))
	sam.tab <- sam.tab[!duplicated(sam.tab[,s.name.c]), ]
	cat("\n SAM rows !dup",nrow(sam.tab))
	m.bed.sam <- merge(x=bed, y=sam.tab, by.x=b.name.c, by.y=s.name.c, all=F, sort=T)
	cat("\n BED-SAM  rows",nrow(m.bed.sam))
	if(sam.strand){
		m.bed.sam <- m.bed.sam[,c(1,7:20,6)]
		m.bed.sam <- m.bed.sam[order(m.bed.sam[,s.chr.c],m.bed.sam[,s.pos.c],decreasing=F), ]
	}else{
		m.bed.sam <- m.bed.sam[,c(1,6)]
	}
	return(m.bed.sam)
}

#aa <- bed.strand.2.sam(bed="/home/korschi/1/sources/RNAcentral/pirbase_hsa.hg19.bed",
#		sam.tab=a, b.strand.c=6, b.name.c=4, s.name.c=1, s.chr.c=3, s.pos.c=4, sam.strand=T)
# 1 min




#### create BED file from intermediate SAM table from custom.sam01()

bed.from.strand.sam <- function(outF, sam.tab, s.name.c, s.chr.c, s.pos.c, s.strand.c, s.len.c){
	# create BED and export based on custom.sam01()
	# outF: BED path/filename,  sam.tab: existing SAM table
	# s.name.c, s.chr.c, s.pos.c, s.strand.c, s.len.c: column numbers
	bed.end <- sam.tab[,s.pos.c] + sam.tab[,s.len.c] - 1
	bed <- cbind(sam.tab[,c(s.chr.c,s.pos.c)], bed.end, sam.tab[,s.name.c], ".", sam.tab[,s.strand.c])
	write.table(bed, file=outF, append=F, quote=F, sep="\t", row.names=F, col.names=F)
	return()
}

#bed.from.strand.sam(outF="/home/korschi/1/sources/RNAcentral/r_star_pirbase_gold.bed",
#		sam.tab=a2, s.name.c=1, s.chr=3, s.pos=4, s.strand.c=16, s.len.c=17)




#### create FASTA file from intermediate SAM table from custom.sam01()

fasta.from.strand.sam <- function(outF, sam.tab, s.name.c, s.chr.c, s.pos.c, s.seq.c){
	# create FASTA and export based on custom.sam01()
	# outF: FASTA path/filename,  sam.tab: existing SAM table
	# s.name.c, s.chr.c, s.pos.c, s.seq.c: column numbers
	nr <- nrow(sam.tab)
	fasta.name <- paste(">",sam.tab[,s.name.c]," ",sam.tab[,s.chr.c],":", sam.tab[,s.pos.c], sep="")
	fasta.seq <- sam.tab[,s.seq.c]
	erg <- matrix("",(nr*2),1)
	k <- 1
	for(i in 1:nr){
		erg[k,] <- fasta.name[i]
		k <- k+1
		erg[k,] <- fasta.seq[i]
		k <- k+1
	}
	write.table(erg, file=outF, append=F, quote=F, sep="", row.names=F, col.names=F)
	return()
}

#fasta.from.strand.sam(outF="/home/korschi/1/sources/RNAcentral/r_star_pirbase_gold.fasta",
#		sam.tab=a2, s.name.c=1, s.chr=3, s.pos=4, s.seq.c=10)





