Calculate and plot statistics on an R matrix

My goal is to iterate over a datamatrix, calculate mean, standard error of the mean, and then have a bar plot representing the 2 groups of my data matrix in R.

My code below works as a for loop. However, my initial goal was to output each 8 plots in one .png file. To that end, I posted the question here.

One of the 'comments' was to use an 'apply' and 'function'. I never wrote a code in R using apply or function.

Here's my input file:

TranscriptID    GeneID  Biotype TranscriptName  CommonNAme  GeneName    TSS-ID  Locus-ID    DNp63D-DMECs-1  DNp63D-DMECs-2  DNp63D-DMECs-3  DNp63WTMECs-1   DNp63WTMECs-2   Fold    2-tailedtest
Test1   TestA   protein_coding  Fun1    Ex1 Ex1 ExA1    ExA1    1.15E-08    2.68E-12    0.005077929 4.99E-07    6.38E-08    6.02E+03    0.495089687
Test2   TestB   protein_coding  Fun2    Ex2 Ex2 ExA2    ExA2    3.69E-08    0.014129129 0.075213367 0.121370367 0.404553833 1.13E-01    0.123434776
Test3   TestC   protein_coding  Fun3    Ex3 Ex3 ExA3    ExA3    4.89E-05    0   0   6.58E-05    1.64E-34    4.96E-01    0.643007583
Test4   TestA   protein_coding  Fun4    Ex4 Ex4 ExA4    ExA4    0.058629449 0   0   0.056200966 0.253314667 1.26E-01    0.180082201
Test5   TestB   protein_coding  Fun5    Ex5 Ex5 ExA5    ExA5    7.80E-06    0   0   1.42E-11    4.20E-36    3.66E+05    0.495026427
Test6   TestC   protein_coding  Fun6    Ex6 Ex6 ExA6    ExA6    0   0   0   0   2.41E-101   0.00E+00    0.272228401
Test7   TestA   protein_coding  Fun7    Ex7 Ex7 ExA7    ExA7    3.77E-08    0.023945749 0.077103517 0.262936167 0.2940195   1.21E-01    0.004479038
Test8   TestB   protein_coding  Fun8    Ex8 Ex8 ExA8    ExA8    9.30E-09    4.82E-14    0.000827853 8.19E-07    7.47E-07    3.52E+02    0.496141526


Here is my for loop code:

input <- read.delim(file="MECs-DNp63IsoformLevels.txt", header=TRUE, sep="\t")
input<-as.matrix(input)

for (i in 1:nrow(input)) {
mean1 <- mean(as.numeric(input[i,12:13]))
mean2 <- mean(as.numeric(input[i,9:11]))
sd1 <- sd(as.numeric(input[i,12:13]))
sd2 <- sd(as.numeric(input[i,9:11]))
sem1 <- sd2/sqrt(length(input[i,12:13]))
sem2 <- sd1/sqrt(length(input[i,9:11]))

mean_sem <- data.frame(mean=c(mean1, mean2), sem=c(sem1, sem2), group=c("WT", "DNp63D-D"))
mean_sem$group<-factor(mean_sem$group, levels=mean_sem\$group, ordered=TRUE) #this prevents ggplot from ordering the x-axis alphabaetically and keeps the order as the input dataframe
theme_set(theme_gray(base_size = 20))
print(i)
p<- ggplot(mean_sem, aes(x=group, y=mean)) +
geom_bar(stat='identity', width=.3, colour="black", fill=c("blue", "red")) +
geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem),
width=.2) +
geom_line(aes(colour=group)) +
scale_colour_manual(values=c("blue", "red")) +
xlab('Genotype of MECs') +
ylab('Quantile Norm FPKM')

q = p +ggtitle(input[i,5])

ggsave(filename=paste(input[i,5],'.png', sep=""), plot=q)

}


Here is my code with the apply function:

input <- read.delim(file="MECs-DNp63IsoformLevels.txt", header=TRUE, sep="\t")
input<-as.matrix(input)
apply(input, 1, function(input) { mean1=mean(as.numeric(input[9:11]))
mean2=mean(as.numeric(input[12:13]))
sd1= sd(as.numeric(input[9:11]))
sd2 = sd(as.numeric(input[12:13]))

sem1= sd1/sqrt(length(input[9:11]))
sem2= sd2/sqrt(length(input[12:13]))

mean_sem = data.frame(mean=c(mean1, sem1), sem=c(sd1, sem2), group=c("WT", "DNp63D-D"))

pdf("Test.pdf")
p=ggplot(mean_sem, aes(x=group, y=mean))+ geom_bar(stat='identity', width=.3, colour="black", fill=c("blue", "red"))+ geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem), width=.2)+geom_line(aes(colour=group)) + scale_colour_manual(values=c("blue", "red")) + xlab('Genotype of MECs') + ylab('Quantile Norm FPKM')
q = p +ggtitle(input[1])
ggsave(filename=paste(input[1],'.png', sep=""), plot=q) #plots each figure and gives it a name similar to column 1.
## I need it to plot each 8 figures in 1 png file

})


Can my for loop or function code be optimized? This is literally my first function/apply code I wrote in R and it was more trial and error on the syntax as I saw some people using 'c' with others using multiple {{}{} in each function.

• Can you elaborate on ## I need it to plot each 8 figures in 1 png file? Can you provide 16 lines of data so we see which column can be used to group rows 8 by 8? – flodel Feb 3 '16 at 12:29
• @floder. sure. The way the script above is written, it iterates over each row, then plots one barplot for each for row and outputs each in 1 png file. In fact my data consists of 100's or rows. hence, instead of having 100's of .png files, i was hoping to output 8 plots/png file. I listed an example of the expected outputs here . Does this answer your question? Also, do you think my current formatting/write up make sense? Thanks – BioProgram Feb 3 '16 at 16:20
• It you collect your data into a data.frame instead of plotting in the loop/*apply, you can use facet_wrap when plotting afterwards, which would likely give you what you're looking for. – alistaire Feb 10 '16 at 17:22
• @alistaire can you advise on how to collect it as dataframe? it keeps overwriting. – BioProgram Feb 10 '16 at 21:17
• If you're in an *apply function, you'll have to use <<- to break out of its environment so you can append your computations to something; with a for loop, you can just use <-. Note that pre-allocating a data structure of the correct dimensions can speed up your loop, if necessary. You could probably also refactor your code with dplyr, which is usually convenient for these kinds of operations. – alistaire Feb 10 '16 at 21:24

This is the first time I have answered anything on SO, and I am not really much of an authority on R or programming in general, so it this may not be optimal, but I thought I'd try and contribute for once.

Functions

Firstly I would create a function for the SEM, just because it is handy. I also created a function that wrapped that function and mean() to use in apply()

sem <- function(x){
sd(x) / sqrt(length(x))
}

sum.stats <- function(x){
m <- mean(x)
e <- sem(x)
c(mean = m, sem = e)
}


As you said, the apply() function overwrites the previous iterations. I think this is because functions return only the last instruction evaluated, so you have to include everything you want back in a return()

Create df

I applied sem() to each group of data (df1 is the original dataset) separately - you do both groups together as you did above, but this is more general

grp.a <- apply(df1[, 12:13], 1, sum.stats)
grp.b <- apply(df1[, 9:11], 1, sum.stats)


This produced a dataframe for each group, so these need to be combined for ggplot (Note the reversed expand.grid - by happy accident, I found that this ordered the groups non-alphabetically as you said you wanted)

df2 <- data.frame(rev(expand.grid(obs = 1:8, group = c('WT', 'DN'))),
rbind(t(grp.a), t(grp.b)))


Plot

For the plot, I am not sure about this, but I don't think apply() works with ggplot() because ggplot() needs a dataframe - apply() provides a vector.

Also I chose to do away with the legend since the groups are named in the axis labels anyway, and because I think geom_point() is more appropriate with those wide error bars - but of course, the choice is yours :)

library(ggplot2)

ggplot(df2, aes(group, mean)) +
geom_point() +
geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem), width = 0.2) +
xlab('Genotype of MECs') +
ylab('Quantile Norm FPKM') +
facet_wrap(~ obs, scales = 'free')

ggsave(filename='plot.png', width = 12, height = 8, units = 'in', dpi = 300)


Hope this is useful