My Computational Journal Summer 2012 - Revision history

Katrina Sherbina: /* July 12, 2012 */ Corrected a grammatical error.

2012-10-08T22:33:22Z

July 12, 2012: Corrected a grammatical error.

←Older revision		Revision as of 22:33, 8 October 2012
Line 1,017:		Line 1,017:
	Before the 20,000 iteration reverse stochastic simulation could finish, the entire simulation needed to be restarted setting the fix_P to 1 (it was 0) so that the prodution rates are not estimated.		Before the 20,000 iteration reverse stochastic simulation could finish, the entire simulation needed to be restarted setting the fix_P to 1 (it was 0) so that the prodution rates are not estimated.

-	While the simulation was running during the work day (it has not finished yet), we performed within array normalization and between array normalization for S. paradoxus microarray data gathered during a cold shock experiment. In looking over the normalized data (after both normalization procedures), we noticed some very drastic log fold change values (ex. around -5 and downwards to -15). As this caused us some concern, we made boxplots of the data before normalization, after within array normalization, and between array normalization. The log fold change values after the between array normalization for t120 flask 1 (there was only one flask for t120 for analysis) were mostly negative. There is a very long tail for the box plot below 0 reaching all the way to -15. We then decided to also produce MA plots for before normalization and after within array normalization to corroborate what we saw with t120 flask 1. For t120 flask 1, there was a prominent tail with a negative slope in the before and after within array normalization MA plots. The log fold change values in this tail had increasingly negative log fold changes and increasing intensity values. What was reassuring to see was the robustness of loess normalization as this procedure did not try to somehow correct this negative tail by improperly shifting the values in the tail. For the other chips, the normalization seemed to work well getting rid of hooks. After looking at the graphs, we went back to the pictures of the microarray slides with Dr. Dahlquist. We ~~say~~ that for the t120 flask 1 microarray, a majority of the spots were very green and there were some blocks that fluoresced very little and if so the spots were green. Tomorrow, we will try to extract from the MA plots which genes consist of the negative slope tail by selecting the cutoffs for the M values and A values where the tail seems to begin and finding the genes that have M values and A values within the cutoff.	+	While the simulation was running during the work day (it has not finished yet), we performed within array normalization and between array normalization for S. paradoxus microarray data gathered during a cold shock experiment. In looking over the normalized data (after both normalization procedures), we noticed some very drastic log fold change values (ex. around -5 and downwards to -15). As this caused us some concern, we made boxplots of the data before normalization, after within array normalization, and between array normalization. The log fold change values after the between array normalization for t120 flask 1 (there was only one flask for t120 for analysis) were mostly negative. There is a very long tail for the box plot below 0 reaching all the way to -15. We then decided to also produce MA plots for before normalization and after within array normalization to corroborate what we saw with t120 flask 1. For t120 flask 1, there was a prominent tail with a negative slope in the before and after within array normalization MA plots. The log fold change values in this tail had increasingly negative log fold changes and increasing intensity values. What was reassuring to see was the robustness of loess normalization as this procedure did not try to somehow correct this negative tail by improperly shifting the values in the tail. For the other chips, the normalization seemed to work well getting rid of hooks. After looking at the graphs, we went back to the pictures of the microarray slides with Dr. Dahlquist. We saw that for the t120 flask 1 microarray, a majority of the spots were very green and there were some blocks that fluoresced very little and if so the spots were green. Tomorrow, we will try to extract from the MA plots which genes consist of the negative slope tail by selecting the cutoffs for the M values and A values where the tail seems to begin and finding the genes that have M values and A values within the cutoff.

	We also performed an ANOVA test with the normalized data to find which genes had significant differential expression during at least one time point. Over 1000 genes had significant differential expression according to the ANOVA (p value<0.05) but there were only 4 significantly differentially expressed genes after performing the Benjamini & Hochberg p value correction (p value<0.05). Looking at the graphs of the data and the full model from the ANOVA test for those 4 genes with B&H p values < 0.05, we saw that the t120 flask 1 data really skewed the model. The results of this ANOVA test are prone to suspicion however considering the normalized log fold changes for t120 flask 1 and the fact that there were at most 2 replicates for each timepoint.		We also performed an ANOVA test with the normalized data to find which genes had significant differential expression during at least one time point. Over 1000 genes had significant differential expression according to the ANOVA (p value<0.05) but there were only 4 significantly differentially expressed genes after performing the Benjamini & Hochberg p value correction (p value<0.05). Looking at the graphs of the data and the full model from the ANOVA test for those 4 genes with B&H p values < 0.05, we saw that the t120 flask 1 data really skewed the model. The results of this ANOVA test are prone to suspicion however considering the normalized log fold changes for t120 flask 1 and the fact that there were at most 2 replicates for each timepoint.

Katrina Sherbina: /* June 12, 2012 */ Corrected some logical and grammatical errors.

2012-10-08T22:27:24Z

June 12, 2012: Corrected some logical and grammatical errors.

←Older revision		Revision as of 22:27, 8 October 2012
Line 790:		Line 790:
	===June 12, 2012===		===June 12, 2012===

-	I focused today mainly on analyzing the most recent microarray data from the first mRNA decay experiment that Dr. Dahlquist performed. There were three sets of data differing in how the PMT was set: manually for each ~~arrray~~, automatically by the program for each array, or one PMT for each array based on an average. We decided ~~not~~ to hold off on the within array normalization and between array normalization (the protocol that was created last summer) because these kind of normalization procedures are based on the assumption that a majority of the genes will not have a significant change in expression. However, we do expect to see a change in the mRNA decay experiment.	+	I focused today mainly on analyzing the most recent microarray data from the first mRNA decay experiment that Dr. Dahlquist performed. There were three sets of data differing in how the PMT was set: manually for each array, automatically by the program for each array, or one PMT for each array based on an average. We decided to hold off on the within array normalization and between array normalization (the protocol that was created last summer) because these kind of normalization procedures are based on the assumption that a majority of the genes will not have a significant change in expression. However, we do expect to see a change in the mRNA decay experiment.

-	I began by first making MA plots for the non-normalized data (one plot for each GPR file). I used a script similar to the one that was used last summer where the M values and A values were calculated by hand from the RG file that was created when reading in the targets. I did have to change the M and A matrices (after the tapply function was performed) to data frames so that I could ~~deleted~~ the controls Arabidopsis and 3XSSC. To deleted the controls, I first specified ~~the~~ the rownames (gene names) within R rather an importing a sv file of the gene names as was done previously:	+	I began by first making MA plots for the non-normalized data (one plot for each GPR file). I used a script similar to the one that was used last summer where the M values and A values were calculated by hand from the RG file that was created when reading in the targets. I did have to change the M and A matrices (after the tapply function was performed) to data frames so that I could delete the controls Arabidopsis and 3XSSC. To deleted the controls, I first specified the rownames (gene names) within R rather an importing a csv file of the gene names as was done previously:

	RGt<-tapply(RG$R[,1],as.factor(GeneList),mean)		RGt<-tapply(RG$R[,1],as.factor(GeneList),mean)

Katrina Sherbina: /* Week 10 */ Added entry for July 19, 2012.

2012-07-20T00:07:00Z

Week 10: Added entry for July 19, 2012.

←Older revision		Revision as of 00:07, 20 July 2012
Line 1,038:		Line 1,038:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:52, 18 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:52, 18 July 2012 (EDT)
		+
		+	===July 19, 2012===
		+
		+	Today, we finished analyzing the growth curve data from 2009 to this summer to determine if to determine whether or not a difference in growth rate in 13 degrees between the wildtype and another strain (deletion strain or S. paradoxus) is significant when compared to the growth rate at 30 degrees. The following strains were compared to the wildtype: dCIN5, dGLN3, dRPN4, dZAP1, dHMO1, dINO2, dOPI1, dRIM101, dGAT1, dGZF3, S. paradoxus, dURE2, dGAT1, dGZF3, dHAP3, and dHAP5. There were two sets of data (from two different growth curve experiments) for dGAT1 and dGZF3. To test the hypothesis, a general linear model was constructed for all 16 comparisons in SPSS to obtain the hypothesized model residuals and the full model residuals. In the case of the hypothesized residuals, only 7 of the 8 parameters were used in the calculation (all except the slope for the deletion strain or S. paradoxus, depending on the comparison, for 30 degrees Celsius). For the full model residuals, all 8 parameters were used in the calculation. For all 16 comparisons, very large F statistics were calculated from the hypothesized residuals and full model residuals (all above 100 and some into the 1000s). As was expected from the F statistics, very small p values were obtained for all 16 comparisons. The only p value that SPSS could compute was that for the comparison between wild type and dZAP1 (on the order of 10^-14). For the other comparisons, Graph Pad was used to try to obtain a p value. The best that it could do was calculate the p value to be less than <0.0001 for all of the comparisons. Therefore, the difference in growth rate between the 13 degrees and 30 degrees for the wildtype and the other strains is significant.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:06, 19 July 2012 (EDT)

Katrina Sherbina: /* Week 9 */ Added entry for July 18, 2012.

2012-07-19T00:52:10Z

Week 9: Added entry for July 18, 2012.

←Older revision		Revision as of 00:52, 19 July 2012
Line 1,028:		Line 1,028:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:58, 13 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:58, 13 July 2012 (EDT)
		+
		+	==Week 10==
		+
		+	===July 18, 2012===
		+
		+	Sometime last week, Dr. Dahlquist sent revised growth curve data for several deletion strains and S. paradoxus. Since we spent Monday working on posters for the Mathematical Biology Conference and the presentation for the SURP Symposium (which was yesterday), I returned to the statistical analysis of the growth curve data. I already had matrices set up for GLM testing for wildtype versus URE2, GAT1, and GZF3. However, I redid these matrices with the revised growth curve data that Dr. Dahlquist sent. In addition, I worked on setting up the matrices for the other deletion strains as well as S. paradoxus. Each matrix has a dependent variable column and eight other columns corresponding to eight parameters: slope for wildtype during cold shock, slope for wildtype during optimal temperature (which is more like a delta representing the change in the growth rate between cold and optimal temperatures), intercept for wildtype during cold shock, intercept for wildtype during optimal temperature, and those corresponding slopes and intercepts for the deletion strain. Each matrix is set up to compare the growth rate of one deletion strain at a time to the growth rate of the wildtype. Tomorrow, I will be finishing the matrices for the remaining deletion strains so that i can perform the GLM testing before research ends for the summer.
		+
		+	I also spent a part of today continuing to work on the determinstic model poster for the Mathematical Biology conference. Namely, I redid the "AND" gate, "OR" gate, and sigmoid model representations in 3D. From what Dr. Fitzpatrick said today, the graph we called the sigmoid model actually was a representation of the production of a target gene in the case of an activator and repressor with an "OR" gate type. The original representation of the "OR" gate type actually corresponded to the production of a target gene in the case of two activators with an "OR" gate. I kept these two graphs and substituted the "AND" gate graph with two representations of the production of a target gene in the case of two activators with the sigmoid model and in the case of an activator and a repressor with the sigmoid model.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:52, 18 July 2012 (EDT)

Katrina Sherbina: /* Week 9 */ Added entry for July 13, 2012.

2012-07-14T00:58:27Z

Week 9: Added entry for July 13, 2012.

←Older revision		Revision as of 00:58, 14 July 2012
Line 1,022:		Line 1,022:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 21:07, 12 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 21:07, 12 July 2012 (EDT)
		+
		+	===July 13, 2012===
		+
		+	Today I spent the day working on the PowerPoint presentation for the SURP Symposium and the poster about the deterministic model for the Mathematical Biology Conference. Intially, the the gene expression graphs in the PowerPoint came from running the MATLAB simulation to solve the differential equation and optimize parameters with just wildtype data and with just dCIN5 data. Today, I outputted expression graphs for the 21 genes based on the widltype model and dCIN5 model produced from the simulation using data from all strains. In the intial graphs, there seemed to be a differences in the expression profile of SWI6 (a no input gene) between the wildtype model and the dCIN5 model. However, none of the input genes had a difference in their expression profiles in the models produced for the wildtype and dCIN5 strains from the simulation run with data from all strains. This is what we expected because by definition the expression of the genes with no inputs in the model should not change if a gene encoding a transcription factor in that same model is deleted. These graphs provide evidence that using the data from all strains produces a more robust model than using data from just one strain.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:58, 13 July 2012 (EDT)

Katrina Sherbina: /* Week 9 */ Added entry for July 12, 2012.

2012-07-13T01:07:36Z

Week 9: Added entry for July 12, 2012.

←Older revision		Revision as of 01:07, 13 July 2012
Line 1,012:		Line 1,012:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:44, 11 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:44, 11 July 2012 (EDT)
		+
		+	===July 12, 2012===
		+
		+	Before the 20,000 iteration reverse stochastic simulation could finish, the entire simulation needed to be restarted setting the fix_P to 1 (it was 0) so that the prodution rates are not estimated.
		+
		+	While the simulation was running during the work day (it has not finished yet), we performed within array normalization and between array normalization for S. paradoxus microarray data gathered during a cold shock experiment. In looking over the normalized data (after both normalization procedures), we noticed some very drastic log fold change values (ex. around -5 and downwards to -15). As this caused us some concern, we made boxplots of the data before normalization, after within array normalization, and between array normalization. The log fold change values after the between array normalization for t120 flask 1 (there was only one flask for t120 for analysis) were mostly negative. There is a very long tail for the box plot below 0 reaching all the way to -15. We then decided to also produce MA plots for before normalization and after within array normalization to corroborate what we saw with t120 flask 1. For t120 flask 1, there was a prominent tail with a negative slope in the before and after within array normalization MA plots. The log fold change values in this tail had increasingly negative log fold changes and increasing intensity values. What was reassuring to see was the robustness of loess normalization as this procedure did not try to somehow correct this negative tail by improperly shifting the values in the tail. For the other chips, the normalization seemed to work well getting rid of hooks. After looking at the graphs, we went back to the pictures of the microarray slides with Dr. Dahlquist. We say that for the t120 flask 1 microarray, a majority of the spots were very green and there were some blocks that fluoresced very little and if so the spots were green. Tomorrow, we will try to extract from the MA plots which genes consist of the negative slope tail by selecting the cutoffs for the M values and A values where the tail seems to begin and finding the genes that have M values and A values within the cutoff.
		+
		+	We also performed an ANOVA test with the normalized data to find which genes had significant differential expression during at least one time point. Over 1000 genes had significant differential expression according to the ANOVA (p value<0.05) but there were only 4 significantly differentially expressed genes after performing the Benjamini & Hochberg p value correction (p value<0.05). Looking at the graphs of the data and the full model from the ANOVA test for those 4 genes with B&H p values < 0.05, we saw that the t120 flask 1 data really skewed the model. The results of this ANOVA test are prone to suspicion however considering the normalized log fold changes for t120 flask 1 and the fact that there were at most 2 replicates for each timepoint.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 21:07, 12 July 2012 (EDT)

Katrina Sherbina: /* Week 9 */ Added entry for July 11, 2012.

2012-07-12T00:44:45Z

Week 9: Added entry for July 11, 2012.

←Older revision		Revision as of 00:44, 12 July 2012
Line 1,003:		Line 1,003:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:25, 10 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:25, 10 July 2012 (EDT)
		+
		+	===July 11, 2012===
		+
		+	The reverse simulation for the wildtype strain using the stochastic model finished some time this morning. When looking at the expression graphs generated by the simulation, it was evident that the model did not fit all of the genes very well. For some genes, the model curve was above that of the data. In future graphs that we generate, it may be useful to include the upper and lower error bounds as they were once shown in the expression graphs from the original deterministic model simulation to make a better judgement of the goodness of fit of the model. To try to find a model that fits the data better, we repeated the simulation with different values of a (calculated in the fstochmin function) using n = 1000 for the stochastic optimization. After outputting the expresson graph for each gene in the network for each value of a on the same set of axes, we decided to use a = 0.05. Nick started a simulation with a = 0.05 and n=20000 which will hopefully be done by tomorrow.
		+
		+	Another comparison we made today the percent difference between the input network weights and the network weights optimized by the stochastic simulation that finished this morning. It seems that for a majority of the genes the optimized network weights were not much different from the input weights. However, for some genes (particularly YAP6) the weights changed sign.
		+	Along the lines of improving the model, Dr. Fitzpatrick also suggested changing the inital guesses we use for the network weights from the optimized network weights from the determinstic model to those weights multiplied by the production rate. However, this is a simulation that we will probably consider tomorrow. In waiting for simulations to finish, we worked on the presentation for the SURP symposium.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:44, 11 July 2012 (EDT)

Katrina Sherbina: /* Week 9 */ Added entry for July 10, 2012.

2012-07-11T00:25:31Z

Week 9: Added entry for July 10, 2012.

←Older revision		Revision as of 00:25, 11 July 2012
Line 991:		Line 991:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:07, 9 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:07, 9 July 2012 (EDT)
		+
		+	===July 10, 2012===
		+
		+	The line that computes the model in the stochastic sim fun is correct. With Dr. Fitzpatrick, we did a series of debugging procedures on the forward simulation. The forward simulation calls only the paramters, stoch LSE, stochastic sim fn and the network rates function. The problem that stumped us for some time was that Nlist (which is computed in the stochastic sim fu) was not being zeroed out after one iteration of the for loop from 1 to the number of Monte Carlo's (nmc). This did not seem to be a problem when the original simple stochastic simulation and network rates function were run. The problem arose in taking out the loop from 1 to the nmc and making it a separate function. In order to zero out Nlist, we needed to add a line of code before the for loop (x0 = N0;) and replace the assignment Nlist(1,:) (in the for loop) so that it is x0' not N0'.
		+
		+	After running through the forward simulation for just the wildtype strain, it we saw that Nlist was being zeroed out because the final model outputted had 1's in the first row, as we expected. The entire forward simulation on my computer took over an hour to finish. In looking over the results, the log2 optimized concernations of the three input genes did not change much over the time course of the simulation (0 to 60 in increments of 5 minutes). The log2 optimized fold changes for SWI6 were for the most part negative and the same magnitude throughout a majority of the simulation. Those for FHL1 and SKO1 were positive and of the same magnitude through the simulation. In looking at the expression graphs that were outputted, specfically those for FHL1, SKO1, and SWI6, Dr. Fitzpatrick suggested that we make sure that there is some basal level of transcription accounted in the code. In order to do so, the network rates function was altered so that jump_rates(ii,1) and jump_rates(ii,2) were calculated twice (for when upj and downj are calculated and for when pro and deg are calculated). In addition, Dr. Fitzpatrick suggested increasing the number of iterations n that occur of the for loop in the fstochmin function to see if we can get a better fit of the model to the data. Thus, we increased n from 1000 to 5000.
		+
		+	We began the reverse simulation for just the wildtype strain on two computers. There is a marked difference in the speed of the calculations between the two computers as determined by observing the counter change as the fstochmin function goes through multiple iterations calculating theta1. By now, Nick's computer is over halfway done with the fstochmin iterations and my computer is still below 1000. When we begun running the simulation for all the strain data, we will most likely have to work on the computers in the math department.
		+
		+	Toward the end of the day, I also began setting up the matrices to use GLM tests to determine whether or not a difference in growth rate in 13 degrees between the wildtype and a deletion strain is significant when comparing to the growth rate at 30 degrees. I have thus far completed the matrices for wt versus GAT1, wt versus GZF3, and wt versus URE2. Each matrix has 8 parameters (4 slopes and 4 intercepts). Before I move on to the other matrix, I will check with Dr. Fitzpatrick to make sure that I am understanding the hypothesis for the GLM correctly.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:25, 10 July 2012 (EDT)

Katrina Sherbina: /* Week 8 */ Added entry July 9, 2012.

2012-07-10T00:07:03Z

Week 8: Added entry July 9, 2012.

←Older revision		Revision as of 00:07, 10 July 2012
Line 977:		Line 977:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:41, 6 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:41, 6 July 2012 (EDT)
		+
		+	==Week 9==
		+
		+	===July 9, 2012===
		+
		+	The simulation that I started running last Friday ended with an error message yesterday night that the penalty_out term was not defined. After putting the penalty_out in the global statement, we reran the simulation with just the wildtype strain. The simulation began before lunch and when we came back after lunch there was yet another error message. Our primary concern today was why it was taking so long for the stochastic code to finish running. To debug, we changed iestimate to 0 so that we only run the forward simulation. We also changed the time that was being used in the stochastic simulation function from the 101 times from 0 to 60 min to simtime (0 to 60 min by increments of 5 minutes). The forward simulation finished on my computer within a matter of minutes. However, the model that was outputted was not correct. This was particularly evident in the expression graphs outputted as few if any actually began at 0 as they should when taking the log2 of the model in which the first row should be 1's corresponding to time 0 minutes. At this point, there is an error occuring in the stochastic sim fun when the model is being computed. When running the code I compiled, the first model outputted through one iteration of the for loop in the simulation does have 1's in the first row for the model. However, by the time the last model is outputted, the first row has changed to values other than 1's. When running the code Nick compiled (which seems the same with the exception of a few of the variables being changed), the first row in the first model outputted consists of all 0.1's. In the last model that is outputted, the first row has changed from the first model.
		+
		+	I have some questions about the line of code in the stochastic sim fun that calculates the model:
		+
		+	model = model + interp1(tlist,Nlist,time)/nmc;
		+
		+	With the way this line of code is set up, it seems like the rows that have already been calculated for the previous model change with the new model that is computed through the next iteration of the for loop. However, in taking out the summation of the interpolation and the previous model, the final model is still not what we are looking for it to be.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:07, 9 July 2012 (EDT)

Katrina Sherbina: /* Week 8 */ Added entry for July 6, 2012.

2012-07-07T00:41:12Z

Week 8: Added entry for July 6, 2012.

←Older revision		Revision as of 00:41, 7 July 2012
Line 968:		Line 968:

	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:45, 5 July 2012 (EDT)		[[User:Katrina Sherbina\|Katrina Sherbina]] 20:45, 5 July 2012 (EDT)
		+
		+
		+	===July 6, 2012===
		+
		+	The dimensions of several variables in the fstochmin function were changed to make sure that what needed to be an array was an array and that what needed to be a scalar was a scalar. In addition, a line of code was added at the end of the for loop that just states theta1 so that we can see the progress of the optimization. In addition, some small changes were added to the way L and penalty_out are calculated at the end of the general_stoch_ls function.
		+
		+	We began the entire stochastic simulation from the beginning before lunch. Now at the end of the work day, the simulation has not yet stopped. Interestingly, a similar simulation that is still running on Nick's computer is showing the same kind of subplot that was generated during the deterministic simulation to see how the model is fitting the day with each step of the optimization. Howvever, even though I have the same line of code to specify this subplot, no figure has appeared over the course of the simulation. Instead, I can only see how the theta1 values are changing. At the moment, there all seem to be within -1.5 and 1.5. I will keep the simulation running through out the weekend and check on it when I come back to work Monday.
		+
		+	[[User:Katrina Sherbina\|Katrina Sherbina]] 20:41, 6 July 2012 (EDT)