#############################
## R PROGRAMMING EXERCISES ##
#############################
## Workshop: Programming in R
## May 14, 2011 

##########################################
## (1) Loops: for/while and apply loops ##
##########################################
## Create a sample matrix

myMA <- matrix(rnorm(500), 100, 5, dimnames=list(1:100, paste("C", 1:5, sep="")))

## (1.1) Compute the mean of each row in myMA by applying the mean function in a for loop

myve_for <- NULL
for(i in seq(along=myMA[,1])) {
	myve_for <- c(myve_for, mean(as.numeric(myMA[i, ])))
}
myResult <- cbind(myMA, mean_for=myve_for)
myResult[1:4, ]

## (1.2) Compute the mean of each row in myMA by applying the mean function in a while loop

z <- 1
myve_while <- NULL
while(z <= length(myMA[,1])) {
	myve_while <- c(myve_while, mean(as.numeric(myMA[z, ])))
	z <- z + 1
}
myResult <- cbind(myMA, mean_for=myve_for, mean_while=myve_while)
myResult[1:4, ]

## (1.3) How can one confirm that the results from both mean calculations are identical?

sum(myResult[,6] != myResult[,7])

## (1.4) Compute the mean of each row in myMA by applying the mean function in an apply loop

myve_apply <- apply(myMA, 1, mean)
myResult <- cbind(myMA, mean_for=myve_for, mean_while=myve_while, mean_apply=myve_apply)
myResult[1:4, ]

## (1.5) Compute the row means of myMA by taking advantage of the built-in looping behavior 
## of the mean() function. This will require the transpose of myMA and its conversion into 
## a data frame.

mymean <- mean(as.data.frame(t(myMA)))
myResult <- cbind(myMA, mean_for=myve_for, mean_while=myve_while, mean_apply=myve_apply, mean_int=mymean)
myResult[1:4, ]

## (1.6) How can one confirm that the values from all four mean calculations are identical?
## The unique() and length() functions can be used for this step.

myq <- apply(myResult[,6:9], 1, function(x) length(unique(x)))
sum(myq != 1)

## (1.7) Replace one mean value in the last column by '0' and repeat the identity check in step (1.6) 

myResult[1,8] <- 0
myq <- apply(myResult[,6:9], 1, function(x) length(unique(x)))
sum(myq != 1)

###################
## (2) Functions ##
###################

## (2.1) Save the above code to a script file named "myMeans.R". Then convert the code to a 
## function (myMeans) that returns by default the final data frame from step (1.5) with all 
## four mean calculations. To test and debug the function, source your script and then
## execute your myMeans function like this:

	## source("myMeans.R")
	## myMeans() 

## (2.2) Use the following code as basis to implement a function that allows the user 
## to compute the mean for any combination of columns in a matrix or data frame. 
## One argument of this function should specify the input data set and another should 
## allow the selection of the columns by providing a grouping vector (factor). 

myMA <- matrix(rnorm(100000), 10000, 10, dimnames=list(1:10000, paste("C", 1:10, sep="")))

myList <- tapply(colnames(myMA), c(1,1,1,2,2,2,3,3,4,4), list) 
names(myList) <- sapply(myList, paste, collapse="_")
myMAmean <- sapply(myList, function(x) apply(myMA[,x], 1, mean))
myMAmean[1:4,] 

myMAcomp <- function(myMA=myMA, group=c(1,1,1,2,2,2,3,3,4,4), myfct=mean) {
	myList <- tapply(colnames(myMA), group, list)
	names(myList) <- sapply(myList, paste, collapse="_")
	myMAmean <- sapply(myList, function(x) apply(myMA[,x], 1, myfct))
	return(myMAmean)
}
myMAcomp(myMA=myMA, group=c(1,1,1,2,2,2,3,3,4,4), myfct=mean)[1:4,] 

#######################################################
## (3) Improving Speed Performance by Avoiding Loops ##
#######################################################

## (3.1) Implement a speed improved version of your previous function from step (2.2).
## Use the following code for its implementation. In this code the R loop over the
## data intensive dimension is avoided by using the rowMeans() function. 

myList <- tapply(colnames(myMA), c(1,1,1,2,2,2,3,3,4,4), list)
myMAmean <- sapply(myList, function(x) rowMeans(myMA[,x]))
colnames(myMAmean) <- sapply(myList, paste, collapse="_")
myMAmean[1:4,]

## (3.2) Compare the speed performance of the implementation from step 2.2 with the 
## speed improved version from step (3.1) using the system.time() function. 

###############################################
## (4) Looping over lists with lapply/sapply ##
###############################################
## Create a sample list populated with character vectors of different lengths

setlist <- lapply(11:30, function(x) sample(letters, x, replace=TRUE))
names(setlist) <- paste("S", seq(along=setlist), sep="") 

## (4.1) Compute the length for all pairwise intersects of the vectors stored in 'setlist'.
## The intersects can be determined with the '%in%' function like this: 
##              sum(setlist[[1]] %in% setlist[[2]])

setlist <- sapply(setlist, unique)
olMA <- sapply(names(setlist), function(x) sapply(names(setlist), function(y) sum(setlist[[x]] %in% setlist[[y]])))
olMA 

## (4.2) Plot the resulting intersect matrix as heat map
## The heatmap.2() function from the gplots library can be used for this.
library("gplots")
heatmap.2(olMA, trace="none", Colv="none", Rowv="none", dendrogram="none", col=colorpanel(40, "darkred", "orange", "yellow")) 

######################################################################################
## (5) Batch operations on many files: import, export, calculations, renaming, etc. ##
######################################################################################
## Follow the instructions of this manual section:
## http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/R_Programming.html

########################################################
## (6) Basic Utilities for DNA Sequence Manipulations ##
########################################################

## (6.1) Write a RevComp() function that returns the reverse and complement of a DNA sequence.
## Include an argument that will allow to return only the reversed sequence, the complemented 
## sequence or the reversed and complemented sequence. The following R functions will be useful 
## for the implementation:
	## x <- c("ATGCATTGGACGTTAG") 
	## substring(x, 1:nchar(x), 1:nchar(x)) 
	## rev(x) 
	## paste(x, collapse="") 
	## chartr("ATGC", "TACG", x)

revComp <- function(myseq, type="RC") {
	seqmod <- myseq[1] # The substring function below works only on a single string.
	if(type=="R") {
		seqmod <- substring(seqmod, 1:nchar(seqmod), 1:nchar(seqmod))
		seqmod <- rev(seqmod)
		seqmod <- paste(seqmod, collapse="")  	
	}
        if(type=="C") { 
		seqmod <- chartr("ATGC", "TACG", seqmod) 
        }
	if(type=="RC") {
		seqmod <- chartr("ATGC", "TACG", seqmod) 
		seqmod <- substring(seqmod, 1:nchar(seqmod), 1:nchar(seqmod))
		seqmod <- rev(seqmod)
		seqmod <- paste(seqmod, collapse="")  	
	}
	return(seqmod)
}
cat("\n## (1) REVERSE AND COMPLEMENT OF DNA SEQUENCES")
cat("\n## (1.1) Creates random sample DNA sequences\n\t myDNA <- sapply(1:100, function(x) paste(sample(c(\"A\",\"T\",\"G\",\"C\"), 20, replace=T), collapse=\"\"))\n")
cat("\t\t # Creates 100 sequences each with 20 residues.\n")
cat("\n## (1.2) How to run revComp function\n\t revComp(myseq=myDNA[1], type=\"RC\")\n")
cat("\t\t # Option R: reverse \n\t\t # Option C: complement \n\t\t # Option RC: reverse and complement\n")

## (6.2) Write a function that applies the RevComp() function to many sequences stored in a vector.

cat("\n## (1.3) Apply function to many sequences with sapply/lapply \n\t sapply(myDNA, function(x) revComp(myseq=x, type=\"RC\"))\n")

## (6.3) Write a function that will translate one or many DNA sequences in all three reading frames into proteins. 
## The following commands will simplify this task:

	## AAdf <- read.table(file="AA.txt", header=T, sep="\t")  
        ## AAv <- AAdf[,2]; names(AAv) <- AAdf[,1] 
        ## y <- gsub("(...)", "\\1_", y) 
        ## y <- unlist(strsplit(y, "_"))
        ## y <- y[grep("∧...$", y)] 
        ## AAv[y]

## Import genetic code and store it in named vector AAv
cat("\n## (2) TRANSLATION OF DNA INTO AA SEQUENCES")
AAdf <- read.table(file="http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/AA.txt", header=T, sep="\t")
AAv <- AAdf[,2]; names(AAv) <- AAdf[,1]
cat("\n## (2.1) The genetic code was imported and is stored in named vector AAv\n")

## Define translate function
translateDNA <- function(myseq, frame=1, code=AAv) {
	seqmod <- myseq[1] # Function assumes only a single string
	tripSplit <- function(seqmod=seqmod) { # Internal function for converting DNA string to vector of triplets 
		seqmod <- gsub("(...)", "\\1_", seqmod) 
		seqmod <- unlist(strsplit(seqmod, "_")) 
		seqmod <- seqmod[grep("^...$", seqmod)] 
		return(seqmod)	
	}
	if(frame==1) {
		seqmod <- tripSplit(seqmod) 
	}
	if(frame==2) {
		seqmod <- gsub("^.", "", seqmod)
		seqmod <- tripSplit(seqmod) 
	}
	if(frame==3) {
		seqmod <- gsub("^..", "", seqmod)
		seqmod <- tripSplit(seqmod) 
	}
	myAA <- AAv[seqmod]
	myAA <- paste(myAA, collapse="")
	return(myAA)
}
cat("\n## (2.2) How to run translateDNA function \n\t translateDNA(myseq=myDNA[1], frame=1, code=AAv) \n")
cat("\t\t # The three possible reading frames (1, 2 or 3) can be selected under the frame option.\n")

## (3.3) Apply function to many sequences.
cat("\n## (2.3) Apply function to many sequences with sapply/lapply \n\t mypep <- sapply(myDNA, function(x) translateDNA(myseq=x, frame=1))\n")
cat("\n## (2.4) Return result in data frame \n\t data.frame(DNA=names(mypep), PEP=as.vector(mypep))[1:12,]\n")

## Non R codefor web upload:
## Convert to html  
##        :runtime! syntax/2html.vim
## Upload
##        scp R_programming_exercises.R tgirke@cache.ucr.edu:~/public_html/Documents/R_BioCond/My_R_Scripts/