# TODO: Add comment # # Author: E.Korsching Feb 11, 2014 - GPL v3 ############################################################################### ### create empty database structure db.create.structure <- function(DB.name="", DB.path=""){ # create a new and empty sqlite database structure with the given name and path # the database structure is hard coded in this function # if no path is given the current R work directory is used # if no name is given -> stop # ini require(RSQLite) work.dir <- getwd() if(DB.name==""){ cat("\n Please give a name - stop\n"); return() } if(DB.path==""){ DB.path.name <- paste(work.dir,DB.name,sep="/") }else{ DB.path.name <- paste(DB.path,DB.name,sep="/") } # If the named database does not yet exist, one is created if(file.exists(DB.path.name)){ # file.access(names, mode={0124}) cat("\n Database file is already existing - stop\n") return() }else{ db <- dbConnect(dbDriver("SQLite"), dbname=DB.path.name) cat("\n Database created in : ", DB.path.name, "\n") } # create database table dbSendQuery(conn=db, "CREATE TABLE sample ( sample_id INTEGER PRIMARY KEY AUTOINCREMENT, sample_description VARCHAR(40), batch INTEGER )" ) dbSendQuery(conn=db, "CREATE TABLE variations ( var_id INTEGER PRIMARY KEY AUTOINCREMENT, sample_id INTEGER not null, mutation INTEGER, chr INTEGER, position INTEGER, var VARCHAR(10), quality REAL, frequency REAL, coverage INTEGER, reads_var INTEGER, balance REAL, genes VARCHAR(40), type VARCHAR(40), known VARCHAR(40), prediction VARCHAR(40), maf VARCHAR(40), pcranno VARCHAR(40), FOREIGN KEY(sample_id) REFERENCES sample(sample_id) )" ) # The next three commands show: cat("\n\n Tables \n") print(dbListTables(db)) # The tables in the database cat("\n Table : sample : fields : \n") print(dbListFields(db, "sample")) # The columns in a table cat("\n Table : variations : fields : \n") print(dbListFields(db, "variations")) # The columns in a table # dbReadTable(db, "School") # The data in a table cat("\n") # close connection cat(" Database connection closed: ", dbDisconnect(db),"\n") return() } ### delete database file db.remove.total <- function(path.name=""){ if(path.name==""){ cat("\n please give a path-name of the database to be removed - stop\n"); return() } file.remove(path.name) } ### add data set db.add.data <- function(DB.path.name="", data=NULL, df.name="", group.type=0){ # append data to a sqlite database with the given name and path # and table variations -- this also effects the table sample # the database structure which is assumed to be resulted from function: db.create.structure() # data: one data.frame # df.name: if data should get another name # if no name / no data is given -> stop # group.type: can be a signed integer : -1, 0, 1, ... # at the moment: -1: uncertain, 0: normal, 1: mutation to be tested # dbClearResult(res, ...) for freeing memory # ini append <- T require(RSQLite) if(DB.path.name==""){ cat("\n Please give a name - stop\n"); return() } if(is.null(data)){ cat("\n Please provide some data - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } ## write data frame into relational database # get number of rows nr <- nrow(data) cat("\n nr ", nr, "\n") # get col one description if(df.name==""){ tmp <- deparse(substitute(data)) }else{ tmp <- df.name } # populate table sample res <- dbSendQuery(conn=db, paste("INSERT INTO sample (sample_description,batch) VALUES ('",tmp,"','",group.type,"')",sep="")) # get max index sample res <- dbSendQuery(conn=db, "SELECT max(sample_id) FROM sample") tmp <- fetch(res, n=-1) # returned data.frame with one col # add index in col one of data.frame tmpN <- rep(x=NA,times=nr) tmp <- rep(x=tmp[1,1],times=nr) data <- cbind(id=tmpN, sampleid=tmp, data) # append or overwrite variations table dbWriteTable(conn=db, name="variations", value=data, row.names=F, append=append) # close connection dbDisconnect(db) return() } ### insert column (at the end) db.insert.column <- function(DB.path.name="", table.name="", col.name="", data.type=""){ # insert column into the database # get relative database path/name if(DB.path.name==""){ cat("\n Please give a DB.path.name - stop\n"); return() } if(table.name==""){ cat("\n Please give a table.name - stop\n"); return() } if(col.name==""){ cat("\n Please give a col.name - stop\n"); return() } if(data.type==""){ cat("\n Please give a data.type - stop\n"); return() } # If the named database does exist, open it require(RSQLite) if(file.exists(paste(getwd(),DB.path.name,sep=""))){ db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # insert column res <- dbSendQuery(conn=db, paste( "ALTER TABLE ",table.name," ADD ",col.name," ",data.type, sep="") ) # close connection dbDisconnect(db) return() } ### insert bed PCR information during import insert.PCR.anno <- function(df.anno="", chr.loop="", position.loop="", genes.loop=""){ # insert PCR information during import # get chr, position and gene_symbol(s) and return pcr fragment name(s) # if more than one gene_symbol in genes.loop : separator : space char # filter first gene name pos.end <- regexpr(pattern=" ", text=genes.loop, fixed=T) if(pos.end>0){ #not found: -1 genes.loop <- substr(genes.loop, 1, pos.end-1) } # find in df.anno tmp <- get(df.anno, pos=1, inherits=F) tmp1 <- grepl( pattern=genes.loop, x=tmp[,"pcr.name"] ) # filter on gene tmp <- tmp[tmp1,,drop=F] tmp1 <- grepl( pattern=chr.loop, x=tmp[,"chr"] ) # filter on chr tmp2 <- tmp[tmp1,,drop=F] if(nrow(tmp2)!=0){ tmp3 <- tmp2[tmp2[,"start"]<=position.loop & tmp2[,"end"]>=position.loop , ,drop=F] # filter on hit range(s) tmp3.nr <- nrow(tmp3) if(tmp3.nr>0){ if(tmp3.nr==1){ # found one return(tmp3[,"pcr.name"]) }else{ # found more (ambiguous) # line up all found, separated by space tmpname <- "" for(i in 1:tmp3.nr){ if(i==tmp3.nr){ tmpname <- paste(tmpname,tmp3[i,"pcr.name"],sep="") }else{ tmpname <- paste(tmpname,tmp3[i,"pcr.name"]," ",sep="") } } return(tmpname) } } } # nothing found return("") } ### insert bed PCR information in the db - fill db completely new db.insert.PCR.anno <- function(DB.path.name="", df.anno=""){ # insert PCR information in the db - fill db completely (fill new) # database column must exist, anno must exist if(DB.path.name==""){ cat("\n give a relative path to the database \n"); return() } if(df.anno==""){ cat("\n give a data.frame name with annotation information \n"); return() } # we look for global variables only in the global environment if(!exists(x=df.anno, where=1, mode="any", inherits=F)){ cat("\n data.frame object: ",df.anno," not found in global workspace \n") return() } if(DB.path.name==""){ cat("\n give a relative path and name to the database file \n"); return() } DB.path.name <- paste(getwd(),DB.path.name,sep="") # If the named database does exist, open it require(RSQLite) if(file.exists(DB.path.name)){ db <- dbConnect(dbDriver("SQLite"), dbname=DB.path.name) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # get all columns of interest from database db.set <- dbGetQuery(conn=db, "SELECT var_id,chr,position,genes FROM variations" ) cat("\n var names: ",names(db.set),"\n mode: ",sapply(db.set, class),"\n class: ",class(db.set),"\n") # change names for R names(db.set) <- c("var.id","chr","position","genes") cat("\n var names new: ",names(db.set),"\n mode: ",sapply(db.set, class),"\n") db.set.nr <- nrow(db.set) cat("\n db.set.nr: ",db.set.nr) # process entry by entry k1 <- 0 # found counter k2 <- 0 # ambiguous counter k3 <- 0 # no counter for(i in 1:db.set.nr){ id.loop <- db.set[i,"var.id"] # ???? Error in res[i, "var.id"] : object of type 'S4' is not subsettable #not defined the variable res, before trying to subset it. res might also be a function or something else # in our case res was a return var in the following database queries chr.loop <- db.set[i,"chr"] position.loop <- db.set[i,"position"] genes.loop <- db.set[i,"genes"] # filter first gene name pos.end <- regexpr(pattern=" ", text=genes.loop, fixed=T) if(pos.end>0){ #not found: -1 genes.loop <- substr(genes.loop, 1, pos.end-1) } # find in df.anno tmp <- get(df.anno, pos=1, inherits=F) tmp1 <- grepl( pattern=genes.loop, x=tmp[,"pcr.name"] ) # filter on gene tmp <- tmp[tmp1,,drop=F] tmp1 <- grepl( pattern=chr.loop, x=tmp[,"chr"] ) # filter on chr tmp2 <- tmp[tmp1,,drop=F] if(nrow(tmp2)!=0){ tmp3 <- tmp2[tmp2[,"start"]<=position.loop & tmp2[,"end"]>=position.loop , ,drop=F] # filter on hit range(s) tmp3.nr <- nrow(tmp3) if(tmp3.nr>0){ if(tmp3.nr==1){ # write back into database - value db.out <- dbSendQuery(conn=db, paste("UPDATE variations SET pcranno='",tmp3[,"pcr.name"],"' WHERE var_id=",id.loop ,sep="")) k1 <- k1+1 }else{ # line up all found, separated by space tmpname <- "" for(i in 1:tmp3.nr){ if(i==tmp3.nr){ tmpname <- paste(tmpname,tmp3[i,"pcr.name"],sep="") }else{ tmpname <- paste(tmpname,tmp3[i,"pcr.name"]," ",sep="") } } # write back into database - ambiguous db.out <- dbSendQuery(conn=db, paste("UPDATE variations SET pcranno='",tmpname,"' WHERE var_id=",id.loop ,sep="")) k2 <- k2+1 } }else{ k3 <- k3+1 # no result } }else{ k3 <- k3+1 # no result } } # close connection dbDisconnect(db) # cat("\n processed rows: ",db.set.nr," pcr name exact: ",k1," pcr 1:n : ",k2," no pcr: ",k3,"\n\n") return() } ### get import number get.batch.import.number <- function(DB.path.name=""){ # get the last batch number # which is indicating a batch import process # and increase the batch by +1 # if it is the first time set the number to '0' if(DB.path.name==""){ cat("\n Please give a name - stop\n"); return() } # If the named database does exist, open it require(RSQLite) if(file.exists(paste(getwd(),DB.path.name,sep=""))){ db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # get max index sample and go back res <- dbSendQuery(conn=db, "SELECT max(batch) FROM sample") tmp <- unlist(fetch(res, n=-1)) # returned data.frame with one col # close connection dbDisconnect(db) if(is.na(tmp)){ return(1) # start with 1 (database empty) }else{ return(tmp+1) # increment by 1 } } ### get unique gene names in database get.unique.genes <- function(DB.path.name="", batch=NULL, first=T, out.path.name="/rel/path/name"){ # get all unique genes of a specific database # if batch=NULL : complete database is used, otherwise one or a vector of batch numbers is given # if out.path.name is NULL: names are returned to the terminal # assume: multiple entries in the gene field are separated by a --space-- character # look into: table: variations , variable: genes : # first=T: first gene in this field is used, F: all genes in this field are used if(DB.path.name==""){ cat("\n Please give a name - stop\n"); return() } # If the named database does exist, open it require(RSQLite) db.file <- paste(getwd(),DB.path.name,sep="") out.file <- paste(getwd(),out.path.name,sep="") if(file.exists(db.file)){ db <- dbConnect(dbDriver("SQLite"), dbname=db.file) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # get all gene entries if(is.null(batch)){ sql.1 <- paste("SELECT genes FROM variations" ,sep="") }else{ sql.1 <- paste("SELECT variations.genes FROM sample INNER JOIN variations ON sample.sample_id=variations.sample_id WHERE sample.batch IN (",paste(batch,collapse=","),")" ,sep="") } res <- dbGetQuery(conn=db, statement=sql.1) # close connection dbDisconnect(db) # select # filter on unique rows (reduce data amount) in the one col data.frame and create a vector res <- unlist(unique(res)) res <- gsub(pattern="\"", replacement="", x=res) #remove quotes res.len <- length(res) if(first){ # pick the first gene of each entry res2 <- vector(mode="character",length=res.len) for(i in 1:res.len){ res2[i] <- unlist(strsplit(res[i], split=" ", fixed=T))[1] } }else{ # pick all genes of each entry res2 <- vector(mode="character",length=1) for(i in 1:res.len){ res2 <- c(res2, unlist(strsplit(res[i], split=" ", fixed=T)) ) } res2 <- res2[-1] # remove first dummy entry } # filter again on unique entries (be specific) res2 <- unique(res2) # export or return values if(is.null(out.path.name)){ return(res2) }else{ write.table(x=data.frame(res2,stringsAsFactors=F), file=out.file, append=F, quote=F, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=F) return() } } ### set mutation(true)/polymorphism(SNP) flags - see coding_schemes.txt set.flags.variations <- function(DB.path.name="", table.field="", flag.value="", work.list=NULL){ # set flags/values for a certain field according to a work.list # DB.path.name: sqlite database with the given name and relative path # database table is a constant: 'variations' # table.field: a field from 'variations' # flag.value: the value which is inserted # work.list: a tab separated table with named columns: at least: chr(omosome), position, var(iation), further columns are ignored # give a relative path and a file name # ini require(RSQLite) if(DB.path.name==""){ cat("\n Please give a realtive database path/name - stop\n"); return() } if(table.field==""){ cat("\n Please give a database field - stop\n"); return() } if(flag.value==""){ cat("\n Please give a value for the flag - stop\n"); return() } if(is.null(work.list)){ cat("\n Please give a realtive path/name for the TAB separated text file - stop\n"); return() } db.file <- paste(getwd(),DB.path.name,sep="") work.file <- paste(getwd(),work.list,sep="") # read work.list workL <- read.table(file=work.file, header=T, sep="\t", quote="", dec=".", stringsAsFactors=F) workL.nr <- nrow(workL) # If the named database does exist, open it if(file.exists(db.file)){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=db.file) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } #work n1 <- 0 # counter success n2 <- 0 # counter entry not found n3 <- 0 # counter total found for(i in 1:workL.nr){ sql.st1 <- paste("SELECT * FROM variations WHERE chr = ",paste(workL[i,"chr"],sep=""), " AND position = ",paste(workL[i,"position"],sep=""), " AND var LIKE '",paste(workL[i,"var"],sep=""),"'", sep="") id.0 <- dbGetQuery(conn=db, sql.st1 ) id.0.len <- nrow(id.0) n3 <- n3+id.0.len cat("\n i ",i," ",paste(workL[i,"chr"],sep="")," ",paste(workL[i,"position"],sep="")," ",paste(workL[i,"var"],sep="")," found : ",id.0.len) if(id.0.len>0){ sql.st2 <- paste("UPDATE variations SET ",table.field," = ",paste(flag.value,sep=""), " WHERE chr = ",paste(workL[i,"chr"],sep=""), " AND position = ",paste(workL[i,"position"],sep=""), " AND var LIKE '",paste(workL[i,"var"],sep=""),"'", sep="") dbSendQuery(conn=db, sql.st2 ) n1 <- n1 +1 }else{ cat("\n i ",i," ",paste(workL[i,"chr"],sep="")," ",paste(workL[i,"pos"],sep="")," ",paste(workL[i,"var"],sep="")," missed") n2 <- n2 +1 } } # close connection dbDisconnect(db) # feedback cat("\n\n number of rows processed: ", workL.nr) cat("\n successful: ", n1) cat("\n missed: ", n2) cat("\n total found in database: ", n3,"\n") return() } ### get specific position db.get.all.specific.position <- function(DB.path.name="", batch="", position=NULL, var=NULL, var.show=F, out.path.name=""){ # get all occurances of a certain position / variant out of the # sqlite database with the given name and relative path # the database structure which is assumed resulted from function: db.create.structure() # if no name is given -> stop # dbClearResult(res, ...) for freeing memory # batch: selects certain sample group(s) (c("1","3")) or all ("") # position: position is mandatory, var: is optional, var list can be created first by: var.show=T # out.path.name: relativ path and prefix name # ini 1 require(RSQLite) if(DB.path.name==""){ cat("\n Please give a name - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } ## give insight into var of database if(var.show){ res <- dbGetQuery(conn=db, "SELECT var FROM variations" ) res <- unlist(unique(res)) cat("\n var types in database: \n") return(res) } # ini 2 if(is.null(position)){ cat("\n Please provide one genomic position - stop\n"); return() } if(is.null(var)){ cat("\n No var type given - only position is used\n") }else{ cat("\n var type given: ",var,"\n") } ## search for position in relational database # get data.frame with all selected database lines if(batch==""){ if(is.null(var)){ res <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE variations.position LIKE '",position,"'",sep="")) }else{ res <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE position LIKE '",position,"' AND var LIKE '",var,"'" ,sep="")) } }else{ if(is.null(var)){ res <- dbGetQuery(conn=db, paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id=variations.sample_id WHERE variations.position LIKE '",position,"' AND sample.batch IN (",paste(batch,collapse=","),")" ,sep="")) }else{ res <- dbGetQuery(conn=db, paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id=variations.sample_id WHERE variations.position LIKE '",position,"' AND variations.var LIKE '",var,"' AND sample.batch IN (",paste(batch,collapse=","),")" ,sep="")) } } res.nr <- nrow(res) if(res.nr!=0){ ## replace the sample_id by the sample_description links <- unique(res[,"sample_id"]) len.links <- length(links) res <- transform(res, sample_id=as.character(sample_id)) for(i in 1:len.links){ res2 <- dbGetQuery(conn=db, paste("SELECT * FROM sample WHERE sample_id LIKE '",links[i],"'",sep="")) res[res[,"sample_id"]==as.character(links[i]),"sample_id"] <- res2[,"sample_description"] } # close connection dbDisconnect(db) if(out.path.name!=""){ write.table(x=res, file=paste(getwd(), out.path.name,".",position,".",format(Sys.time(), "%Y%m%d%H%M"),".txt",sep=""), append=F, quote=T, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) cat("\n number of rows: ",res.nr,"\n") return() }else{ return(res) } }else{ cat("\n no results - stop \n") return() } } ### export complete table(s) db.export.sample <- function(DB.path.name="", export.path="", orderby="batch", unip=F, batch=""){ # get complete patient table from # sqlite database with the given name and path # export path with relativ path and name # and export in a csv format to disk # orderby: sample_id, sample_description, batch # unip: get all unique patients of a database # batch: restrict to certain batch parts e.g. c("1","3") or all ("") # ini require(RSQLite) if(DB.path.name==""){ cat("\n Please give a relative path and name - stop\n"); return() } if(export.path==""){ cat("\n Please give a relative export path and name - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) cat("\n Database open \n") }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } ## search for position in relational database # get data.frame with all selected database lines if(batch==""){ res <- dbGetQuery(conn=db, paste("SELECT * FROM sample ORDER BY ",orderby," COLLATE NOCASE ASC", sep="")) }else{ res <- dbGetQuery(conn=db, paste("SELECT * FROM sample WHERE sample.batch IN (",paste(batch,collapse=","),") ORDER BY ",orderby," COLLATE NOCASE ASC", sep="")) } # close connection cat("\n Database connection closed: ", dbDisconnect(db),"\n") if(unip){ # filter by patient # split description by dot in two parts (last and rest) res.nr <- nrow(res) res.tmp <- vector(mode="character",length=res.nr) for(i in 1:res.nr){ res.out <- unlist(strsplit(x=res[i,2], split=".", fixed=T)) res.out.len <- length(res.out) res.tmp[i] <- res.out[res.out.len] # take last element } # order res.tmp <- res.tmp[order(res.tmp)] # filter res <- res.tmp[!duplicated(res.tmp)] # res <- res.tmp } # export in TAB separated table write.table(x=res, file=paste(getwd(),export.path,sep=""), append=F, quote=F, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) # cat("\n Table exported to /..",export.path, "\n") return() } db.export.sample.variations <- function(DB.path.name="", export.path.name=""){ # get complete patient table from # sqlite database with the given name and path # export path with relativ path and name # and export in a csv format to disk # ini require(RSQLite) if(DB.path.name==""){ cat("\n Please give a relative path and name - stop\n"); return() } if(export.path.name==""){ cat("\n Please give a relative export path and name - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) cat("\n Database open ") }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } ## search for position in relational database # get data.frame with all selected database lines res <- dbGetQuery(conn=db, "SELECT * FROM sample INNER JOIN variations ON sample.sample_id = variations.sample_id ORDER BY sample_id COLLATE NOCASE ASC") # close connection cat("\n Database connection closed: ", dbDisconnect(db),"\n") # export in TAB separated table write.table(x=res, file=paste(getwd(),export.path.name,sep=""), append=F, quote=T, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) # cat("\n Table exported to /..",export.path.name, "\n") return() } ### import / overwrite database table (! only useful after export and modification of such a table) db.overwrite.table.sample <- function(DB.path.name="", data.path.name=""){ # overwrite database table sample with TAB separated data file # first export database table into file, then modify # then import - without ! header row and row numbers # NA: empty cells ! # after import check data.formats in the database etc. # data: path and file name with TAB separated data # ini require(RSQLite) if(DB.path.name==""){ cat("\n Please give relative path and name - stop\n"); return() } if(data.path.name==""){ cat("\n Please provide relative path/file name with data - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) cat("\n Database open \n") }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } if(!file.exists(paste(getwd(),data.path.name,sep=""))){ cat("\n File does not exist / or wrong name/path - stop \n") return() } ## write data frame into relational database # delete content res <- dbSendQuery(conn=db, "DELETE FROM sample") # import complete csv tables via R into the database dbWriteTable(conn=db, name="sample", value=paste(getwd(),data.path.name,sep=""), sep="\t", row.names=F, header=F, overwrite=F, append=T) # close connection cat("\n Database connection closed: ", dbDisconnect(db),"\n") return() } db.overwrite.table.variations <- function(DB.path.name="", data.path.name=""){ # overwrite database table sample with TAB separated data file # first export database table into file, then modify # then import - without ! header row and row numbers # NA: empty cells ! # after import check data.formats in the database etc. # data: path and file name with TAB separated data # ini require(RSQLite) if(DB.path.name==""){ cat("\n Please give a relative path and name - stop\n"); return() } if(data.path.name==""){ cat("\n Please provide relative path/file name with data - stop\n"); return() } # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ # file.access(names, mode={0124}) db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) cat("\n Database open \n") }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } if(!file.exists(paste(getwd(),data.path.name,sep=""))){ cat("\n File does not exist / or wrong name/path - stop \n") return() } ## write data frame into relational database # delete content res <- dbSendQuery(conn=db, "DELETE FROM variations") # import complete csv tables via R into the database dbWriteTable(conn=db, name="variations", value=paste(getwd(),data.path.name,sep=""), sep="\t", row.names=F, header=F, overwrite=F, append=T) # close connection cat("\n Database connection closed: ", dbDisconnect(db),"\n") return() } ### statistics on database # function to generate statistics on frequency db.stat.frequency <- function(DB.path.name="", batch=NULL, sampleID=NULL, what1="", geneToWork=NULL, out.path.name="", workspace.suffix=""){ # statistics on frequency in 'good' group (what1="0": no mutation in variations table) # DB.path.name: sqlite database with the given name and relative path # batch : NULL : complete database will be analyzed or vector of batch numbers # sampleID : NULL : no sample_id used, else also sample_id s used # what1 : group to be used (mutation score) : variations.mutation # if not given the whole table will be used according to var type # geneToWork: gene name (or vector of names) to work on # out.path.name: relative path and filename (prefix) for PDF output # workspace.suffix: name of the resulting workspace variables stat. , freq. , bala. # ini require(RSQLite) require(colorspace) source("/home/korschi/eclipseR/0functions/0general/hist.plot.simple.R") source("/home/korschi/eclipseR/0functions/0general/trim.distribution.R") if(DB.path.name==""){ cat("\n Please give a database name and rel. path - stop\n"); return() } if(out.path.name==""){ cat("\n Please give a report name and rel. path - stop\n"); return() } if(workspace.suffix==""){ cat("\n Please give a suffix for the result variables - stop\n"); return() } if(is.null(geneToWork)){ cat("\n No gene given \n") return() } # set trim parameter (coherent to db.test.dist.basic.dist.frequency() !) trim.d <- "upper" trim.v <- 0.05 # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # get number of rows in variations t.rows <- dbGetQuery(conn=db, "SELECT Count(*) FROM variations" ) ## make groups, extract ids, create results # establish group if(is.null(batch)){ # we take all batch groups if(what1==""){ tmp.out <- paste("\n aa. we take the whole ",sep="") sql.st <- "SELECT * FROM variations" }else{ tmp.out <- paste("\n ab. we take mutation : ",what1,sep="") sql.st <- paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE variations.mutation LIKE '",what1,"' AND sample.batch = 0" , sep="") } }else{ # we take some batch groups and/or sample_id s if(what1==""){ if(!is.null(sampleID)){ tmp.out <- paste("\n ba. we take the whole ",sep="") sql.st <- paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE sample.batch IN (",paste(batch,collapse=","),") OR sample.sample_id IN (",paste(sampleID,collapse=","),")", sep="") }else{ tmp.out <- paste("\n bb. we take the whole ",sep="") sql.st <- paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE sample.batch IN (",paste(batch,collapse=","),")" , sep="") } }else{ if(!is.null(sampleID)){ tmp.out <- paste("\n bc. we take mutation : ",what1,sep="") sql.st <- paste("SELECT variations.* INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE variations.mutation LIKE '",what1,"' AND sample.batch IN (",paste(batch,collapse=","),")" , sep="") }else{ tmp.out <- paste("\n bd. we take mutation : ",what1,sep="") sql.st <- paste("SELECT variations.* INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE variations.mutation LIKE '",what1,"' AND ( sample.batch IN (",paste(batch,collapse=","),") OR sample.sample_id IN (",paste(sampleID,collapse=","),") )", sep="") } } } id.0 <- dbGetQuery(conn=db, sql.st ) id.0.nr <- nrow(id.0) cat("\n nr ",id.0.nr) # close db dbDisconnect(db) # results table stat <- data.frame(matrix(NA,1,16)) # with one dummy row names(stat) <- c("chr","position","var","mutation","genes","n","freq.mean","freq.median","freq.sd","freq.t.mean","freq.t.sd","balance.mean","freq.max","freq.min","bal.max","bal.min") # calculations # output to PDF pdf(file=paste(getwd(),out.path.name,".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title="statistics", pointsize=12) par(mfrow=c(2,2)) k <- 0 # unique row counter kk <- 0 # all row counter # create two separate lists containing all frequency distributions & names AND balance distributions & names freq.l <- list(NULL) freq.n <- NULL bala.l <- list(NULL) bala.n <- NULL for(i in geneToWork){ tmp <- grepl(pattern=i, x=id.0[,"genes"]) # caveat: check if grep is creating non overlapping sets id.1 <- id.0[tmp,] # sub set cat("\n geneToWork ",i," nr : ",nrow(id.1)) # unique var var.uni <- unique(id.1[,"var"]) var.uni <- var.uni[order(var.uni)] # cat("\n var.uni ",var.uni) # strings from database always own "quotation marks" for(j in var.uni){ cat("\n j ", j) tmp1 <- grepl(pattern=paste("^",j,"$",sep=""), x=id.1[,"var"]) # find exact (but check) id.2 <- id.1[tmp1,] # sub set # unique pos pos.uni <- unique(id.2[,"position"]) pos.uni <- pos.uni[order(pos.uni)] # cat("\n pos.uni ",pos.uni) for(l in pos.uni){ tmp2 <- id.2[id.2[,"position"]==l, ] # save frequency freq.l[[length(freq.l)+1]] <- tmp2[,"frequency"] freq.n <- c(freq.n, paste(i, gsub(pattern="\"", replacement="", x=j), l, sep=".")) # j: strip quotation marks # save balance bala.l[[length(bala.l)+1]] <- tmp2[,"balance"] bala.n <- c(bala.n, paste(i, gsub(pattern="\"", replacement="", x=j), l, sep=".")) # create information for stat table and graphics if(sum(tmp2[,"mutation"])>0){ mut <- 1 }else{ mut <- 0 } # check if one entry is manually flagged as mutation tmp2.nr <- nrow(tmp2) k <- k+1 kk <- kk + tmp2.nr cat("\n n : ", tmp2.nr," unique row counter : ",k," all row counter : ",kk) if(tmp2.nr==0){ next } # if no result skip output # save some information tmp.df <- stat[1,,drop=F] # create data.frame tmp.df[1, "chr"] <- tmp2[1,"chr"] tmp.df[1, "position"] <- tmp2[1,"position"] tmp.df[1, "var"] <- tmp2[1,"var"] tmp.df[1, "mutation"] <- mut tmp.df[1, "genes"] <- tmp2[1,"genes"] tmp.df[1, "n"] <- tmp2.nr # calculate some values tmp.df[1, "freq.mean"] <- mean( tmp2[,"frequency"] ) tmp.df[1, "freq.sd"] <- sd( tmp2[,"frequency"] ) tmp.df[1, "freq.t.mean"] <- mean( trim.distribution(x=tmp2[,"frequency"], fraction=trim.v, where=trim.d) ) # means no trimming in small sets ! tmp.df[1, "freq.t.sd"] <- sd( trim.distribution(x=tmp2[,"frequency"], fraction=trim.v, where=trim.d) ) tmp.df[1, "freq.median"] <- median( tmp2[,"frequency"] ) tmp.df[1, "balance.mean"] <- mean( tmp2[,"balance"] ) tmp.df[1, "freq.max"] <- max( tmp2[,"frequency"] ) tmp.df[1, "freq.min"] <- min( tmp2[,"frequency"] ) tmp.df[1, "bal.max"] <- max( tmp2[,"balance"] ) tmp.df[1, "bal.min"] <- min( tmp2[,"balance"] ) stat <- rbind(stat,tmp.df) # draw 2 text blocks & 2 graphs plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") mtext(text=paste("Korsching",format(Sys.time(), "%a %d %b %Y"),sep=", "), side=3, line=1, adj=0, cex=0.5) text(x=0.5, y=8, labels=paste("frequency sd: ",round(tmp.df[1, "freq.sd"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=7, labels=paste("frequency mean (trim ",trim.d," ",trim.v,"): ",round(tmp.df[1, "freq.t.mean"],3)," sd: ",round(tmp.df[1, "freq.t.sd"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=6, labels=paste("frequency median: ",round(tmp.df[1, "freq.median"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=2, labels=paste("frequency min/max : ",round(tmp.df[1, "freq.min"],3)," - ",round(tmp.df[1, "freq.max"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=1, labels=paste("balance min/max : ",round(tmp.df[1, "bal.min"],3)," - ",round(tmp.df[1, "bal.max"],3),sep=""), adj=c(0,0.5), col="blue") mtext(text="histogram", side=1, line=1) plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") text(x=0.5, y=9, labels=paste("row #",k,sep=" "), adj=c(0,0.5)) text(x=0.5, y=5, labels=paste("var: ",tmp.df[1, "var"],", chr",tmp.df[1, "chr"],"-",tmp.df[1, "position"], sep=" "), adj=c(0,0.5), col="blue") text(x=0.5, y=4, labels=paste(" gene: ",tmp.df[1, "genes"],sep=" "), adj=c(0,0.5), col="blue") text(x=0.5, y=3, labels=paste(" mean: frequency: ",round(tmp.df[1, "freq.mean"],3)," balance: ",round(tmp.df[1, "balance.mean"],3),sep=""), adj=c(0,0.5), col="blue") mtext(text="histogram", side=1, line=1) # mtext(text="density", side=1, line=1) # if(k==933){ dev.off(); return(pos.uni) } hist.plot.simple(x=tmp2[,"frequency"], xcolor=T, bin.num=F, bin.size=F, bin.fix=NULL, bin.fix.open=F, norm=F, n.factor=1, curve=F, xlab="frequency", ylab="counts/bin", xlim=c(0,1), y.max=NULL, x.at=c(0,0.25,0.5,0.75,1), y.at=NULL, y.num.format=F, ylog=F, bar.width=1, offset=0, digits=1, col=c("blue","red","green"), cex=1, smooth.f=0, subfolder.filename=NULL, h.title="", s.title=T, add=F) # if(tmp2.nr>=20){ # density distribution # plot( density(x=tmp2[,"frequency"], bw="SJ", kernel="gaussian"), main="" ) # }else{ # plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") # empty plot # text(x=0.5, y=5, labels="no density plot : count < 20", adj=c(0,0.5), cex=1.2) # } hist.plot.simple(x=tmp2[,"balance"], y=tmp2[,"frequency"], bin.num=F, bin.size=F, bin.fix=NULL, bin.fix.open=F, norm=F, n.factor=1, curve=F, xlab="balance", ylab="counts/bin", xlim=c(0,1), y.max=NULL, x.at=c(0,0.25,0.5,0.75,1), y.at=NULL, y.num.format=F, ylog=F, bar.width=1, offset=0, digits=1, col=c("blue","red","green"), cex=1, smooth.f=0, subfolder.filename=NULL, h.title="", s.title=T, add=F) } } } # remove first dummy row nrs <- nrow(stat) if(nrs>1){ stat <- stat[-1,] } # remove first list element and join list vectors and names # freq.l <- freq.l[!sapply(freq.l, is.null)] # more elemnts are NULL ! freq.l <- freq.l[-1] names(freq.l) <- freq.n bala.l <- bala.l[-1] names(bala.l) <- bala.n # sort according to chr , position , variant stat <- stat[order(stat[,1],stat[,2],stat[,3]), ] nrs <- nrow(stat) tmp.path <- paste(getwd(), out.path.name,".",what1,".",format(Sys.time(), "%Y%m%d%H%M"),".txt",sep="") # save results write.table(x=stat, file=tmp.path, append=F, quote=F, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) assign(paste("freq",workspace.suffix,sep="."), value=freq.l, pos=1) assign(paste("bala",workspace.suffix,sep="."), value=bala.l, pos=1) assign(paste("stat",workspace.suffix,sep="."), value=stat, pos=1) # report to pdf par(mfrow=c(1,1)) plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") text(x=0.5, y=9, labels=paste( "\n total number of lines in table : ", unlist(t.rows), "\n number of lines in selection : ", id.0.nr, "\n statistics - number of rows : ", nrs,"\n", sep=""), adj=c(0,0.5)) dev.off() # report to console cat("\n out path matrix : ", tmp.path, "\n",tmp.out, "\n",sql.st, "\n total number of lines in table : ", unlist(t.rows), "\n number of lines in selection : ", id.0.nr, "\n statistics - number of rows : ", nrs, "\n" ) return() } # test how many times a row appears how.many.dupli <- function(x, ref=NULL){ # calculating how many occurrences of a row (or part of a row) in x are appearing in a reference set (ref) # add a column with counts to x # if only x is given the dublicates in x will be calculated if(is.null(ref)){ ref <- x x <- x[!duplicated(x), ] } ncx <- ncol(x) ncr <- ncol(ref) nrx <- nrow(x) nrr <- nrow(ref) if(ncx!=ncr){ cat("\n we need the same number of columns in both parameters \n"); return() } cat("\n x - number of rows ",nrx,"\n ref - number of rows ",nrr,"\n") # make a string of all if more than one column included if(ncx>1){ x.collapse <- apply(x,1,paste,collapse="") ref.collapse <- apply(ref,1,paste,collapse="") }else{ x.collapse <- x ref.collapse <- ref } # work tmp <- vector(mode="numeric", length=nrx) for(i in 1:nrx){ tmp[i] <- sum( grepl(pattern=paste("^",x.collapse[i],"$",sep=""), x=ref.collapse, fixed=F) ) # find exact } res <- cbind(x,counts=tmp) return(res) } #aa <- how.many.dupli(x=test.data01.b[,c("chr","pos","var")], ref=test.data01[,c("chr","pos","var")]) #aa <- how.many.dupli(x=test.data01[,c("chr","pos","var")]) #x - number of rows 336 #ref - number of rows 5991 ### select new mutations # a priori do in terminal: sudo R # Load the BSgenome package and download the hg19 reference sequence # For documentation see http://www.bioconductor.org/packages/release/bioc/html/BSgenome.html # update R core packages via apt and all other packages: update.packages(ask = FALSE, dependencies = c('Suggests')) # bioconductor via: biocLite("BiocUpgrade") #source("http://www.bioconductor.org/biocLite.R") #biocLite("BSgenome") #library #biocLite("BSgenome.Hsapiens.UCSC.hg19") # human genome (~850 MB). #biocLite("BSgenome.Hsapiens.NCBI.GRCh38") # more recent (~807 MB) # function to fetch flanking sequence -- see also below get.flank <- function(position=102741896, chr="chr12", alleles="[C/A]", offset=10){ # by GGD Team Stephen Turner # http://gettinggeneticsdone.blogspot.de/2011/04/using-rstats-bioconductor-to-get.html # offset: length of flanking sequence left/right # alleles: default : [N/N] or what you prefer require(BSgenome) require(BSgenome.Hsapiens.UCSC.hg19) # require(BSgenome.Hsapiens.NCBI.GRCh38) leftflank <- getSeq(Hsapiens,chr,position-offset,position-1) rightflank <- getSeq(Hsapiens,chr,position+1,position+offset) return(paste(leftflank,alleles,rightflank,sep="")) } # get.flank(position=102741896, chr="chr12", alleles="[C/A]", offset=10) is returning "TGGCAACTCC[C/A]TTCCATTTGC" # get.flank(position=102741896, chr='12', alleles="[C/A]", offset=10) is returning "TAAGCCGTCA[C/A]ATCTAAGTTA" # which is exactly what you'd see if you searched db SNP: # http://www.ncbi.nlm.nih.gov/sites/entrez?db=snp&cmd=search&term=rs1520218 #getSeq(Hsapiens,"chr1",100000,100020) # detach("package:BSgenome.Hsapiens.NCBI.GRCh38", unload=TRUE) # CACTAAGCACACAGAGAATAA #getSeq(Hsapiens,"1",100000,100020) # CACTAAGCACACAGAGAATAA get.flank2 <- function(position=102741896, chr="12", alleles="[C/A]", offset=10){ ################ not functional -maybe never # by biomaRt getSequence function # offset: length of flanking sequence left/right # alleles: default : [N/N] or what you prefer # getSequence(): wrapper of getBM() require(biomaRt) # set data connection bioM.ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl", host="www.ensembl.org", verbose=F) f_value <- paste(chr, ":", position-offset, ":", position-1, sep="") leftflank <- getBM(attributes=c("cdna"), filters="chromosomal_region", values=f_value, mart=bioM.ensembl) # one or more rows f_value <- paste(chr, ":", position+1, ":", position+offset, sep="") rightflank <- getBM(attributes=c("cdna"), filters="chromosomal_region", values=f_value, mart=bioM.ensembl) return(paste(leftflank,alleles,rightflank,sep="")) } # get.flank2(position=102741896, chr="12", alleles="[C/A]", offset=10) is returning "TGGCAACTCC[C/A]TTCCATTTGC" # function to get corresponding chromosome identifier (very slow, if not changed to EBI server! see below) get.chr.name <- function(provider="", db.name="", gene.symbol=""){ # function to get corresponding chromosome identifier from gene name # provider: name of database provider, db.name: name of database, gene.symbol: to look for: one or vector (max. 500) # use: # get.chr.name(provider="ensembl", db.name="hsapiens_gene_ensembl", gene.symbol="BRCA1") # biomaRt query require(biomaRt) # set data connection bioM.ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl", host="www.ensembl.org", verbose=F) # send query (max. 500) & retrieve results result1 <- getBM(attributes=c("chromosome_name","hgnc_symbol"), filters="hgnc_symbol", values=gene.symbol, mart=bioM.ensembl) # reduce to unique gene_symbols (check if useful!) result2 <- result1[!duplicated(result1[,"hgnc_symbol"]), ] # return(result2) } # get list of all known database sites #bioM.db.present <- listMarts() # set data connection #bioM.ensembl <- useMart("ensembl") # show all available datasets #bioM.ensembl.dset <- listDatasets(mart=bioM.ensembl) # inform on filters #bioM.ensembl.f <- listFilters(bioM.ensembl) # -- the data and type we give to get something e.g. hgnc_symbol # inform on attributes #bioM.ensembl.a <- listAttributes(bioM.ensembl) # -- the corresponding data we want to get back e.g. entrezgene, hgnc_symbol # fragment: ## if window not 0 then extend the basic table by all those positions which exist +- window #if(window.size!=0){ # not functional : now array - adjust ! # vec.pos <- vector(mode="numeric",length=(2*window.size+1)*nr.uniPos) # # get one position and return this position plus all extended positions # f1 <- function(x, wsize){ # is working:wsize=0 # vec.len <- 2*wsize+1 # vec.pos <- vector(mode="numeric",length=vec.len) # k <- -wsize # for(i in 1:vec.len){ # vec.pos[i] <- x + k # k <- k+1 # } # return(vec.pos) # } # vec.pos <- as.vector(sapply(uniPositions, f1, wsize=window.size)) # # in the case that there is an overlap # uniPositions <- unique(vec.pos) #} # sub-function fill.ref.set <- function(db=db, basic.set=basic.set, id.1=id.1){ # sub-function to db.test.on.mutations - create the ref set for(i in id.1[,"sample_id"]){ id.1.one <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE sample_id LIKE '",i,"'" ,sep="") ) basic.set <- rbind(basic.set, id.1.one) } # remove first dummy row if(nrow(basic.set)>1){ basic.set <- basic.set[-1,] } # return(basic.set) } # sub-function set.mutation <- function(db, var.id, mcode.in){ # sub-function to db.test.on.mutations - set mutation # codes: 0: no mutation (manually curated reference data set sample.batch=0) or not tested if mutation or not # (sample.batch: 0: reference, and increasing interger numbers for each batch import) # 1: mutation: manually set or manually curated ( 0->1 or 2->1 ) # 2: mutation: automatically detected (here and in check.test.ref [mutation "2" in report] ) res <- dbSendQuery(conn=db, paste("UPDATE variations SET mutation=",mcode.in," WHERE var_id=",var.id ,sep="") ) return() } # sub-function brca 1+2 check.test.ref.brca12 <- function(db, geneToWork, basic.set, result.set, id.2, mcode, mcode.in, flagIt, ref.stat){ # sub-function to db.test.on.mutations - check each line in test block against reference block (basic.set) # get a per sample test list # warning w.1 <- NULL n1 <- 0 # counter set flag events mut.set <- result.set # suspected mutation set neg.set <- result.set # negative observation set ntr <- 0 # number of total rows read (not acontrol on how many lines are finally processed) listgnf <- 0 # gene in given list not found id.2.len <- length(id.2[,"sample_id"]) # for each patient (sample_id) for(i in 1:id.2.len){ # get sample_id i.id <- id.2[i,"sample_id"] # get name i.name <- id.2[i,"sample_description"] # select test set id.2.one <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE sample_id LIKE '",i.id,"'" ,sep="") ) nr <- nrow(id.2.one) cat("\n nr ",nr) # for each test row of patient (sample_id) for(j in 1:nr){ var.id <- id.2.one[j, "var_id"] # read id gene <- id.2.one[j, "genes"] # read gene gene <- unlist(strsplit(x=gene, split=" ", fixed=T))[1] # only use first name; assume space as separator chr <- id.2.one[j, "chr"] # read chr position <- id.2.one[j, "position"] # read position var <- id.2.one[j, "var"] # read var freq <- id.2.one[j, "frequency"] # read freq # work on gene if in list geneToWork whichgene <- grepl(pattern=gene, x=geneToWork) # gene in list ? return logical - [check if specific] , unique pattern needed n.whichgene <- sum(whichgene) # exist gene+pos+var in ref? if(n.whichgene==1){ # sub-select all entries in ref - [check if specific] sel.gene <- geneToWork[whichgene] tmp <- grepl(pattern=sel.gene, x=basic.set[,"genes"]) & # select gene - [any sub-string] ## chr would be an alternative grepl(pattern=paste("^",position,"$",sep=""), x=basic.set[,"position"]) & # select position - [exact] grepl(pattern=paste("^",var,"$",sep=""), x=basic.set[,"var"]) # select var - [exact] tmp2 <- basic.set[tmp, ,drop=F] tmp2.nr <- nrow(tmp2) }else if(n.whichgene>1){ stop(paste("gene name ambiguous : ",gene,sep="")) }else if(n.whichgene==0){ listgnf <- listgnf +1 next } # not existing if(tmp2.nr==0){ # get flanking sequence fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # report mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="new entry", id.2.one[j, ])) next } # del/ins handling tmp3 <- grepl(pattern="(del|ins)", x=tmp2[,"var"], ignore.case=T) # select del-ins - [fuzzy] tmp4 <- tmp2[tmp3, ,drop=F] tmp4.nr <- nrow(tmp4) # del/ins exist -- check threshold -- based on the stat table if(tmp4.nr>0){ # select line in stat table tmp.logi <- ref.stat[,"chr"]==chr & ref.stat[,"position"]==position & ref.stat[,"var"]==var tmp.stat <- ref.stat[tmp.logi, , drop=F] # check number of rows (expected: one) tmp.stat.nr <- nrow(tmp.stat) tmp.i <- 0 # counter on solutions used # 1 : case found in stat if(tmp.stat.nr==1){ tmp.mut <- tmp.stat[,"mutation"] fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # cat("\n chr ",chr," pos ",position," var ",var,"\n") # cat("\n tmp.stat.nr ",tmp.stat.nr," w.1 ",w.1,"\n") if(freq>tmp.stat[,"freq.max"] | tmp.mut==mcode){ #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("freq > or mut, max=",tmp.stat[,"freq.max"]," n=",tmp.stat[,"n"]," mut=",tmp.mut,sep=""), id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=mcode.in); n1 <- n1 +1 } }else{ #neg neg.set <- rbind(neg.set, data.frame(description=i.name, flanking=fla, flag=paste("freq <= max=",tmp.stat[,"freq.max"]," n=",tmp.stat[,"n"]," mut=",tmp.mut,sep=""), id.2.one[j, ])) } tmp.i <- tmp.i +1 } # 0 : case -not- found in stat if(tmp.stat.nr==0){ w.1 <- paste(w.1," j : ",j,", l ",l,", tmp.stat.nr : ",tmp.stat.nr,"\n",sep="") fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("-not- in stat table - no flag",sep=""), id.2.one[j, ])) tmp.i <- tmp.i +2 } # >1 : multiple cases in stat if(tmp.stat.nr>1){ w.1 <- paste(w.1," j : ",j,", l ",l,", tmp.stat.nr : ",tmp.stat.nr,"\n",sep="") fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("-multiple- entries in stat table - no flag",sep=""), id.2.one[j, ])) tmp.i <- tmp.i +4 } # error handling if( !(tmp.i %in% c(0,1,2,4)) ){ stop(paste("wrong del/ins handling : ",tmp.i,sep="")) } next } # SNP handling fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # 2 old entries are existing, one entry is new (total 3) if(tmp2.nr<=2){ if( any(grepl(pattern=mcode, x=tmp2[,"mutation"])) ){ # any old one flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="<=2 old obs, flag mut", id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=1); n1 <- n1 +1 } }else{ # old one not flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="<=2 old obs, no flag", id.2.one[j, ])) } # >2 existing entries }else{ if( any(grepl(pattern=mcode, x=tmp2[,"mutation"])) ){ # any old one flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=">2 old obs, flag mut", id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=1); n1 <- n1 +1 } }else{ # old one not flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=">2 old obs, no flag", id.2.one[j, ])) } } }#for each row ntr <- ntr + nr }#for patient/sample_id # warning if(!is.null(w.1)){ cat("\n warning : ",w.1) } # report cat("\n test set - fetched lines : ", ntr) cat("\n test set - gene not in work list - lines : ", listgnf) cat("\n mutations flagged in the db (as ",mcode.in,") : ", n1) # remove first dummy row nrs <- nrow(mut.set) if(nrs>1){ mut.set <- mut.set[-1,] } nrs <- nrow(neg.set) if(nrs>1){ neg.set <- neg.set[-1,] } # report nrs <- nrow(mut.set) cat("\n suspected mutation set - number of lines : ", nrs) nrs <- nrow(neg.set) cat("\n negative observation set - number of lines : ", nrs) # return( list(mut.set=mut.set, neg.set=neg.set) ) } # sub-function br8 indel frequency cutoff check.test.ref.br8indel <- function(db, geneToWork, basic.set, result.set, id.2, mcode, mcode.in, flagIt, ref.stat, ref.freq){ # sub-function to db.test.on.mutations - check each line in test block against reference block (basic.set) # get a per sample test list # warning w.1 <- NULL mut.set <- result.set # suspected mutation set neg.set <- result.set # negative observation set ntr <- 0 # number of total rows read listgnf <- 0 # gene not found (in geneToWork) listnoindel <- 0 # test line is no indel listnostat1 <- 0 # test line is not in stat table -first listnostat2 <- 0 # test line is not in stat table -second id.2.len <- length(id.2[,"sample_id"]) for(i in 1:id.2.len){ # get sample_id i.id <- id.2[i,"sample_id"] # get name i.name <- id.2[i,"sample_description"] # select test set id.2.one <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE sample_id LIKE '",i.id,"'" ,sep="") ) nr <- nrow(id.2.one) cat("\n nr ",nr) for(j in 1:nr){ # for each test row var.id <- id.2.one[j, "var_id"] # read id gene <- id.2.one[j, "genes"] # read gene gene <- unlist(strsplit(x=gene, split=" ", fixed=T))[1] # only use first name; assume space as separator chr <- id.2.one[j, "chr"] # read chr position <- id.2.one[j, "position"] # read position var <- id.2.one[j, "var"] # read var freq <- id.2.one[j, "frequency"] # read freq # work on gene if in list geneToWork whichgene <- grepl(pattern=gene, x=geneToWork) # gene in list ? return logical - [check if specific] , unique pattern needed n.whichgene <- sum(whichgene) # is the gene in the search list? if(n.whichgene==1){ # get gene name sel.gene <- geneToWork[whichgene] }else if(n.whichgene>1){ stop(paste("gene name ambiguous : ",gene,sep="")) }else if(n.whichgene==0){ listgnf <- listgnf +1 next } # is the test line an indel? tmp <- grepl(pattern="(del|ins)", x=var, ignore.case=T) if(!tmp){ neg.set <- rbind(neg.set, data.frame(description=i.name, flanking="", flag="var", id.2.one[j, ])) listnoindel <- listnoindel +1 next } # get flanking seqence fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # build search name for freq list l.name <- paste(gene,var,position,sep=".") # search in freq list l.freq <- ref.freq[[l.name]] # freq.br8.2.3.4.28072014[["ATM.insT.108117675"]] if(is.null(l.freq)){ mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("indel - no statistics/freq distribution",sep=""), id.2.one[j, ])) listnostat1 <- listnostat1 +1 next } # is in upper 10% ? l.quant <- quantile(l.freq, .9) # select line in stat table tmp.logi <- ref.stat[,"chr"]==chr & ref.stat[,"position"]==position & ref.stat[,"var"]==var tmp.stat <- ref.stat[tmp.logi, , drop=F] # check number of rows (expected: one) tmp.stat.nr <- nrow(tmp.stat) # one case found in stat if(tmp.stat.nr==1){ if(freq>=l.quant){ #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("indel - freq >= ",l.quant,"(.9), freq=",freq," n=",tmp.stat[,"n"],sep=""), id.2.one[j, ])) }else{ #neg neg.set <- rbind(neg.set, data.frame(description=i.name, flanking=fla, flag=paste("indel - freq < ",l.quant,"(.9), freq=",freq," n=",tmp.stat[,"n"],sep=""), id.2.one[j, ])) } } # multiple cases in stat if(tmp.stat.nr>1){ w.1 <- paste(w.1," j : ",j,", l ",l,", tmp.stat.nr : ",tmp.stat.nr,"\n",sep="") #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("-multiple- entries in stat table",sep=""), id.2.one[j, ])) } # case -not- found in stat - should be already filtered by !is.null(l.freq) if(tmp.stat.nr==0){ listnostat2 <- listnostat2 +1 } } # for j test row ntr <- ntr + nr # block of lines per sample id } # for i sample no. # warning if(!is.null(w.1)){ cat("\n warning : ",w.1) } # report cat("\n test set - fetched lines : ", ntr) cat("\n test set - gene not in work list - lines : ", listgnf) cat("\n test set - no indel - lines : ", listnoindel) cat("\n test set - no stat (1) - lines : ", listnostat1) cat("\n test set - no stat (2) - lines : ", listnostat2) # cat("\n mutations flagged in the db (as ",mcode.in,") : ", n1) # remove first dummy row nrs <- nrow(mut.set) if(nrs>1){ mut.set <- mut.set[-1,] } nrs <- nrow(neg.set) if(nrs>1){ neg.set <- neg.set[-1,] } # report nrs <- nrow(mut.set) cat("\n suspected mutation set - number of lines : ", nrs) nrs <- nrow(neg.set) cat("\n negative observation set - number of lines : ", nrs) # return( list(mut.set=mut.set, neg.set=neg.set) ) } # sub-function br8 advanced 1 check.test.ref.br8adv1 <- function(db, geneToWork, basic.set, result.set, id.2, mcode, mcode.in, flagIt, ref.stat){ # sub-function to db.test.on.mutations - check each line in test block against reference block (basic.set) # get a per sample test list # warning w.1 <- NULL n1 <- 0 # counter set flag events mut.set <- result.set # suspected mutation set neg.set <- result.set # negative observation set ntr <- 0 # number of total rows read (not acontrol on how many lines are finally processed) listgnf <- 0 # gene in given list not found id.2.len <- length(id.2[,"sample_id"]) # for each patient (sample_id) for(i in 1:id.2.len){ # get sample_id i.id <- id.2[i,"sample_id"] # get name i.name <- id.2[i,"sample_description"] # select test set id.2.one <- dbGetQuery(conn=db, paste("SELECT * FROM variations WHERE sample_id LIKE '",i.id,"'" ,sep="") ) nr <- nrow(id.2.one) cat("\n nr ",nr) # for each test row of patient (sample_id) for(j in 1:nr){ var.id <- id.2.one[j, "var_id"] # read id gene <- id.2.one[j, "genes"] # read gene gene <- unlist(strsplit(x=gene, split=" ", fixed=T))[1] # only use first name; assume space as separator chr <- id.2.one[j, "chr"] # read chr position <- id.2.one[j, "position"] # read position var <- id.2.one[j, "var"] # read var freq <- id.2.one[j, "frequency"] # read freq # work on gene if in list geneToWork whichgene <- grepl(pattern=gene, x=geneToWork) # gene in list ? return logical - [check if specific] , unique pattern needed n.whichgene <- sum(whichgene) # exist gene+pos+var in ref? if(n.whichgene==1){ # sub-select all entries in ref - [check if specific] sel.gene <- geneToWork[whichgene] tmp <- grepl(pattern=sel.gene, x=basic.set[,"genes"]) & # select gene - [any sub-string] grepl(pattern=paste("^",chr,"$",sep=""), x=basic.set[,"chr"]) & # select chr - [exact] grepl(pattern=paste("^",position,"$",sep=""), x=basic.set[,"position"]) & # select position - [exact] grepl(pattern=paste("^",var,"$",sep=""), x=basic.set[,"var"]) # select var - [exact] tmp2 <- basic.set[tmp, ,drop=F] tmp2.nr <- nrow(tmp2) }else if(n.whichgene>1){ stop(paste("gene name ambiguous : ",gene,sep="")) }else if(n.whichgene==0){ listgnf <- listgnf +1 next } # not existing if(tmp2.nr==0){ # get flanking sequence fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # report mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="new entry", id.2.one[j, ])) next } # del/ins handling tmp3 <- grepl(pattern="(del|ins)", x=tmp2[,"var"], ignore.case=T) # select del-ins - [fuzzy] tmp4 <- tmp2[tmp3, ,drop=F] tmp4.nr <- nrow(tmp4) # del/ins exist -- check threshold -- based on the stat table if(tmp4.nr>0){ # select line in stat table tmp.logi <- ref.stat[,"chr"]==chr & ref.stat[,"position"]==position & ref.stat[,"var"]==var tmp.stat <- ref.stat[tmp.logi, , drop=F] # check number of rows (expected: one) tmp.stat.nr <- nrow(tmp.stat) tmp.i <- 0 # counter on solutions used # 1 : case found in stat if(tmp.stat.nr==1){ tmp.mut <- tmp.stat[,"mutation"] fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # cat("\n chr ",chr," pos ",position," var ",var,"\n") # cat("\n tmp.stat.nr ",tmp.stat.nr," w.1 ",w.1,"\n") if(freq>tmp.stat[,"freq.max"] | tmp.mut==mcode){ #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("freq > or mut, max=",tmp.stat[,"freq.max"]," n=",tmp.stat[,"n"]," mut=",tmp.mut,sep=""), id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=mcode.in); n1 <- n1 +1 } }else{ #neg neg.set <- rbind(neg.set, data.frame(description=i.name, flanking=fla, flag=paste("freq <= max=",tmp.stat[,"freq.max"]," n=",tmp.stat[,"n"]," mut=",tmp.mut,sep=""), id.2.one[j, ])) } tmp.i <- tmp.i +1 } # 0 : case -not- found in stat if(tmp.stat.nr==0){ w.1 <- paste(w.1," j : ",j,", l ",l,", tmp.stat.nr : ",tmp.stat.nr,"\n",sep="") fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("-not- in stat table - no flag",sep=""), id.2.one[j, ])) tmp.i <- tmp.i +2 } # >1 : multiple cases in stat if(tmp.stat.nr>1){ w.1 <- paste(w.1," j : ",j,", l ",l,", tmp.stat.nr : ",tmp.stat.nr,"\n",sep="") fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) #pos mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=paste("-multiple- entries in stat table - no flag",sep=""), id.2.one[j, ])) tmp.i <- tmp.i +4 } # error handling if( !(tmp.i %in% c(0,1,2,4)) ){ stop(paste("wrong del/ins handling : ",tmp.i,sep="")) } next } # SNP handling fla <- get.flank(position=position, chr=paste("chr",chr,sep=""), alleles="[.]", offset=10) # 2 old entries are existing, one entry is new (total 3) if(tmp2.nr<=2){ if( any(grepl(pattern=mcode, x=tmp2[,"mutation"])) ){ # any old one flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="<=2 old obs, flag mut", id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=1); n1 <- n1 +1 } }else{ # old one not flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag="<=2 old obs, no flag", id.2.one[j, ])) } # >2 existing entries }else{ if( any(grepl(pattern=mcode, x=tmp2[,"mutation"])) ){ # any old one flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=">2 old obs, flag mut", id.2.one[j, ])) # set mut. flag in database if(flagIt){ set.mutation(db=db, var.id=var.id, mcode.in=1); n1 <- n1 +1 } }else{ # old one not flagged as mutation mut.set <- rbind(mut.set, data.frame(description=i.name, flanking=fla, flag=">2 old obs, no flag", id.2.one[j, ])) } } }#for each row ntr <- ntr + nr }#for patient/sample_id # warning if(!is.null(w.1)){ cat("\n warning : ",w.1) } # report cat("\n test set - fetched lines : ", ntr) cat("\n test set - gene not in work list - lines : ", listgnf) cat("\n mutations flagged in the db (as ",mcode.in,") : ", n1) # remove first dummy row nrs <- nrow(mut.set) if(nrs>1){ mut.set <- mut.set[-1,] } nrs <- nrow(neg.set) if(nrs>1){ neg.set <- neg.set[-1,] } # report nrs <- nrow(mut.set) cat("\n suspected mutation set - number of lines : ", nrs) nrs <- nrow(neg.set) cat("\n negative observation set - number of lines : ", nrs) # return( list(mut.set=mut.set, neg.set=neg.set) ) } # main - detect mutations/polymorphisms function db.test.on.mutations <- function(DB.path.name="", ref1="", test2="", geneToWork=NULL, flagIt=F, ref.stat=NULL, ref.freq=NULL, mcode=1, mcode.in=2, out.path.name="", method="brca12"){ # get all occurances of a known or new mutation in new imported cases # and report those in a table; also mark new mutations with 1 in variants.mutation column # new cases: group in sample table batch has e.g. 1 # so reference group has ref1= e.g. "0" or c("0","1") and test group has test2 = e.g. 1 or "" etc. # geneToWork: HGNC gene name to work on (can be a vector of strings) # flagIt: T: flag new mutations in database, F: only produce output lists # ref.stat: data.frame holding statistics of existing database entries # ref.freq: corresponding frequency distributions # mcode: mutation code to be used - see coding_scheme.txt # mcode.in: mutation code to be set # output will be grouped by sample # DB.path.name: sqlite database with the given name and relative path # out.path.name: relative output path and file name # method: different algorithms, see available subfunctions # ini require(RSQLite) DB.path.name <- paste(getwd(),DB.path.name,sep="") if(DB.path.name==""){ cat("\n Please give a name - stop\n"); return() } if(is.null(ref.stat)){ cat("\n Please give the statistics data.frame - stop\n"); return() } if(is.null(ref.freq) & method=="br8simple"){ cat("\n Please give the statistics frequency list - stop\n"); return() } if(is.null(geneToWork)){ cat("\n Please give at least one gene name - stop\n"); return() } cat("\n check (!) : ref1 = ", ref1, " test2 = ", test2, "\n") # If the named database does exist, open it if(!file.exists(DB.path.name)){ # file.access(names, mode={0124}) cat("\n Database does not exist / or wrong name/path - stop \n") return() } ## open database db <- dbConnect(dbDriver("SQLite"), dbname=DB.path.name) ## make groups, extract ids, create results # establish groups len.ref1 <- length(ref1) len.test2 <- length(test2) id.1 <- NULL id.2 <- NULL for(i in 1:len.ref1){ tmp <- dbGetQuery(conn=db, paste("SELECT * FROM sample WHERE batch LIKE '",ref1[i],"'",sep="")) if(nrow(tmp)>0 & !is.null(id.1)){ # order of if clauses matters id.1 <- rbind(id.1,tmp) } if(nrow(tmp)>0 & is.null(id.1)){ id.1 <- tmp } } for(i in 1:len.test2){ tmp <- dbGetQuery(conn=db, paste("SELECT * FROM sample WHERE batch LIKE '",test2[i],"'",sep="")) if(nrow(tmp)>0 & !is.null(id.2)){ # order of if clauses matters id.2 <- rbind(id.2,tmp) } if(nrow(tmp)>0 & is.null(id.2)){ id.2 <- tmp } } ## extract basic set # determine length of record tmp <- dbGetQuery(conn=db, "SELECT * FROM variations LIMIT 2") # tmp.n <- names(tmp) tmp.nc <- ncol(tmp) # create data.frame basic.set <- data.frame(matrix(NA,1,tmp.nc)) names(basic.set) <- tmp.n # use this template also for the results snd submit the template to the subfunction result.set <- basic.set result.set <- cbind("description","flanking","flag",result.set) # extend by additional columns names(result.set) <- c("description","flanking","flag",tmp.n) # fill basic.set basic.set <- fill.ref.set(db=db, basic.set=basic.set, id.1=id.1) ### deleted : variation around position - make function, if needed (see fragment) ## create list with mutation samples if(method=="brca12"){ result.set <- check.test.ref.brca12(db=db, geneToWork=geneToWork, basic.set=basic.set, result.set=result.set, id.2=id.2, mcode=mcode, mcode.in=mcode.in, flagIt=flagIt, ref.stat=ref.stat) } if(method=="br8simple"){ result.set <- check.test.ref.br8simple(db=db, geneToWork=geneToWork, basic.set=basic.set, result.set=result.set, id.2=id.2, mcode=mcode, flagIt=flagIt, ref.stat=ref.stat, ref.freq=ref.freq) } if(method=="br8indel"){ result.set <- check.test.ref.br8indel(db=db, geneToWork=geneToWork, basic.set=basic.set, result.set=result.set, id.2=id.2, mcode=mcode, flagIt=flagIt, ref.stat=ref.stat, ref.freq=ref.freq) } if(method=="br8adv1"){ result.set <- check.test.ref.br8adv1(db=db, geneToWork=geneToWork, basic.set=basic.set, result.set=result.set, id.2=id.2, mcode=mcode, mcode.in=mcode.in, flagIt=flagIt, ref.stat=ref.stat) } # report cat("\n reference set - number of lines : ", nrow(basic.set)) # output to TAB text file pfname1 <- paste(getwd(), out.path.name,".mut.",paste("genes",length(geneToWork),sep=""),".",format(Sys.time(), "%Y%m%d%H%M"),".txt",sep="") cat("\n output suspected mutations: ",pfname1) write.table(x=result.set[["mut.set"]], file=pfname1, append=F, quote=T, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) pfname2 <- paste(getwd(), out.path.name,".neg.",paste("genes",length(geneToWork),sep=""),".",format(Sys.time(), "%Y%m%d%H%M"),".txt",sep="") cat("\n output negative observations: ",pfname2,"\n") write.table(x=result.set[["neg.set"]], file=pfname2, append=F, quote=T, sep="\t", eol="\n", na="NA", dec=".", row.names=F, col.names=T) # close dbDisconnect(db) cat("\n end \n") # return(result.set) return() } # analyse test distribution versus basic distribution db.test.dist.basic.dist.frequency <- function(DB.path.name="", batch=NULL, ref.stat=NULL, ref.freq=NULL, ref.bala=NULL, out.path.name=""){ # generate statistics on frequency e.g. of a new batch and does compare to an existing statistics table (reference) # based on data from db.stat.frequency() # DB.path.name: sqlite database with the given name and relative path # out.path.name: relative path and filename (prefix) for PDF output # batch : number or vector of batch numbers to be analyzed (data blocks of the database) # ref.stat : reference table to look into, ref.freq: list of frequencies according to stat table, ref.bala: list of balances according to stat table # ini require(RSQLite) source("/home/korschi/eclipseR/0functions/0general/hist.plot.simple.R") source("/home/korschi/eclipseR/0functions/0general/trim.distribution.R") if(DB.path.name==""){ cat("\n Please give a database name and rel. path - stop\n"); return() } if(out.path.name==""){ cat("\n Please give a name and rel. path for PDF output file - stop\n"); return() } if(is.null(batch)){ cat("\n Please give a batch number (or a vector of numbers) to be analyzed - stop\n"); return() } if(is.null(ref.stat)){ cat("\n Please give a statistics data.frame name (R variable) - stop\n"); return() } if(is.null(ref.freq)){ cat("\n Please give a frequency list name (R variable) - stop\n"); return() } if(is.null(ref.bala)){ cat("\n Please give a balance list name (R variable) - stop\n"); return() } # set trim parameter (coherent to db.stat.frequency() !) trim.d <- "upper" trim.v <- 0.05 # If the named database does exist, open it if(file.exists(paste(getwd(),DB.path.name,sep=""))){ db <- dbConnect(dbDriver("SQLite"), dbname=paste(getwd(),DB.path.name,sep="")) }else{ cat("\n Database does not exist / or wrong name/path - stop \n") return() } # get number of rows in variations t.rows <- dbGetQuery(conn=db, "SELECT Count(*) FROM variations" ) ## make groups, extract ids, create results # we take the new batch group(s) tmp.out <- paste("\n we take batch ",batch,sep="") sql.st <- paste("SELECT variations.* FROM sample INNER JOIN variations ON sample.sample_id = variations.sample_id WHERE sample.batch IN (",paste(batch,collapse=","),")" , sep="") id.0 <- dbGetQuery(conn=db, sql.st ) id.0.nr <- nrow(id.0) # close dbDisconnect(db) # col names : # var_id sample_id mutation chr position var quality frequency coverage reads_var balance genes type known prediction maf pcranno # stat.. table # c("chr","position","var","mutation","genes","n","freq.mean","freq.median","freq.sd","freq.t.mean","freq.t.sd","balance.mean","freq.max","freq.min","bal.max","bal.min") # freq.. list # format of element name(s): "gene name . variant . position" e.g. "RAD51C.T>A.56798139" content: vector of frequency values # calculations # output to PDF pdf(file=paste(getwd(),out.path.name,".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title="test&location", pointsize=12) par(mfrow=c(2,2)) # work on each batch row k1 <- 0 # counter k2 <- 0 # counter for(i in 1:id.0.nr){ # for each test row tmp.rs <- NULL tmp.rf <- NULL tmp.rb <- NULL var.id <- id.0[i, "var_id"] # read var_id gene <- id.0[i, "genes"] # read gene(s) chr <- id.0[i, "chr"] # read chr position <- id.0[i, "position"] # read position var <- id.0[i, "var"] # read var frequency <- id.0[i, "frequency"] # read frequency balance <- id.0[i, "balance"] # read balance cat("\n counter ",i," gene ",gene," var ",var," position ",position) # find in ref.stat and ref.freq # extract first gene -- check if intended ! first.gene <- sub(pattern="^\"", replacement="", x=gene) # rep.: \\1 first blank ^\"(\\w*)\\s.* extract up to the last blank ^\"(.*)\\s.* first.gene <- sub(pattern="\\s.*", replacement="", x=first.gene) first.gene <- sub(pattern="\"$", replacement="", x=first.gene) #cat("\n first.gene ",first.gene) first.ref.genes <- sub(pattern="^\"", replacement="", x=ref.stat[,"genes"]) # filter on first genes first.ref.genes <- sub(pattern="\\s.*", replacement="", x=first.ref.genes) first.ref.genes <- sub(pattern="\"$", replacement="", x=first.ref.genes) # select stat tmp.rs <- ref.stat[first.ref.genes==first.gene & ref.stat[,"var"]==var & ref.stat[,"position"]==position , , drop=F] tmp.rs.nr <- nrow(tmp.rs) #cat("\n nrow(tmp.rs) ",tmp.rs.nr," tmp.rs ",paste(tmp.rs,collapse=";")) if(tmp.rs.nr>1){ cat("\n multiple rows in ref.stat -- check -- stop ") return() } # select freq var <- gsub(pattern="\"", replacement="", x=var) #unlist(strsplit(var, split='"', fixed=T))[2] : not optimal tmp.n <- paste(first.gene,var,position,sep=".") #cat("\n tmp.n ",tmp.n) tmp.rf <- ref.freq[[tmp.n]] #cat("\n tmp.rf ",tmp.rf) # if(){} # check if multiple entries ? # select balance tmp.rb <- ref.bala[[tmp.n]] #cat("\n tmp.rb ",tmp.rb) #if(i==116){ return(first.ref.genes) } if(tmp.rs.nr>0 & !is.null(tmp.rf)){ # plot distribution etc. k1 <- k1 +1 # draw 2 text blocks & 2 graphs plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") mtext(text=paste("Korsching",format(Sys.time(), "%a %d %b %Y"),sep=", "), side=3, line=1, adj=0, cex=0.5) text(x=0.5, y=8, labels=paste("row #",i,sep=" "), adj=c(0,0.5), col="red") text(x=0.5, y=7, labels="", adj=c(0,0.5)) text(x=0.5, y=6, labels=paste("frequency sd",round(tmp.rs[1, "freq.sd"],3),sep=": "), adj=c(0,0.5), col="blue") text(x=0.5, y=5, labels=paste("frequency median",round(tmp.rs[1, "freq.median"],3),sep=": "), adj=c(0,0.5), col="blue") text(x=0.5, y=4, labels=paste("frequency (trim ",trim.d," ",trim.v,") mean: ",round(tmp.rs[1, "freq.t.mean"],3)," sd: ",round(tmp.rs[1, "freq.t.sd"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=3, labels="", adj=c(0,0.5)) text(x=0.5, y=2, labels=paste("frequency min/max : ",round(tmp.rs[1, "freq.min"],3)," - ",round(tmp.rs[1, "freq.max"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=1, labels=paste("balance min/max : ",round(tmp.rs[1, "bal.min"],3)," - ",round(tmp.rs[1, "bal.max"],3),sep=""), adj=c(0,0.5), col="blue") mtext(text="histogram", side=1, line=1) plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") text(x=0.5, y=9, labels="", adj=c(0,0.5)) text(x=0.5, y=8, labels="", adj=c(0,0.5)) text(x=0.5, y=7, labels=paste("REFERENCE: ",tmp.rs[1, "var"],", chr",tmp.rs[1, "chr"],"-",tmp.rs[1, "position"], sep=" "), adj=c(0,0.5), col="blue") text(x=0.5, y=6, labels=paste(" gene: ",tmp.rs[1, "genes"],sep=" "), adj=c(0,0.5), col="blue") text(x=0.5, y=5, labels=paste(" mean: frequency: ",round(tmp.rs[1, "freq.mean"],3)," balance: ",round(tmp.rs[1, "balance.mean"],3),sep=""), adj=c(0,0.5), col="blue") text(x=0.5, y=4, labels="", adj=c(0,0.5)) text(x=0.5, y=3, labels=paste("TEST: ",var," , chr ",chr," - ",position,sep=""), adj=c(0,0.5), col="red") text(x=0.5, y=2, labels=paste(" gene: ",gene,sep=""), adj=c(0,0.5), col="red") text(x=0.5, y=1, labels=paste(" frequency: ",frequency," balance: ",balance,sep=""), adj=c(0,0.5), col="red") mtext(text="histogram", side=1, line=1) # mtext(text="density", side=1, line=1) hist.plot.simple(x=tmp.rf, xcolor=T, bin.num=F, bin.size=F, bin.fix=NULL, bin.fix.open=F, norm=F, n.factor=1, curve=F, xlab="frequency", ylab="counts/bin", xlim=c(0,1), y.max=NULL, x.at=c(0,0.25,0.5,0.75,1), y.at=NULL, y.num.format=F, ylog=F, bar.width=1, offset=0, digits=1, col=c("blue","red","green"), cex=1, smooth.f=0, subfolder.filename=NULL, h.title="", s.title=T, verbose=F, add=F) # add a line for the presently tested new variant abline(v=frequency, col="green", lwd=3) # if(length(tmp.rf)>=20){ # plot( density(x=tmp.rf, bw="SJ", kernel="gaussian"), main="" ) # }else{ # plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") # empty plot # text(x=0.5, y=5, labels="no density plot : count < 20", adj=c(0,0.5), cex=1.2) # } hist.plot.simple(x=tmp.rb, y=tmp.rf, bin.num=F, bin.size=F, bin.fix=NULL, bin.fix.open=F, norm=F, n.factor=1, curve=F, xlab="balance", ylab="counts/bin", xlim=c(0,1), y.max=NULL, x.at=c(0,0.25,0.5,0.75,1), y.at=NULL, y.num.format=F, ylog=F, bar.width=1, offset=0, digits=1, col=c("blue","red","green"), cex=1, smooth.f=0, subfolder.filename=NULL, h.title="", s.title=T, verbose=F, add=F) # add a line for the presently tested new variant abline(v=balance, col="green", lwd=3) }else{ # no plot if no statistics exist (and frequency (or balance) list should have no entry too) if(tmp.rs.nr==0 & is.null(tmp.rf)){ k2 <- k2 +1 # but count those events # could be due to that # a new event is observed # not all genes of the basic set are analysed # no corresponding statistics table is used # register event and export in file #### open } } } par(mfrow=c(1,1)) plot(x=c(0:10), y=c(0:10), type="n", axes=F, xlab="", ylab="") text(x=0.5, y=9, labels=paste( "\n batch: ",paste(batch,collapse=","), "\n ref.stat: ",deparse(substitute(ref.stat)), "\n ref.freq: ",deparse(substitute(ref.freq)), "\n ref.bala: ",deparse(substitute(ref.bala)), "\n number of lines in batch : ", id.0.nr, "\n ref statistics & frequency(balance) list : ",k1," (number of observations)", "\n no ref statistics : ",k2, "\n due to -a new event is observed, -not all genes are analysed, -wrong stat table used", sep=""), adj=c(0,0.5)) dev.off() # report some run details cat( "\n\n batch: ",batch," ref.stat: ",deparse(substitute(ref.stat)), " ref.freq: ",deparse(substitute(ref.freq)), " ref.bala: ",deparse(substitute(ref.bala)), "\n number of lines in batch : ", id.0.nr, "\n ref statistics & frequency(balance) list : ",k1," (number of observations)", "\n no ref statistics : ",k2, "\n due to -a new event is observed, -not all genes are analysed, -wrong stat table used\n" ) # close return() } # get SNP annotations get.snp.anno.ebi <- function(chr=NULL,pos.start=NULL,pos.end=NULL){ # get SNP annotations from EBI # give chr+postion get 'rs'-number and other information if(is.null(chr)){ cat("\n Please give a chromosome number - stop\n"); return() } if(is.null(pos.start)){ cat("\n Please give a genomic (start) position - stop\n"); return() } if(is.null(pos.end)){ pos.end <- pos.start+1 } require("biomaRt") # marts <- listMarts() # snpmart <- useMart("snp") # snpsets <- listDatasets(snpmart) snpmart.hs <- useMart("snp", dataset="hsapiens_snp") # filters <- listFilters(snpmart.hs) # attributes <- listAttributes(snpmart.hs) res <- getBM(attributes=c('refsnp_id','allele','chrom_start','chrom_strand','minor_allele','minor_allele_freq','clinical_significance','synonym_name'), filters=c('chr_name','chrom_start','chrom_end'), values=list(chr,pos.start,pos.end), mart=snpmart.hs) cat("\n end ") return(res) } #aa <- date() #get.snp.anno.ebi(chr=c(11,8,22),pos.start=c(108150208,90967512,29091788)) #aa;date() # 8 s to 2h 12 and failed #var: delT genes: ATM #save.image() #var: delT genes: NBN #var: T>C genes: CHEK2 get.rs.loc.tab <- function(DB.path.name="", db.table="", chr=NULL, pos=NULL){ # get rs number from local sqlite database # give chr+postion get 'rs'-number (and maybe further information) # DB.path.name : relativ path & name # X_CHROM POS ID REF ALT QUAL FILTER INFO # 2 min if(DB.path.name==""){ cat("\n Please give a database path and name - stop\n"); return() } if(is.null(chr)){ cat("\n Please give a chromosome number - stop\n"); return() } if(is.null(pos)){ cat("\n Please give a genomic position - stop\n"); return() } #ini require(RSQLite) require(sqldf) DB.path.name <- paste(getwd(),DB.path.name,sep="") # open database if(!file.exists(DB.path.name)){ cat("\n Database file is not existing - stop\n") return() }else{ db <- dbConnect(dbDriver("SQLite"), dbname=DB.path.name) } sql.c <- paste("SELECT ID,REF,ALT FROM ",db.table, " WHERE X_CHROM like '",chr,"' AND POS=",pos,sep="") cat("\n sql ",sql.c) # get rs number e.rs <- dbGetQuery(conn=db, statement=sql.c) cat("\n rows ",nrow(e.rs)) dbDisconnect(db) # Close connection return(e.rs) } #date() #get.rs.loc.tab(DB.path.name="/All", db.table="all_rs", chr=8, pos=90945548) #date() #stat.br8.2.3.4.28072014[,c(1,2),drop=F]