comparison region_motif_compare.r @ 37:76e108a19e64 draft

Uploaded

36:b868314320e2

# 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: [e-mail]
# Date: 15/02/11
# Note: This script is meant to be invoked with the following command
# R --slave --vanilla -f ./region_motif_compare.r --args <workingdir> <pwm_file> 
#        <intab1> <intab2> <enriched_tab> <depleted_tab> <plots_png>
# <workingdir> is working directory of galaxy installation
# Dependencies: region_motif_data_manager, plotting.r

# Auxiliary function to concatenate multiple strings
concat <- function(...) {
	input_list <- list(...)
	return(paste(input_list, sep="", collapse=""))

}

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

37:76e108a19e64

# 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: [e-mail]
# Date: 15/02/11
# Note: This script can be invoked with the following command
# R --slave --vanilla -f ./region_motif_compare.r --args <workingdir> <pwm_file> 
#        <intab1> <intab2> <enriched_tab> <depleted_tab> <plots_png>
# <workingdir> is the directory where plotting.r is saved
# Dependencies: region_motif_data_manager, plotting.r, 

# Auxiliary function to concatenate multiple strings
concat <- function(...) {
	input_list <- list(...)
	return(paste(input_list, sep="", collapse=""))

}

startTime = Sys.time()
cat("Running ... Started at:", format(startTime, "%a %b %d %X %Y"), "...\n")

# Loading motif position weight matrix (pwm) file 
cat("Loading motif postion weight matrices...\n")
lines = scan(pwmFile, what="character", sep="\n", quiet=TRUE)
indices = which(grepl("MOTIF", lines))
names(indices) = lapply(indices, function(i) {
    nameline = lines[i]
    name = substr(nameline, 7, nchar(nameline))
    })

pwms = sapply(indices, function(i) {
    infoline = unlist(strsplit(lines[i+1], " "))
    alength = as.numeric(infoline[4])
    width = as.numeric(infoline[6])
    subset = lines[(i+2):(i+2+width-1)]
    motiflines = strsplit(subset, " ")
    motif = t(do.call(rbind, motiflines))
    motif = apply(motif, 2, as.numeric)
    }, simplify=FALSE, USE.NAMES=TRUE)

# Loading input tab files
cat("Loading and reading input region motif count files...\n")
region1DF = read_tsv(inTab1)
region2DF = read_tsv(inTab2)
region1Counts = region1DF$counts
region2Counts = region2DF$counts
names(region1Counts) = region1DF$motif