# TODO: Add comment # # Author: E.Korsching, adjusted again 02-2024 ############################################################################### # import quantifier data of salmon, kallisto, idxstat, qualimap quant.import <- function(pathf="", pathn="", type="s", rmlast=F, pattern=NULL){ # import the result files of the read quantifier # (!) assumption: all results have the same ids, the same number of rows # pathf: path to the master folder, one above of the results folders # pathn: path to the ordered sample name file # type: output: k: kallisto, s: salmon, m: sam idxstat, q: qualimap count # rmlast: T,F: sam idxstat remove last line (other or non mapper) # pattern: filter on result folder names including this pattern (could be a regex term) # read folder list d.set <- list.dirs(path=pathf, full.names=F, recursive=F) if(!is.null(pattern)){ d.set <- grep(pattern, x=d.set, value=T) } nflist <- length(d.set) # number of samples if(nflist==0){ stop("\n Folder is empty") } # adjust names for R d.set.c <- gsub("_", ".", d.set, fixed=T) # note: the list.dirs() sample order (d.set, d.set.c) does not correpond # to the order of pathn (namelist.txt) file (sample.n) # therefore a mapping file needs to be created and applied sample.n <- scan(file=pathn, what=" ") # read names separated by space from file sample.n.l <- length(sample.n) if(sample.n.l==nflist){ # create mapping table o.str <- regmatches(d.set.c, regexpr("(\\d+)(?!.*\\d)", text=d.set.c, perl=T)) # extract last number block substr(x,14,18) #print(o.str); stop("o.str") o.vec <- order(as.integer(o.str), decreasing=F) #print(o.vec); stop("o.vec") d.set.co <- d.set.c[o.vec] # order d.set.c increasing s2n.map <- cbind(sample=d.set.co, name=sample.n) import.map <- cbind(sample=d.set.c, rn=d.set.c) f.map <- merge(x=import.map, y=s2n.map, by.x="rn", by.y="sample") # check if merge() is creating x order # col 3 is the new now corresponding name order }else{ stop(paste("sample number in namelist.txt ",sample.n.l," != name number in list.dirs ",nflist,sep="")) } assign("a.map", list(d.set.c=d.set.c, o.str=o.str, sample.n=sample.n, s2n.map=s2n.map, f.map=f.map), envir=.GlobalEnv) # cf. with namelist.txt #stop("see a.map") ## import types - check all h.names <- c("name","seq.length","count") if(type=="s"){ fname <- "quant.sf" c.names.o <- c(1,2,5) # salmon: "Name" "Length" "EffectiveLength" "TPM" "NumReads" get 1,2,5 h.l <- T }else if(type=="k"){ fname <- "abundance.tsv" c.names.o <- c(1,2,4) # kallisto: "target_id" "length" "eff_length" "est_counts" "tpm" get 1,2,4 h.l <- T }else if(type=="m"){ fname <- "mcounts.tsv" c.names.o <- c(1,2,3) # sam idxstat: "name" "length" "counts" "not_counted" get 1,2,3 h.l <- F }else if(type=="q"){ fname <- "qm.111HUmt.tsv" c.names.o <- c(1,2) # qualimap count: "name" "counts" get 1,2 h.l <- F } if(type=="s"|type=="k"){ for(i in 1:nflist){ pftmp <- paste(pathf, d.set[i], fname, sep="/") if(file.exists(pftmp) & file.size(pftmp)>2){ tmp <- read.table(file=pftmp, header=h.l, sep="\t", dec=".", stringsAsFactors=F) nr <- nrow(tmp) tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] # (!) ,1 r.names <- tmp[,1] r.names <- sub(pattern="\\s.*","",r.names,perl=T) # space group separated, keep first element seq.length <- tmp[,2] break }else{ warning("salmon/kallisto file might not exist or is empty") # quant.sf } } if(!exists("nr")){ cat("\nsalmon/kallisto all input files are not there or useless - so nr is NULL\n if you want to have an empty result structure use an adequate zero template file\n") return(NULL) } counts1 <- data.frame(matrix(0,nr,nflist)) for(i in 1:nflist){ pftmp <- paste(pathf, d.set[i], fname, sep="/") if(!file.exists(pftmp) | file.size(pftmp)<=2){ counts1[,i] <- rep(0,nr) warning(paste("salmon/kallisto file might not exist or is empty - created a 0 vector : ",pftmp,sep="")) }else{ tmp <- read.table(file=pftmp, header=h.l, sep="\t", dec=".", stringsAsFactors=F) print(nrow(tmp)) tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] counts1[,i] <- tmp[,3] } } dimnames(counts1)[[1]] <- r.names dimnames(counts1)[[2]] <- f.map[,3] return(list(counts=counts1, seq.length=seq.length)) }else if(type=="m"){ for(i in 1:nflist){ pftmp <- paste(pathf, d.set[i], fname, sep="/") if(file.exists(pftmp) & file.size(pftmp)>2){ tmp <- read.table(file=pftmp, header=h.l, sep="\t", dec=".", stringsAsFactors=F) nr <- nrow(tmp) if(rmlast){ tmp <- tmp[-nr,] } # remove last line (other or non mapper) tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] # (!) ,1 r.names <- tmp[,1] r.names <- sub(pattern="\\s.*","",r.names,perl=T) # space group separated, keep first element seq.length <- tmp[,2] break }else{ warning("samtools idxstats file might not exist or is empty") # mcount } } if(!exists("nr")){ cat("\nmcounts all input files are not there or useless - so nr is NULL\n if you want to have an empty result structure use an adequate zero template file\n") return(NULL) } counts1 <- data.frame(matrix(0, (nr-1), nflist)) for(i in 1:nflist){ pftmp <- paste(pathf, d.set[i], fname, sep="/") if(!file.exists(pftmp) | file.size(pftmp)<=2){ counts1[,i] <- rep(0,(nr-1)) warning(paste("samtools idxstats file might not exist or is empty - created a 0 vector : ",pftmp,sep="")) }else{ tmp <- read.table(file=paste(pathf, d.set[i], fname, sep="/"), header=h.l, sep="\t", dec=".", stringsAsFactors=F) if(rmlast){ tmp <- tmp[-nr,] } # remove last line (other or non mapper) print(nrow(tmp)) tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] counts1[,i] <- tmp[,3] } } #assign("a",list(counts1=counts1,rnames=r.names,samplen=sample.n),envir=.GlobalEnv) dimnames(counts1)[[1]] <- r.names dimnames(counts1)[[2]] <- f.map[,3] return(list(counts=counts1, seq.length=seq.length)) }else if(type=="q"){ # compensate for a qualimap empty bam error k <- 1 while (T) { tmp <- try( read.table(file=paste(pathf, d.set[k], fname, sep="/"), header=h.l, sep="\t", dec=".", stringsAsFactors=F) ) if(!inherits(tmp, 'try-error')){ nr <- nrow(tmp) break } k <- k+1 if(k>nflist){ stop("all qualimap associated bam files have : 0 rows") } } tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] # (!) ,1 r.names <- tmp[,1] r.names <- sub(pattern="\\s.*","",r.names,perl=T) # space group separated, keep first element counts1 <- data.frame(matrix(0,nr,nflist)) for(i in 1:nflist){ pftmp <- paste(pathf, d.set[i], fname, sep="/") tmp <- try( read.table(file=pftmp, header=h.l, sep="\t", dec=".", stringsAsFactors=F) ) if(!inherits(tmp, 'try-error')){ print(nrow(tmp)) tmp <- tmp[order(tmp[,1], decreasing=F), c.names.o] counts1[,i] <- tmp[,2] }else{ warning(paste("qualimap case 0 rows in ",pftmp,sep="")) counts1[,i] <- rep(0, nr) } } dimnames(counts1)[[1]] <- r.names dimnames(counts1)[[2]] <- f.map[,3] cat("\ncheck for case 0 rows...\n") return(list(counts=counts1)) }else{ stop("import type not found") } } #a <- quant.import(pathf="/home/korschi/1/13/a_0/03star1out", pathn="/home/korschi/1/13/a_0/namelist.txt", type="s", rmlast=F, pattern=NULL) #a <- quant.import(pathf="/home/korschi/1/13/a_0/03star1out", pathn="/home/korschi/1/13/a_0/namelist.txt", type="m", rmlast=T, pattern=NULL) #a <- quant.import(pathf="/home/korschi/on2/karo_redo2024/work_sr/03star1out", # pathn="/home/korschi/on2/karo_redo2024/work_sr/namelist.txt", # type="s", # rmlast=F, # pattern=NULL) #### translate our custom ENST/G numbers to HGNC # dimnames(mapTable.mir.trna)[[2]] trans.ENSTG.hgnc <- function(x, y, ens.y=1, hgnc.y=3){ # translate custom ENST/G to hgnc ids # x: ENST/G vector, y: map table (presumably: "ensembl_transcript_id","ensembl_gene_id","hgnc_symbol") # ens.y: column number corresponding to x, hgnc.y: result column number, check all xlen <- length(x) out <- vector("character",xlen) for(i in 1:xlen){ out[i] <- y[which(y[,ens.y] %in% x[i], arr.ind=T), hgnc.y] } return(out) } #aaa <- trans.ENSTG.hgnc(x=row.names(a1$counts), y=mapTable.mir.trna) # needs to be there: mapTable.mir.trna #write.table(mapTable.mir.trna,file="mapTable.mir.trna",quote=F,sep="\t",row.names=F,col.names=T) #### translate standard qualimap count table based on ENSG numbers by HGNC names etc. trans.ENSG2.hgnc <- function(x, ens="ensembl_gene_id", out=c("ensembl_gene_id","hgnc_symbol"), host=1){ # translate ENST/G to hgnc ids # x: ENST/G vector # ens: filter, "ensembl_gene_id" or "ensembl_transcript_id" corresponding to x # out: attributes xlen <- length(x) a1 <- f.biomart(x, filter=ens, attri=out, host=host) ui <- unique(a1[,1]) ui.l <- length(ui) # all genes returned? if(xlen!=ui.l){ cat("\n xlen,ui.l",xlen,ui.l) } a2 <- merge(x=data.frame(a=x), y=a1, by.x="a", by.y=ens) return(a2) } #a <- trans.ENSG2.hgnc(x=row.names(aa$counts)[-1], ens="ensembl_transcript_id", out=c("ensembl_transcript_id","hgnc_symbol"), host=1) #### there are 4 wrapper functions for quant.import() --- see also import.R script # # for 32_se_mir.sh # # k.raw.branch() miR/tRNA salmon quant (quant.sf) + idxstat quant (mcounts) # # k.raw.branch.2() premature_miR, pmm Ensembl 111 set, idxstat quant (mcounts) # # k.raw.branch.3() star mte (all RNA), bowtie chr view, bowtie piRNA, idxstat quant (mcounts) # # k.raw.branch.4() star mte (all RNA), salmon quant (quant.sf) #### import parameter for the isubM data frame [only for k.raw.branch()] # miR/tRNA # # folder eval fla column # "03star1out" "s" 1 1 # "03star1out" "m" 1 1 # "04bt21out" "m" 0 0 # "04bt22out" "q" 1 2 # "05bbm1out" "m" 0 0 # # 33_se_mir_ref_comparison.sh # modified shell script [based on 33_se_mir.sh] # # folder eval fla column # "03star2out" "s" 1 1 # "03star2out" "m" 1 1 ## example for 33_se_mir_ref_comparison.sh # #test.mt <- k.raw.branch(dpath="/home/korschi/on2/test2/MAFttt", # isubM=data.frame(folder=c("03star1out","03star1out","03star2out","03star2out","04bt21out","05bbm1out"), # eval=c("s","m","s","m","m","m"), # fla=c(1,1,1,1,0,0), # column=c(1,1,1,1,0,0)), # outp="res20260617test", outname="test_mt", trta=mapTable.mir.trna) # different miR output variants k.raw.branch <- function(dpath="", isubM=NULL, outp="", outname="", trta, test=F){ # import data and arrange in a custom way # export raw unsorted and sorted in a tsv file # dpath: data path # isubM: build on 3 evaluation categories: s (salmon), m (idxstat), q (QualiMap) # and on folders # provide a string data frame with 4 columns: folder, evaluation type. # outp: tsv session save path # outname: identifier for analysis # trta: mapping table custom ENST to hgnc (like) names # test: export all single tables and a joint version # initialization trna.r <- c(1:260) # used for split nl <- paste(dpath,"/namelist.txt", sep="") # work isubM.l <- nrow(isubM) m.l <- NULL # master list for(i in 1:isubM.l){ st.i <- paste(isubM[i, 1], isubM[i, 2], sep="_") cat("\n",st.i,"\n") xt <- quant.import(pathf=paste(dpath, isubM[i,1], sep="/"), pathn=nl, type=isubM[i,2], rmlast=if(isubM[i,2]=="m"){T}else{F}, pattern=NULL) if(!is.null(xt)){ m.l[[st.i]] <- k.raw.branch.sub1(x=xt, outp=outp, outname=outname, trna.r=trna.r, trta=trta, pp=st.i, fla=isubM[i, 3], ens.y=isubM[i, 4]) } } cat("\n") return(m.l) } k.raw.branch.sub1 <- function(x, outp, outname, trna.r, trta, pp, fla=0, ens.y=0){ # fla=1 for sts stm btq , fla=0 for btm bbm # focus on counts x <- x$counts # adjust col names #dimnames(x)[[2]] <- paste(pp,".", dimnames(x)[[2]], sep="") # translate names if(fla==1){ row.names(x) <- trans.ENSTG.hgnc(x=row.names(x), y=trta, ens.y=ens.y, hgnc.y=3) } # for fla=1 -> ens.y: sts stm: 1 btq: 2 # order all for split x <- x[order(row.names(x),decreasing=F),] # find a column which is not complete zero sortc <- fc.nonzero(x) # order & export # split mir and trna x.trna <- x[trna.r,] x.mir <- x[-trna.r,] # sort with first column, decreasing x.trna <- x.trna[order(x.trna[,sortc], decreasing=T),] x.mir <- x.mir[order(x.mir[,sortc], decreasing=T),] write.table(x.trna, file=paste(outp,"/",outname,"_",pp,"_tRNA.tsv",sep=""), append=F, quote=F, sep="\t", row.names=T,col.names=NA) write.table(x.mir, file=paste(outp,"/",outname,"_",pp,"_miR.tsv",sep=""), append=F, quote=F, sep="\t", row.names=T,col.names=NA) return(list(mir=x.mir,trna=x.trna)) } # k.raw.branch(dpath="/home/korschi/1/13/a_0", isubM=c("sts","stm"), outp="res20240331", outname="aa", trta=mapTable.mir.trna, test=F) #a1 <- "hsa-miR-23a-3p" #a$sts$counts[which(row.names(a$sts$counts)==a1,arr.ind=T),] #a$stm$counts[which(row.names(a$stm$counts)==a1,arr.ind=T),] #a$btm$counts[which(row.names(a$btm$counts)==a1,arr.ind=T),] #a$btq$counts[which(row.names(a$btq$counts)==a1,arr.ind=T),] #a$bbm$counts[which(row.names(a$bbm$counts)==a1,arr.ind=T),] #a$all[which(row.names(a$all)==a1,arr.ind=T),] fc.nonzero <- function(x){ # find a column which is not complete zero nc.c <- ncol(x) for(i in 1:nc.c){ tmp <- sum(x[,i]) if(tmp!=0){ sortc <- i break } if(i==nc.c){ sortc <- 1 } } return(sortc) } # premature_mir pmm k.raw.branch.2 <- function(dpath="", isubf="", outp="", outname="", host=1){ # for the imports of premature_mir pmm # isubf: here only one folder name source("../0functions/NGS/biom.R") source("../0functions/NGS/quant.import.R") nl <- paste(dpath,"/namelist.txt", sep="") #file.copy(paste(dpath,"/04bt22out", sep=""), paste(getwd(),"/",outp,"/",prefix,"_pmm.tsv",sep=""), overwrite=F) stpm <- quant.import(pathf=paste(dpath,"/",isubf, sep=""), pathn=nl, type="m", rmlast=T, pattern=NULL) stpm <- stpm$counts # find a column which is not complete zero sortc <- fc.nonzero(stpm) a <- trans.ENSG2.hgnc(x=row.names(stpm), ens="ensembl_transcript_id", out=c("ensembl_transcript_id","hgnc_symbol"), host=host) row.names(a) <- a[,1] a <- merge(a, stpm, by=0, all=T) a <- a[!duplicated(a[,3]),] row.names(a) <- a[,3] a <- a[,-c(1,2,3)] a <- a[order(a[,sortc],decreasing=T),] write.table(a,file=paste(outp,"/",outname,".tsv",sep=""),append=F,quote=F,sep="\t",row.names=T,col.names=NA) return(a) } # k.raw.branch.2(dpath="/home/korschi/1/13/a_0", isubf="03star2out", outp="res20240331", outname="b0_pmm") # star mte, bowtie chr view, bowtie piRNA mcounts k.raw.branch.3 <- function(dpath="", isubf="", outp="", outname=""){ # for the imports of star mte, bowtie chr view, bowtie piRNA # isubf: here only one folder source("../0functions/NGS/biom.R") source("../0functions/NGS/quant.import.R") nl <- paste(dpath,"/namelist.txt", sep="") a <- quant.import(pathf=paste(dpath,"/",isubf, sep=""), pathn=nl, type="m", rmlast=T, pattern=NULL) a <- a$counts # find a column which is not complete zero sortc <- fc.nonzero(a) a <- a[order(a[,sortc],decreasing=T),] write.table(a,file=paste(outp,"/",outname,".tsv",sep=""),append=F,quote=F,sep="\t",row.names=T,col.names=NA) return(a) } # k.raw.branch.3(dpath="/home/korschi/1/13/a_0", isubf="04bt22out", outp="res20240331", outname="b0_btch") # k.raw.branch.3(dpath="/home/korschi/1/13/b_0", isubf="04bt23out", outp="res20240331", outname="c0_btpi") # star mte salmon quant [not used at the moment, but an option] k.raw.branch.4 <- function(dpath="", isubf="", outp="", outname=""){ # for the imports of star mte (salmon quant) # isubf: here only one folder source("../0functions/NGS/biom.R") source("../0functions/NGS/quant.import.R") nl <- paste(dpath,"/namelist.txt", sep="") a <- quant.import(pathf=paste(dpath,"/",isubf, sep=""), pathn=nl, type="s", rmlast=T, pattern=NULL) a <- a$counts # find a column which is not complete zero sortc <- fc.nonzero(a) a <- a[order(a[,sortc],decreasing=T),] write.table(a,file=paste(outp,"/",outname,".tsv",sep=""),append=F,quote=F,sep="\t",row.names=T,col.names=NA) return(a) } # merge some imported and comparable data sets by common intersection merge.intersec.name <- function(x, col.num, rm=F){ # merge by common intersection, and ev. remove rows # x: list of data.frame names, col.num: column number for the merging procedure (in all elements identical) 0:row.names # rm: F/T: remove all zero rows xl <- length(x) if(xl==1){ stop("more than one data frame") } # xl=2 ### eval(parse(text="b.t1.mt$mir$stm")) # solution for get() list sub-elements a <- merge(x=eval(parse(text=x[1])), y=eval(parse(text=x[2])), by.x=col.num, by.y=col.num) row.names(a) <- a[,1] a <- a[,-1] if(xl>=3){ for(i in 3:xl){ a <- merge(x=a, y=eval(parse(text=x[i])), by.x=col.num, by.y=col.num) row.names(a) <- a[,1] a <- a[,-1] } } if(rm){ a <- a[apply(a,1,sum)>0, ] } a <- a[order(a[,1],decreasing=T), ] return(a) } #a6 <- merge.intersec.name(x=c("b.0.btpi","b.1.btpi","b.2.btpi"), col.num=0, rm=T) # check function p.quotient.n.c <- function(x){ # pairwise quotient of neighbouring columns # matrix, data.frame source("../0functions/0general/is.even.R") nr <- nrow(x) nc <- ncol(x) if(nc<2){ stop("number of columns <2") } if(is.even(nc)){ k <- 1 a <- matrix(0,nr,(nc/2)) snc <- seq(1,nc,2) for(i in snc){ a[,k] <- x[,i]/x[,(i+1)] k <- k+1 } }else{ stop("uneven number of columns") } return(a) } #p.quotient.n.c(cbind(c(1,2,2,1), c(4,5,4,5))) #p.quotient.n.c(cbind(c(1,2,2,1), c(4,5,4,5), c(4,5,4,5), c(5,5,5,5)))