Mercurial > repos > jeremyjliu > region_motif_enrichment
view region_motif_compare.r @ 15:81bd3a218c51 draft
Deleted test motif database
author | jeremyjliu |
---|---|
date | Wed, 12 Nov 2014 15:10:51 -0500 |
parents | 19d2cffb8db3 |
children |
line wrap: on
line source
# Name: region_motif_compare.r # Description: Reads in two count files and determines enriched and depleted # motifs (or any location based feature) based on poisson tests and gc # corrections. All enrichment ratios relative to overall count / gc ratios. # Author: Jeremy liu # Email: jeremy.liu@yale.edu # Date: 14/07/03 # Note: This script is meant to be invoked with the following command # R --slave --vanilla -f ./region_motif_compare.r --args <workingdir> <db> <intab1> <intab2> # <enriched_tab> <depleted_tab> <plots_png> # <workingdir> is working directory of galaxy installation # <db> types: "t" test, "p" pouya, "j" jaspar jolma, "m" mouse, "c" combined # Dependencies: plotting.r # Auxiliary function to concatenate multiple strings concat <- function(...) { input_list <- list(...) return(paste(input_list, sep="", collapse="")) } # Supress all warning messages to prevent Galaxy treating warnings as errors options(warn=-1) # Set common and data directories args <- commandArgs() workingDir = args[7] commonDir = concat(workingDir, "/tools/my_tools") dbCode = args[8] # dbCode "c" implemented when pwmFile is loaded if (dbCode == "t" | dbCode == "p") { pwmFile = concat(commonDir, "/region_motif_db/pwms/pouya.pwms.from.seq.RData") } else if (dbCode == "j") { pwmFile = concat(commonDir, "/region_motif_db/pwms/jaspar.jolma.pwms.from.seq.RData") } else if (dbCode == "m") { pwmFile = concat(commonDir, "/region_motif_db/pwms/mm9.pwms.from.seq.RData") } else if (dbCode == "c") { # rest of dbCode "c" implemeted when pwmFile loaded pwmFile = concat(commonDir, "/region_motif_db/pwms/pouya.pwms.from.seq.RData") pwmFile2 = concat(commonDir, "/region_motif_db/pwms/jaspar.jolma.pwms.from.seq.RData") } else { pwmFile = concat(commonDir, "/region_motif_db/pwms/pouya.pwms.from.seq.RData") } # Set input and reference files inTab1 = args[9] inTab2 = args[10] enrichTab = args[11] depleteTab = args[12] plotsPng = args[13] # Load dependencies source(concat(commonDir, "/region_motif_lib/plotting.r")) # Auxiliary function to read in tab file and prepare the data read_tsv <- function(file) { data = read.table(file, sep="\t", stringsAsFactors=FALSE) names(data)[names(data) == "V1"] = "motif" names(data)[names(data) == "V2"] = "counts" return(data) } startTime = Sys.time() cat("Running ... Started at:", format(startTime, "%a %b %d %X %Y"), "...\n") # Loading motif position weight matrix (pwm) file and input tab file #cat("Loading and reading input region motif count files...\n") load(pwmFile) # pwms data structure if (dbCode == "c") { # Remaining implementation of dbCode "c" combined temp = pwms load(pwmFile2) pwms = append(temp, pwms) } region1DF = read_tsv(inTab1) region2DF = read_tsv(inTab2) region1Counts = region1DF$counts region2Counts = region2DF$counts names(region1Counts) = region1DF$motif names(region2Counts) = region2DF$motif # Processing count vectors to account for missing 0 count motifs, then sorting #cat("Performing 0 count correction and sorting...\n") allNames = union(names(region1Counts), names(region2Counts)) region1Diff = setdiff(allNames, names(region1Counts)) region2Diff = setdiff(allNames, names(region2Counts)) addCounts1 = rep(0, length(region1Diff)) addCounts2 = rep(0, length(region2Diff)) names(addCounts1) = region1Diff names(addCounts2) = region2Diff newCounts1 = append(region1Counts, addCounts1) newCounts2 = append(region2Counts, addCounts2) region1Counts = newCounts1[sort.int(names(newCounts1), index.return=TRUE)$ix] region2Counts = newCounts2[sort.int(names(newCounts2), index.return=TRUE)$ix] # Generate gc content matrix gc = sapply(pwms, function(i) mean(i[2:3,3:18])) # Apply poisson test, calculate p and q values, and filter significant results #cat("Applying poisson test...\n") rValue = sum(region2Counts) / sum(region1Counts) pValue = sapply(seq(along=region1Counts), function(i) { poisson.test(c(region1Counts[i], region2Counts[i]), r=1/rValue)$p.value }) qValue = p.adjust(pValue, "fdr") indices = which(qValue<0.1 & abs(log2(region1Counts/region2Counts/rValue))>log2(1.5)) # Setting up output diagnostic plots, 4 in 1 png image png(plotsPng, width=800, height=800) xlab = "region1_count" ylab = "region2_count" lim = c(0.5, 5000) layout(matrix(1:4, ncol=2)) par(mar=c(5, 5, 5, 1)) # Plot all motif counts along the linear correlation coefficient plot.scatter(region1Counts+0.5, region2Counts+0.5, log="xy", xlab=xlab, ylab=ylab, cex.lab=2.2, cex.axis=1.8, xlim=lim, ylim=lim*rValue) abline(0, rValue, untf=T) abline(0, rValue*2, untf=T, lty=2) abline(0, rValue/2, untf=T, lty=2) # Plot enriched and depleted motifs in red, housed in second plot plot.scatter(region1Counts+0.5, region2Counts+0.5, log="xy", xlab=xlab, ylab=ylab, cex.lab=2.2, cex.axis=1.8, xlim=lim, ylim=lim*rValue) points(region1Counts[indices]+0.5, region2Counts[indices]+0.5, col="red") abline(0, rValue, untf=T) abline(0, rValue*2, untf=T, lty=2) abline(0, rValue/2, untf=T, lty=2) # Apply and plot gc correction and loess curve #cat("Applying gc correction, rerunning poisson test...\n") ind = which(region1Counts>5) gc = gc[names(region2Counts)] # Reorder the indices of pwms to match input data lo = plot.scatter(gc,log2(region2Counts/region1Counts),draw.loess=T, xlab="gc content of motif",ylab=paste("log2(",ylab,"/",xlab,")"), cex.lab=2.2,cex.axis=1.8,ind=ind) # This function is in plotting.r gcCorrection = 2^approx(lo$loess,xout=gc,rule=2)$y save(gc, file="gc.RData") # Recalculate p and q values, and filter for significant entries pValueGC = sapply(seq(along=region1Counts),function(i) { poisson.test(c(region1Counts[i],region2Counts[i]),r=1/gcCorrection[i])$p.value }) qValueGC=p.adjust(pValueGC,"fdr") indicesGC = which(qValueGC<0.1 & abs(log2(region1Counts/region2Counts*gcCorrection))>log2(1.5)) # Plot gc corrected motif counts plot.scatter(region1Counts+0.5, (region2Counts+0.5)/gcCorrection, log="xy", xlab=xlab, ylab=paste(ylab,"(normalized)"), cex.lab=2.2, cex.axis=1.8, xlim=lim, ylim=lim) points(region1Counts[indicesGC]+0.5, (region2Counts[indicesGC]+0.5)/gcCorrection[indicesGC], col="red") abline(0,1) abline(0,1*2,untf=T,lty=2) abline(0,1/2,untf=T,lty=2) # Trim results, compile statistics and output to file # Only does so if significant results are computed if(length(indicesGC) > 0) { # Calculate expected counts and enrichment ratios #cat("Calculating statistics...\n") nullExpect = region1Counts * gcCorrection enrichment = region2Counts / nullExpect # Reorder selected indices in ascending pvalue #cat("Reordering by ascending pvalue...\n") indicesReorder = indicesGC[order(pValueGC[indicesGC])] # Combine data into one data frame and output to two files #cat("Splitting and outputting data...\n") outDF = data.frame(motif=names(pValueGC), p=as.numeric(pValueGC), q=qValueGC, stringsAsFactors=F, region_1_count=region1Counts, null_expectation=round(nullExpect,2), region_2_count=region2Counts, enrichment=enrichment)[indicesReorder,] names(outDF)[which(names(outDF)=="region_1_count")]=xlab names(outDF)[which(names(outDF)=="region_2_count")]=ylab indicesEnrich = which(outDF$enrichment>1) indicesDeplete = which(outDF$enrichment<1) outDF$enrichment = ifelse(outDF$enrichment>1, round(outDF$enrichment,3), paste("1/",round(1/outDF$enrichment,3))) write.table(outDF[indicesEnrich,], file=enrichTab, quote=FALSE, sep="\t", append=FALSE, row.names=FALSE, col.names=TRUE) write.table(outDF[indicesDeplete,], file=depleteTab, quote=FALSE, sep="\t", append=FALSE, row.names=FALSE, col.names=TRUE) } # Catch display messages and output timing information catchMessage = dev.off() cat("Done. Job started at:", format(startTime, "%a %b %d %X %Y."), "Job ended at:", format(Sys.time(), "%a %b %d %X %Y."), "\n")