Conversation
This reverts commit da8f33c.
…anges before starting new branch.
…o make scatter plots and calculate correlation in MATLAB_CODE_me. Saved outputs constrain_flux_regulation in VariablesSaved. Renamed min variable to min_model to avoid conflict with min function.
… values. Saved resulting figures.
…into structures called sigRxn... while long program ran on desktop.
…n PC, one on desktop.
…alue. compared 1E-1 FBA and 1E-1 fBA and saved variables.�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D�[D
…uxes_allrxns for all 6 weights.
… positive or negative.
| offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); | ||
|
|
||
| %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | ||
| % set the glucose uptake based on media |
There was a problem hiding this comment.
These lower bounds only change glucose abundance for 4 medium conditions. The code matlab/scripts/nutrient_sensitivty/medium_LB_constraints.m in master can improve your drug sensitivity predictions.
| load supplementary_software_code hdacexpfcs hdacexpallgeneids | ||
| %contains gene expression data after treatment with hdac inhibitors | ||
|
|
||
| % Converting arrays to tables and renaming some Variables to make |
There was a problem hiding this comment.
Once you prettify these variables, you should save them in vars and just load them rather than formatting them during every iteration within this script.
ScottCampit
left a comment
There was a problem hiding this comment.
@fanela, there are several changes that need to be made before pushing this to the master branch. General stuff that I saw:
- Variable names are still super confusing - please help make the code readable by changing the variables names to something more understandable (ie instead of using
v1, useKMTi_aucdirectly, which is much more clear. - There's a ton of code that is repetitive between the three new scripts in the
MeCorrdirectory. Since you're reusing these variables, save them as.matfiles and load them at the beginning of each script. That way, it's easier to figure out what this script is actually doing. - I found a few bugs that should be addressed and that should change your results and potentially your interpretation of the results.
Also, I moved some of your old matlab scripts in the drug_sensitivity folder. I though MeCorr was too generic, because there are several modules that compute correlation values between expression and Me flux. If you can, move these scripts there and save the new matlab variables in the vars directory.
| fluxes_allrxns= NaN(height(exptidcelllinemediamatch),length(model2.rxns)); | ||
| g_rate= NaN(height(exptidcelllinemediamatch), length(model2.rxns)); | ||
|
|
||
| for i = 1:height(exptidcelllinemediamatch) |
There was a problem hiding this comment.
Can you change i, ii, and iii into better variable names that are more accurate? It becomes really confusing reading the rest of the code
There was a problem hiding this comment.
Repeat this for a lot of the variables below. Our goal is to have this code more readable so that it is easy to figure out what's going on without comments.
| @@ -0,0 +1,419 @@ | |||
| %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |||
There was a problem hiding this comment.
Rename the three files with MATLAB_CODE_* files - it's not clear from the get-go what the difference between allrxns, me, and methyl are without opening and carefully reading the file contents.
| ii = find(ismember(ctd2celllineidname_id, exptidcelllinemediamatch(i,2))); | ||
| iii = find(ismember(celllinenames_ccle1, ctd2celllineidname(ii,1))); | ||
| if ~isempty(iii) | ||
| iii = iii(1); |
There was a problem hiding this comment.
Why is iii set to the first value if it's not empty? If length(iii) > 1, then you may want to average ccle_expression_metz for those columns, then get ongenes and offgenes. This gets the average gene expression values for the duplicated cell line.
| rho_p= NaN(height(hmei_list),length(model2.rxns)); | ||
| for j = 1:height(hmei_list) | ||
| fx = find(ismember(ctd2compoundidname_name, hmei_list(j,1))); | ||
| ix = ismember(drug_auc_expt(:,3), ctd2compoundidname_id(fx,1)); sum(ix); |
There was a problem hiding this comment.
separate code after semicolon to new line to improve readability
| %% Create Table of Significant Reactions by Correlation value | ||
| % Workflow: Change threshold. Change struct field name (e.g. above3) | ||
| % accordingly | ||
| sigTF= (abs(rho) > 0.3); |
There was a problem hiding this comment.
Why are you taking the absolute value of rho? Negative values are important and have a statistical interpretation.
There was a problem hiding this comment.
Also if you want to get significant reactions, you would want to filter using the p-values as well as R value. That is, you want the reaction to have an R-value >= 0.3 and a p-value < 0.05.
There was a problem hiding this comment.
Line 195: I coded abs(rho) > 0.3 so that the negative correlations below -0.3 would be selected as well. I am including those correlations in downstream analysis.
There was a problem hiding this comment.
I see, you should create two variables have rho > 0.3 and rho < -0.3, unless you wrote in a way to distinguish between the negatively and positively correlated entries by saving their respective indices first. That way, you can combine them later to make a plot while keeping the respective scalar directions.
|
|
||
|
|
||
| [ix, aapos] = ismember( mediareactions1(20:37),acetylation_model.rxns); | ||
| posgluc = 1385; % glucose uptake reaction in recon1. |
There was a problem hiding this comment.
Don't hard-code this - use find(ismember(model.rxns, 'EX__glc_D_e')) to find the glucose exchange reaction. Same for glutamine and linoleic acid.
|
|
||
|
|
||
| model3.c(rxnpos) = epsilon_acetylation; | ||
| model3.c(rxnpos) = epsilon_methylation; |
|
|
||
|
|
||
| model3.c(rxnpos) = epsilon_acetylation; | ||
| model3.c(rxnpos) = epsilon_methylation; |
There was a problem hiding this comment.
Also repetitive - model3.c should not change the way your code is written.
| acetgenessorted = unqgenes(spos); | ||
| [sx spos] = sort(acet_screen_rpmi(:,1)); | ||
| acetgenessorted = unqgenes(spos); |
There was a problem hiding this comment.
change variable names to something more generic since this code looks at methylation now
|
|
||
| [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); | ||
| if ~isempty(solf.x) && ~isnan(solf.x) | ||
| acet_screen_rpmi(i,2) = solf.x(objpos); % impact on growth |
There was a problem hiding this comment.
Change variable name so it's more obvious that it's the growth rate from the constrained metabolic model. Also, are you intentionally using the acetylation model? It is confusing since you're storing the objective coefficient value into the variable name epsilon_methylation.
No description provided.