Egem2

matlab/MeCorr/constrain_flux_regulation.m

@@ -1,4 +1,5 @@
-function [fluxstate_gurobi,grate, solverobj] =  constrain_flux_regulation(model1,onreactions,offreactions,kappa,rho,epsilon,mode,epsilon2,minfluxflag)
+function [fluxstate_gurobi, grate, solverobj] =  constrain_flux_regulation(model1,...
+    onreactions,offreactions,kappa,rho,epsilon,mode,epsilon2,minfluxflag)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if (~exist('mode','var')) || (isempty(mode))

matlab/MeCorr/methylVariables.m

@@ -1,8 +1,10 @@
 % Plot correlation values between AUC of growth inhibition of cell 
 % lines and methylation flux
 
-% Workflow: 1)Run make_eGEM 2)Run this script 3)Assign epsilon_methylation
-% 3)Run last module of MATLAB_CODE_methyl (3 sections)
+% Workflow: 1)Run make_eGEM  
+cd '.\..\MeCorr';
+% 2)Run methylVariables 3)Assign epsilon_methylation 4)Run last module of 
+...MATLAB_CODE_methyl (3 sections)
 
 % initCobraToolbox;
 % changeCobraSolver('gurobi');
@@ -21,16 +23,21 @@
 %index_DNMTi = [2 , 78, 181, 402, 444];     % DNMT = DNA MT
 
 % Get all 3 rows of s3 for drug_auc_me. Replaces drug_auc_expt
-drug_auc_me= zeros(664*8, 3);
+% 7/1: drug_auc_expt has all drugs tested. drug_auc_me only has methyl drugs 
+% advantage: less to search. matlab_code already extracts interested drugs
+% (variable hmei_auc_expt_me).
+drug_auc_me= s3;
+
+% drug_auc_me= zeros(664*8, 3);
+% drug_auc_me(1:597, :)= s3(1293:1889, :); % drug 3
+% drug_auc_me(665:1291, :)= s3(96133:96759, :); %drug 178
+% drug_auc_me(1328:1766, :)= s3(149039:149477, :); %drug 280
+% drug_auc_me(1992:2435, :)= s3(201699:202142, :); %drug 374
+% drug_auc_me(2656:3302, :)= s3(204605:205251, :); %drug 380
+% drug_auc_me(3320:3630, :)= s3(228209:228519, :); %drug 421
+% drug_auc_me(3984:4606, :)= s3(232799:233421, :); %drug 431
+% drug_auc_me(4648:4718, :)= s3(257167:257237, :); %drug 475
 
-drug_auc_me(1:597, :)= s3(1293:1889, :); % drug 3
-drug_auc_me(665:1291, :)= s3(96133:96759, :); %drug 178
-drug_auc_me(1328:1766, :)= s3(149039:149477, :); %drug 280
-drug_auc_me(1992:2435, :)= s3(201699:202142, :); %drug 374
-drug_auc_me(2656:3302, :)= s3(204605:205251, :); %drug 380
-drug_auc_me(3320:3630, :)= s3(228209:228519, :); %drug 421
-drug_auc_me(3984:4606, :)= s3(232799:233421, :); %drug 431
-drug_auc_me(4648:4718, :)= s3(257167:257237, :); %drug 475
 % Col 1 is drug id #, col 2 is cell line id #, col 3 is auc. Switch col 2 
 ...and 3 to match drug_auc_expt
 celllineid_tmp= drug_auc_me(:, 2);
@@ -44,14 +51,14 @@
 drug_auc_me.Properties.VariableNames{'drug_auc_me3'}='index_ccl';
 
 % Other variables to replace with methylation data
-ctd2celllineidname_id_me= s2(:, 1);
+ctd2clidname_id_me= s2(:, 1);
 
-ctd2celllineidname_me= [s2(:,2), s2(:,4), s2(:,5)]; 
-ctd2celllineidname_me.Properties.VariableNames{'cell_line_name'}='celllinenames_ccle1';
+ctd2clidname_me= [s2(:,2), s2(:,4), s2(:,5)]; 
+ctd2clidname_me.Properties.VariableNames{'cell_line_name'}='celllinenames_ccle1';
 
 s1= readtable('Data\Ludlow2015_SmallMolecInformer.xlsx','Sheet','s1','Range','A:B');
-ctd2compoundidname_id_me= s1(:, 1); 
-ctd2compoundidname_name_me= s1(:, 2);
+ctd2cpdidname_id_me= s1(:, 1); 
+ctd2cpdidname_name_me= s1(:, 2);
 
 % Below: written into MATALB_CODE_methyl
 % exptidcelllinemediamatch= array2table(exptidcelllinemediamatch);

matlab/scripts/make_eGEM.m

@@ -2,8 +2,12 @@
 % @author: Scott Campit & Lauren Fane
 
 % Initialize parameters
-initCobraToolbox;
-changeCobraSolver('gurobi');
+x=input('Initialize Cobra Toolbox and change solver? "yes"/"no": ');
+if x == "yes"
+    initCobraToolbox;
+    changeCobraSolver('gurobi');
+elseif x == "no"
+end
 
 % Load AcGEM model (Shen et al., 2019)
 load ./../shen-et-al/supplementary_software_code acetylation_model;
@@ -248,6 +252,20 @@
 %     'checkDuplicate', 'true');
 %same chemical rxn as H3K9 methylation -> should i add it? How do i
 %distinguish between locations?
-%% save
+
+% Fill in empty reactions
+min = addReaction(min,'SLC17A1','reactionName','Sodium/phosphate cotransporter 1',...
+    'reactionFormula','po4[e] + 3 na1[e] -> po4[c] + 3 na1[c]',...
+    'subSystem','Transport, Extracellular',...
+    'geneRule','(SLC17A1) or (SLC17A2) or (SLC17A3) or (SLC17A4)');
+min = addReaction(min,'SLC17A3','reactionName','Sodium/phosphate cotransporter 4',...
+    'reactionFormula','po4[e] + 3 na1[e] -> po4[c] + 3 na1[c]',...
+    'subSystem','Transport, Extracellular',...
+    'geneRule','(SLC17A1) or (SLC17A2) or (SLC17A3) or (SLC17A4)');
+%% save 
+% rename so that variable does not conflict with the function min
 min = checkDuplicateRxn(min,'S',1,1);
-save('./../models/min.mat', 'min');
+%save('./../models/min.mat', 'min');
+%save('./../models/min2SLC.mat', 'min');
+min_model=min;
+clearvars min