annotate Rscripts/ridb-regression.R @ 24:385d21a8d0a0

New tool to Query multiple public repositories for elemental compositions from accurate mass values detected by high-resolution mass spectrometers

##
#
# Performs regression analysis using either 3rd degree polynomial- or linear-method
#
##

# Commandline arguments
args  <- commandArgs(TRUE)
if (length(args) < 7)
	stop(cat("Missing arguments, usage:\n\tRscript ridb-regression.R RI-database ",
             "ouput_file logfile min_residuals range_mod pvalue rsquared method ",
			 "plot(yes/no) plot_archive"))

ridb <- args[1]
out_file <- args[2]
logfile <- args[3]
min_residuals <- as.integer(args[4])
range_mod <- as.integer(args[5])
pvalue <- as.double(args[6])
rsquared <- as.double(args[7])
method <- args[8]
plot <- tolower(args[9])
if (plot == 'true')
	plot_archive = args[10]

# Do not show warnings etc.
sink(file='/dev/null')

progress <- c()
logger <- function(logdata) {
	## Logs progress, adds a timestamp for each event
	#cat(paste(Sys.time(), "\t", logdata, "\n", sep="")) ## DEBUG
	progress <<- c(progress, paste(Sys.time(), "\t", logdata, sep=""))
}

logger("Reading Retention Index Database..")

# Read Retention Index Database
ridb <- read.csv(ridb, header=TRUE, sep="\t")
logger(paste("\t", nrow(ridb), "records read.."))
# Get a unique list 
gc_columns <- unique(as.vector(as.matrix(ridb['Column.name'])[,1]))
cas_numbers <- unique(as.vector(as.matrix(ridb['CAS'])[,1]))

add_poly_fit <- function(fit, gc1_index, gc2_index, range) {
	pval = anova.lm(fit)$Pr
	r.squared = summary(fit)$r.squared

	data = rep(NA, 11)
	# Append results to matrix
	data[1] = gc_columns[gc1_index] # Column 1
	data[2] = gc_columns[gc2_index] # Column 2
	data[3] = coefficients(fit)[1]  # The 4 coefficients
	data[4] = coefficients(fit)[2]
	data[5] = coefficients(fit)[3]
	data[6] = coefficients(fit)[4]
	data[7] = range[1]              # Left limit
	data[8] = range[2]              # Right limit
	data[9] = length(fit$residuals) # Number of datapoints analysed
	data[10] = pval[1]              # p-value for resulting fitting
	data[11] = r.squared            # R-squared
	return(data)
}


add_linear_fit <- function(fit, gc1_index, gc2_index, range) {
	pval = anova.lm(fit)$Pr
	r.squared = summary(fit)$r.squared
	
	data = rep(NA, 7)
	# Append results to matrix
	data[1] = gc_columns[gc1_index] # Column 1
	data[2] = gc_columns[gc2_index] # Column 2
	data[3] = coefficients(fit)[1]  # The 4 coefficients
	data[4] = coefficients(fit)[2]
	data[7] = length(fit$residuals) # Number of datapoints analysed
	data[8] = pval[1]               # p-value for resulting fitting
	data[9] = r.squared             # R-squared
	return(data)
}


add_fit <- function(fit, gc1_index, gc2_index, range, method) {
	if (method == 'poly')
		return(add_poly_fit(fit, gc1_index, gc2_index, range))
	else
		return(add_linear_fit(fit, gc1_index, gc2_index, range))
}


plot_fit <- function(ri1, ri2, gc1_index, gc2_index, coeff, range, method) {
    if (method == 'poly')
        pol <- function(x) coeff[4]*x^3 + coeff[3]*x^2 + coeff[2]*x + coeff[1]
    else
        pol <- function(x) coeff[2]*x + coeff[1]
    pdf(paste('regression_model_',
              make.names(gc_columns[gc1_index]), '_vs_', 
              make.names(gc_columns[gc2_index]), '.pdf', sep=''))
    curve(pol, 250:3750, col="red", lwd=2.5, main='Regression Model', xlab=gc_columns[gc1_index],
          ylab=gc_columns[gc2_index], xlim=c(250, 3750), ylim=c(250, 3750))
    points(ri1, ri2, lwd=0.4)
    # Add vertical lines showing left- and right limits when using poly method
    if (method == 'poly')
        abline(v=range, col="grey", lwd=1.5)
    dev.off()
}

# Initialize output dataframe
if (method == 'poly') {
    m <- data.frame(matrix(ncol = 11, nrow = 10))
} else {
    m <- data.frame(matrix(ncol = 9, nrow = 10))
}


get_fit <- function(gc1, gc2, method) {
	if (method == 'poly')
		return(lm(gc1 ~ poly(gc2, 3, raw=TRUE)))
	else
		return(lm(gc1 ~ gc2))
}

# Permutate
k <- 1
logger(paste("Permutating (with ", length(gc_columns), " GC-columns)..", sep=""))

for (i in 1:(length(gc_columns)-1)) {
	logger(paste("\tCalculating model for ", gc_columns[i], "..", sep=""))
	breaks <- 0
	for (j in (i+1):length(gc_columns)) {
		col1 = ridb[which(ridb['Column.name'][,1] == gc_columns[i]),]
		col2 = ridb[which(ridb['Column.name'][,1] == gc_columns[j]),]
		
		# Find CAS numbers for which both columns have data (intersect)
		cas_intersect = intersect(col1[['CAS']], col2[['CAS']])
		
		# Skip if number of shared CAS entries is < cutoff
		if (length(cas_intersect) < min_residuals) {
			breaks = breaks + 1
			next
		}
		# Gather Retention Indices
		col1_data = col1[['RI']][match(cas_intersect, col1[['CAS']])]
		col2_data = col2[['RI']][match(cas_intersect, col2[['CAS']])]
		
		# Calculate the range within which regression is possible (and move if 'range_mod' != 0)
		range = c(min(c(min(col1_data), min(col2_data))), max(c(max(col1_data), max(col2_data))))
		if (range_mod != 0) {
			# Calculate percentage and add/subtract from range
			perc = diff(range) / 100
			perc_cutoff = range_mod * perc
			range = as.integer(range + c(perc_cutoff, -perc_cutoff))
		}
		
		# Calculate model for column1 vs column2 and plot if requested
		fit = get_fit(col1_data, col2_data, method)
		m[k,] = add_fit(fit, i, j, range, method)
		
		if (plot == 'true')
			plot_fit(col1_data, col2_data, i, j, coefficients(fit), range, method)
		
		# Calculate model for column2 vs column1 and plot if requested
		fit = get_fit(col2_data, col1_data, method)
		m[k + 1,] = add_fit(fit, j, i, range, method)
		
		if (plot == 'true')
			plot_fit(col2_data, col1_data, j, i, coefficients(fit), range, method)
		
		k = k + 2
	}
	logger(paste("\t\t", breaks, " comparisons have been skipped due to nr. of datapoints < cutoff", sep=""))
}

# Filter on pvalue and R-squared
logger("Filtering on pvalue and R-squared..")
if (method == 'poly') {
    pval_index <- which(m[,10] < pvalue)
    rsquared_index <- which(m[,11] > rsquared)
} else {
    pval_index <- which(m[,8] < pvalue)
    rsquared_index <- which(m[,9] > rsquared)
}
logger(paste(nrow(m) - length(pval_index), " models discarded due to pvalue > ", pvalue, sep=""))

logger(paste(nrow(m) - length(rsquared_index), " models discarded due to R-squared < ", rsquared, sep=""))

# Remaining rows
index = unique(c(pval_index, rsquared_index))

# Reduce dataset
m = m[index,]
sink()

# Place plots in the history as a ZIP file
if (plot == 'true') {
    logger("Creating archive with model graphics..")
    system(paste("zip -9 -r models.zip *.pdf > /dev/null", sep=""))
    system(paste("cp models.zip ", plot_archive, sep=""))
}

# Save dataframe as tab separated file
logger("All done, saving data..")
header = c("Column1", "Column2", "Coefficient1", "Coefficient2", "Coefficient3", "Coefficient4", 
           "LeftLimit", "RightLimit", "Residuals", "pvalue", "Rsquared")
if (method != 'poly')
	header = header[c(1:4, 7:11)]
write(progress, logfile)
write.table(m, file=out_file, sep="\t", quote=FALSE, col.names=header, row.names=FALSE)