My Computational Journal Summer 2012
Week 1
May 14, 2012
Today code for a simple ANOVA and a comparison of two strains was tested to see if the results produced by our adviser could be reproduced. The output without modifying the code for the comparison of two strains using wildtype and dGLN3 for the comparison could be reproduced. A few difficulties have been faced when trying to modify the simple ANOVA code. It was attempted to match the results produced for dGLN3 using the two strain comparison code with the results of the simple ANOVA code with the input specified as the log fold concentrations for the dGLN3 data. After being initially unsuccessful, the code was revisited and looked at more closely. It was found that two matrices in the matrix division had an unequal number of rows (matrix X and matrix Y). The indices (which specify columns of the input file) were modified after which matrix X had the same number of rows as matrix Y. The following code was used:
ind15 = ind15-indx(1)+1; ind30 = ind30-indx(1)+1; ind60 = ind60-indx(1)+1; ind90 = ind90-indx(1)+1; ind120 = ind120-indx(1)+1;
However, the results from the simple ANOVA for dGLN3 still did not match the results for dGLN3 from the comparison of two strains code.
Katrina Sherbina 21:31, 14 May 2012 (EDT)
May 15, 2012
A separate set of indices were made to designate the rows of the X matrices for the reduced and the full model. This separate set of indices (one for each timpeoint) removed rows of zeros that were previously added when substituting the select 0's for 1's to create the two X matrices. As a result, it was possible to call upon the original indices to extract the log fold changes corresponding to a specific deletion to create the Y matrix. After these corrections, the two strain comparisons were performed between wildtype and each one of the deletion strains.
Further modifications were made to the two strain comparison code to be able to make a comparison of all of the strains simulatenously. To do so, indices were added to designate the columns in the input corresponding to each of the deletion strains. A separate set of indices (one index per timepoint) for each strain as in the case of the two strain comparison. The parameters p were increased to 25 (five timepoints for each deletion strain) and the constraints q were increased to 20 to yield an X and Xh matrices with the appropriate numbers of columns. Changes were made to the out_data(ii,[number]) lines to account for the increase in the number of strains being compared. In the output, the data for each of the timepoints for each of the deletion strains matched the corresponding data in the two strain comparisons.
The next task is figure out how to handle missing data. Since both GCAT and Ontario chips were used, log fold change concentrations are missing for some genes for some timepoints for the wildtype (these cells show up as NaN). One possibility is to modify the X matrices so that they take into account the NaN's for each gene. On this note, a for loop was begun to try to exclude timepoints for a gene for which there is an NaN.
Katrina Sherbina 20:38, 15 May 2012 (EDT)
May 16, 2012
NaN's were removed from the Y matrix for each gene by using the command YY = Y(~isnan(Y(:,1)),1). This truncated the Y matrix from 43 rows to 34 rows. However, in doing so, it was not possible to solve for beta (the solution to the linear system XB=Y). Another method to deal with the NaN's was to replace them with a very large number using the command Y(isnan(Y(:,1)),1)=1e6 in order to be able to recognize the timepoints and flasks for which there were NaN's in the output. However, this messed up the B&H corrections.
The next step will build on trucating the Y matrix by getting rid of the rows with the NaN's altogether. To keep this step, the X and Xh matrix pairs must be formulated so that they are unique to each gene. A possibilty is to use a for loop to define each of the indices for each gene by the columns that do not contain an NaN for that gene. The Y matrix calculations would have to be within this loop.
Katrina Sherbina 21:33, 16 May 2012 (EDT)
May 17, 2012
A for loop was experimented with to generate different X and Xh matrix pairs for each gene to take into account those genes that have an NaN for a a flask for a timepoint in the wildtype. This loop contained a majority of the same coding that was used to set up the indices in the original two strain comparison code. The different lines of code are
for ii=1:n I = find(~isnan(a(ii,1:23))); ind15 = I(find(I>=1&I<=4)); ind30 = I(find(I>=5&I<=9)); ind60 = I(find(I>=10&I<=13)); ind90 = I(find(I>=14&I<=18)); ind120 = I(find(I>=19&I<=23));
indX = find(indx); end
These new lines of code retain the original column numbers of matrix a (the matrix containing the log fold change data) in the indices. The NaN's were then also removed from the corresponding Y vectors. At first, the loop to compute the Y matrix was put into the loop to compute the X and Xh pairs. However, the code ran for too long at that point. Then the Y loop was separated from the loop to compute the X and Xh matrices. An output was generated. However, for any of the genes that had at least one NaN, the entire row for that gene was blank. In dissecting this new code several problems were found. First, the X matrices for genes that had NaN's were not being computed correctly. Since the indices still had the original column numbers of matrix a, the ones were not put into the X matrix in the correct places. Briefly, cell arrays were looked at as a solution instead of numerical arrays because there can be blank cells in a cell array but numerical arrays do not contain blanks. However, this potential solution proved to complicate the calculations that need to be made.
May 18, 2012
(This entry is being written on May 20th regarding the work that was done on May 18th.)
The first issue that was tackled was the improper formation of the X and Xh matrices. For any of the genes that had at least one NaN, the 1's were put in the wrong positions in the X and Xh matrices. This is because the indices that were used to replace the 0's in the X and Xh matrices with 1 did not start at 1. This was necessary to later exlude the cells that had NaN's in the Y matrix. Rather than getting rid of the method to select the cells that do not have NaN's, the indices that were used to replace the 0's with 1's in the X and Xh matrices were renumbered if indx (the matrix of the indices for each time point) was less than 23 (the maximum number of columns correponding to the wildtype if none of the cells have an NaN for that gene). The following added bit of code was put into the for loop for the X and Xh matrices after the indx, indX, i1, and i2 lines of code:
if indx<23
ind15 = find(ind15);
ind30 = ind30-ind15(end)+1;
ind60 = ind60-ind30(1)-1;
ind90 = ind90-ind60(1);
ind120 = ind120-ind90(1);
end
However, this may have fixed the X and Xh matrices. However, in the final output, the rows corresponding to the genes with NaN's were still completely blank.
At this point, the wrong direction was taken to find a solution to this. Some experimenting was done to produce X and Xh matrices all of which had 43 rows but some of which had something other than a 0 or 1 in the cell to designate the NaN in the original a matrix. Along this line, the NaN's were kept in the Y matrix. The goal was to be able to somehow ignore the NaN's then when calculating beta and betah.
It was considered to keep the NaN's because it was thought that the problems with the output had to do with some X, Xh, and Y matrices having less than 43 rows while others had 43 rows. Then looking back at how linear systems of equations are solved from what I learned in my linear algebra class, I realized that the number of rows could not be the reason why the some rows of genes in the output were completely blank. Even though the number of rows changed, the number of columns corresponding to the different type points never did. All the X matrices had 10 columns and all the Xh matrices had 5 columns. Therefore, there should be 5 betas being calculated (one for each timepoint) for the wildtype and for dGLN3 for all of the genes regardless of whether or not a gene had NaN's.
The next step in troubleshooting the code will be testing it with just a select number of genes that have NaN's for some timpeoints and flasks.
Katrina Sherbina 23:44, 20 May 2012 (EDT)
Week 2
May 21, 2012
The two strain comparison code was tested with an output that contained only seven genes (six of which had an NaN). At this point, there were no blank rows. To fix the blank rows that appear for genes with NaN's, the computation for the y loop was put in the same loop as the computations for the X and Xh matrices. This ensured that when the betas and betahswere being caluclated that the X and Xh matrices changed accordingly (depending on the gene for which the computation was being done). Previously, the code was set up so that he betas and betahs were calculated with the last X and Xh matrix that was computed in the loop.
Several changes were made to this code to compare all strains (wildtype, dCIN5, dGLN3, dHMO1, and dZAP1). First, indices were created for dCIN5, dGLN3, dHMO1, and dZAP1 (the two strain comparison code only had matrices for the widlytpe and dGLN3). The indices for the deletion strains were constructed for each timepoint similar to how they were constructed for dGLN3 making sure to specifiy the appropriate columns in the original data for each strain.
So that the X and Xh matrices were of the correct size and had 1's and 0's positions in the correct columns and rows, the indices for each timepoint for each strain were adjusted in the for loop.
i3 = ind1202m(end); i4 = indx3(1);
ind153m = ind153-i4+1+i3; ind303m = ind303-i4+1+i3; ind603m = ind603-i4+1+i3; ind903m = ind903-i4+1+i3; ind1203m = ind1203-i4+1+i3;
i5 = ind1203m(end); i6 = indx4(1);
ind154m = ind154-i6+1+i5; ind304m = ind304-i6+1+i5; ind604m = ind604-i6+1+i5; ind904m = ind904-i6+1+i5; ind1204m = ind1204-i6+1+i5;
i7 = ind1204m(end); i8 = indx5(1);
ind155m = ind155-i8+1+i7; ind305m = ind305-i8+1+i7; ind605m = ind605-i8+1+i7; ind905m = ind905-i8+1+i7; ind1205m = ind1205-i8+1+i7;
In addition, the the calculations for t, t2, and N were expanded to include all of the indices for all of the strains. The X and Xh matrix for genes with no NaN's had dimensions 103x25 and 103x5, respectively. The X and Xh matrices for the genes with NaN's had dimensions 94x25 and 94x5, respectively.
The coputation for the Y matrix was also altered to include all of the strains:
Y = [a(ii,indx)';a(ii,indx2)';a(ii,indx3)';a(ii,indx4)';a(ii,indx5)'];
The lines of code corresponding to the dimensions of the out_data were also altered to include all of the betas and betahs.
In addition, a few plots from the all strain comparison of the time versus expression genes with both low and high p values were saved. To do so, the if statement for the plots was removed. The save statement used to save each plot generated separately was
saveas(gcf,['Fig' num2str(ii) '.jpg'])
It would be beneficial to change this line of code so that the name of the figure corresponding to the gene for which the expression was being plotted.
Also, an avi file was made of the plots of time versus expression for genes with a p value of less than 0.05 and then those with less than 0.04 (this p value was observed to be the cut off for the B&H corrected p values in the final output for the all strian comparison). The plots for the genes with a p value less than 0.04 have to be redone because a banner that popped up in the figure window was also recorded (it blocks the gene name).
Before the for loop, the following line of code was added
aviobj = avifile('allstrain.avi','fps',15,'compression','Indeo3');
Some fiddling may still need to be done with the number of frames per second.
The following lines of code were added at the end of the for loop (after the if statement for to generate the plots):
%To get the handles from the figure open. Saves the entire figure window that is open. frame = getframe(gcf); %Adds the currently open plot to the avi file. aviobj = addframe(aviobj,frame);
Right after the for loop, a command line was added to close the avi file:
aviobj = close(aviobj);
The avi file had to be opened in QuickTime player. It would not open in windows media player when either the compression format was Indeo3 or Indeo5.
May 22, 2012
Some changes were made to the code to write an avi file. However, despite all of the revisions, the avi file was not being made properly. For the first few seconds the plots that showed up for some of the genes looked nothing like they should. To be able to even play the file, the compression had to be set to 'None'. This is what probably caused the file to be 4 GB. They appearance of the strange graphs at the beginning of the movie may be attributed to the large size of the file.
Initially, the p values showing up in the plots in the avi file did not match the p values for the corresponding genes in the output. However, this error was fixed by adjusting the if statement in the loop that creates the plots for genes with a p value of less than 0.05.
if pval<0.05
figure(1),plot(t,Y,'bo'),hold on,plot(t,Ym,'r+','LineWidth',3),plot(t,Yh,'k^','LineWidth',3),
title([b{ii+1,1} ' W = ',num2str(W),' pval = ',num2str(pval)]),...
xlabel('time, minutes'),ylabel('log fold change'),legend('data','model'),drawnow
hold off
nsig = nsig + 1;
[ii,nsig]
end
The rest of today was spent dissecting the code we currently have for the deterministic model. The goal is to split the code for the general parameter estimation function into multiple parts, namely, reading the files into matlab, the least squares error function, plotting the graphs, and writing the output.
Katrina Sherbina 20:38, 22 May 2012 (EDT)
May 23, 2012
The lines of code to create an avi file from the plots from the all strains comparison were fixed (thanks to Dr. Fitzpatrick). The changes include the following line of code before the for loop:
figure(1) set(gcf, 'DoubleBuffer','on'); yourMovieFileName = ['all_strain_comparison_p_value_less_than_0.05_Cinepak.avi'];
The frame = getframe(gcf) and aviobj = addframe(aviobj,frame) inside the if statement to generate the plots was removed. The aviobj=close(aviobj) line of code immediately after the for loop was replaced with
movie2avi(movieArray,yourMovieFileName,'compression','Cinepak','fps',5,'quality',100);
This time the movie generated with Cinepak compression code be played on Windows Media Player on an XP operating system. The size of this movie was significantly smaller than the size of the previous movie made without compression that could play with Windows Media Player.
There are three different scripts for the statistical analysis in MATLAB:
- all_strains_compare_movie.m to do a comparison for all strains
- wildtype_deletion_strain_compare.m to do a comparison between the wildtype and another strain
- two_deletion_strain_compare.m to do a comparison between two deletion strains
The two strain comparison code was divided into two (the scripts in the last two bullets above) because the for loop designed to remove NaN's is not necessary for the deletion strain data. The principle differences between these two scripts is the calculations to renumber the indices. The following is the part of the script for the wildtype to deletion strain comparisons that contains differences with the other two strain comparison script:
%These indices are changed depending on what deletion strain is being analyzed ind152 = [33,34,35,36]; ind302 = [37,38,39,40]; ind602 = [41,42,43,44]; ind902 = [45,46,47,48]; ind1202 = [49,50,51,52];
p = 10; q = 5;
alpha = 0.05; % significance level for the tests
n152 = length(ind152); n302 = length(ind302); n602 = length(ind602); n902 = length(ind902); n1202 = length(ind1202);
indx2 = [ind152,ind302,ind602,ind902,ind1202];
n = length(a(:,1));
out_data = zeros(n,12); nsig = 0;
for ii=1:n
I = find(~isnan(a(ii,1:23)));
ind15 = I(find(I>=1&I<=4));
ind30 = I(find(I>=5&I<=9));
ind60 = I(find(I>=10&I<=13));
ind90 = I(find(I>=14&I<=18));
ind120 = I(find(I>=19&I<=23));
indx = [ind15,ind30,ind60,ind90,ind120];
indX = find(indx);
if length(indx)<23
ind15 = find(ind15);
ind30 = ind30 - ind15(end);
ind60 = ind60 - ind30(1) - 1;
ind90 = ind90 - ind60(1);
ind120 = ind120-ind90(1);
end
n15 = length(ind15);
n30 = length(ind30);
n60 = length(ind60);
n90 = length(ind90);
n120 = length(ind120);
i1 = indX(end);
i2 = indx2(1);
ind152m = ind152-i2+1+i1; ind302m = ind302-i2+1+i1; ind602m = ind602-i2+1+i1; ind902m = ind902-i2+1+i1; ind1202m = ind1202-i2+1+i1; t = [ones(n15,1)*15;ones(n30,1)*30;ones(n60,1)*60;ones(n90,1)*90;ones(n120,1)*120;... ones(n152,1)*15;ones(n302,1)*30;ones(n602,1)*60;ones(n902,1)*90;ones(n1202,1)*120];
t2 = [min(t):5:max(t)]'; N = length(indx)+length(indx2);
X = zeros(N,p); Xh = zeros(N,p-q);
Xh(ind15,1) = 1; Xh(ind30,2) = 1; Xh(ind60,3) = 1; Xh(ind90,4) = 1; Xh(ind120,5) = 1; Xh(ind152m,1) = 1; Xh(ind302m,2) = 1; Xh(ind602m,3) = 1; Xh(ind902m,4) = 1; Xh(ind1202m,5) = 1;
X(ind15,1) = 1; X(ind30,2) = 1; X(ind60,3) = 1; X(ind90,4) = 1; X(ind120,5) = 1; X(ind152m,6) = 1; X(ind302m,7) = 1; X(ind602m,8) = 1; X(ind902m,9) = 1; X(ind1202m,10) = 1; end
The calculations for the Y matrices, beta's, betah's, fstats, pvals, and plots are within the for loop. These are the same calculations in the deletion strain to deletion strain comparison script. These caluclations are the only ones in the for loop for the deletion strain to deletion strain comparison script. The following is a part of the deletion strain to deletion strain comparison script that contains differences in comparison to the previous script.
%These indices are changed depending on what deletion strain is being analyzed ind15 = [91,92,93,94]; ind30 = [95,96,97,98]; ind60 = [99,100,101,102]; ind90 = [103,104,105,106]; ind120 = [107,108,109,110];
ind152 = [120,121,122,123]; ind302 = [124,125,126,127]; ind602 = [128,129,130,131]; ind902 = [132,133,134,135]; ind1202 = [136,137,138,139];
p = 10; q = 5;
alpha = 0.05; % significance level for the tests
n15 = length(ind15); n30 = length(ind30); n60 = length(ind60); n90 = length(ind90); n120 = length(ind120); n152 = length(ind152); n302 = length(ind302); n602 = length(ind602); n902 = length(ind902); n1202 = length(ind1202);
indx = [ind15,ind30,ind60,ind90,ind120]; indX = find(indx); indx2 = [ind152,ind302,ind602,ind902,ind1202];
i0 = ind15(1); i1 = indX(end); i2 = indx2(1);
ind15m = find(ind15); ind30m = ind30-i0+1; ind60m = ind60-i0+1; ind90m = ind90-i0+1; ind120m = ind120-i0+1;
ind152m = ind152 -i2+1+i1; ind302m = ind302 -i2+1+i1; ind602m = ind602 -i2+1+i1; ind902m = ind902 -i2+1+i1; ind1202m = ind1202 -i2+1+i1;
t = [ones(n15,1)*15;ones(n30,1)*30;ones(n60,1)*60;ones(n90,1)*90;ones(n120,1)*120;... ones(n152,1)*15;ones(n302,1)*30;ones(n602,1)*60;ones(n902,1)*90;ones(n1202,1)*120];
t2 = [min(t):5:max(t)]';
N = length(indx)+length(indx2);
X = zeros(N,p); Xh = zeros(N,p-q);
Xh(ind15m,1) = 1; Xh(ind30m,2) = 1; Xh(ind60m,3) = 1; Xh(ind90m,4) = 1; Xh(ind120m,5) = 1; Xh(ind152m,1) = 1; Xh(ind302m,2) = 1; Xh(ind602m,3) = 1; Xh(ind902m,4) = 1; Xh(ind1202m,5) = 1;
X(ind15m,1) = 1; X(ind30m,2) = 1; X(ind60m,3) = 1; X(ind90m,4) = 1; X(ind120m,5) = 1; X(ind152m,6) = 1; X(ind302m,7) = 1; X(ind602m,8) = 1; X(ind902m,9) = 1; X(ind1202m,10) = 1;
n = length(a(:,1));
out_data = zeros(n,12);
nsig = 0;
The dimensions of the X and Xh matrices in this case stay the same for all genes.
The rest of today was spent continuing to make changes to the MATLAB script for the deterministic model. Nick succeeded in separating the general parameter estimation function into four parts:
- Input and labeling of the master sheet containing the network, log2 fold change concentrations for each strain, optimization parameters, production rates, degradation rates, network b's, and simulation times.
- The call for the general least squares error and optimization
- Generating graphs
- Creating the output
The input of the mastersheet in this script is different from that in the previous general parameter estimation function. In the actual input sheet, the log2 concentrations were divided into several sheets (one for each deletion strain) and inputted sheet by sheet into MATLAB. In addition, the log2 concentrations for each flask for each timepoint were inputted rather than the average log2 concentrations for each timepoint. Indicators for the timepoints, strain name, and strain index were added to the optimization parameters sheet. This required modifying the for loop for inputting the optimization parameters.
The next steps with the deterministic model include making some changes to the general least squares error script where the error = ((log2(x1) - d1).^2) line of code needs to be modified to take into account the replicates in the log fold change data for each timepoint for each strain.
Katrina Sherbina 21:22, 23 May 2012 (EDT)
May 24, 2012
In order to modify the error calculation in the general least squares error function (in MATLAB) to take into account the change from using average log fold changes for each timpepoint to using all of the data, several modifications were made to other scripts. (To make any of these changes, I first dissected the various functions to understand the various parts more.)
The log fold change data had to be imported separately for each strain (this modification was made earlier in the week):
[log2FC(1).d,TX1] = xlsread(input_file,'wt_log2FC'); [log2FC(2).d,TX1] = xlsread(input_file,'dCIN5_log2FC'); [log2FC(3).d,TX1] = xlsread(input_file,'dGLN3_log2FC'); [log2FC(4).d,TX1] = xlsread(input_file,'dHMO1_log2FC'); [log2FC(5).d,TX1] = xlsread(input_file,'dZAP1_log2FC');
Nevertheless, averages for each timepoint were later calculated for the purposes of later plotting standard deviations on the time versus log fold change plots.
The driver to detect active genes will be taken out. All of the genes in the network are considered active genes. The assignment of a number to each active gene was done with the detect active genes script. Now the active changes are "detected" in the matrices with the log fold change data:
active = find(log2FC(1).d(2:end,1));
This array is called for in the general network dynamic function. The ode45 function in matlab approximates a solution to the differential equation defined in the general network dynamic function.
Previously, the error was calculated using the model and the average log fold changes for each timepoint from the data. However, this error calculation needs to be changed so that the unaveraged log fold changes can be used. A preliminary equation that was experiment with was
myerror = ((log2(x1) - d1([1 5 10],:)).^2)+((log2(x1) - d1([2 6 11],:)).^2)+ ((log2(x1) - d1([3 7 12],:)).^2) +((log2(x1) - d1([4 8 13],:)).^2);
The problem here is that the equation does not take into account the last replicate for the t30 timepoint (because t30 has five replicates while the other timpeoints have four replicates).
Another concern is the actual data d and d1 that is being fed into the general least squares error function and other functions. For now, both d and d1 correspond to just the wildtype data. However, we will want to use the deletion strain data for the model as well.
Katrina Sherbina 21:29, 24 May 2012 (EDT)
May 25, 2012
The functions written in MATLAB for the deterministic model were modified to include all of the data (not just the average log fold changes for each timepoint) in the various calculations.
Modifications were made in the parameters script and the error calculation (now called errormat) in the general least quares error function. The following lines of code were added to the script:
% % The first row of the log2FC data indicating all of the replicate timepoints
td = (log2FC(1).d(1,:));
% % Finds the indices in td that correspond to each timepoint in tspan.
tindex(1).indx = find(td == tspan(1));
tindex(2).indx = find(td == tspan(2));
tindex(3).indx = find(td == tspan(3));
n_times = length(tspan);
% % Log2FC data for the wildtype
d = (log2FC(1).d(2:end,:))';
% % The average log2FC for each timepoint for each gene.
means = zeros(length(log2FC(1).d(2:end,1)),n_times);
for iT = 1:n_times
means(:,iT) = mean(log2FC(1).d(2:end,tindex(iT).indx(1):tindex(iT).indx(end)),2);
end
% % Standard deviations for each timepoint for each gene. Necessary for the graphs
% % that will be produced.
sigmas = zeros(length(log2FC(1).d(2:end,1)),n_times);
for iT = 1:n_times
sigmas(:,iT) = std(log2FC(1).d(2:end,tindex(iT).indx(1):tindex(iT).indx(end)),0,2);
end
It is necessary to still find the standard deviations and the means to be able to produce graphs for each gene of the log fold change versus time with the upper and lower error bounds for each timepoint.
The error calculation in the general least squares error function was modified with the following lines of code:
errormat = 0;
for iT = 1:length(tspan)
for iF = 1:length(tindex(iT).indx)
errormat = errormat+(log2(x1(iT,:))-d(tindex(iT).indx(iF),:)).^2;
end
end
% This is the sum of all the errors
L = sum(errormat(:));
In the general least squares error function, a modification was also made to the first figure produced (with the counter and lse_out). At first, all of the data (not the average log fold changes) was being used to specify the second subplot in the figure which fits the model to the data. When viewing the figure, the model was being fit to only three of the thirteen sets of data points in the second subplot. There were thriteen sets of data points (in columns) since there are thirteen columns in the data due to the replicates in each timepoint. Instead of using all of the data, the code to produce the subplots was modified to instead use the average log fold change for each timepoint:
if(100*round(counter/100)==counter) figure(1),subplot(211),plot(theta,'d'), title(['counter = ' num2str(counter) ' lseout = ' num2str(lse_out)]) subplot(212),plot(means','*'),hold on,plot(log2(x1)), hold off,pause(.1) end
The script specifying the plots for each gene of log fold change versus time was modified. The if else statement having to do with cssigflg was removed. This designation had to do with the concentration sigmas (standard deviations) inputted into MATLAB when they were still in the input excel spreadsheet. However, since the standard deviations are now being calculated in MATLAB, there is no need to determine the kind of the plot to be produced depending on the value of cssigflg. The following is the new code for the plots:
if igraph == 1
concd = 2.^d(:,active);
error_up = (means(active,:)+1.96*sigmas);
error_dn = (means(active,:)-1.96*sigmas);
for ii=1:n_active_genes
figure(ii+2),hold on
plot(time,log2(model(:,ii)),td,d(:,active(ii)),'o',tspan,error_up(ii,:),'r+',tspan,
error_dn(ii,:),'r+','LineWidth',3),axis([tmin tmax -3 3])
title(TX1{1+(ii),2},'FontSize',16)
xlabel('Time (minutes)','FontSize',16)
ylabel('Expression (log2 fold change)','FontSize',16)
end
end
The global statement at the top of each script or function was also changed due to the aforementioned changes:
global A prorate degrate active inactive tspan d wts b alpha lse_out penalty_out counter parms parmnames tindex means
Further work will need ot be done on the plots. In the plots for MAL33 and ZAP1, the upper error bounds are not shown. The y axis may not to be adjusted to be larger than [-3,3]. Also, the model for HAP5, the model goes beyond the error bounds for the t60 timepoint.
In the output excel spreadsheet, there are no log2 optimized concentrations nor network optimized B's for FHL1, SKO1, or SWI6.
In addition to dealing with the aforementioned concerns, we need to explore how to include the deletion strain data in the model.
Katrina Sherbina 20:23, 25 May 2012 (EDT)
