# TODO: Add comment # # Author: E.Korsching 01-2024 ############################################################################### ## fetch the quantification from salmon and kallisto ## import tab separated columns #a <- kal.sal.import(path.n="/home/korschi/1/10quantRAW/salmonraw", type="s") #b <- kal.sal.import(path.n="/home/korschi/1/10quantRAW/kallistoraw", type="k") kal.sal.import <- function(path.n="", type="s", pattern="mat_s"){ # Kallisto/Salmon quantification import function # path.n: path name to the master folder, one above of the results folders # type: k: kallisto, s: salmon quantifier data # pattern: file name search element or NULL f1 <- function(x, nr, olen, ncc, fns, nflist, h.names){ # kal.sal.import - split table out <- vector("list",olen) names(out) <- h.names for(i in 1:olen){ tmp <- x[,seq(i,ncc,olen),drop=F] dimnames(tmp)[[2]] <- fns out[[i]] <- tmp } return(out) } # read folder list d.set <- dir(path=path.n, pattern=pattern) d.set.c <- gsub("_", ".", d.set, fixed=T) #f.set <- list.files(path="/home/korschi/1/10quantRAW/salmonraw", pattern="\\.R", all.files=F, full.names=F, recursive=F, ignore.case=F, include.dirs=F, no.. =F) nflist <- length(d.set) if(nflist==0){ stop("\n Folder is empty - stop") } ## import # check first tmp <- read.table(file=paste(path.n,d.set[1],if(type=="s"){"quant.sf"}else{"abundance.tsv"},sep="/"), header=T, sep="\t", dec=".", stringsAsFactors=F) # salmon: "Name" "Length" "EffectiveLength" "TPM" "NumReads" # kallisto: "target_id" "length" "eff_length" "est_counts" "tpm" # get column names h.names <- names(tmp) # prune column names h.names <- gsub("_", ".", h.names, fixed=T) # get row names and order names hsa.. Homo.. ENS.. # assure the right order for joining the columns -- assumption all results have the same ids and the same number of rows r.names.o <- order(tmp[,1], decreasing=T) nr <- nrow(tmp) nc <- ncol(tmp) # which cols # order columns if(type=="s"){ # salmon c.names.o <- c(2,3,5,4) }else{ # kallisto c.names.o <- c(2,3,4,5) } olen <- length(c.names.o) ncc <- nflist*olen collect1 <- data.frame(matrix(0,nr,ncc)) dimnames(collect1)[[1]] <- tmp[r.names.o, 1] k <- 1 for(i in 1:nflist){ tmp <- read.table(file=paste(path.n,d.set[i],if(type=="s"){"quant.sf"}else{"abundance.tsv"},sep="/"), header=T, sep="\t", dec=".", stringsAsFactors=F) tmp <- tmp[r.names.o,] # order rows collect1[,k:(k+olen-1)] <- tmp[,c.names.o] # oder columns k <- k+olen } out <- f1(collect1, nr=nr, olen=olen, ncc=ncc, fns=d.set.c, nflist=nflist, h.names=h.names[c.names.o]) return(out) } f.biomart <- function(x, IDT=1, filT=NULL, attrib_out=NULL){ # wrapper - biomart # x: identifier (see EBI) # IDT: simple i/o, give 1: ENST, 2: gene symbol, 3: ENSG, # only 4: multi i/o # filT: hgnc_symbol, external_gene_name, external_synonym, uniprot_gn_symbol, wikigene_name # attrib_out: ensembl_gene_id, entrezgene_id, hgnc_symbol [not more than app. 2 external] # # require('biomaRt'); mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl")) # attributePages(mart) listAttributes() # keytypes(mart) # to filter: ,pattern="" - ensembl_transcript_id entrezgene_id # columns(mart) if(IDT==4 & (is.null(filT) | is.null(attrib_out))){ stop("filT & attrib_out need one or more db string(s)") } require('biomaRt') mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl")) Sys.sleep(1) if(IDT==1){ filT <- "ensembl_transcript_id" attR <- c("ensembl_gene_id","ensembl_transcript_id","transcript_biotype","mirbase_id","hgnc_symbol") } if(IDT==2){ filT <- "hgnc_symbol" # external_gene_name external_synonym attR <- c("ensembl_gene_id","transcript_biotype","mirbase_id","hgnc_symbol") } if(IDT==3){ filT <- "ensembl_gene_id" attR <- c("ensembl_gene_id","hgnc_symbol") } if(IDT==4){ bm.list <- vector("list",1) k <- 1 for(i in filT){ bm.list[[k]] <- getBM(filters=i, attributes=attrib_out, values=x, mart=mart) k <- k+1 Sys.sleep(2) } return(bm.list) } bm.list <- getBM(filters=filT, attributes=attR, values=x, mart=mart) return(bm.list) } # bm.list <- bm.list[!duplicated(bm.list[,1]),] # unique ensembl_transcript_id # x3 <- merge(x2, bm.list, by.x="ensembl_transcript_id1", by.y="ensembl_transcript_id") # join # x3 <- x3[order(x3[,"padj"], decreasing=F), ] # order by padj # return(x3) #res01.anno <- f.biomart(row.names(s.nik.data), IDT=2) #res01.anno <- f.biomart(row.names(s.nik.data), IDT=4, # filT=c("hgnc_symbol", "external_gene_name", "external_synonym", "uniprot_gn_symbol", "wikigene_name"), # attrib_out=c("ensembl_gene_id", "hgnc_symbol")) expr.pattern.plot <- function(x, lim=NULL, z=NULL, out="", row.names=F, group.lines=NULL, avg=F, grlab=NULL){ # plot - expression pattern # x: data frame / matrix expression values, z: annotation data.frame with two columns ENSG,HGNC # lim: lower and upper limit of separation between the three expression ranges # group.lines: end position of each(!) group, row left to right, except last group (needed also for group based calculation) # avg: average the groups, grlab: if avg:T add group labels vector # out: (path)+file name, row.names: T/F: x dimnames if(avg & (is.null(group.lines) | is.null(grlab)) ){ stop("add between group.lines and group labels") } xcn <- dimnames(x)[[2]] xrn <- dimnames(x)[[1]] nr <- nrow(x) nc <- ncol(x) group.lines1 <- c(group.lines,nc) # add last group end y <- x # create output template insert.start <- function(x){ # expr.pattern.plot xlen <- length(x); x2len <- 2*xlen x.r <- vector("integer",x2len) k <- 1 for(i in 1:xlen){ x.r[k] <- if(i==1){ 1 }else{ x[(i-1)]+1 } x.r[k+1] <- x[i] k <- k+2 } return(x.r) } group.avg <- function(x, gr){ # expr.pattern.plot glen <- length(gr); glen2 <- glen/2 res <- vector("numeric",glen2) ir <- seq(1,glen,2) k <- 1 for(i in ir){ res[k] <- mean(x[gr[i]:gr[i+1]]) k <- k+1 } return(res) } limits.r <- function(x, lim){ # expr.pattern.plot ra <- range(x) radd <- ra[2]-ra[1] ou <- c(lo=(lim[1]*radd)+ra[1], up=(lim[2]*radd)+ra[1]) return(ou) } if(avg){ g1 <- insert.start(group.lines1) x1 <- t(apply(x, 1, group.avg, gr=g1)) dimnames(x1)[[2]] <- grlab }else{ x1 <- x } tmp1 <- t(apply(x1, 1, limits.r, lim=lim))# t(apply(x 1, quantile, probs=lim)) tmp1 <- tmp1[xrn,] if(!is.null(z)){ dimnames(z)[[1]] <- z[,1]; z <- z[xrn,]; z <- z[,2] } ll <- -0.15; mm <- 0; hh <- 0.15 for(i in 1:nr){ y[i , x1[i,]>tmp1[i,2] ] <- hh y[i , x1[i,]>=tmp1[i,1] & x1[i,]<=tmp1[i,2] ] <- mm y[i , x1[i,]