Hi Ungku, it's really difficult for us to take a huge block of code and understand where an error happened. There's several things that can help us help you:
1) First and foremost, what is the error message or undesired behavior you're experiencing? 2) Second, please pare down the code to the place where you're experiencing a problem. Maybe just generate some simple data and try making the plot from there without all the calculations and formatting options, and see if it works then. If it still doesn't, post that simplified code and someone here will be much more able to help you. The posting guide, http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>, has some further tips. Best, Ethan On Sun, Jun 19, 2011 at 7:05 PM, Ungku Akashah <kasla...@yahoo.com> wrote: > HELLO, anybody... could you help me to check the below coding for volcano. > what is the mistake? > what the plot could not display? > > > > > > > > # volcano_plot.r > # > # Author: Amsha Nahid, Jairus Bowne, Gerard Murray > # Purpose: Produces a volcano plot > # > # Input: Data matrix as specified in Data-matrix-format.pdf > # Output: Plots log2(fold change) vs log10(t-test P-value) > # > # Notes: Group value for control must be alphanumerically first > # Script will return an error if there are more than 2 groups > > # > # Load the data matrix > # > # Read in the .csv file > data<-read.csv("input5.csv", sep=",", row.names=1, header=TRUE) > # Get groups information > groups<-data[,1] > # Get levels for groups > grp_levs<-levels(groups) > if (length(levels(groups)) > 2) > print("Number of groups is greater than 2!") else { > > # > # Split the matrix by group > # > new_mats<-c() > for (ii in 1:length(grp_levs)) > new_mats[ii]<-list(data[which(groups==levels(groups)[ii]),]) > > # > # Calculate the means > # > # For each matrix, calculate the averages per column > submeans<-c() > # Preallocate a matrix for the means > means<-matrix( > nrow = 2, > ncol = length(colnames(data[,-1])), > dimnames = list(grp_levs,colnames(data[,-1])) > ) > # Calculate the means for each variable per sample > for (ii in 1:length(new_mats)) > {submeans[ii]<-list(apply(new_mats[[ii]][,-1],2,mean,na.rm=TRUE)) > means[ii,]<-submeans[[ii]]} > > # > # Calculate the fold change > # > folds<-matrix( > nrow=length(means[,1]), > ncol=length(means[1,]), > dimnames=list(rownames(means),colnames(means)) > ) > for (ii in 1:length(means[,1])) > for (jj in 1:length(means[1,])) > folds[ii,jj]<-means[ii,jj]/means[1,jj] > > # > # t-test P value data > # > > > pvals<-matrix(nrow=ncol(data[,-1]),ncol=1,dimnames=list(colnames(data[-1]),"P-Value")) > > > # > # Perform the t-Test > # > for(ii in 1:nrow(pvals)) { > > pvals[ii,1]<-t.test(new_mats[[1]][,ii+1],new_mats[[2]][,ii+1])$p.value > } > > m<-length(pvals) > x_range<-range(c( > min( > range(log2(folds[2,])), > range(c(-1.5,1.5)) > ), > max( > range(log2(folds[2,])), > range(c(-1.5,1.5)) > ) > )) > y_range<-range(c( > min(range(-log10(pvals)), > range(c(0,2)) > ), > max(range(-log10(pvals)), > range(c(0,2)) > ) > )) > > # > # Plot data > # > # Define a function, since it's rather involved > volcano_plot<-function(fold, pval) > {plot(x_range, # x-dim > y_range, # y-dim > type="n", # empty plot > xlab="log2 Fold Change", # x-axis title > ylab="-log10 t-Test P-value", # y-axis title > main="Volcano Plot", # plot title > ) > abline(h=-log10(0.05),col="green",lty="44")# horizontal line at > P=0.05 > abline(v=c(-1,1),col="violet",lty="1343") # vertical lines at > 2-fold > # Plot points based on their values: > for (ii in 1:m) > # If it's below 0.05, we're not overly interested: purple. > if (-log10(pvals[ii])>(-log10(0.05))) { > # Otherwise, more checks; > # if it's greater than 2-fold decrease: blue > if (log2(folds[2,][ii])>(-1)) { > # If it's significant but didn't change much: orange > if (log2(folds[2,][ii])<1) { > points( > log2(folds[2,][ii]), > -log10(pvals[ii]), > col="orange", > pch=20 > ) > # Otherwise, greater than 2-fold increase: red > } else { > points( > log2(folds[2,][ii]), > -log10(pvals[ii]), > col="red", > pch=20 > ) > } > } else { > points( > log2(folds[2,][ii]), > -log10(pvals[ii]), > col="blue", > pch=20 > ) > } > } else { > points( > log2(folds[2,][ii]), > -log10(pvals[ii]), > col="purple", > pch=20 > ) > } > } > # Plot onscreen via function > x11() > volcano_plot(folds,pvals) > > # Return table to analyse results > > # > # Generate figures as image files > # > # (Uncomment blocks as necessary) > > ##### jpeg ##### > # pic_jpg<-function(filename, fold, pval) > # {# Start jpeg device with basic settings > # jpeg(filename, > # quality=100, # image quality > (percent) > # bg="white", # background colour > # res=300, # image resolution > (dpi) > # units="in", width=8.3, height=5.8) # image dimensions > (inches) > # par(mgp=c(5,2,0), # axis margins > # # (title, labels, line) > # mar=c(6,6,4,2), # plot margins > (b,l,t,r) > # las=1 # horizontal labels > # ) > # # Draw the plot > # volcano_plot(folds, pvals) > # dev.off() > # } > # pic_jpg("volcano_plot.jpg") > ##### end jpeg ##### > > > #### png ##### > # pic_png<-function(filename, fold, pval) > # {# Start png device with basic settings > # png(filename, > # bg="white", # background colour > # res=300, # image resolution > (dpi) > # units="in", width=8.3, height=5.8) # image dimensions > (inches) > # par(mgp=c(5,2,0), # axis margins > # # (title, labels, line) > # mar=c(6,6,4,2), # plot margins > (b,l,t,r) > # las=1 # horizontal labels > # ) > # # Draw the plot > # volcano_plot(folds, pvals) > # dev.off() > # } > # pic_png("volcano_plot.png") > #### end png ##### > > > # #### tiff ##### > # pic_tiff<-function(filename, fold, pval) > # {# Start tiff device with basic settings > # tiff(filename, > # bg="white", # background colour > # res=300, # image resolution > (dpi) > # units="in", width=8.3, height=5.8) # image dimensions > (inches) > # compression="none" # image compression > # # (one of none, lzw, > zip) > # par(mgp=c(5,2,0), # axis margins > # # (title, labels, line) > # mar=c(6,6,4,2), # plot margins > (b,l,t,r) > # las=1 # horizontal labels > # ) > # # Draw the plot > # volcano_plot(folds, pvals) > # dev.off() > # } > # pic_tiff("volcano_plot.tif") > # #### end tiff ##### > > > > > # > # Legacy code which allows for blue/red to be independent of 0.05 > # (purple is limited to the middle lower region) > # > ##### > # for (ii in 1:m) > # if (log2(folds[2,][ii])<1) { > # if (log2(folds[2,][ii])>-1) { > # if (-log10(pvals[ii])<(-log10(0.05))) { > # points( > # log2(folds[2,][ii]), > # -log10(pvals[ii]), > # col="purple", > # pch=20 > # ) > # } else { > # points( > # log2(folds[2,][ii]), > # -log10(pvals[ii]), > # col="orange", > # pch=20 > # ) > # } > # } else { > # points( > # log2(folds[2,][ii]), > # -log10(pvals[ii]), > # col="blue", > # pch=20 > # ) > # } > # } else { > # points( > # log2(folds[2,][ii]), > # -log10(pvals[ii]), > # col="red", > # pch=20 > # ) > # } > > # If function from above needs to be closed > } > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
