My Computational Journal

Week 1

May 20, 2011

The R statistical package was used to read nonnormalized data in the form of a GPR file for all of the microarray chips prepared in the lab. It was found that the script used yesterday to read the nonnormalized data and then to normalize the data was incorrect. Without clearly specifying whether the mean or median values for each labelling dye should be used, R would with some parts of the script use mean values and with others the median values. In going back to the data today, we specified that we wanted the median values to be used in calculating log ratios. Specifically, we used R to calculate the log base 2 ratio of red to green dye was calculated for the Dahlquist wildtype flask 3 t15/t0 microarray chip . The ratio was calculated by subtracting the median value for the red background from the median value of the red foreground and then this difference was divided by the difference between the median value for the green foreground and the median value for the green background. The negative log ratio was calculated because the chip being analyzed was dye swapped. The first five log ratios calculated by R (two of which were duplicates of two of the three genes looked at) were compared to the first five log ratios in the GPR file corresponding to the chip being analyzed. It was observed that the log ratios of R were the negative values of the log ratios in the GPR file which was expected because of the dye swap done on the chip. However, we should ask Dr. Dalhquist why she used a red to green ratio instead of a green to red ratio for the dye swaps. Also, the log ratios generated by R were then compared to the log ratios within an Excel file of the log ratios for all of the chips. There were no duplicate genes in the Excel file. Therefore, an average of the log ratios calculated in R for each of the two duplicated genes was found and compared to the log ratios within the Excel file. It was observed that the former log ratios did not correspond to the log ratios for the three genes on the Dahlquist wildtype flask 3 t15/t0 chip. We need to talk to Dr. Dahlquist regarding how she did her calculations to maybe shed some light on why the log ratios calculated by R did not match the log ratios that were provided in Excel.

Katrina Sherbina 16:38, 20 May 2011 (EDT)

Week 2

May 25, 2011

Scatter plots were generate for normalized log fold change ratios versus nonnormalized fold change ratios for Dahlquist wildtype flask 3 t15/t0 chip, Dalhquist wildtype flask 3 t90/t0 chip, delGLN3 flask 1 t60/t0, delGLN3 flask 1 t30/t0, and delGLN3 flask 2 t15/t0 as indicated in GPR files. While the scatter plot for the wildtype flask showed a linear relationship between nonnormalized and normalized log fold change ratios with an R^2 value of 1, the scatter plots for each of the delGLN3 flasks did not display a very linear relationship between nonnormalized and normalized log fold change ratios. In the scatter plots for the deletion strain, the data points were very dispersed around the origin. The plot for delGLN3 flask 1 showed much more spread and had a lower r^2 value than the plot for the delGLN3 flask 2. After observing the nonlinear nature of the scatter plots for the deletion strain, scatter plots were generated for the normalized values for the red channel versus the nonnormalized values for the red channel and normalized values for the green channel versus the nonnormalized values for the green channel for delGLN3 flask 1 t60/t0, delGLN3 flask 1 t30/t0, and delGLN3 flask 2 t15/t0. In all of the scatter plots, the data points represented something like a diagonal arrow facing down toward the origin. The data points were widely dispersed close to the origin but then began to display a linear relationship as the nonnormalized value for the color increased.

In addition, data from GPR files for the wildtype strain was normalized using the R statistical package. A source code was written that performed normalization within arrays and normalization between arrays and also created MA plots before and after within array normalization and a box plot of log fold change ratios (M) after within array normalization. While both normalizations were performed and the MA plots were generated, the box plot could not be generated. In addition, we were not able to write the results of the normalizations to a cvs.file using the write.table command.

Katrina Sherbina 20:42, 25 May 2011 (EDT)

May 26, 2011

Using the data generated by completing the within array normalization yesterday, all of the M values for replicate genes were average using the tapply code in R [tapply(MA$M,as.factor(MA$genes[,5]),mean)], which Dr. Fitzpatrick supplied. In writing this data to a table, it was realized that the way in which the code was written told R to only average for the replicate genes in only one microarray chip (flask 3 t30/t0). Then, a matrix of the M values was created to which a for loop was applied so that the M values for the replicate genes were average for all microarray chips. Then, a scatter plot was create in Excel comparing these normalized values to nonnormalized values provided in an Excel sheet by Dr. Dahlquist. There was no linear relationship between the nonnormalized and normalized values.

Then, only one GPR file corresponding to one chip (flask 4 t30/t0) was looked at. This GPR file was read by R to create an RG file. Then a code was written to get rid of the data for the Arabadopsis and 3XSSC genes. Then the log fold changes for all of the other genes was calculated in R by finding the log 2 ratio of the difference of the green and green background divided by the difference of the red and red background. The green to red ratio was found instead of the red to green ratio because flask 4 t30/t0 is a dye swap. Then the log ratios for the replicate genes were average. Using a scatter plot, these log ratios were compared to the log fold change ratios averaged for replicates generated by the within array normalization. This scatter plot showed that there was a strong linear relationship between the nonnormalized data (the log fold change generated in R) and the within array normalized data.

targets<-readTargets(file.choose())
f <- function(x) as.numeric(x$Flags > -99)
RG <- read.maimages(targets, source="genepix.median", wt.fun=f)
RG2<-RG[RG$genes$Name!="Arabidopsis",]
RG3<-RG2[RG2$genes$Name!="3XSSC",]
GeneList<-RG3$genes$Name
lr<-log2((RG3$G-RG3$Gb)/(RG3$R-RG3$Rb))
lravg<-tapply(lr,as.factor(GeneList),mean)
write.table(lravg,"Log_Ratios_RG_Avg_Rep_Minus_Controls.csv",sep=",")

Writing the within array normalized data using the write.table function was still unsuccessful.

Katrina Sherbina 20:24, 26 May 2011 (EDT)

May 27, 2011

The tapply function in R was applied to the log fold change ratios from the nonnormalized GPR files for Flask 3 t15/t0 and Flask 4 t15/t30 to average all of the replicates. These log fold change ratios were all multiplied by -1 because both Flask 3 and Flask 4 were dye swapped. The corrected log fold change ratios were compared to the log fold change ratios in a Master Data sheet provided by Dr. Dahlquist. A scatter plot was created comparing the corrected log fold change ratios versus the log fold change ratios from the Master Data sheet. The scatter plot displayed what looked to be like a linear relationship with the exception that the data points were more dispersed around 0. The tapply function in R was also applied to the log fold change ratios from the normalized GPR files for Flask 3 t15/t0 and Flask 4 t15/t30 to average all of the replicates. These log fold change ratios were all multiplied by -1 because both Flask 3 and Flask 4 were dye swapped. Due to the way in which the Master Data was compiled, it was predicted that these corrected log fold changes would be of the same value and magnitude as the log fold changes in the Master Data. In comparing the corrected log fold changes versus the log fold changes in the Master Data, this was found to be case.

Also, the code developed yesterday to analyze the GPR files for the wild type strains was built upon. The matrix generated by creating a for loop for the within array normalization data was able to be saved to a csv.file and viewed within Excel. A matrix was also created for the between array normalization data using a for loop. This matrix was also successfully written to a .csv file.

targets<-readTargets(file.choose())
f <- function(x) as.numeric(x$Flags > -99)
RG <- read.maimages(targets, source="genepix.median", wt.fun=f)
par(mfrow=c(1,2))
plotMA(RG, main="Before Within Array Normalization")
MA<-normalizeWithinArrays(RG,method="loess",bc.method="normexp")
plotMA(MA, main="After Within Array Normalization")
M1<-tapply(MA$M[,1],as.factor(MA$genes[,5]),mean)
n1<-length(M1)
n0<-length(MA$M[1,])
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for (i in 1:14) {MM[,i]<-tapply(MA$M[,i],as.factor(MA$genes[,5]),mean)}
write.table(MM,"WANorm_Avg_Reps_Log_Fold_Change_Matrix.csv",sep=",")
MAScale<-normalizeBetweenArrays(MA, method="scale")
M2<-tapply(MAScale$M[,1],as.factor(MAScale$genes[,5]),mean)
n2<-length(M2)
n0<-length(MAScale$M[1,])
MM2<-matrix(nrow=n1,ncol=n0)
MM2[,1]=M2
for (i in 1:14) {MM2[,i]<-tapply(MAScale$M[,i],as.factor(MAScale$genes[,5]),mean)}
write.table(MM2,"BANorm_Avg_Reps_Log_Fold_Change_Matrix.csv",sep=",")

In addition, box plots were generated for the data before between array normalization, which was also after within array normalization, and for after between normalization.

par(mfrow=c(3,6))
for(i in 1:14){boxplot(MM[,i], ylim=c(-4,10))}  (Before within array normalization)
for(i in 1:14){boxplot(MM2[,i], ylim=c(-5,12))} (After within array normalization)

Then this code involving within array normalization, between array normalization, for loops for both types of normalization, MA plots, and box plots was modified to analyze the data from GPR files for the dGLN3 strain and the dCIN5 strain.

So far, only the GPR files from Ontario microarray chips have been normalized. The next step is to find a way in R to normalize the GPR files from the GCAT microarray chips.

Katrina Sherbina 19:19, 27 May 2011 (EDT)

Week 3

June 1, 2011

One Master Excel workbook was created with a spreadsheet of the GCAT chip normalizations that a coworker completed today, a spreadsheet of all of the within array normalizations already done last week for all of the Ontario chips, a spreadsheet with the between array normalizations for the Ontario chips, and also a spreadsheet that integrated both the GCAT chip normalizations and the within array normalizations for the Onatario chips.

To create the integrated spreadsheet mentioned above, Microsoft Access was used to merge the normalized data from the GCAT and the Ontario chips eliminating any genes in the GCAT chips that were not also in the Ontario chips. The steps are as follow.

1. Save the Excel files of the data for the GCAT and Ontario chiops as tab delimited files.
2. Create a new database on Access.
3. Import the data (File->Get External Data->Import)
4. Go through the import Wizard: specify the data as delimited, keep the delimiter as tab and the text qualifier as none and indicate that the first row contains field names, and choose my primary key as the ID names (the genes). The import wizard is gone through twice once with the GCAT data and once with the Ontario data.
5. In the window for the current database, go to queries and select "Create query in Design View".
6. Add both imported tables (the GCAT and Ontario).
7. In the "Select Query" window, join GCAT ID and Ontario ID with a line. Right click on the line and press the "Join Properties" option and then third option (to include all records from the Ontario data and only those of the GCAT data that are also within the Ontario data).
8. Select all of the fields in the GCAT query window and drag into the first box in the "Field" row in the table below.
9. Select all of the fields in the Ontario query window and drag into the next free box in the "Field" row in the table below.
10. Create a new table for this joined data (Query->Make-Table Query).
11. Copy and pastet the new table into a new Excel spreadsheet.

Katrina Sherbina 20:27, 1 June 2011 (EDT)

June 2, 2011

It was found that the MA plots generated before and after within array normalization only corresponded to the first microarray chip (first GPR file) in the targets file imported into R. A new code was written to generate MA plots before and after within array normalizations for all of the GPR files in the targets file. The number of rows and columns of MA plots, the number of iterations for the for loop, and the limits of the y-axis were altered for each strain.

par(mfrow=c(3,5))
for (i in 1:14) {plotMA(RG[,i],ylim=c(-4,4))}
for (i in 1:14) {plotMA(MA[,i])}

Originally, individual graphs of boxplots were generated side by side for before and after between array normalization for all the GPR files in the targets file. A new code was written to generate all of the boxplots for all of the GPR files in one graph. The limits for the y-axis were changed for each strain.

x<-as.matrix(MA$M)
boxplot(x[,1],x[,2],x[,3],x[,4],x[,5],x[,6],x[,7],x[,8],x[,9],x[,10],x[,11],x[,12],x[,13],x[,14],ylim=c(-6,6))
y<-as.matrix(MAScale$M)
boxplot(y[,1],y[,2],y[,3],y[,4],y[,5],y[,6],y[,7],y[,8],y[,9],y[,10],y[,11],y[,12],y[,13],y[,14],ylim=c(-6,6))

In addition, box plots were generated in one graph for all the GPR files in the targets file before any kind of normalization. The code was similar to the code for the boxplot above with the exception that in place of the x or y values the log base 2 ratio of red-red background to green-green background was found for each individual microarray chip. For chips that word dye swapped, the negative log base 2 ratio was found.

Katrina Sherbina 19:34, 2 June 2011 (EDT)

Weel 4

June 8, 2011

Using the multtest package that comes with the R program, we began calculating p values to analyze whether or not genes that were differentially expressed through experimentation were differentially expressed by chance or not. Only the between array normalized data from the wildtype strains from time point t15 were used. This data was extracted from its master file and written to a new .csv file. This file was then imported into R. In multtest, the method used was the step-down minP and the test used was the one sample t test. When first running the code, few numerical p values were generated. Mostly the output was a series of NaN's. In modifying the code to include standardize=FALSE, multtest was able to successfully produce p-values.

data1<-read.csv("wt_t15_normalized.csv",header=TRUE,sep=",")
data2<-as.matrix(data1[,2:4])
data3<-rep(0,length(data1[,2:4]))
seed<-99
exp<-MTP(X=data2,Y=data3,test="t.onesamp",B=100,method="sd.minP",seed=seed,standardize=FALSE)
write.table(as.matrix(exp@rawp),"wt_t15_pvalues.csv",sep=",")

The next step is to compare the raw p-values to p-values generated by hand with the same data in Excel using Dr. Dahlquist's method.

Katrina Sherbina 20:07, 8 June 2011 (EDT)

June 9, 2011

The T statistic and P values for the data collected for the wild type strains for all replicates at t15 were calculated in Excel using the formulas AVERAGE((range of cells)/(STDEV(range of cells)/SQRT(number of replicates)) and TDIST(ABS(cell containing T statistic),degrees of freedom,2) where degrees of freedom is the number of replicates minus 1, respectively. The Bonferroni correction was also calculated for each data point.

Using the multtest package in R, the raw p values and the adjusted p values were collected for the same data mentioned above using multiple testing procedues such as the single-step maxT (ss.maxT), single-step minP (ss.minP), step-down maxT (sd.maxT), and step-down minP (sd.minP). Also, raw p and adjusted p values were calculated by controlling the false discover rate.

data1<-read.csv("wt_t15_normalized.csv",header=TRUE,sep=",")
data2<-as.matrix(data1[,2:4])
data3<-rep(0,length(data1[,2:4]))
seed<-99
exp<-MTP(X=data2,Y=data3,test="t.onesamp",standardize=FALSE,typeone="fdr",fdr.method="conservative",B=3000,seed=seed)
exp1<-MTP(X=data2,Y=data3,test="t.onesamp",standardize=FALSE,method="sd.minP",B=3000,seed=seed)
exp2<-MTP(X=data2,Y=data3,test="t.onesamp",standardize=FALSE,method="sd.maxT",B=3000,seed=seed)
exp3<-MTP(X=data2,Y=data3,test="t.onesamp",standardize=FALSE,method="ss.minP",B=3000,seed=seed)
exp4<-MTP(X=data2,Y=data3,test="t.onesamp",standardize=FALSE,method="ss.maxT",B=3000,seed=seed)

Also four graphical summaries were created for the results of the ss.minP: number of rejected hypotheses vs. Type I error rate, sorted adjusted p-values vs. number of rejected hypotheses, adjusted p-values vs. test statistics, and adjusted p-values versus index.

The next step will be to compare the raw and adjusted p values calculated by the R multiple testing procedures with the p value and Bonferroni corrections calculated in Excel, respectively.

Katrina Sherbina 19:43, 9 June 2011 (EDT)

June 10, 2011

For wildtype t15 data, a scatter plot was created comparing the t-statistic derived raw p values with the raw p values obtained from performing a multtest using FDR and found that the t-stat derived raw p values were more conservative than the FDR raw p values. Also, for the dCIN5 data for all time point, f-statistic derived raw p values were calculated. These raw p values were comparied to the raw p values from a multtest using FDR, ss.minP, and sd.minP in three separate scatter plots. Each of these three scatter plots showed no correlation between the f-statistic derived raw p values and any of the p values calculated by the different multtests performed. Also, a benjamini hochberg correction wqas applied to the f-statistic derived p values. These were then compared to the adjusted p-values calculated using the aformentioned multtests. Again, it was not possible to discern any relationship between the values.

In addition, Dr. Fitzpatrick suggested comparing t-statistic derived p values for the dCIN5 data to the f-statistic derived p calues for the same data. However, if the same t-statistic calculations performed on the wildtype strain t15 data are performed on the dCIN5 data, this comparison cannot be made because while the f-statistic derived p values take into account all time points, at least as was calculated today, the t-statistic derived p values are for each individual time point.

Katrina Sherbina 19:48, 10 June 2011 (EDT)

Week 5

June 13, 2011

Several errors were found in the master excel spread sheet, the compilation of all of the normalizing done with all of the GPR files in R. To remedy these errors, within array and between array normalization was performed on the GCAT chips and the new numbers were inputted into the master excel spread sheet. In addition, the within array and between array normalized data for dZAP1 was corrected within the master excel spread sheet.

In preparation for adjusting the current model for the gene regulatory network, means and standard deviations were calculated with the between array normalized data for fifteen genes in the wildtype, dCIN5, dGLN3, dHMO1, and dZAP1 strains. An original spread sheet for the modeling was modified to create five separate Excel spread sheets for each of the five strains inputting the means calculated as aforementioned for the log 2 concentrations and the standard deviations calculated as aforementioned for the concentration sigmas.

Katrina Sherbina 19:40, 13 June 2011 (EDT)

June 14, 2011

Each of the five excel spread sheets created yesterday (with the log 2 concentrations and concentration sigmas) were used in running already created estimation codes for MatLab. For each strain, a graph was created modeling the the log fold changes versus time for each of the fifteen yeast genes focused on. While the graphs of the wildtype strain showed marked differences in the log fold changes with time, many of the graphs for the deletion strains did not.

We began to look at changing the alpha values for the wildtype and dGLN3 data to find an alpha value that gives a middle ground between a large least squares error but a small penalty and a small least squares error but a large penalty. We ran the estimation codes for MatLab again for the two strains for six different alpha values (0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001). An alpha value of 0.005 was used the first time the estimation codes were run for all strains.

Katrina Sherbina 21:13, 14 June 2011 (EDT)

June 15, 2011

In analyzing the least squares error and penalty data computed for six different alpha values yesterday, we decided to try three intermediate alpha values between 0.005 and 0.001, namely 0.002, 0.003, and 0.004, for both the wildtype and dGLN3 data. A scatter plot was generated plotting the least squares error versus the penalty values for all ten alpha values. The optimal alpha value was found to be 0.002.

Next, we computed the relative error for the production rates, network weights, and network thresholds for the alpha values aforementioned for the wildtype and dGLN3 data. This error was calculated by comparing the values generated by an alpha value of 0.002 with the values generated by each of the other alpha values. The error was scaled by the sum of the squared sum of the values generated by an alpha value of 0.002 and he squared sum of the values generated by the other alpha values being compared. There was a general trend between the relative errors calculated for the production rates, network weights, and network thresholds. The relative error would decrease as the alpha value initially decreased but then it began to increase with decreasing alpha values beginning at the middle alpha value. The goal of calculating these relative error values is to determine whether the production rates or the network weights or the network thresholds will be the most difficult to estimate.

Katrina Sherbina 20:08, 15 June 2011 (EDT)

June 17, 2011

For the dGLN3 and wiltype data, the production rates, network weights, and network thresholds were graphed versus the ten alpha values that have been used with the MatLab estimation codes. It was found that the network weights were the most dispersed. In addition, the production rates, network weights, and network thresholds generated by the alpha values of 0.001, 0.002, and 0.005 were each graphed in separate scatter plots against a numerical index. Within each scatter plots, there were similar trends between the data generated by the alpha values.

Yesterday, it was found that we did not have the GPR's from flask 3 for dGLN3 when performing the normalization. Therefore, within array and between array normalization was performed within R again with the additional data. The between array normalized data was used to compute the average log fold change and standard deviation for the fifteen genes constituting the current predicted regulatory network. The comparisons and plots mentioned in the previous paragraph were not generated for the corrected between array normalized data because both the wildtype and previous dGLN3 data showed simliar trends in these plots and comparisons. Therefore, it was decided in conjuction with Dr. Fitzpatrick that it was not necessary to generate these plots again.

YEASTRACT was used to generate a new network between the fifteen transcription factors and target genes already being studies. In addition a new network was generated in YEASTRACT with these transcrption factors and genes in addition to GLN3, HMO1, and ZAP1. New excel data sheets were created to use with MatLab code using the network determined by YEASTRACT. In addition, degradation rates for GLN3, HMO1, and ZAP1 were inputed from another research papaer and production rates were calculated by doubling the degradationr rates. In addition, the average log fold change and standard deviation was calulated for these three genes added to the network for the wildtype and each of the deletion strains.

We began running the estimation codes in MATLAB again with these new Excel sheets. We will be running further simulations to test which alpha value will be the best to use in additional simulations with the new network.

Katrina Sherbina 19:45, 17 June 2011 (EDT)

Week 6

June 20, 2011

The first thing done today was changing all of the RSC1 labels (RSC1 being a synonym for AFT1) to AFT1. The 15 gene and 18 gene networks previously generated with YEASTRACT were done so using documented regulations only from direct evidence. Today we generated new 15 gene and 18 gene networks in YEASTRACT using documented regulations from both direct and indirect evidence. In changing the RSC1 labels to AFT 1, the networks generated for the 15 and 18 gene networks using direct evidence were different than those previously generated. Therefore, the MatLab simulations that were run last Friday have to be redone. In addition, we generated two networks with a set of 21 genes (this network contains all the genes from the 15 gene network in addition to 6 more genes). Both generated using YEASTRACT, one was generated with documented regulations from only direct evidence and the other was generated with documented regulations from both direct and indirect evidence.

New excel spreadsheets were generated for each of the new networks generated today in preparation for running new MatLab simulations. For the 21 gene networks, the average log fold changes and the standard deviations were calculated for the 6 genes added to the previous 15 genes to create the 21 gene network. For any of the genes whose transcripts did not have a half life value in the source being used to calculate degradation rates, a half life of 25.5 minutes was used. We used the degradation rate calcualted using a half life of 25.5 minutes to also calculated the production rate.

We also experimented with what changes are made to the network generated in YEASTRACT for the original 15 genes being used when different specifications are choosen (such as documented regulations or potential regulations and the specifications within them). We are still trying to figure out where exactly YEASTRACT is getting its direct and indirect evidence from.

Katrina Sherbina 20:46, 20 June 2011 (EDT)

June 21, 2011

The objective was to determine the most appropriate alpha value to use for the new 15 gene network. New MatLab simulations were done with the new 15 gene network obtained by using direct evidence only and the wildtype and dGLN3 data for five alpha values (0.01, 0.005, 0.002, 0.001, and 0.0005). The least squares error was graphed versus the penalty for both the widltype and dGLN3 strains. Upon looking at the scatter plot for the dGLN3 data, it was decided to also try alpha values of 0.003, 0.0015, and 0.0001. When comparing the scatter plot for the wildtype data with that for the dGLN3 data, the two scatter plots seem to suggest that two different alpha values are most appropriate for running further simulations. It will need to be decided whether further simulations must be run to determine the most appropriate alpha value to use for the new 15 gene network generated by using both direct and indirect evidence.

Then the next step is to determine appropriate alpha values by running similar simulations with the 18 gene and 21 gene networks.

Katrina Sherbina 20:47, 21 June 2011 (EDT)

June 22, 2011

We run simulations in MatLab for all strains for the two new 15 gene networks (direct evidence only and direct and indirect evidence). We began running simulations for the two new 18 gene networks (direct evidence only and direct and indirect evidence). All of these simulations were performed using an alpha value of 0.002.

Replicate simulations were performed for wildtype and dGLN3 for the two new 15 gene networks (direct evidence only and direct and indirect evidence) using the alpha value 0.002 to confirm that the optimized network weights and network thresholds are the same within the replicate experiments within the individual networks.

Katrina Sherbina 21:33, 22 June 2011 (EDT)

June 24, 2011

First, we revisited and edited the comparisons we began two days ago for replicate simulations. In contrast to the intial comparisons we made, we made error calcuations and other calculations within Excel to compare the optimized network weights and network thresholds between the networks within the same trial simulation.

1. Calcualted the absolute value of each entry within each of the four matrices (direct network trial 1, direct & indirect network trial 1, direct network trial 2, direct & indirect trial 2) for the wildtype and dGLN3. The sum of each of these matrices was found.
2. The sum A was found of corresponding networks within the same trial (i.e. sum of direct network trial 1 and direct & indirect network trial 1).
3. A new matrix consisting of relative errors was calculated for each trial by finding the difference of corresponding network matrix entries (i.e. difference between first entry in direct network trial 1 matrix and the first entry in direct & indirect network trial 2) and dividing by the sum A.
4. Cells that corresponded to a network connection that's only in one of the networks (direct vs. direct plus indirect) were colored with a yellow fill. Cells that were common to both networks were colored with an orange fill.
5. The signs of the weights and the thresholds were compared between trials within the same network and also between two networks within the trial simulation. A new matrix was generated in which each cell represented a comparison of the signs using the Excel function =if(sign(cell 1)=sign(cell 2),0,1). The sum of this matrix was found to determine how many sign changes occurred.

We ran MatLab simulations for the wildtype for the 15 gene direct and direct plus indirect networks using alpha values of 0.005 and 0.001. Then, we repeated the same calculations above for this data. Due to time constraints, we could not do the same for dGLN3 today. However, we plan to run these simulations and perform the calcuations for the dGLN3 data next week.

We also redid normalization in R for both dHMO1 and dZAP1 because we observed sign changes at some time points when looking at regulatory pathways in GeneMapp modeled with Dr. Dahlquist's original data versus those modeled by with the normalized data in R. Unfortunately we came upon a disturbing discovery when comparing the normalization we did today with the normalized data in the Excel spreadsheet. We found that the normalization changed when the controls 3XSSC and Arabidopsis on the microarray were removed before normalization in comparison to when the controls were kept for the normalization. We found indiscrepencies in the sign of the log fold changes between two sets of normalized data both of which were presumably normalized with the controls. We are currently consulting with Dr. Fitzpatrick and Dr. Dahlquist to determine whether or not the controls should be kept when normalizing the data. We also have to determine why we are getting different results with normalization when performing replicate processes. Consequently, within and between array normalization will be redone next week for all strains and compared to the data we currently have to see if there are any inconsistencies. This also means that the MatLab simulations that have been performed this week may need to be redone. Despite this, the comparisons between replicate experiments that have been mentioned above will be kept because these comparisons look at the output and are thus indpendent of the input. The same trends should be observed in the comparisons if new data is inputted into MatLab to perform replicate simulations.

Week 7

June 27, 2011

My coworker and I went through all of the data for all five strains and went through the normalization process in R all over again individually.

The first set of code was used to create an MA plot of the data for each strain before any normalization was done.

par(mfrow=c(4,4))

This designates how many rows and columns of graphs should be fitted onto one window. The numbers changing depending on how many target files are being imported into R.

targets<-readTargets(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RG<-read.maimages(targets, source="genepix.median", wt.fun=f)
dyeswap<-c(1,-1,-1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1)

The last line of code is necessary to designate which gpr files corresponding to microarry chips that have been dyeswapped. The "-1" represents a dye-swapped chip. This line of code needs to be changed for each strain in order to accurately reflect which gpr files correspond to chips that are dyeswapped for specific deletion strains.

GeneList<-RG$genes$Name
for(i in 1:14) {lr<-dyeswap[i]*(log2((RG$R-RG$Rb)/(RG$G-RG$Gb)))}
r0<-length(lr[1,])
RX<-tapply(lr[,1],as.factor(GeneList),mean)
r1<-length(RX)
M<-matrix(nrow=r1,ncol=r0)
M[,1]=RX
for(i in 1:14) {M[,i]<-tapply(lr[,i],as.factor(GeneList),mean)}

The "1:14" in the for loop designiates how many target files there are. This also has to be altered for data from each strain.

la<-(1/2*log2((RG$R-RG$Rb)*(RG$G-RG$Gb)))
r3<-length(la[1,])
RQ<-tapply(la[,1],as.factor(GeneList),mean)
r4<-length(RQ)
A<-matrix(nrow=r4,ncol=r3)
A[,1]=RQ
for(i in 1:14) {A[,i]<-tapply(la[,i],as.factor(GeneList),mean)}
for(i in 1:14) {plot(A[,i],M[,i],xlab="A",ylab="M",ylim=c(-5,5),xlim=c(0,15))}

The "1:14" in the four loop designiates how many target files there are. This also has to be altered for data from each strain. The y and x limits to the plots have been kept the same for all of the strains as decided upon looking at the data. However, these limitations may need to be changed with new data. The first line in the above set of code was used to find the average log fold change for each spot. All of the duplicate spots (spots on the microarray that corresponding to the same gene) were averaged using the second to last line of code.

This next set of code was used to create an MA plot of the data for each strain after within normalization was applied. For within array normalization, the controls were retained.

par(mfrow=c(4,4))
MA<-normalizeWithinArrays(RG, method="loess", bc.method="normexp")
M1<-tapply(MA$M[,1],as.factor(MA$genes[,5]),mean)
n0<-length(MA$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for(i in 1:14) {MM[,i]<-tapply(MA$M[,i],as.factor(MA$genes[,5]),mean)}
write.table(MM,"New_wt_WANorm_Matrix.csv",sep=",")
X1<-tapply(MA$A[,1],as.factor(MA$genes[,5]),mean)
y0<-length(MA$A[1,])
y1<-length(X1)
AA<-matrix(nrow=y1,ncol=y0)
AA[,1]=X1
for(i in 1:14) {AA[,i]<-tapply(MA$A[,i],as.factor(MA$genes[,5]),mean)}
for(i in 1:14) {plot(AA[,i],MM[,i],ylab="M",xlab="A",ylim=c(-5,5),xlim=c(0,15))}

Before between array normalization could be applied, the .csv file that was written had to be accessed. The first two rows, corresponding to the controls, were deleted. In addition the Master Index was also deleted.

MA2<-read.csv("New_wt_WANorm_Matrix.csv",header=TRUE,sep=",")
MA3<-as.matrix(MA2)
MAB<-normalizeBetweenArrays(MA3,method="scale",targets=NULL)
MAC<-as.matrix(MAB)
write.table(MAC,"New_wt_BANorm_Matrix.csv",sep=",")

The first line of the above set of code reinputted the .csv file that was edited to remove the controls. The second line of code transformed this data into a matrix so that between array normalization could be performed.

In addition to normalizating and creating the MA plots, boxplots were also created for each strain graphing the average log fold change.

par(mfrow=c(1,3))
targets<-readTargets(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RG<-read.maimages(targets, source="genepix.median", wt.fun=f)
dyeswap<-c(1,-1,-1,1,-1,-1,1,-1,1,-1,-1,1,-1,-1)
GeneList<-RG$genes$Name
for(i in 1:14) {lr<-dyeswap[i]*(log2((RG$R-RG$Rb)/(RG$G-RG$Gb)))}
r0<-length(lr[1,])
RX<-tapply(lr[,1],as.factor(GeneList),mean)
r1<-length(RX)
RR<-matrix(nrow=r1,ncol=r0)
RR[,1]=RX
for(i in 1:14) {RR[,i]<-tapply(lr[,i],as.factor(GeneList),mean)}
boxplot(RR[,1],RR[,2],RR[,3],RR[,4],RR[,5],RR[,6],RR[,7],RR[,8],RR[,9],RR[,10],RR[,11],RR[,12],RR[,13],RR[,14],ylim=c(-5,5))
MA<-normalizeWithinArrays(RG, method="loess", bc.method="normexp")
M1<-tapply(MA$M[,1],as.factor(MA$genes[,5]),mean)
n0<-length(MA$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for(i in 1:20) {MM[,i]<-tapply(MA$M[,i],as.factor(MA$genes[,5]),mean)}
boxplot(MM[,1],MM[,2],MM[,3],MM[,4],MM[,5],MM[,6],MM[,7],MM[,8],MM[,9],MM[,10],MM[,11],MM[,12],MM[,13],MM[,14],ylim=c(-5,5))
MA2<-read.csv("New_wt_WANorm_Matrix.csv",header=TRUE,sep=",")
MA3<-as.matrix(MA2)
MAB<-normalizeBetweenArrays(MA3,method="scale",targets=NULL)
boxplot(MAB[,1],MAB[,2],MAB[,3],MAB[,4],MAB[,5],MAB[,6],MAB[,7],MAB[,8],MAB[,9],MAB[,10],MAB[,11],MAB[,12],MAB[,13],MAB[,14],ylim=c(-5,5))

Again, before between array normalization could occur, the .csv file containing the within array normalized data minus the controls had to be imported into R.

All of the data generated after within array normalization and between array normalization for all strains was compiled into one Excel file. For each time point for each strain, the average log fold change was found by averaging the log fold change for all of the flasks within the time point. The standard deviation was found for each time point for each strain. The t statistic was calculated for each time point for each strain using the formula

=AVERAGE(range of cells)/(STDEV(range of cells)/SQRT(number of replicates)).

The p value was calculated for each time point for each strain using the formula

=TDIST(ABS(cell containing T statistic),degrees of freedom,2).

June 28, 2011

Using all of the between array normalized data generated yesterday along with the statistics calculated yesterday, new Excel spread sheets were created for the deletion strains to input into MATLAB. These spread sheets followed the same format as discussed during week 5. While the production rates and degradation rates remained the same, the log fold changes and the concentration sigmas had to be changed in accordance with the new normalized data. The alpha value remained 0.002. These spread sheets have yet to be created for the wildtype because we are having issues performing between array normalization for the GCAT chips for the wildtype.

MATLAB simualtions were begun for dZAP1 and dHMO1 for the 15 gene direct network generated by YEASTRACT. Unexpectedly, the simulation for dHMO1 took much longer than the simulation for dZAP1.

Katrina Sherbina 20:46, 28 June 2011 (EDT)

Week 8

July 5, 2011

To minimize chances for human error, the procedure for normalizing all of the GPRs was altered. Rather than normalizing the wildtype and the four deletion strains separately, all of the Ontario chip GPRs were normalized at once. In order to do so, first a one target file had to be created listing all of the GPRs corresponding to all of the Ontario chips from the wildtype and deletion strains. The GCAT chips were still normalized separately and the target files previously created for the top and bottom of each of the chips were used. Also, a .csv file consiting of a singla column of 1's and -1's was created along with the targets file for the Ontario chps to indicate which chips were dye swapped. The dye-swapped chips were multipled by -1. In addition, after performing within array normalization on the Ontario and GCAT chips separately, a new .csv file was created so that the Ontario and GCAT chips could be normalized between arrays in one step. Also, rather than dyeswapping all of the appropirate chips by hand in Excel, a few lines of code were added so that R could perform the dyeswapping.

The code dealing with the Ontario chips up through within array normalization is as follows.

targets<-readTargets(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RG<-read.maimages(targets, source="genepix.median", wt.fun=f)
MA<-normalizeWithinArrays(RG, method="loess", bc.method="normexp")

#Create a blank matrix which has as many rows as there are unique genes and as many columns as there are GPR files.
#Average all of the replicate spots by using the tapply function. Fill in the matrix previously created with these values.
M1<-tapply(MA$M[,1],as.factor(MA$genes[,5]),mean)
n0<-length(MA$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for(i in 1:94) {MM[,i]<-tapply(MA$M[,i],as.factor(MA$genes[,5]),mean)}

#Import the .csv file denoting which chips are dyeswapped.
ds<-read.csv("dyeswap_matrix.csv",header=TRUE,sep=",")
#Create a new matrix.
MN<-matrix(nrow=6191,ncol=94)
#Multiply the values in the matrix containing the averaged replicates log fold changes by the dyeswap list created.
for (i in 1:94) {MN[,i]<-ds[i,]*MM[,i]}

#Write the data to a table.
write.table(MN,"ont_WANorm_dyeswapped_Matrix.csv",sep=",")

Open up the exported table in Excel and delete the first two rows of values corresponding to the controls.

To within array normalize the GCAT chips,

#Read the top and bottom chips into R separately and rbind them 
#to model the data.frame after ontario chips, 14000 spots in each column.
targets<-readTargets(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RT<-read.maimages(targets, source="genepix.median", wt.fun=f)
targets<-readTargets(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RB<-read.maimages(targets, source="genepix.median", wt.fun=f)
RG<-rbind(RT,RB)

#Normalize within each array.
MA<-normalizeWithinArrays(RG,method="loess", bc.method="normexp")

#Use the tapply function to average duplicate spots so that only unique
#spots remain in the data.
M<-matrix(nrow=6403,ncol=9)
for(i in 1:9) {M[,i]<-tapply(MA$M[,i],as.factor(MA$genes[,4]),mean)}

#Write the table.
write.table(M,"GCAT_WANorm_rep.csv",sep=",")

Before between array normalization is performed, the exported table with the within array normalized GCAT chips must be put into Microsoft Access to eliminate any entries (genes) that are not within the Ontario chips. Follow the procedure outlined in the June 1, 2011 entry within this computational journal. Rather then importing all of the wildtype Ontario chip data into Microsoft Access, it is sufficient ot import just a list of the Ontario chip ID's in the form of a list in a .csv file. Once it has been edited in Access, the within array normalized data for the GCAT chips must be combined with the within array normalized data for the Ontario chips into one .csv file. Then, between array normalization is performed as follows.

#Read the .csv file containing the within array normalized data for all of the chips into R.
MA2<-read.csv("ALL_WANorm_dyeswapped_Matrix.csv",header=TRUE,sep=",")
#Force the .csv file into a matrix.
MA3<-as.matrix(MA2)
#Perform between array normalization.
MAB<-normalizeBetweenArrays(MA3,method="scale",targets=NULL)
MAC<-as.matrix(MAB)
#Write the data to a table.
write.table(MAC,"ALL_BANorm_dyeswapped_Matrix.csv",sep=",")

All of this code was performed by both myself and my coworker individually. There were some discrepencies in the data we exported after within array normalization and between array normalization. However, our values differed by at most a number to of the magnitude E-10. However, when we individually repeated the process, there was no difference between the replicate in comparison to the intial output.

We still need to regenerate all of the MA plots and all of the boxplots by making some additions to the aforementioned lines of code.

Katrina Sherbina 20:55, 5 July 2011 (EDT)

Week 9

July 11, 2011

A protocol was written up to perform within array normalization and between array normalization as previously with the exception that all of the ontario chip GPR files are treated at once rather than separately by deletion strain. This protocol is currently up on the Dahlquist Lab OpenWetWare page. However, in further looking into between array normalization, it was questioned whether or not this was the proper way to scale the data. It is not clear whether R finds the MAD of all of the chips or if it finds the MAD of individual chips. As a result, we tried to program R to scale each GPR by its own MAD rather than the MAD of all of the GPRs. The following line of code was used.

for (i in 1:15) {MS[,i]<-MX2[,i]/mad(MX2[,i],center=median(MX2[,i]),na.rm=FALSE,low=FALSE,high=FALSE)}

The scaled data was plotted against the data scaled using between array normalization. It was found that the difference between the data scaled with the above line of code and the data scaled by between array normalization is a scale factor of approximately 3.

Another problem that we began to address today was how to keep the genes as row names and the strain, flask, and timepoint as column names within R so that the row names and the column names do not have to be added after outputting the normalized data from R. First, we created a .csv file with three columns the first one which had a sample of the GPR files, the second one which had the strain, flask, and timepoint corresponding to the GPR, and the third which had a "1" or a "-1" depending on whether the chip was dyeswapped. After importing the file, we used the following code to assign column and row names through within array normalization and averaging replicates.

#Read in file with GPR names, "user-friendly" names, and dyeswap.
targets<-read.csv("GPR_and_column_headers.csv")
f<-function(x) as.numeric(x$Flags > -99)
RGO<-read.maimages(targets, source="genepix.median", wt.fun=f)

#Set from what part of the file the row names are taken.
row<-read.csv(file.choose())
row2<-as.matrix(row)
#Set from what part of the file the column names are taken.
col<-targets[,2]
#Designate where the dyeswap list comes from.
ds<-targets[,3]

#Normalize Within Arrays
MAO<-normalizeWithinArrays(RGO, method="loess", bc.method="normexp")

#Set up a for loop to average replicates.
M1<-tapply(MAO$M[,1],as.factor(MAO$genes[,5]),mean)
n0<-length(MAO$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for(i in 1:15) {MM[,i]<-tapply(MAO$M[,i],as.factor(MAO$genes[,5]),mean)}

#Set up a for loop to dyeswap the necessary GPRs.
M1<-tapply(MAO$M[,1],as.factor(MAO$genes[,5]),mean)
n0<-length(MAO$M[1,])
n1<-length(M1)
MN<-matrix(nrow=n1,ncol=n0)
MN[,1]=M1
for (i in 1:15) {MN[,i]<-ds[i]*MM[,i]}

#Turn the matrix into a data frame.
#Label the rows and columns.
MX<-as.data.frame.matrix(MN)
rownames(MX)<-row2
colnames(MX)<-col

#Delete the rows with the controls.
MX2<-MX[c(3:6191),]

Thus, in turning the output of within array normalization into a data frame, we could assign column names and row names within R. However, we cannot, at least for the moment, find a way to use similar code to retain column and row names after the data has been scaled by the MAD. We are concerned that if we try the using a data frame we still cannot be confident that we can just insert the same column names and row names after scaling has occurred because we are not sure if in doing this scaling the gene order might be rearranged.

Katrina Sherbina 20:17, 11 July 2011 (EDT)

July 12, 2011

The code in the previous entry was further refined today. The controls were taken out earlier, namely right after replicates were averaged. Then, the dyeswaps were performed using a more simplified code.

#Delete the rows with the controls.
MM2<-MM[c(3:6191),]

#Set up a for loop to dyeswap the necessary GPRs.
rows<-length(MM2[,1])
cols<-length(MM2[1,])
MN<-matrix(nrow=rows,ncol=cols)
for (i in 1:15) {MN[,i]<-ds[i]*MM2[,i]/mad(MM2[,i])}

In addition, the code to turn the resulting matrix after dyeswapping into a data frame was slightly altered. A "row3" was specified using the code

row3<-row2[c(3:6191)]

and this "row 3" was used to define the row names of data frame of dyeswapped within array normalized data.

Also, we succeeded in preventing the column headers to be shifted over to the left in outputting the normalized data from R as a .csv file. The following code was used

write.table(MX,"MAD_scaling_rep.csv",sep=",",col.names=NA,row.names=TRUE)

In addition, we began to look into how to get rid of GCAT data that does not correspond to genes on the Ontario chips in R rather than Access as we have done thus far. First, within array normalization and averaging of replicates was performed as recorded in the journal entry for July 5th with the exception that the matrix for the tapply function was not hard coded. Rather the following matrix was created first and then the tapply function was used.

m1<-tapply(MAG$M[,1],as.factor(MAG$genes[,4]),mean)
g0<-length(MAG$M[1,])
g1<-length(m1)
M<-matrix(nrow=g1,ncol=g0)
M[,1]=m1
for(i in 1:9) {M[,i]<-tapply(MAG$M[,i],as.factor(MAG$genes[,4]),mean)}

Then the resulting matrix was converted to a data frame and the row and column names were labeled with the GCAT gene names and the GPR file names, respectively.

#Read GCAT ID's (GCAT_ID.csv)
GCATID<-read.csv(file.choose(),header=TRUE,sep=",")
#Define columns.
GCATNames<-GCATID[,1]
r<-as.matrix(GCATNames)
GPR<-as.matrix(RGG$targets)
col<-GPR
#Convert to data frame
M2<-as.data.frame.matrix(M)
rownames(M2)<-r
colnames(M2)<-col

Then, the GCAT data was filtered so that the data for the genes on the GCAT chips that were not located on the Ontario chips was taken out. In order to do so, first a list of Ontario gene names was imported and this list was used to filter the data for the GCAT chips.

#Read ontario ID (Ont_Chip_GeneID.csv)
ontID<-as.data.frame(read.csv(file.choose(),header=TRUE,sep=","))

#Filter GCATs.
data<-M2
rownames(data)<-GCATID[,1]
names.to.keep<-ontID[,1]
GO<-subset(data,rownames(data) %in% names.to.keep)

The data for the GCAT genes not on the Ontario chips was deleted. However, rows corresponding to Ontario genes that the GCAT chips did not have data for were also deleted from the Therefore, instead of producing the desired output with 6189 rows (as has been produced in Access because it keeps the rows where the data for the genes unique to the GCAT chips was), the output contained only 6136 rows. As a result, we experimented with R finding the Ontario chip genes that were not in the GCAT chips. We used the following lines of code.

#Filter Ontario chips.
ont1<-subset(ontID,Ontario_ID!="Arabidopsis")
ont2<-subset(ont1,Ontario_ID!="3XSSC")
OntNames<-ont2[c(1:6403),]
ONT<-subset(OntNames,!OntNames %in% GCATNames) 

Then, we tried to bind this list of Ontario gene names not on the GCAT chips ("ONT") with the filtered within array normalized GCAT data ("GO"). We ran into errors using the rbind function to do this. Rather than adding the set of rows with the Ontario genes not on the GCAT chips, a single row was added to the end of the filtered within array normalized GCAT data with numbers. We tried to manipulate "ONT" to have as many rows and as columns as GO and to label the columns of ONT after rows and columns were added with the names of the GCAT GPR files to see if this would enable "ONT" and "GO" to be merged using rbind. However, this was also not successful.

Katrina Sherbina 20:22, 12 July 2011 (EDT)

July 13, 2011

The code for the GCAT chips listed in the previous two entries was altered.

First, rather than importing two separate text files (one for the top of the GCAT chips and one for the bottom of the GCAT chips), one .csv file was created and imported with a column indicating whether the GPR file corresponds to the top or bottom of the chip; a column listing all of the GPR files; and a column naming the GPR files with the strain, timepoint, and flask. This file was separated into the tops of the GCAT chips and the bottom of the GCAT chips. Then, the two were combined using the rbind function.

targets<-read.csv(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RT<-read.maimages(targets[1:9,2], source="genepix.median", wt.fun=f)
f<-function(x) as.numeric(x$Flags > -99)
RB<-read.maimages(targets[10:18,2], source="genepix.median", wt.fun=f)
RGG<-rbind(RT,RB)

Within array normalization was performed using the same code as before. Duplicate spots were averaged using the same lines of code was before with the exception that the line M[,1]=m1 was taken out.

Then, the GCAT gene names were read in and defined with a variable to later use to set row names. Also, the GPR files names were given a variable so as to be used later to set column names.

GCATID<-read.csv(file.choose(),header=TRUE,sep=",")
GCATNames<-as.matrix(GCATID[,1])
GPR<-as.matrix(targets[1:9,3])

The matrix M (corresponding to the data with averaged replicates) was converted to a data frame and the the GCAT gene names were set to the row names while the GCAT GPR file names were set as the column names.

M2<-as.data.frame.matrix(M)
rownames(M2)<-GCATNames
colnames(M2)<-GPR[,1]

As written in the previous journal entry, the .csv file with the Ontario gene names was read in. Then this list of names was filtered to get rid of the controls. Then, a subset of the filtered Ontario gene names was taken corresponding to the Ontario genes that the GCAT chips do not have data for.

The code to filter the within array normalized GCAT data so that it only contains Ontario gene names was slightly altered to reduce unnecessary lines of code.

rownames(M2)<-GCATNames
names.to.keep<-ontID[,1]
GO<-subset(M2,row.names(M2) %in% names.to.keep)

Then, the aforementioned subset of the filtered Ontario gene names was successfully merged with the filtered within array normalized GCAT data. First a matrix was created with as many rows as the number of Ontario gene names that the GCAT chips did not have data for and with as many columns as GCAT GPR files. This matrix was converted to a data frame and the rows were labeled with the Ontario genes not on the GCAT chips while the columns were labeled with the GPR file names. Then, this data frame was merged to the filtered within array normalized GCAT data frame using rbind.

new<-matrix(nrow=length(ONT[,1]),ncol=length(GPR[,1]))
new2<-as.data.frame.matrix(new)
rownames(new2)<-ONT$Ontario_ID
colnames(new2)<-GPR[,1]
final<-rbind(GO,new2)

The rows in the data frame were sorted so that the gene names would be listed in alphabetical order.

final.sort<-final[order(row.names(final)),]

Then, the data was scalaed by the MAD. First, a matrix had to be created with the appropriate number of columns and rows. Then this matrix had to be converted to a data frame, the rows of which were labeled with the row names from the final.sort data frame and the columns of which were labeled by the GPR file names. Only afterward was the data scaled with the MAD.

r<-length(ont2$Ontario_ID)
col<-length(GPR[,1])
MADM<-matrix(nrow=r,ncol=col)
for (i in 1:9) {MADM[,i]<-final.sort[,i]/mad(final.sort[,i],na.rm=TRUE)}
GMAD<-as.data.frame.matrix(MADM)
rownames(GMAD)<-row.names(final.sort)
colnames(GMAD)<-GPR[,1]

A new .csv file was created for all the Ontario chips listing the GPR file names in one column; the GPR files designated by strain, timepoint, and flask in another column; and whether or not the chip was dyeswapped in another column. This .csv file was imported into R. Then similar code was used as in the entries for July 11th and July 12th to within array normalized and scale normalize the data.

The following code was used to read in the .csv file aforementioned.

targets<-read.csv(file.choose())
f<-function(x) as.numeric(x$Flags > -99)
RGO<-read.maimages(targets[,1], source="genepix.median", wt.fun=f) 

Then, a .csv file with the Ontario gene names was read into to later use for row names. Also, the portion of the targets file with the strain, timepoint, and flask was defined by a different variable to later use for column names. The portion fo the targets with the dyeswap list was also designated with its own variable.

row<-read.csv(file.choose())
column<-targets[,2]
ds<-targets[,3]

Then, the data was within array normalized with the same code as before.

MAO<-normalizeWithinArrays(RGO, method="loess", bc.method="normexp")

Then the replicates were averaged using very similar lines of code with the exception that the for loop was designated for 1:94, corresponding to the 94 Ontario GPR files.

M1<-tapply(MAO$M[,1],as.factor(MAO$genes[,5]),mean)
n0<-length(MAO$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
MM[,1]=M1
for(i in 1:94) {MM[,i]<-tapply(MAO$M[,i],as.factor(MAO$genes[,5]),mean)}

The above matrix (MM) was converted to a data frame, the rows of which were labeled with the Ontario gene names and the columns of which were labeled with the strain, timepoint, and flask designations of the GPR files.

MM2<-as.data.frame.matrix(MM)
rownames(MM2)<-row$Ontario_ID
colnames(MM2)<-column

Then the controls were removed.

MM3<-subset(MM2,row.names(MM2)!="Arabidopsis")
MM4<-subset(MM3,row.names(MM3)!="3XSSC")

Then the GPR files that needed to be dyeswapped were dyeswapped. Also, simultaneously, each GPR file was scaled by its own MAD.

rows<-length(MM4[,1])
cols<-length(MM4[1,])
MN<-matrix(nrow=rows,ncol=cols)
for (i in 1:94) {MN[,i]<-ds[i]*MM4[,i]/mad(MM4[,i])}

Then this dyeswapped and scaled matrix was turned into a data frame. The rows were labeled using the same row names as in the data frame (MM4) with the controls removed. The columns were labeled with the strain, timepoint, and flask designations of the GPR files.

MX<-as.data.frame.matrix(MN)
rownames(MX)<-row.names(MM4)
colnames(MX)<-column

The resulting data frames from the within array normalization and scale normalization of the Ontario chips (MX) and the from the within array normalization and scale normalization of the GCAT chips (GMAD).

merged<-cbind(GMAD,MX)

Katrina Sherbina 20:25, 13 July 2011 (EDT)

July 14, 2011

The code to perform scale normalization for the GCAT and Ontario chips once both sets of data are within array normalized has been changed. Instead of scale normalizing each data set separately, first the GCAT and Ontario chip within array normalized data is merged. Then, each column within the merged data is scaled using the MAD of that column.

The revised code is shown below:

#Set the directory to where all the Ontario chip GPR files are located.
#Read in file with GPR names, "user-friendly" names, and dyeswap 
#for the Ontario chips.
Otargets<-read.csv(file.choose())
#Read in the GPR files by designating the column within the .csv file 
#that has the list of GPR files.
f<-function(x) as.numeric(x$Flags > -99)
RGO<-read.maimages(Otargets[,1], source="genepix.median", wt.fun=f)

#Read in the .csv file with the Ontario gene names (Ontario_Gene_ID.csv)
#genes=List of genes on the Ontario chips.
Ogenes<-read.csv(file.choose())
#row=names of the genes on the Ontario chips.
row<-Ogenes[,1]
#column=list of all of the GPR files for the Ontario chips.
#Designate which column of the Otargets file contains the names of the GPR files.
column<-Otargets[,2]
#Extract the dyeswap list (ds) from the Otargets file.
ds<-Otargets[,6]

#Normalize Within Arrays
MAO<-normalizeWithinArrays(RGO, method="loess", bc.method="normexp")

#Create a blank matrix MM.
M1<-tapply(MAO$M[,1],as.factor(MAO$genes[,5]),mean)
n0<-length(MAO$M[1,])
n1<-length(M1)
MM<-matrix(nrow=n1,ncol=n0)
#Set up a for loop to average replicates.
#The output of the loop is written to the matrix MM.
for(i in 1:94) {MM[,i]<-tapply(MAO$M[,i],as.factor(MAO$genes[,5]),mean)}

#Create a new blank matrix MN.
M2<-tapply(MAO$M[,1],as.factor(MAO$genes[,5]),mean)
n2<-length(MAO$M[1,])
n3<-length(M2)
MN<-matrix(nrow=n3,ncol=n2)
#Set up a for loop to dyeswap the chips that need to be dyeswapped.
#The output of the loop is written to the matrix MN.
for(i in 1:94) {MN[,i]<-ds[i]*MM[,i]}

#Convert matrix MN with the averaged replicates to a data frame.
MN2<-as.data.frame.matrix(MN)
#Insert the Ontario chip gene names as row names.
rownames(MN2)<-row
#Insert the GPR file names as the column names.
colnames(MN2)<-column
#Get rid of the controls (Arabidopsis and 3XSSC).
MN3<-subset(MN2,row.names(MN2)!="Arabidopsis")
MN4<-subset(MN3,row.names(MN3)!="3XSSC")

#Change directory to where all the GCAT chip GPR files are located.

#Read in .csv file with all GPR file names (designated as top or bottom 
#in a column of the .csv file),and the GPR files each labeled with 
#strain, timepoint, and flask.
Gtargets<-read.csv(file.choose())
#Read in the GPR files for the top of the GCAT chips by designating the 
#column within the .csv file that has that list of GPR files.
f<-function(x) as.numeric(x$Flags > -99)
RT<-read.maimages(Gtargets[1:9,2], source="genepix.median", wt.fun=f)
#Read in the GPR files for the bottom of the GCAT chips by designating the 
#column within the .csv file that has that list of GPR files.
f<-function(x) as.numeric(x$Flags > -99)
RB<-read.maimages(Gtargets[10:18,2], source="genepix.median", wt.fun=f)

#Bind RT and RB so that for each timepoint/flask the top and bottom of the
#chip is in one column.
RGG<-rbind(RT,RB)

#Normalize within each array.
MAG<-normalizeWithinArrays(RGG,method="loess", bc.method="normexp")

#Create a blank matrix MM.
m1<-tapply(MAG$M[,1],as.factor(MAG$genes[,4]),mean)
y0<-length(MAG$M[1,])
x1<-length(m1)
M<-matrix(nrow=x1,ncol=y0)
#Set up a for loop to average replicates.
#The output of the loop is written to the matrix MM.
for(i in 1:9) {M[,i]<-tapply(MAG$M[,i],as.factor(MAG$genes[,4]),mean)}

#Read the .csv file containing all of the GCAT gene names (GCAT_Gene_ID.csv).
GCATID<-read.csv(file.choose(),header=TRUE,sep=",")
#Set the column with the GCAT gene names as the variable GCATNames.
GCATNames<-GCATID[,1]
#Set the GPR files as the variable GPR.
GPR<-Gtargets[1:9,3]

#Convert matrix M with the averaged replicates to a data frame.
M2<-as.data.frame.matrix(M)
#Label the rows of M2 with the GCAT gene names.
rownames(M2)<-GCATNames
#Label the columns of M2 with the GCAT GPR files.
colnames(M2)<-GPR

#Read the .csv file with the Ontario gene names (Ontario_Gene_ID.csv)
ontID<-as.data.frame(read.csv(file.choose(),header=TRUE,sep=","))

#Get rid of the GCAT genes that are not on the Ontario chips.
names.to.keep<-ontID[,1]
GO<-subset(M2,row.names(M2) %in% names.to.keep)

#Take out the controls from the Ontario chips.
ont1<-subset(ontID,Ontario_ID!="Arabidopsis")
ont2<-subset(ont1,Ontario_ID!="3XSSC")
#Find the Ontario genes that the GCAT chips do not have data for.
ONT<-subset(ont2,!Ontario_ID %in% GCATNames)

#Create a new matrix for the subset of Ontario genes that are not on 
#the GCAT chips. There should be as many rows as there Ontario genes in the
#subset. There should be as many columns as there are GCAT GPR files.
subO<-matrix(nrow=length(ONT[,1]),ncol=length(GPR))

#Convert the matrix to a data frame.
subO2<-as.data.frame.matrix(subO)
#Label the rows with the names of the Ontario genes in the subset.
rownames(subO2)<-ONT$Ontario_ID
#Label the columns with the names of the GPR files.
colnames(subO2)<-GPR

#Bind the filtered GCAT data (GO) to the data frame with the Ontario genes
#not on the GCAT chips (subO2).
G<-rbind(GO,subO2)

#Sort the data frame final so that the genes are in alphabetical order.
G.sort<-G[order(row.names(G)),]

#Write to a table.
write.table(G.sort,"GCAT_and_Ontario_WAnorm.csv",sep=",",col.names=NA,row.names=TRUE)

#Bind GCAT and Ontario within array normalized data.
merged<-cbind(G.sort,MN4)

#Create a blank matrix the as many rows (r) and as many columns (col) as there are 
#in the merged data frame.
r<-length(merged[,1])
col<-length(merged[1,])
MADM<-matrix(nrow=r,ncol=col)
#Scale by the MAD.
for (i in 1:103) {MADM[,i]<-merged[,i]/mad(merged[,i],na.rm=TRUE)}
#Convert the matrix into a data frame.
merged.MAD<-as.data.frame.matrix(MADM)
#Label the rows with the row names from the merged data frame.
rownames(merged.MAD)<-row.names(merged)
#Label the columns with the column names from the merged data frame.
colnames(merged.MAD)<-colnames(merged)

At this point, the columns in the data frame are ordered so that the GCAT chips go first and then the Ontario chips follow. To remedy this, first a .csv file is created with a list of the GPR files in the correct order: strain (first wildtype and then the deletion strains), timepoint, and flask. This file is used as a reference to reorder the columns.

#Read .csv file with the list of all the GPR files in the correct order.
MasterList<-read.csv(file.choose())
#Rearrange the columns of the merged Ontario and GCAT scaled data so that they 
#are in the order that the GPR files are listed in the MasterList.
MMS<-merged.MAD[,match(MasterList[,1],colnames(merged.MAD))]

#Write final to a table.
write.table(merged.MAD,"ONT_and_GCAT_final_scaled_data.csv",sep=",",col.names=NA,row.names=TRUE)

Katrina Sherbina 20:30, 14 July 2011 (EDT)

July 15, 2011

The code listed in the entry for July 14, 2011 to get rid of extraneous lines of code and to sort the columns by strain (wildtype and then deletion strains in alphabetical order), timepoint, and flask a little bit earlier.

The line of code for the GCAT chips which was used to import the list of Ontario genes was taken out because a similar line of code was used for the Ontario chips. The variable designated in the lines of code for the Ontario chips for the list of Ontario genes was changed from Ogenes to ontID to maintain consistency throughout the entire code (both the Ontario and GCAT code).

First, the match function to rearrange the columns was first put in after the merged GCAT and Ontario data was scaled by the MAD. This was altered so that the columns are rearragned after the data is merged but before it is scaled by the MAD. The new code after merging the GCAT and Ontario within array normalized data looks as follows.

#Read .csv file with the list of all the GPR files in the correct order.
MasterList<-read.csv(file.choose())
#Rearrange the columns so that they are ordered by strain (wildtype then deletion strains
#in alphabetical order), timepoint, and then flask.
merged.sort<-merged[,match(MasterList[,1],colnames(merged))]

#Write to a table.
write.table(merged.sort,"GCAT_and_Ontario_WAnorm.csv",sep=",",col.names=NA,row.names=TRUE)

#Create a blank matrix the as many rows (r) and as many columns (col) as there are 
#in the merged data frame.
r<-length(merged.sort[,1])
col<-length(merged.sort[1,])
MADM<-matrix(nrow=r,ncol=col)
#Scale by the MAD.
for (i in 1:103) {MADM[,i]<-merged.sort[,i]/mad(merged.sort[,i],na.rm=TRUE)}
#Convert the matrix into a data frame.
merged.MAD<-as.data.frame.matrix(MADM)
#Label the rows with the row names from the merged and sorted (merged.sort) data frame.
rownames(merged.MAD)<-row.names(merged.sort)
#Label the columns with the column names from the merged and sorted (merged.sort) data frame.
colnames(merged.MAD)<-colnames(merged.sort)

#Write final to a table.
write.table(merged.MAD,"ONT_and_GCAT_final_scaled_data.csv",sep=",",col.names=NA,row.names=TRUE)

The within array and scale normalized data for the Ontario and GCAT chips was analyzed statistically within Excel. First, the average log fold change for each gene was found for each time point for each strain. Second, the standard deviation for each gene for each time point was found for each strain. Third, the t statistic for each gene was calculated for each time point for each strain using the following formula:

=AVERAGE(range of cells)/(STDEV(range of cells)/SQRT(number of replicates))

Fourth, the p value for each gene was calculated for each time point for each strain using the following formula where the degrees of freedom are number of flasks for the timepoint minus 1:

=TDIST(ABS(cell containing T statistic),degrees of freedom,2)

The Bonferroni and Benjamini Hochberg corrections were performed for the aforementioned calculated p values. The Bonferroni correction was done for each timepoint for each strain by multiplying the t-statistic derived p value by the total number of genes (i.e. 6189). To perform the Benjamini Hochberg correction, first the p values for each timepoint for each strain needed to be organized from the least p value to the greatest p value. Second, a numerical index had to be created (i.e. 1, 2, 3, ... , 6189). Then the correction was performed for each timepoint for each deletion strain by multiplying the each p value by the total number of genes (i.e. 6189) and dividing by the respective value in the numerical index.

In addition, the f statistics were found for each gene for each time point for each strain. First, the squared variance for each gene was found for each timepoint for each strain. Second, the average for each gene was found across all the flasks and times for each strain. Third, the MST for each gene was calculated for each strain using the following formula

MST=(sum(average for time i-average across all times and flasks))/T-1

where T is the total number of timepoints, in this case 5.

Fourth, the MSE for each gene was calculated for each strain using the following formula

MSE=(sum((number of flasks for time i-1)*the squared variance for timepoint i)))/N-T

where N is the total number of data points.

The F statistic was calculated for each gene for each strain by dividing the MST by the MSE. The f statistic derived p value was calculated using the following formula:

=FDIST(F statistic, T-1, N-T)

Katrina Sherbina 18:27, 15 July 2011 (EDT)

Week 10

July 18, 2011

MA plots were created before and after normalization for the wildtype and all the deletion strains using code that has been used before with some modifications. In addition, box plots were also created before normalization, after within array normalization, and after scale normalization for the wildtype and all the deletion strains using code that has been used before with some modifications. The following is the code for the MA plots and the box plots that was used today. Some of the variables come from the normalization code listed in the entries for July 14, 2011 and July 15, 2011.

#Ontario graphs 

#MA plots before normalization.Change dimensions to reflect how many graphs
#need to be fit into the window.
par(mfrow=c(4,5))

#Set a variable (GeneList) for all of the Ontario gene IDs before the control have 
#been taken out and before replicates have been averaged.
genelist<-RGO$genes$Name

#Calculate the log fold changes (M values) before normalization has occurred.
#Multiply all of the chips by the dyeswap list (ds[,i]).
#In the for loop, alter the numbers to reflect the number of GPR files in RGO.
for(i in 1:94) {lfm<-ds[i]*(log2((RGO$R-RGO$Rb)/(RGO$G-RGO$Gb)))}

#Create a blank matrix with as many columns as there are columns in lfm and as many 
#rows as there are in lfm after replicates have been averaged.
z0<-length(lfm[1,])
ZX<-tapply(lfm[,1],as.factor(genelist),mean)
z1<-length(ZX)
MO<-matrix(nrow=z1,ncol=z0)
#Calculate the log fold changes (M values) after averaging duplicate genes.
#In the for loop, alter the numbers to reflect the order of the GPR files in RGO
#for the current strain you are working on.
for(i in 1:14) {MO[,i]<-tapply(lfm[,i],as.factor(genelist),mean)}

#Calculate the intensity values (A values) before normalization has occurred.
lfa<-(1/2*log2((RGO$R-RGO$Rb)*(RGO$G-RGO$Gb)))
#Create a blank matrix with as many columns as there are columns in lfa and as many 
#rows as there are in lfa after replicates have been averaged.
z3<-length(lfa[1,])
ZQ<-tapply(lfa[,1],as.factor(genelist),mean)
z4<-length(ZQ)
AO<-matrix(nrow=z4,ncol=z3)
#Calculate the intensity values (A values) after averaging duplicate genes.
#In the for loop, alter the numbers to reflect the order of the GPR files in RGO
#for the current strain you are working on.
for(i in 1:14) {AO[,i]<-tapply(lfa[,i],as.factor(genelist),mean)}

#Plot onto graph.
#In the for loop, alter the numbers to reflect the order of the GPR files in RGO
#for the current strain you are working on.
for(i in 1:14) {plot(AO[,i],MO[,i],xlab="A",ylab="M",ylim=c(-5,5),xlim=c(0,15))}

#MA plots after within array normalization.Change dimensions to reflect how many 
#graphs need to be fit into the window.
par(mfrow=c(4,5))

#Already have the M values (MN) calculated from prior normalization.

#Create a blank matrix with as many columns as there are columns in GPR files and as many 
#rows as there are averaged duplicate genes.
v1<-tapply(MAO$A[,1],as.factor(MAO$genes[,5]),mean)
w0<-length(MAO$A[1,])
w1<-length(v1)
AAO<-matrix(nrow=w1,ncol=w0)
#Calculate the intensity values A after normalization has occurred and after duplicate genes
#have been averaged.
#In the for loop, alter the numbers to reflect the order of the GPR files in RGO
#for the current strain you are working on.
for(i in 1:14) {AAO[,i]<-tapply(MAO$A[,i],as.factor(MAO$genes[,5]),mean)}

#Plot onto graph. The M values are already given by MN. 
#In the for loop, alter the numbers to reflect the order of the GPR files in RGO
#for the current strain you are working on.
for(i in 1:14) {plot(AAO[,i],MN[,i],ylab="M",xlab="A",ylim=c(-5,5),xlim=c(0,15))}

#Before moving on to the MA plots for other strains, graph the boxplots for the strain you are
#currently working on. 

#Boxplots
#Change dimensions to reflect how many graphs need to be fit into the window.
par(mfrow=c(1,3))

#Work with one strain first fitting in the before normalization, after between array normalization, and after
#scale normalization boxplots for that strain into one window. When moving on to the next strain, make sure 
#to reinput the line of code below to designate the dimensions of the window in which the graphs will appear.
#After you have finished creating the boxplots for one strain, go back to the code for the MA plots above and
#create MA plots for all of the other strains. The code to create the genelist, the dyeswapped log fold changes 
#before normalization, and calcuating the intensity values before normalization do not have to be reinputted. 
#However, the dimensions for the window in which the graphs will appear to have to be reset. In the rest of the code, 
#the for loop has to be adjusted to designate the correct GPR files for the strain that you will be working on next. 

#Before normalization
#The number within the brackets designates a GPR file.
#wildtype
boxplot(MO[,1],MO[,2],MO[,3],MO[,4],MO[,5],MO[,6],MO[,7],MO[,8],MO[,9],MO[,10],MO[,1],MO[,12],MO[,13],MO[,14],ylim=c(-5,5))
#dCIN5
boxplot(MO[,15],MO[,16],MO[,17],MO[,18],MO[,19],MO[,20],MO[,21],MO[,22],MO[,23],MO[,24],MO[,25],MO[,26],MO[,27],MO[,28],MO[,29],MO[,30],MO[,31],MO[,32],MO[,33],MO[,34],ylim=c(-5,5))
#dGLN3
boxplot(MO[,35],MO[,36],MO[,37],MO[,38],MO[,39],MO[,40],MO[,41],MO[,42],MO[,43],MO[,44],MO[,45],MO[,46],MO[,47],MO[,48],MO[,49],MO[,50],MO[,51],MO[,52],MO[,53],MO[,54],ylim=c(-5,5))
#dHMO1
boxplot(MO[,55],MO[,56],MO[,57],MO[,58],MO[,59],MO[,60],MO[,61],MO[,62],MO[,63],MO[,64],MO[,65],MO[,66],MO[,67],MO[,68],MO[,69],MO[,70],MO[,71],MO[,72],MO[,73],MO[,74],ylim=c(-5,5))
#dZAP1
boxplot(MO[,75],MO[,76],MO[,77],MO[,78],MO[,79],MO[,80],MO[,81],MO[,82],MO[,83],MO[,84],MO[,85],MO[,86],MO[,87],MO[,88],MO[,89],MO[,90],MO[,91],MO[,92],MO[,93],MO[,94],ylim=c(-5,5))

#After within array normalization
#The number within the brackets designates a GPR file.
#wildtype
boxplot(MN[,1],MN[,2],MN[,3],MN[,4],MN[,5],MN[,6],MN[,7],MN[,8],MN[,9],MN[,10],MN[,1],MN[,12],MN[,13],MN[,14],ylim=c(-5,5))
#dCIN5
boxplot(MN[,15],MN[,16],MN[,17],MN[,18],MN[,19],MN[,20],MN[,21],MN[,22],MN[,23],MN[,24],MN[,25],MN[,26],MN[,27],MN[,28],MN[,29],MN[,30],MN[,31],MN[,32],MN[,33],MN[,34],ylim=c(-5,5))
#dGLN3
boxplot(MN[,35],MN[,36],MN[,37],MN[,38],MN[,39],MN[,40],MN[,41],MN[,42],MN[,43],MN[,44],MN[,45],MN[,46],MN[,47],MN[,48],MN[,49],MN[,50],MN[,51],MN[,52],MN[,53],MN[,54],ylim=c(-5,5))
#dHMO1
boxplot(MN[,55],MN[,56],MN[,57],MN[,58],MN[,59],MN[,60],MN[,61],MN[,62],MN[,63],MN[,64],MN[,65],MN[,66],MN[,67],MN[,68],MN[,69],MN[,70],MN[,71],MN[,72],MN[,73],MN[,74],ylim=c(-5,5))
#dZAP1
boxplot(MN[,75],MN[,76],MN[,77],MN[,78],MN[,79],MN[,80],MN[,81],MN[,82],MN[,83],MN[,84],MN[,85],MN[,86],MN[,87],MN[,88],MN[,89],MN[,90],MN[,91],MN[,92],MN[,93],MN[,94],ylim=c(-5,5))

#After scale normalization
#The number within the brackets designates a GPR file.
#wildtype
boxplot(merged.MAD[,2],merged.MAD[,3],merged.MAD[,4],merged.MAD[,7],merged.MAD[,8],merged.MAD[,9],merged.MAD[,12],merged.MAD[,13],merged.MAD[,16],merged.MAD[,17],merged.MAD[,18],
merged.MAD[,21],merged.MAD[,22],merged.MAD[,23],ylim=c(-5,5))
#dCIN5
boxplot(merged.MAD[,24],merged.MAD[,25],merged.MAD[,26],merged.MAD[,27],merged.MAD[,28],merged.MAD[,29],merged.MAD[,30],merged.MAD[,31],merged.MAD[,32],merged.MAD[,33],
merged.MAD[,34],merged.MAD[,35],merged.MAD[,36],merged.MAD[,37],merged.MAD[,38],merged.MAD[,39],merged.MAD[,40],merged.MAD[,41],merged.MAD[,42],merged.MAD[,43],ylim=c(-5,5))
#dGLN3
boxplot(merged.MAD[,44],merged.MAD[,45],merged.MAD[,46],merged.MAD[,47],merged.MAD[,48],merged.MAD[,49],merged.MAD[,50],merged.MAD[,51],merged.MAD[,52],merged.MAD[,53],
merged.MAD[,54],merged.MAD[,55],merged.MAD[,56],merged.MAD[,57],merged.MAD[,58],merged.MAD[,59],merged.MAD[,60],merged.MAD[,61],merged.MAD[,62],merged.MAD[,63],ylim=c(-5,5))
#dHMO1
boxplot(merged.MAD[,64],merged.MAD[,65],merged.MAD[,66],merged.MAD[,67],merged.MAD[,68],merged.MAD[,69],merged.MAD[,70],merged.MAD[,71],merged.MAD[,72],merged.MAD[,73],
merged.MAD[,74],merged.MAD[,75],merged.MAD[,76],merged.MAD[,77],merged.MAD[,78],merged.MAD[,79],merged.MAD[,80],merged.MAD[,81],merged.MAD[,82],merged.MAD[,83],ylim=c(-5,5))
#dZAP1
boxplot(merged.MAD[,84],merged.MAD[,85],merged.MAD[,86],merged.MAD[,87],merged.MAD[,88],merged.MAD[,89],merged.MAD[,90],merged.MAD[,91],merged.MAD[,92],merged.MAD[,93],
merged.MAD[,94],merged.MAD[,95],merged.MAD[,96],merged.MAD[,97],merged.MAD[,98],merged.MAD[,99],merged.MAD[,100],merged.MAD[,101],merged.MAD[,102],merged.MAD[,103],ylim=c(-5,5))

#GCAT graphs

#MA plots before normalization. Change dimensions to reflect how many graphs
#need to be fit into the window.
par(mfrow=c(3,3))

#Set a variable (GeneList) for all of the GCAT gene IDs before the control have 
#been taken out and before replicates have been averaged.
GeneList<-RGG$genes$ID

#Calculate the log fold changes (M values) before normalization has occurred.
lr<-(log2((RGG$R-RGG$Rb)/(RGG$G-RGG$Gb)))

#Create a blank matrix with as many columns as there are columns in lr and as many 
#rows as there are in lr after replicates have been averaged.
r0<-length(lr[1,])
RX<-tapply(lr[,1],as.factor(GeneList),mean)
r1<-length(RX)
MG<-matrix(nrow=r1,ncol=r0)
#Calculate the log fold changes (M values) after averaging duplicate genes.
#In the for loop, alter the numbers to reflect the number of GPR files.
for(i in 1:9) {MG[,i]<-tapply(lr[,i],as.factor(GeneList),mean)}

#Calculate the intensity values (A values) before normalization has occurred.
la<-(1/2*log2((RGG$R-RGG$Rb)*(RGG$G-RGG$Gb)))
#Create a blank matrix with as many columns as there are columns in la and as many 
#rows as there are in la after replicates have been averaged.
r3<-length(la[1,])
RQ<-tapply(la[,1],as.factor(GeneList),mean)
r4<-length(RQ)
AG<-matrix(nrow=r4,ncol=r3)
#Calculate the intensity values (A values) after averaging duplicate genes.
#In the for loop, alter the numbers to reflect the number of GPR files.
for(i in 1:9) {AG[,i]<-tapply(la[,i],as.factor(GeneList),mean)}

#Plot onto graph. In the for loop, alter the numbers to reflect the number of GPR files.
for(i in 1:9) {plot(AG[,i],MG[,i],xlab="A",ylab="M",ylim=c(-5,5),xlim=c(0,15))}

#MA plots after within array normalization. Change dimensions to reflect how many 
#graphs need to be fit into the window.
par(mfrow=c(3,3))

#Already have the M values (M) calculated from prior normalization.

#Create a blank matrix with as many columns as there are columns in GPR files and as many 
#rows as there are averaged duplicate genes.
X1<-tapply(MAG$A[,1],as.factor(MAG$genes[,4]),mean)
y0<-length(MAG$A[1,])
y1<-length(X1)
AAG<-matrix(nrow=y1,ncol=y0)
#Calculate the intensity values A after normalization has occurred and after duplicate genes
#have been averaged.. In the for loop, alter the numbers to reflect the number of GPR files.
for(i in 1:9) {AAG[,i]<-tapply(MAG$A[,i],as.factor(MAG$genes[,4]),mean)}

#Plot onto graph. The M values are already given by M[,i]. In the for loop, 
#alter the numbers to reflect the number of GPR files.
for(i in 1:9) {plot(AAG[,i],M[,i],ylab="M",xlab="A",ylim=c(-5,5),xlim=c(0,15))}

#Boxplots
#Change dimensions to reflect how many graphs need to be fit into the window.
par(mfrow=c(1,3))

#Before normalization.
#The number within the brackets designates a GPR file.
boxplot(MG[,1],MG[,2],MG[,3],MG[,4],MG[,5],MG[,6],MG[,7],MG[,8],MG[,9],ylim=c(-5,5))

#After within array normalization.
#The number within the brackets designates a GPR file.
boxplot(M[,1],M[,2],M[,3],M[,4],M[,5],M[,6],M[,7],M[,8],M[,9],ylim=c(-5,5))

#After scale normalization.
#The number within the brackets designates a GPR file.
boxplot(merged.MAD[,1],merged.MAD[,5],merged.MAD[,6],merged.MAD[,10],merged.MAD[,11],merged.MAD[,14],merged.MAD[,15],merged.MAD[,19],merged.MAD[,20],ylim=c(-5,5))

The Bonferroni and Benjamini Hochberg corrections were also done for the p values generated by the f statistics using the same calculations as described in the journal entry for July 15, 2011.

Katrina Sherbina 20:34, 18 July 2011 (EDT)

July 21, 2011

Input Excel files for MATLAB simulations were created using the most recent within array and scale normalized data generated in R. Input files for the wildtype, dCIN5, dHMO1, dGLN3, and dZAP1 data were created for a 21 gene network based on a Lee et al. (2002) paper, a 21 gene network the connections in which were generated by YEASTRACT using direct evidence, and a 24 gene network (the 21 gene network plus GLN3, HMO1, and ZAP1) the connections in which were generated by YEASTRACT using direct evidence. For each of these input files, the log 2 concentrations for each gene in the network were taken from the average log fold changes for each time point for each strain calculated from the normalized data. The concentration sigmas for each gene in the network were taken from the standard deviations for each timepoint for each strain calculated from the normalized data. The production rates and degradation rates were copied from a previous template input MATLAB file. The alpha value was set to 0.002 as was used in MATLAB simulation previously.

Katrina Sherbina 16:14, 21 July 2011 (EDT)