diff --git a/git b/git new file mode 100644 index 0000000..e69de29 diff --git a/matlab/MeCorr/CorrScatter1E-2/BIX.jpg b/matlab/MeCorr/CorrScatter1E-2/BIX.jpg new file mode 100755 index 0000000..375bf57 Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/BIX.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/BRD.jpg b/matlab/MeCorr/CorrScatter1E-2/BRD.jpg new file mode 100755 index 0000000..6aebcfe Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/BRD.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/CBB.jpg b/matlab/MeCorr/CorrScatter1E-2/CBB.jpg new file mode 100755 index 0000000..104443d Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/CBB.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/GSK.jpg b/matlab/MeCorr/CorrScatter1E-2/GSK.jpg new file mode 100755 index 0000000..416c2ca Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/GSK.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/QW-BI.jpg b/matlab/MeCorr/CorrScatter1E-2/QW-BI.jpg new file mode 100755 index 0000000..ba2de07 Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/QW-BI.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/UNC0321.jpg b/matlab/MeCorr/CorrScatter1E-2/UNC0321.jpg new file mode 100755 index 0000000..76075a0 Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/UNC0321.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/UNC0638.jpg b/matlab/MeCorr/CorrScatter1E-2/UNC0638.jpg new file mode 100755 index 0000000..5050930 Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/UNC0638.jpg differ diff --git a/matlab/MeCorr/CorrScatter1E-2/methylstat.jpg b/matlab/MeCorr/CorrScatter1E-2/methylstat.jpg new file mode 100755 index 0000000..aba638e Binary files /dev/null and b/matlab/MeCorr/CorrScatter1E-2/methylstat.jpg differ diff --git a/matlab/MeCorr/Data/Ludlow2015_SmallMolecInformer.xlsx b/matlab/MeCorr/DataLudlow2015/Ludlow2015_SmallMolecInformer.xlsx similarity index 100% rename from matlab/MeCorr/Data/Ludlow2015_SmallMolecInformer.xlsx rename to matlab/MeCorr/DataLudlow2015/Ludlow2015_SmallMolecInformer.xlsx diff --git a/matlab/MeCorr/Data/Ludlow_SelectDrugs.xlsx b/matlab/MeCorr/DataLudlow2015/Ludlow_SelectDrugs.xlsx similarity index 100% rename from matlab/MeCorr/Data/Ludlow_SelectDrugs.xlsx rename to matlab/MeCorr/DataLudlow2015/Ludlow_SelectDrugs.xlsx diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S10_mastercpdid&PC.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S10_mastercpdid&PC.xlsx new file mode 100755 index 0000000..9fad201 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S10_mastercpdid&PC.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S11_aboutCellLines.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S11_aboutCellLines.xlsx new file mode 100755 index 0000000..b8f97d8 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S11_aboutCellLines.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S12_set11.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S12_set11.xlsx new file mode 100755 index 0000000..7c1fe09 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S12_set11.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S13_GSEA.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S13_GSEA.xlsx new file mode 100755 index 0000000..ac7f472 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S13_GSEA.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S14_lipid.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S14_lipid.xlsx new file mode 100755 index 0000000..bbba266 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S14_lipid.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S15_mapCCLEtoCTD2.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S15_mapCCLEtoCTD2.xlsx new file mode 100755 index 0000000..ea0760e Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S15_mapCCLEtoCTD2.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S2_ccl.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S2_ccl.xlsx new file mode 100755 index 0000000..11c5358 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S2_ccl.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S3_cpd.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S3_cpd.xlsx new file mode 100755 index 0000000..9d1c282 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S3_cpd.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S4_auc.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S4_auc.xlsx new file mode 100755 index 0000000..35a449e Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S4_auc.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S5_set4.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S5_set4.xlsx new file mode 100755 index 0000000..f695009 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S5_set4.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S6_SignifCorr.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S6_SignifCorr.xlsx new file mode 100755 index 0000000..f75747b Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S6_SignifCorr.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S8_cpd&stats.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S8_cpd&stats.xlsx new file mode 100755 index 0000000..a270b94 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S8_cpd&stats.xlsx differ diff --git a/matlab/MeCorr/DataRees2016/nchembio.1986-S9_set8.xlsx b/matlab/MeCorr/DataRees2016/nchembio.1986-S9_set8.xlsx new file mode 100755 index 0000000..f731848 Binary files /dev/null and b/matlab/MeCorr/DataRees2016/nchembio.1986-S9_set8.xlsx differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_boxplot.jpg new file mode 100755 index 0000000..5f01792 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_scatter.jpg new file mode 100755 index 0000000..ac1133e Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/bix_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_boxplot.jpg new file mode 100755 index 0000000..9762ff7 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_scatter.jpg new file mode 100755 index 0000000..058c812 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/brd_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_boxplot.jpg new file mode 100755 index 0000000..126a138 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_scatter.jpg new file mode 100755 index 0000000..ea18c34 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/cbb_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_boxplot.jpg new file mode 100755 index 0000000..9afdab8 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_scatter.jpg new file mode 100755 index 0000000..111bc2b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/gsk_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_boxplot.jpg new file mode 100755 index 0000000..5d7304b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_scatter.jpg new file mode 100755 index 0000000..c8a1226 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/methylstat_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_boxplot.jpg new file mode 100755 index 0000000..86b22b5 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_scatter.jpg new file mode 100755 index 0000000..e954b9c Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/qw-bi_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_boxplot.jpg new file mode 100755 index 0000000..3592382 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_scatter.jpg new file mode 100755 index 0000000..68f8ebd Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0321_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_boxplot.jpg new file mode 100755 index 0000000..0df9b7b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_scatter.jpg new file mode 100755 index 0000000..77cde3b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-1_scatterAndBoxplot/unc0638_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/bix.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/bix.jpg new file mode 100755 index 0000000..6ec0ecd Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/bix.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/brd.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/brd.jpg new file mode 100755 index 0000000..247c389 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/brd.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/cbb.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/cbb.jpg new file mode 100755 index 0000000..11f4cca Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/cbb.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/gsk.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/gsk.jpg new file mode 100755 index 0000000..112c3f8 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/gsk.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/methylstat.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/methylstat.jpg new file mode 100755 index 0000000..34157da Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/methylstat.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/qw-bi.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/qw-bi.jpg new file mode 100755 index 0000000..2cb2997 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/qw-bi.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0321.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0321.jpg new file mode 100755 index 0000000..e2b3269 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0321.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0638.jpg b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0638.jpg new file mode 100755 index 0000000..f5a103d Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-2_boxplot_methylFlux/unc0638.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_boxplot.jpg new file mode 100755 index 0000000..a68076c Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_scatter.jpg new file mode 100755 index 0000000..3856a14 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/bix_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_boxplot.jpg new file mode 100755 index 0000000..545bc23 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_scatter.jpg new file mode 100755 index 0000000..ae3aea0 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/brd_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_boxplot.jpg new file mode 100755 index 0000000..1527a65 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_scatter.jpg new file mode 100755 index 0000000..6f71b02 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/cbb_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_boxplot.jpg new file mode 100755 index 0000000..7a54edc Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_scatter.jpg new file mode 100755 index 0000000..57843d5 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/gsk_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_boxplot.jpg new file mode 100755 index 0000000..536fccb Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_scatter.jpg new file mode 100755 index 0000000..bc30c38 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/methylstat_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_boxplot.jpg new file mode 100755 index 0000000..ab8bdb6 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_scatter.jpg new file mode 100755 index 0000000..7877a79 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/qw-bi_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_boxplot.jpg new file mode 100755 index 0000000..cae153e Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_scatter.jpg new file mode 100755 index 0000000..e31d786 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0321_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_boxplot.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_boxplot.jpg new file mode 100755 index 0000000..0a63d0d Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_boxplot.jpg differ diff --git a/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_scatter.jpg b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_scatter.jpg new file mode 100755 index 0000000..006b6a3 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/1E-4_scatterAndBoxplot/unc0638_scatter.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BIX.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BIX.jpg new file mode 100755 index 0000000..2fd69d0 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BIX.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BRD.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BRD.jpg new file mode 100755 index 0000000..2a05512 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/BRD.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/CBB.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/CBB.jpg new file mode 100755 index 0000000..f6127a2 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/CBB.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GSK.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GSK.jpg new file mode 100755 index 0000000..7bf43aa Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GSK.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GreenCurve.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GreenCurve.jpg new file mode 100755 index 0000000..a0281e0 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/GreenCurve.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/MeFluxDistrib.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/MeFluxDistrib.jpg new file mode 100755 index 0000000..ef6593b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/MeFluxDistrib.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/PredictedvsObservedAUC.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/PredictedvsObservedAUC.jpg new file mode 100755 index 0000000..1b94ebd Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/PredictedvsObservedAUC.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/QW-BI.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/QW-BI.jpg new file mode 100755 index 0000000..5873251 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/QW-BI.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0321.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0321.jpg new file mode 100755 index 0000000..4919643 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0321.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0638.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0638.jpg new file mode 100755 index 0000000..4e00665 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/UNC0638.jpg differ diff --git a/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/methylstat.jpg b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/methylstat.jpg new file mode 100755 index 0000000..a0bf9f8 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/AcetExprData_pairedBoxplots_Rees2016/methylstat.jpg differ diff --git a/matlab/MeCorr/Figures050719/AcetylaFlux.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/AcetylaFlux.fig similarity index 100% rename from matlab/MeCorr/Figures050719/AcetylaFlux.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/AcetylaFlux.fig diff --git a/matlab/MeCorr/Figures050719/BulkH3K9Acetyl.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/BulkH3K9Acetyl.fig similarity index 100% rename from matlab/MeCorr/Figures050719/BulkH3K9Acetyl.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/BulkH3K9Acetyl.fig diff --git a/matlab/MeCorr/Figures050719/acetyl_levels_in_diffgrowthconditions.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/acetyl_levels_in_diffgrowthconditions.fig similarity index 100% rename from matlab/MeCorr/Figures050719/acetyl_levels_in_diffgrowthconditions.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/Figures050719_acet/acetyl_levels_in_diffgrowthconditions.fig diff --git a/matlab/MeCorr/Figures061219/Fig6incomplete.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/Fig6incomplete.fig similarity index 100% rename from matlab/MeCorr/Figures061219/Fig6incomplete.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/Fig6incomplete.fig diff --git a/matlab/MeCorr/Figures061219/LBH-589.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/LBH-589.fig similarity index 100% rename from matlab/MeCorr/Figures061219/LBH-589.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/LBH-589.fig diff --git a/matlab/MeCorr/Figures061219/LastMod_Distrib.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/LastMod_Distrib.fig similarity index 100% rename from matlab/MeCorr/Figures061219/LastMod_Distrib.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/LastMod_Distrib.fig diff --git a/matlab/MeCorr/Figures061219/Mod1.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/Mod1.fig similarity index 100% rename from matlab/MeCorr/Figures061219/Mod1.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/Mod1.fig diff --git a/matlab/MeCorr/Figures061219/belinostat.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/belinostat.fig similarity index 100% rename from matlab/MeCorr/Figures061219/belinostat.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/belinostat.fig diff --git a/matlab/MeCorr/Figures061219/differentialSensitivity.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/differentialSensitivity.fig similarity index 100% rename from matlab/MeCorr/Figures061219/differentialSensitivity.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/differentialSensitivity.fig diff --git a/matlab/MeCorr/Figures061219/entinostat.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/entinostat.fig similarity index 100% rename from matlab/MeCorr/Figures061219/entinostat.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/entinostat.fig diff --git a/matlab/MeCorr/Figures061219/vorinostat.fig b/matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/vorinostat.fig similarity index 100% rename from matlab/MeCorr/Figures061219/vorinostat.fig rename to matlab/MeCorr/Figures-LF/PreliminaryFigures/originalCodeAndData/vorinostat.fig diff --git a/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-1_pFBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-1_pFBA.jpg new file mode 100755 index 0000000..bab7e71 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-1_pFBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-2.jpg b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-2.jpg new file mode 100755 index 0000000..34a84c6 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-2.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-6.jpg b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-6.jpg new file mode 100755 index 0000000..72c0527 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/MethylFluxHist_1E-6.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA.jpg new file mode 100755 index 0000000..8b80c26 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA_magn.jpg b/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA_magn.jpg new file mode 100755 index 0000000..83b6f07 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/bix_flux_1E-1_pFBA_magn.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/methylFluxHist_1E-4_FBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/methylFluxHist_1E-4_FBA.jpg new file mode 100755 index 0000000..8e95601 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/methylFluxHist_1E-4_FBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA.jpg new file mode 100755 index 0000000..77779d0 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA_magn.jpg b/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA_magn.jpg new file mode 100755 index 0000000..bb551ec Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/qw-bi_flux_1E-1_pFBA_magn.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_g_1E-1_pFBA_magn.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_g_1E-1_pFBA_magn.jpg new file mode 100755 index 0000000..8794e57 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_g_1E-1_pFBA_magn.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth)1E-1_pFBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth)1E-1_pFBA.jpg new file mode 100755 index 0000000..d51cd2b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth)1E-1_pFBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-1_pFBA_magn.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-1_pFBA_magn.jpg new file mode 100755 index 0000000..1010fcd Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-1_pFBA_magn.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2.jpg new file mode 100755 index 0000000..32e1265 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2_pFBA_magn.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2_pFBA_magn.jpg new file mode 100755 index 0000000..fe90707 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-2_pFBA_magn.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-3.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-3.jpg new file mode 100755 index 0000000..0e78d91 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-3.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-4_FBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-4_FBA.jpg new file mode 100755 index 0000000..ad911b1 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-4_FBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-6_FBA.jpg b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-6_FBA.jpg new file mode 100755 index 0000000..ad911b1 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/SignifP/unc0638_growth_1E-6_FBA.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.fig b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.fig new file mode 100755 index 0000000..8818a6a Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.fig differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.jpg new file mode 100755 index 0000000..3eaf89f Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/MeFluxDistrib.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/bix-01294.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/bix-01294.jpg new file mode 100755 index 0000000..e2795de Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/bix-01294.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/brd-a02303741.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/brd-a02303741.jpg new file mode 100755 index 0000000..f655e8b Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/brd-a02303741.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/cbb-1007.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/cbb-1007.jpg new file mode 100755 index 0000000..7701a91 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/cbb-1007.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/gsk-j4.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/gsk-j4.jpg new file mode 100755 index 0000000..8351b5a Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/gsk-j4.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/methylstat.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/methylstat.jpg new file mode 100755 index 0000000..a4f5593 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/methylstat.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/qw-bi-011.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/qw-bi-011.jpg new file mode 100755 index 0000000..fe1e600 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/qw-bi-011.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0321.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0321.jpg new file mode 100755 index 0000000..1f13da5 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0321.jpg differ diff --git a/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0638.jpg b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0638.jpg new file mode 100755 index 0000000..457c4b5 Binary files /dev/null and b/matlab/MeCorr/Figures-LF/boxplot_Ludlow2015_min2SlcRxns/unc0638.jpg differ diff --git a/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-1.jpg b/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-1.jpg new file mode 100755 index 0000000..a26e05d Binary files /dev/null and b/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-1.jpg differ diff --git a/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-3.jpg b/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-3.jpg new file mode 100755 index 0000000..b36d23e Binary files /dev/null and b/matlab/MeCorr/Figures-LF/grateVarCorrected_pFBA/MethylFluxHist_1E-3.jpg differ diff --git a/matlab/MeCorr/MATLAB_CODE.asv b/matlab/MeCorr/MATLAB_CODE.asv deleted file mode 100755 index d521983..0000000 --- a/matlab/MeCorr/MATLAB_CODE.asv +++ /dev/null @@ -1,491 +0,0 @@ - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % this code reproduces all the main text figures from the manuscript (Shen et al) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of nutrient sources on acetylation - figure 2A -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -load recon1 -%load supplementary_software_code acetylation_model %contains metabolic model with nuclear acetylation reaction -load supplementary_software_code labels media_exchange1 mediareactions1 %list of nutrient conditions and uptake rates -posgluc = 1385; % glucose uptake reaction in recon1. -% LF changed reaction to find - rxnpos = find(ismember(metabolicmodel.rxns,'LYSMTF1n')); % rxnpos = 2451 -objpos = find(metabolicmodel.c); %biomass objective -minfluxflag = 0; % no PFBA -%epsilon_acetylation = 1E-3; -epsilon_methylation = 1E-2; % or 1E-1 - - for kappatype = 1:2 - if kappatype == 1, kappa = 10; else kappa = 0.01; end - %kappatype=1 means high [nutrient]. kappytype=2 means low - for i = 1:50 - kappa1 = kappa; - if (kappatype == 2) && (ismember(i,[2,3,5:19])) % trace elements - kappa1 = kappa/100; - elseif (kappatype == 1) && (ismember(i,[1;4])) % glucose or glutamine - kappa1 = 3; - end - model2 = metabolicmodel; - % change media.. - [ix, pos] = ismember(mediareactions1(i), model2.rxns); - model2.lb(pos) = -media_exchange1(i,1)*kappa1; - - [solf.x,sol11] = constrain_flux_regulation(model2,[],[],0,0,0,[],[],minfluxflag); - - str = ['media_change_growth_',num2str(kappatype),'(i,1) = solf.x(objpos);']; - if ~isempty(solf.x) && all(~isnan(solf.x)) - eval(str) - end - - j = 1; - model3 = model2; - model3.c(rxnpos) = epsilon_methylation; - [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); - str = ['media_change_histone_acet_nuc_',num2str(kappatype),'(i,j) = solf.x(rxnpos);']; - if ~isempty(solf.x) && all(~isnan(solf.x)) - eval(str) - end - disp(i) - end - - - disp(kappatype) - end - - labels(2) = {'Glutathione'}; - idx = [1:4,20:50]; - figure; - bar([media_change_histone_acet_nuc_1(idx,1) media_change_histone_acet_nuc_2(idx,1) ],1,'edgecolor','w'); - title('Acetylation levels in different growth conditions','fontweight','bold') - set(gca,'xtick',[1:length(mediareactions1(idx))],'xticklabel',labels(idx),'fontsize',8,'fontweight','bold','XTickLabelRotation',45) - set(gca,'TickDir', 'out') - set(gca,'box','off') - set(gca,'linewidth',2) - set(gcf,'color','white') - set(gca,'fontsize',12) - ylabel('Acetyl- Flux') -h = legend({'Excess','Depletion'}) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of gene deletion on acetylation - figure 2B -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - unqgenes = unique(acetylation_model.genes); %all genes in the model - - for i = 2:length(unqgenes) - - modeltemp = acetylation_model ; - modeltemp = deleteModelGenes(modeltemp,unqgenes(i)); - model3 = modeltemp; - - - [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 - end - - model3.c(rxnpos) = epsilon_acetylation; % epsilon is a weight - [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); - if ~isempty(solf.x) && ~isnan(solf.x) - acet_screen_rpmi(i,1) = solf.x(rxnpos); % impact on acetylation flux - end - - disp(i) - end - - acet_screen_rpmi(1,:) = NaN; - - [sx spos] = sort(acet_screen_rpmi(:,1)); - acetgenessorted = unqgenes(spos); - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of culture media components on acetylation - figure 2C -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -glucflag = logical([1 0 1 1 1 1 0 0 0 0]); -aaflag = logical([1 0 1 0 1 0 0 1 0 1]); -pyrflag = logical([1 0 1 1 0 0 0 0 1 1]); -glnflag = logical([1 0 0 1 1 0 1 0 0 1]); -expval = [1 0 1 1 1 1 1 0 1 1]; - -% part 2 -glucflag1 = [ ones(1,6), 0 ,0, 0, 0]; -expval1 = [ ones(1,7), 0 , 1, 0]; - - -[ix, aapos] = ismember( mediareactions1(20:37),acetylation_model.rxns); -posgluc = 1385; % glucose uptake reaction in recon1. -glnpos = 1386; % glutamine -pyrpos = 1509; % pyruvate -acetatepos = 1238; % acetate -fattyacidpos = 1445; % linoleic acid - -model3 = acetylation_model; -model3.c(rxnpos) = epsilon_acetylation ; -[solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); -wild_type_acet = solf.x(rxnpos); % default acetylation flux - -for i = 1:10 - model3 = acetylation_model; - if ~glucflag(i) - model3.lb(posgluc) = 0; - end - - if ~aaflag(i) - model3.lb(aapos) = 0; - end - - if ~glnflag(i) - model3.lb(glnpos) = 0; - end - - if pyrflag(i) - model3.lb(pyrpos) = -10; - end - - - model3.c(rxnpos) = epsilon_acetylation; - [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); - acet_media_screen_dmem(i,1) = solf.x(rxnpos); % impact on acetylation flux - - -end -disp('consistency with 10 different experimental conditions - Part I') -sum((acet_media_screen_dmem > 0.05*wild_type_acet) == expval') % 9 - - - -for i = 1:10 - model3 = acetylation_model; - if ~glucflag1(i) - model3.lb(posgluc) = 0; - end - - - switch i - case 2 % no ca+ - model3.lb(ismember( {'EX_ca2(e)'},acetylation_model.rxns)) = 0; - case 4 % no Phosphate - model3.lb(ismember( {'EX_pi(e)'},acetylation_model.rxns)) = 0; - case 6 % no vitamins in dmem.. compmosition from sigma - model3.lb(ismember( {'EX_thm(e)';'EX_ribflv(e)';'EX_pydx(e)';'EX_ncam(e)';'EX_inost(e)';'EX_5mthf(e)';'EX_pnto_R(e)';'EX_chol(e)'},acetylation_model.rxns)) = 0; - case 7 % add acetate - model3.lb(acetatepos) = -20; - case 9 % add fatty acid - model3.lb(fattyacidpos) = -5; - end - - - model3.c(rxnpos) = epsilon_acetylation; - [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); - acet_media_screen_dmem1(i,1) = solf.x(rxnpos); % impact on acetylation flux -end - -disp('consistency with 10 different experimental conditions - Part II') -sum((acet_media_screen_dmem1 > 0.05*wild_type_acet) == expval1') % 10 - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of basal metabolic state of cell lines on acetylation -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz %contains CCLE cell line names, gene expression data (z-transformed) -load supplementary_software_code acetlevlistmedia acetlevellist acetlevellistval %contains cell line names, growth media , total bulk acetylation -MODE = 1; % reaction (1) or gene list (0) -epsilon = 1E-2; rho = 1; -kappa = 1; -minfluxflag = 0; % no PFBA - - for i = 1:14 - % match cell line in CCLE data - iii = find(ismember(celllinenames_ccle1, acetlevellist(i))); - if ~isempty(iii) - iii = iii(1); - model2 = acetylation_model; - %find up and down-regulated genes in each cell line - ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); - offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); - - % set the glucose uptake based on media - % default glucose is -5 for rpmi - if ismember({'RPMI'} , acetlevlistmedia(i)) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi - elseif ismember({'DMEM'} , acetlevlistmedia(i)) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*4.5/2;% dmem.. - elseif ismember({'L15'} , acetlevlistmedia(i)) % NO glucose.. LOW Galactose - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -0;% L15 - model2.lb(find(ismember(model2.rxns, {'EX_gal(e)'}))) = -0.9;% - elseif ismember({'McCoy 5A'} , acetlevlistmedia(i)) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3/2;% mccoy - elseif ismember({'IMM'} , acetlevlistmedia(i)) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*4.5/2;% IMDM - - end - - %find reactions from differentially expressed genes - [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); - [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); - - disp(i) - - [fluxstate_gurobi,grate_ccle_exp_acetdat(i,1), solverobj_ccle(i,1)] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[], minfluxflag); % impact on growth - - model2.c(rxnpos) = epsilon_acetylation; - [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[],minfluxflag); - grate_ccle_exp_acetdat(i,2) = fluxstate_gurobi(rxnpos); %acetylation flux - - end - end - - figure; - plot(grate_ccle_exp_acetdat(:,2), acetlevellistval(1,:)','o','markerfacecolor',[ 0.9020 0.3804 0.0039],'markeredgecolor','k') -fb = polyfit(grate_ccle_exp_acetdat(:,2), acetlevellistval(1,:)',1); -fb1 = polyval(fb,grate_ccle_exp_acetdat(:,2)); -hold on; plot(grate_ccle_exp_acetdat(:,2),fb1,'r-','linewidth',2); -text(grate_ccle_exp_acetdat(:,2) + 0.2,acetlevellistval(1,:)',acetlevellist)%,'horizontalalignment','right','fontsize',10,'fontname','helvetica','fontweight','bold') -xlabel('Acetylation flux');ylabel('Bulk H3K9 Acetylation') - - [acetlevelcorr acetlevelcorrpv ] = corr(grate_ccle_exp_acetdat(:,2), acetlevellistval(1,:)') %correlation with h3k9 acetylation - % [acetlevelcorr acetlevelcorrpv ] = corr(grate_ccle_exp_acetdat(:,2), acetlevellistval(1,:)','type','spearman') % - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of nutrient environment of cell lines on sensitivity to deacetylase inhibitor - vorinostat -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - load supplementary_software_code recon1biologpm1match biologratio biolognewpm196 %contains BIOLOG phenotype array data - -for i = 1:96 - model2 = acetylation_model; - if ~isnan(recon1biologpm1match(i)) - model2.lb(posgluc) = -0.1; - model2.lb(recon1biologpm1match(i)) = -10; - - [solf.x,sol11] = constrain_flux_regulation(model2,[],[],0,0,0,[],[],minfluxflag); - disp(i) - - growth_basal2(i,1) = solf.x(objpos); % impact on growth - model2.c(rxnpos) = epsilon_acetylation; - [solf.x,sol11] = constrain_flux_regulation(model2,[],[],0,0,0,[],[],minfluxflag); - growth_basal2(i,2) = solf.x(rxnpos); %acetylation flux - - else - growth_basal2(i,:) = NaN; - end -end - -ixxn2 = ~isnan( growth_basal2(:,1)) & biolognewpm196(:,1) > 2; sum(ixxn2) % -xx = biolognewpm196(ixxn2,2)./biolognewpm196(ixxn2,1); -tel = [2:3,5,7:19]; -yy = [growth_basal2(ixxn2,:)]; -[hh pp] = corr([yy(tel,1) xx(tel) ],'type','spearman') % . -0.36 %correlation between predicted growth and vorinostat sensitivity -[hh pp] = corr([yy(tel,2) xx(tel) ],'type','spearman') % -0.67 %correlation between predicted acetylation flux and vorinostat sensitivity - - -figure; plot(xx(tel), yy(tel,2) ,'o','markerfacecolor',[ 0.9020 0.3804 0.0039],'markeredgecolor','k') % 0.5 -xlim([0.5 2]) -fb = polyfit(xx(tel), yy(tel,2),1); -fb1 = polyval(fb,xx(tel)); -hold on; plot(xx(tel),fb1,'r-','linewidth',2); -ylabel('Acetylation flux');xlabel('Vorinostat treatment vs control ratio from Biolog arrays') - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% impact of basal metabolic state of CCLE cell lines on sensitivity to demethylase inhibitors -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -load supplementary_software_code ctd2celllineidname exptidcelllinemediamatch ctd2celllineidname_id r* %data from seashore-ludlow study, -...contains cell line names , growth media -load supplementary_software_code ctd2compoundidname_id drug_auc_expt ctd2compoundidname_name %data from seashore-ludlow study, -...contains drug names , drug sensitivity data -load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz -%contains CCLE cell line names, gene expression data (z-transformed) -load supplementary_software_code hcommon_exptdat hcommon1 %data from seashore-ludlow study, -...contains drug names , drug sensitivity data for cell lines that were screened against all 4 hdac inhibitors -load supplementary_software_code hdacexpfcs hdacexpallgeneids %contains -...gene expression data after treatment with hdac inhibitors -rxnpos = 2451; % loading the above variables changes rxnpos to 3754, which is out of bounds for fluxstate_gurobi - -MODE = 1; % reaction (1) or gene list (0) -epsilon = 1E-2; rho = 1; -kappa = 1; % parameters for integrating transcriptomics data. kappa is the strength of down regulation of genes (Chandrasekaran & Price, PNAS, 2010) -minfluxflag = 0; % no Pfba -hdactransint = 1; -basalflag = 1; - -grate_ccle_exp_soft = NaN(length(exptidcelllinemediamatch),2); %contains acetylation flux based on basal metabolic state -grate_ccle_exp_soft_hdacsign = NaN(length(exptidcelllinemediamatch),8); %contains acetylation flux based on basal metabolic state and impact of each hdac inhibitor teatment - -for i = 1:length(exptidcelllinemediamatch) - % match cell line data in CCLE with CTD2 - ii = find(ismember(ctd2celllineidname_id, exptidcelllinemediamatch(i,2))); - iii = find(ismember(celllinenames_ccle1, ctd2celllineidname(ii,1))); - if ~isempty(iii) - iii = iii(1); - model2 = metabolicmodel; - - %find up and down-regulated genes in each cell line - ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); - offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % set the glucose uptake based on media - % default glucose is -5 for rpmi - if r1(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi - elseif r2(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% dmem with low glucose - elseif r3(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% Emem.. - elseif r4(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3/2;% mccoy - elseif r5(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% mem - elseif r6(i) - model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3.15/2;% dmem:f12 - end - - %find reactions from differentially expressed genes - [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); - [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); -% - - disp(i) - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % get basal metabolic state based on transcriptome - if basalflag - [fluxstate_gurobi,grate_ccle_exp_soft(i,1), solverobj_ccle(i,1)] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth - model2.c(rxnpos) = epsilon_methylation ; - [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); - grate_ccle_exp_soft(i,2) = fluxstate_gurobi(rxnpos); %acetylation flux - end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % hdac inhibitor impact on transcriptome - if hdactransint - for kk = 4:-1:1 % match drug order - ongenesh = hdacexpallgeneids((hdacexpfcs(:,kk) > 1.1) & (hdacexpfcs(:,kk + 4) < 0.01)); - ongenesh = intersect(ongenesh, model2.genes); - offgenesh = hdacexpallgeneids((hdacexpfcs(:,kk) < 0.9) & (hdacexpfcs(:,kk + 4) < 0.01)); - offgenesh = intersect(offgenesh, model2.genes); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - [~,~,onreactionsh,~] = deleteModelGenes(model2, ongenesh); - [~,~,offreactionsh,~] = deleteModelGenes(model2, offgenesh); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % remove conflicting transcripts - onreactionsh0 = setdiff(onreactionsh, offreactions); - offreactionsh0 = setdiff(offreactionsh, onreactions); - onreactions0 = setdiff(onreactions, offreactionsh); - offreactions0 = setdiff(offreactions, onreactionsh); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - onreactions1h = [onreactions0; onreactionsh0]; - offreactions1h = [offreactions0; offreactionsh0]; - epsilon2 = [repmat(0, size(offreactions1h))]; - model2.c(rxnpos) = 0; - % [fluxstate_gurobi,grate_ccle_exp_soft_hdacsign(i,kk)] = constrain_flux_regulation(model2,onreactions1h,offreactions1h,kappa,rho,epsilon,MODE, epsilon2,minfluxflag); % impact on growth - model2.c(rxnpos) = epsilon_methylation ; - [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions1h,offreactions1h,kappa,rho,epsilon,MODE,epsilon2,minfluxflag); - grate_ccle_exp_soft_hdacsign(i,kk + 4) = fluxstate_gurobi(rxnpos); %acetylation flux - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - end - end - - else - grate_ccle_exp_soft(i,:) = NaN; - grate_ccle_exp_soft_hdacsign(i,:) = NaN; - end -end - - -figure; h = histogram(grate_ccle_exp_soft(:,2),70); - xlabel('Predicted acetylation flux') - ylabel('Total cell lines') -title('Distribution of acetylation flux among CCLE cell lines','fontweight','normal') - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%predicting sensitivity to hdac inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -hdaclist = {'LBH-589','vorinostat','entinostat','belinostat'} -%hdmelist = - -for j = 1:4 -fx = find(ismember(ctd2compoundidname_name, hdaclist(j))) -ix = ismember(drug_auc_expt(:,3), ctd2compoundidname_id(fx)); -sum(ix) % 847. ix is a 1D array. Each drug appears multiple times (appears for each medium in drug_auc_expt) -hdac_auc_dat = drug_auc_expt(ix,:); - -[ix pos] = ismember(hdac_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix) % -v1 = hdac_auc_dat(:,2); -v2 = grate_ccle_exp_soft(pos,:); - -ix0 = (v2(:,2) < 0.05);sum(ix0) % 0.05 is the default threshold value -ix01 = (v2(:,2) > 0.05);sum(ix01) -[hh pp_basal(j,3)] = ttest2(v1(ix0), v1(ix01)); % - - groups = NaN(size(ix0)); - groups(ix0) = 1; groups(ix01) = 2; - vv = NaN(length(v1), 2); - vv(1:sum(groups == 1),1) = v1(groups == 1); - vv(1:sum(groups == 2),2) = v1(groups == 2); - - figure; - %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); - %hold on; - bh = boxplot(v1, groups,'symbol','') - set(gca,'xticklabel',{'Low flux','High flux'}) -ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); -xlabel(hdaclist(j)); -ylim([0 20]) - -ix0 = (v2(:,2) <= prctile(v2(:,2), 25)); sum(ix0) -ix01 = (v2(:,2) > prctile(v2(:,2), 25)); sum(ix01) -[hh pp_basal(j,1)] = ttest2(v1(ix0), v1(ix01)) ;% - -ix0 = (v2(:,2) <= prctile(v2(:,2), 50)); sum(ix0) -ix01 = (v2(:,2) > prctile(v2(:,2), 50)); sum(ix01) -[hh pp_basal(j,2)] = ttest2(v1(ix0), v1(ix01)) ;% - -end -disp(pp_basal) %t-test p-values for each drug - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%predicting variation in sensitivity between hdac inhibitors based on basal -%metabolic state and drug impact - comparison with drug sensitivity data from seashore-ludlow study -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - g1 = grate_ccle_exp_soft_hdacsign(:, 5:8); -g1 = g1 - repmat(ignoreNaN(g1, @median,2),1,4); % - -[ix pos] = ismember(hcommon1, exptidcelllinemediamatch(:,1)); sum(ix) % match cell lines -v1 = hcommon_exptdat; -v1 = v1 - repmat(median(v1,2),1,4); -v2 = g1(pos,:); -v1 = v1(:); - - -ix0 = (v2(:) <= prctile(v2(:), 25)); sum(ix0) -ix01 = (v2(:) >= prctile(v2(:), 25)); sum(ix01) -[hh pp] = ttest2(v1(ix0), v1(ix01)) % 2e-52.. - -ix0 = (v2(:) <= prctile(v2(:), 50)); sum(ix0) -ix01 = (v2(:) >= prctile(v2(:), 50)); sum(ix01) -[hh pp] = ttest2(v1(ix0), v1(ix01)) % 3e-44 - -ix0 = (v2(:) <= 0.05); sum(ix0) -ix01 = (v2(:) >0.05); sum(ix01) -[hh pp] = ttest2(v1(ix0), v1(ix01)) %2e -62 -mean(v1(ix0)) % -0.01 -mean(v1(ix01)) % -2.61 - - groups = NaN(size(ix0)); - groups(ix0) = 1; groups(ix01) = 2; -vv = NaN(length(v1), 2); - - vv(1:sum(groups == 1),1) = v1(groups == 1); - vv(1:sum(groups == 2),2) = v1(groups == 2); -figure; -bh = boxplot(vv,'symbol','' ,'orientation','horizontal') - set(gca,'yticklabel',{'No','Large'}) - xlabel('Observed differential sensitivity between drugs in a cell line') -title({'Predicted differential sensitivity between drugs' 'vs Observed differential sensitivity (AUC)'}) - - -figure -h1 = histfit(vv(:,1),20,'kernel')%, -set(h1(1),'facecolor','g','facealpha',.15,'edgecolor','none') -set(h1(2),'color','g') -hold on -h2 = histfit(vv(:,2),20,'kernel') -set(h2(1),'facecolor','m','facealpha',.15,'edgecolor','none') -set(h2(2),'color','m') -hh = legend([h1(2),h2(2)],'No difference','Large difference') -title(hh,'Predicted differential acetylation') diff --git a/matlab/MeCorr/MATLAB_CODE_allrxns.asv b/matlab/MeCorr/MATLAB_CODE_allrxns.asv new file mode 100644 index 0000000..b3a6588 --- /dev/null +++ b/matlab/MeCorr/MATLAB_CODE_allrxns.asv @@ -0,0 +1,404 @@ + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % this code reproduces all the main text figures from the manuscript (Shen et al) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Bulk methylation model: +% 1) cd ./../scripts. Run make_eGEM to create bulk methylation model +% 2) cd ./../MeCorr. Run section 1 and 3 of this script. Only run Section 2 if you want to +% run the long for-loop +% 3) Parameters of interest to vary: epsilon_methylation, model2 +%clearvars -except min_model % If don't want to run make_eGEM +model2 = min_model; +epsilon_methylation = 1E-1; +%rxnpos = find(ismember(min_model.rxns,'LYSMTF1n')); +minfluxflag = 0; % 0: no Pfba, 1: Pfba +% 4) Change name of files to which variables are saved. End of section 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % impact of basal metabolic state of CCLE cell lines on sensitivity to demethylase inhibitors +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load supplementary_software_code ctd2celllineidname exptidcelllinemediamatch ctd2celllineidname_id r* +%data from seashore-ludlow study, contains cell line names, growth media +load supplementary_software_code ctd2compoundidname_id drug_auc_expt ctd2compoundidname_name +%data from seashore-ludlow study, contains drug names, drug sensitivity data +load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz +%contains CCLE cell line names, gene expression data (z-transformed) +load supplementary_software_code hcommon_exptdat hcommon1 +%data from seashore-ludlow study, contains drug names, drug sensitivity data for cell lines that were screened against all 8 hme inhibitors +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 +% caompatible with function ismember +ctd2celllineidname= array2table(ctd2celllineidname); +ctd2celllineidname.Properties.VariableNames{'ctd2celllineidname1'}= 'ccl_name'; +ctd2celllineidname_id= array2table(ctd2celllineidname_id); +ctd2celllineidname_id.Properties.VariableNames{'ctd2celllineidname_id'}= 'master_ccl_id'; +ctd2compoundidname_name= array2table(ctd2compoundidname_name); +ctd2compoundidname_name.Properties.VariableNames{'ctd2compoundidname_name'}= 'cpd_name'; +ctd2compoundidname_id= array2table(ctd2compoundidname_id); +ctd2compoundidname_id.Properties.VariableNames{'ctd2compoundidname_id'}= 'cpd_id'; +celllinenames_ccle1= array2table(celllinenames_ccle1); +celllinenames_ccle1.Properties.VariableNames{'celllinenames_ccle1'}= 'ccl_name'; +drug_auc_expt= array2table(drug_auc_expt); +drug_auc_expt.Properties.VariableNames{'drug_auc_expt3'}= 'cpd_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt1'}= 'experiment_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt2'}= 'auc'; +exptidcelllinemediamatch=exptidcelllinemediamatch(2:end, :); +exptidcelllinemediamatch= array2table(exptidcelllinemediamatch); %s4 +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch1'}= 'experiment_id'; +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch2'}= 'master_ccl_id'; +hcommon1= array2table(hcommon1); +hcommon1.Properties.VariableNames{'hcommon1'}= 'experiment_id'; + +MODE = 1; % reaction (1) or gene list (0) +epsilon = 1E-2; rho = 1; +kappa = 1; % parameters for integrating transcriptomics data. kappa is the strength of down regulation of genes (Chandrasekaran & Price, PNAS, 2010) +hmetransint = 1; +basalflag = 1; + +hmei_list = {'BRD-A02303741';'BIX-01294';'methylstat';'QW-BI-011';... + 'UNC0321';'CBB-1007';'UNC0638';'GSK-J4'}; +hmei_list= cell2table(hmei_list); +hmei_list.Properties.VariableNames{'hmei_list'}='cpd_name'; +cpd_info= readtable('DataRees2016/nchembio.1986-S3_cpd.xlsx'); +cpd_infoMe= cpd_info(ismember(cpd_info(:, 1), hmei_list), [1,5,6,7,8]); +cpd_infoMe= table2cell(cpd_infoMe); + +% If I don't want to run next for-loop: +%load('VariablesSaved\fluxstate_gurobi'); +%load('VariablesSaved\grate_ccle_exp_soft'); +%% long for-loops +weights= [1E-1, 1E-2, 1E-3, 1E-4, 1E-5, 1E-6]; +for i_weight = 2:length(weights) + epsilon_methylation= weights(i_weight); + grate_ccle_exp_soft = NaN(height(exptidcelllinemediamatch),3); %contains acetylation flux based on basal metabolic state + fluxes_allrxns= NaN(height(exptidcelllinemediamatch),length(model2.rxns)); + g_rate= NaN(height(exptidcelllinemediamatch), length(model2.rxns)); + + for i = 1:height(exptidcelllinemediamatch) + % match cell line data in CCLE with CTD2 + ii = find(ismember(ctd2celllineidname_id, exptidcelllinemediamatch(i,2))); + iii = find(ismember(celllinenames_ccle1, ctd2celllineidname(ii,1))); + if ~isempty(iii) + iii = iii(1); + % model2 = min_model; + % model2 = metabolicmodel; + + %find up and down-regulated genes in each cell line + ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); + offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % set the glucose uptake based on media + % default glucose is -5 for rpmi + if r1(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi + elseif r2(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% dmem with low glucose + elseif r3(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% Emem.. + elseif r4(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3/2;% mccoy + elseif r5(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% mem + elseif r6(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3.15/2;% dmem:f12 + end + + %find reactions from differentially expressed genes + onreactions= findRxnsFromGenes(model2, ongenes); + onreactions= struct2cell(onreactions); + for jj=1:length(onreactions) + onreactions(jj)= onreactions{jj}(1); + end + offreactions= findRxnsFromGenes(model2, offgenes); + offreactions= struct2cell(offreactions); + for jj=1:length(offreactions) + offreactions(jj)= offreactions{jj}(1); + end + % Below 2 lines work for acetyl model, not min methyl model. + % [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); + % [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); + + disp(i) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % get basal metabolic state based on transcriptome + % if basalflag + % [fluxes, grate, solverobj_ccle] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth + % grate_ccle_exp_soft(i,1:2)= grate; % first 2 columns are basal met flux. + % ...3rd column will be met flux with a rxn maximized + % fluxes_allrxns(i,:) = fluxes; % correlate each row of fluxes_allreactions with auc + % model2.c(rxnpos) = epsilon_methylation ; + % [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + % grate_ccle_exp_soft(i,3) = fluxstate_gurobi(rxnpos); %methylation flux + % end + + for rxncount = 1:length(model2.rxns) + if basalflag + %minfluxflag = 1; + %[~, grate, solverobj_ccle] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth + model2.c(rxncount) = epsilon_methylation ; + [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + g_rate(i,rxncount) = fluxstate_gurobi(rxncount); %methylation flux + end + end + else + grate_ccle_exp_soft(i,:) = NaN; + end + end + % Save in local github repo and another local folder. wd is MeCorr + name_e= strcat('grate_', epsilon_methylation); + path1= strcat('VariablesSaved/', name_e); + save(path1, name_e); + path2= strcat('./../../../', name_e); + save(path2, name_e); +end + +% figure; h = histogram(grate_ccle_exp_soft(:,3),70); +% xlabel('Predicted methylation flux') +% ylabel('Total cell lines') +% title('Distribution of methylation flux among CCLE cell lines','fontweight','normal') + +% save('VariablesSaved\grate_3FBA_allRxns', 'g_rate'); +% save('VariablesSaved\fluxesAll_1FBA', 'fluxes_allrxns'); +% save('VariablesSaved\grate_1E-6', 'grate_ccle_exp_soft'); +% save('VariablesSaved\fluxstate_1E-6', 'fluxstate_gurobi'); +%% Correlation between flux and auc for a reaction +% flux_allrexns: 1031 cell lines by 3777 reactions +% drug_auc_expt: extract auc values for the 8 methyl drugs +% 8 x 3000 rhos +% fluxes_allrxns= g_rate; +rho= NaN(height(hmei_list),length(model2.rxns)); +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); + hmei_auc_dat = drug_auc_expt(ix,:); + + % Only uses the experiments of fluxes_allrxns that also are in hmei_auc_dat + [ix, pos] = ismember(hmei_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix); + v1 = hmei_auc_dat(:,2); v1= table2array(v1); + for nCol=1:size(fluxes_allrxns,2) + v4= fluxes_allrxns(pos,nCol); + [rho(j,nCol), rho_p(j,nCol)]= corr(v4, v1, 'rows', 'complete'); + end +end +disp('rho & p-value calculated') +%% Create Table of Significant Reactions by Correlation value +% Workflow: Change threshold. Change struct field name (e.g. above3) +% accordingly +sigRhoTF= (abs(rho) > 0.3); +nSigExpt= sum(sigRhoTF); % sum number of signif expts per rxn (sum each column) +sigIndRxn= (nSigExpt >= 1); +sigRxn= model2.rxns(sigIndRxn); disp(sigRxn) + +[iDrug,iRxn]= find(sigRhoTF); n= length(iDrug); disp(n) +sigRxnS.above3(1:n,1)= num2cell(rho(sigRhoTF)); +sigRxnS.above3(1:n,2)= num2cell(rho_p(sigRhoTF)); +sigRxnS.above3(1:n,3)= model2.rxns(iRxn); +sigRxnS.above3(1:n,4)= model2.rxnNames(iRxn); +sigRxnS.above3(1:n,5)= model2.subSystems(iRxn); +sigRxnS.above3(1:n,6)= cpd_infoMe(iDrug, 1); +sigRxnS.above3(1:n,7)= cpd_infoMe(iDrug, 3); +sigRxnS.above3(1:n,8)= cpd_infoMe(iDrug, 2); + +sigRxnT= sortrows(sigRxnS.above3, 1); % sort by correlation +sigRxnT= cell2table(sigRxnT); +sigRxnT.Properties.VariableNames{'sigRxnT1'}='Correlation'; +sigRxnT.Properties.VariableNames{'sigRxnT2'}='Pvalue'; +sigRxnT.Properties.VariableNames{'sigRxnT3'}='Rxn'; +sigRxnT.Properties.VariableNames{'sigRxnT4'}='RxnName'; +sigRxnT.Properties.VariableNames{'sigRxnT5'}='Subsystem'; +sigRxnT.Properties.VariableNames{'sigRxnT6'}='Compound'; +sigRxnT.Properties.VariableNames{'sigRxnT7'}='CpdActivity'; +sigRxnT.Properties.VariableNames{'sigRxnT8'}='CpdGeneTarget'; + +% Rename variable and save +% sigRxnT_1pFBA_0717=sigRxnT; +% save('sigRxnT_1pFBA_0717', 'sigRxnT_1pFBA_0717'); +% writetable(sigRxnT,('sigRxnT_1pFBA_0717.xlsx')); + +% sigRxnT_maxAllRxns_1FBA=sigRxnT; +% save('sigRxnT_maxAllRxns_1FBA', 'sigRxnT_maxAllRxns_1FBA'); +% writetable(sigRxnT_maxAllRxns_1FBA,('sigRxnT_maxAllRxns_1FBA.xlsx')); +%% Significat reactions by correlation or p-value +sigTF= (abs(rho) > 0.3); +% sigTF= (abs(rho_p) < 0.001); +nSigExpt= sum(sigTF); % sum number of signif expts per rxn (sum each column) +sigIndRxn= (nSigExpt >= 1); +sigRxn= model2.rxns(sigIndRxn); disp(sigRxn) + +[iDrug,iRxn]= find(sigTF); n= length(iDrug); disp(n) +sigRxnS.above3(1:n,1)= num2cell(rho(sigTF)); +sigRxnS.above3(1:n,2)= num2cell(rho_p(sigTF)); +sigRxnS.above3(1:n,3)= model2.rxns(iRxn); +sigRxnS.above3(1:n,4)= model2.rxnNames(iRxn); +sigRxnS.above3(1:n,5)= model2.subSystems(iRxn); +sigRxnS.above3(1:n,6)= cpd_infoMe(iDrug, 1); +sigRxnS.above3(1:n,7)= cpd_infoMe(iDrug, 3); +sigRxnS.above3(1:n,8)= cpd_infoMe(iDrug, 2); + +sigRxnT= sortrows(sigRxnS.above3, 1); % sort by correlation +sigRxnT= cell2table(sigRxnT); +sigRxnT.Properties.VariableNames{'sigRxnT1'}='Correlation'; +sigRxnT.Properties.VariableNames{'sigRxnT2'}='Pvalue'; +sigRxnT.Properties.VariableNames{'sigRxnT3'}='Rxn'; +sigRxnT.Properties.VariableNames{'sigRxnT4'}='RxnName'; +sigRxnT.Properties.VariableNames{'sigRxnT5'}='Subsystem'; +sigRxnT.Properties.VariableNames{'sigRxnT6'}='Compound'; +sigRxnT.Properties.VariableNames{'sigRxnT7'}='CpdActivity'; +sigRxnT.Properties.VariableNames{'sigRxnT8'}='CpdGeneTarget'; +% Rename variable and save +% sigRxnT_1pFBA= sigRxnT; +% save('sigRxnT_1pFBA', 'sigRxnT_1pFBA'); +% writetable(sigRxnT,('sigRxnT_1pFBA.xlsx')); + +% sigRxnT_maxAllRxns_1FBA= sigRxnT; +% save('sigRxnT_maxAllRxns_1FBA', 'sigRxnT_maxAllRxns_1FBA'); +% writetable(sigRxnT_maxAllRxns_1FBA,('sigRxnT_maxAllRxns_1FBA.xlsx')); + +% sigPrxnT_1FBA= sigRxnT; +%% Calculate correlation between flux and auc. Each reaction maximized. +rho3= NaN(height(hmei_list),length(model2.rxns)); +for j= 1:height(hmei_list) + % Extract the auc for a drug in drug_auc_expt + fx = find(ismember(ctd2compoundidname_name, hmei_list(j,1))); + ix = ismember(drug_auc_expt(:,3), ctd2compoundidname_id(fx,1)); + hmei_auc_dat = drug_auc_expt(ix,:); + % Only uses the experiments of fluxes_allrxns that also are in hmei_auc_dat + [ix, pos] = ismember(hmei_auc_dat(:,1), exptidcelllinemediamatch(:,1)); + v1 = hmei_auc_dat(:,2); + v1= table2array(v1); + + for nCol= 1:size(g_rate,2) + v5= g_rate(pos, nCol); + rho3(j,nCol)= corr(v5, v1, 'rows', 'complete'); + end +end + +sigRho3Ind= false(1,length(model2.rxns)); +for j2= 13:13%length(model2.rxns) + ix= ~isnan(rho3(:,j2)); + sigRho3Ind(1,j2)= (sum(ix)>0); +end +sigRxn3= model2.rxns(sigRho3Ind); +disp(sigRxn3) + +save('VariablesSaved\rho3_1E-1', rho3); +save('VariablesSabed\sigRxn3_1E-1', sigRxn3); +%% predicting sensitivity to hme inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%load('grate_pFBA_1E-2'); load('fluxstate_pFBA_1E-2'); +clearvars r_basal +r_basal(height(hmei_list))= struct('rho',0, 'p',0, 'rlowerbound',0, 'rupperbound',0); +pp_flux= zeros(height(hmei_list), 3); +pp_grate= zeros(height(hmei_list), 3); +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) % 847. ix is a 1D array. Each drug appears multiple times (appears for each medium in drug_auc_expt) + hme_auc_dat = drug_auc_expt(ix,:); + + [ix, pos] = ismember(hme_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix); + v1 = hme_auc_dat(:,2); v1= table2array(v1); + v2 = grate_ccle_exp_soft(pos,:); % growth rate + v3 = fluxstate_gurobi(pos,:); + +% % Scatter plots of methylation (metabolic) flux & calculate correlation +% figure; +% scatter(v1, v3) +% ylabel({'Methylation Flux'},'fontname','helvetica'); %'fontweight','bold') +% xlabel({'Sensitivity (AUC)'; hmei_list.cpd_name{j}},'fontname','helvetica'); +% ylim([0, 0.06]) +% [R, P, RL, RU]= corrcoef(v1, v3); +% r_basal(j).rho=R; r_basal(j).p=P; r_basal(j).rlowerbound=RL; r_basal(j).rupperbound=RU; +% str=sprintf('r= %1.3f',r_basal(j).rho(1,2)); +% T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); +% set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + % Boxplot pairs of methylation flux (Actually graphing metabolic flux, + ...but we hypothesize met flux approximates methyl flux) + ix0 = (v3 < 0.05);sum(ix0); % 0.05 is the default threshold value + ix01 = (v3 > 0.05);sum(ix01); + [hh, pp_flux(j,3)] = ttest2(v1(ix0), v1(ix01)); % + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; + vv = NaN(length(v1), 2); + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); + + figure; + %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); + %hold on; + bh = boxplot(v1, groups,'symbol',''); + set(gca,'xticklabel',{'Low flux','High flux'}) + ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); + xlabel(hmei_list{j,1}); + %ylim([0 20]) + str=sprintf('p= %1.4f',pp_flux(j,3)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + ix0 = (v3 <= prctile(v3, 25)); sum(ix0); + ix01 = (v3 > prctile(v3, 25)); sum(ix01); + [hh, pp_flux(j,1)] = ttest2(v1(ix0), v1(ix01)); + + ix0 = (v3 <= prctile(v3, 50)); sum(ix0); + ix01 = (v3 > prctile(v3, 50)); sum(ix01); + [hh, pp_flux(j,2)] = ttest2(v1(ix0), v1(ix01)); + + % Boxplot pairs of growth rate %%%%%% + ix0 = (v2(:,2) < 0.05);sum(ix0); % 0.05 is the default threshold value + ix01 = (v2(:,2) > 0.05);sum(ix01); + [hh, pp_grate(j,3)] = ttest2(v1(ix0), v1(ix01)); % + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; + vv = NaN(length(v1), 2); + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); + + figure; + %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); + %hold on; + bh = boxplot(v1, groups,'symbol',''); + set(gca,'xticklabel',{'Low growth','High growth'}) + ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); + xlabel(hmei_list{j,1}); + %ylim([0 20]) + str=sprintf('p= %1.4f',pp_grate(j,3)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + ix0 = (v2(:,3) <= prctile(v2(:,3), 25)); sum(ix0); + ix01 = (v2(:,3) > prctile(v2(:,3), 25)); sum(ix01); + [hh, pp_grate(j,1)] = ttest2(v1(ix0), v1(ix01)); + + ix0 = (v2(:,3) <= prctile(v2(:,3), 50)); sum(ix0); + ix01 = (v2(:,3) > prctile(v2(:,3), 50)); sum(ix01); + [hh, pp_grate(j,2)] = ttest2(v1(ix0), v1(ix01)); + +% ix0 = (v2(:,2) <= prctile(v2(:,2), 25)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 25)); sum(ix01); +% [hh, pp_grate(j,1)] = ttest2(v1(ix0), v1(ix01)); +% +% ix0 = (v2(:,2) <= prctile(v2(:,2), 50)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 50)); sum(ix01); +% [hh, pp_grate(j,2)] = ttest2(v1(ix0), v1(ix01)); +end +r_basal= struct2table(r_basal); +r_basal.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +disp(r_basal) +pp_flux= array2table(pp_flux); +pp_flux.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +pp_flux.Properties.VariableNames{'pp_flux3'}= 'pf_hilo'; +pp_flux.Properties.VariableNames{'pp_flux1'}= 'pf_25prctile'; +pp_flux.Properties.VariableNames{'pp_flux2'}= 'pf_50prctile'; +disp(pp_flux) %t-test p-values for each drug +pp_grate= array2table(pp_grate); +pp_grate.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +pp_grate.Properties.VariableNames{'pp_grate3'}= 'pg_hilo'; +pp_grate.Properties.VariableNames{'pp_grate1'}= 'pg_25prctile'; +pp_grate.Properties.VariableNames{'pp_grate2'}= 'pg_50prctile'; +disp(pp_grate) +%save('VariablesSaved\pp_grate_1E-6_pFBA', 'pp_grate'); \ No newline at end of file diff --git a/matlab/MeCorr/MATLAB_CODE_allrxns.m b/matlab/MeCorr/MATLAB_CODE_allrxns.m new file mode 100755 index 0000000..7f6a618 --- /dev/null +++ b/matlab/MeCorr/MATLAB_CODE_allrxns.m @@ -0,0 +1,419 @@ + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % this code reproduces all the main text figures from the manuscript (Shen et al) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Bulk methylation model: +% 1) cd ./../scripts. Run make_eGEM to create bulk methylation model +% 2) cd ./../MeCorr. Run section 1 and 3 of this script. Only run Section 2 if you want to +% run the long for-loop +% 3) Parameters of interest to vary: epsilon_methylation, model2 +cd ./../scripts +make_eGEM +cd ./../MeCorr +%clearvars -except min_model % If don't want to run make_eGEM +model2 = min_model; +epsilon_methylation = 1E-1; +%rxnpos = find(ismember(min_model.rxns,'LYSMTF1n')); +minfluxflag = 1; % 0: no Pfba, 1: Pfba +% 4) Change name of files to which variables are saved. End of section 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % impact of basal metabolic state of CCLE cell lines on sensitivity to demethylase inhibitors +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load supplementary_software_code ctd2celllineidname exptidcelllinemediamatch ctd2celllineidname_id r* +%data from seashore-ludlow study, contains cell line names, growth media +load supplementary_software_code ctd2compoundidname_id drug_auc_expt ctd2compoundidname_name +%data from seashore-ludlow study, contains drug names, drug sensitivity data +load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz +%contains CCLE cell line names, gene expression data (z-transformed) +load supplementary_software_code hcommon_exptdat hcommon1 +%data from seashore-ludlow study, contains drug names, drug sensitivity data for cell lines that were screened against all 8 hme inhibitors +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 +% caompatible with function ismember +ctd2celllineidname= array2table(ctd2celllineidname); +ctd2celllineidname.Properties.VariableNames{'ctd2celllineidname1'}= 'ccl_name'; +ctd2celllineidname_id= array2table(ctd2celllineidname_id); +ctd2celllineidname_id.Properties.VariableNames{'ctd2celllineidname_id'}= 'master_ccl_id'; +ctd2compoundidname_name= array2table(ctd2compoundidname_name); +ctd2compoundidname_name.Properties.VariableNames{'ctd2compoundidname_name'}= 'cpd_name'; +ctd2compoundidname_id= array2table(ctd2compoundidname_id); +ctd2compoundidname_id.Properties.VariableNames{'ctd2compoundidname_id'}= 'cpd_id'; +celllinenames_ccle1= array2table(celllinenames_ccle1); +celllinenames_ccle1.Properties.VariableNames{'celllinenames_ccle1'}= 'ccl_name'; +drug_auc_expt= array2table(drug_auc_expt); +drug_auc_expt.Properties.VariableNames{'drug_auc_expt3'}= 'cpd_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt1'}= 'experiment_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt2'}= 'auc'; +exptidcelllinemediamatch=exptidcelllinemediamatch(2:end, :); +exptidcelllinemediamatch= array2table(exptidcelllinemediamatch); %s4 +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch1'}= 'experiment_id'; +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch2'}= 'master_ccl_id'; +hcommon1= array2table(hcommon1); +hcommon1.Properties.VariableNames{'hcommon1'}= 'experiment_id'; + +MODE = 1; % reaction (1) or gene list (0) +epsilon = 1E-2; rho = 1; +kappa = 1; % parameters for integrating transcriptomics data. kappa is the strength of down regulation of genes (Chandrasekaran & Price, PNAS, 2010) +hmetransint = 1; +basalflag = 1; + +hmei_list = {'BRD-A02303741';'BIX-01294';'methylstat';'QW-BI-011';... + 'UNC0321';'CBB-1007';'UNC0638';'GSK-J4'}; +hmei_list= cell2table(hmei_list); +hmei_list.Properties.VariableNames{'hmei_list'}='cpd_name'; +cpd_info= readtable('./../DataRees2016/nchembio.1986-S3_cpd.xlsx'); +cpd_infoMe= cpd_info(ismember(cpd_info(:, 1), hmei_list), [1,5,6,7,8]); +cpd_infoMe= table2cell(cpd_infoMe); + +% If I don't want to run next for-loop: +%load('VariablesSaved\fluxstate_gurobi'); +%load('VariablesSaved\grate_ccle_exp_soft'); +%% long for-loops +weights= [1E-1, 1E-2, 1E-3, 1E-4, 1E-5, 1E-6]; +name= [1, 2, 3, 4, 5, 6]; +for i_weight = 2:length(weights) + epsilon_methylation= weights(i_weight); + grate_ccle_exp_soft = NaN(height(exptidcelllinemediamatch),3); %contains acetylation flux based on basal metabolic state + fluxes_allrxns= NaN(height(exptidcelllinemediamatch),length(model2.rxns)); + g_rate= NaN(height(exptidcelllinemediamatch), length(model2.rxns)); + + for i = 1:height(exptidcelllinemediamatch) + % match cell line data in CCLE with CTD2 + ii = find(ismember(ctd2celllineidname_id, exptidcelllinemediamatch(i,2))); + iii = find(ismember(celllinenames_ccle1, ctd2celllineidname(ii,1))); + if ~isempty(iii) + iii = iii(1); + % model2 = min_model; + % model2 = metabolicmodel; + + %find up and down-regulated genes in each cell line + ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); + offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % set the glucose uptake based on media + % default glucose is -5 for rpmi + if r1(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi + elseif r2(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% dmem with low glucose + elseif r3(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% Emem.. + elseif r4(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3/2;% mccoy + elseif r5(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% mem + elseif r6(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3.15/2;% dmem:f12 + end + + %find reactions from differentially expressed genes + onreactions= findRxnsFromGenes(model2, ongenes); + onreactions= struct2cell(onreactions); + for jj=1:length(onreactions) + onreactions(jj)= onreactions{jj}(1); + end + offreactions= findRxnsFromGenes(model2, offgenes); + offreactions= struct2cell(offreactions); + for jj=1:length(offreactions) + offreactions(jj)= offreactions{jj}(1); + end + % Below 2 lines work for acetyl model, not min methyl model. + % [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); + % [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); + + disp(i) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % get basal metabolic state based on transcriptome + % if basalflag + % [fluxes, grate, solverobj_ccle] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth + % grate_ccle_exp_soft(i,1:2)= grate; % first 2 columns are basal met flux. + % ...3rd column will be met flux with a rxn maximized + % fluxes_allrxns(i,:) = fluxes; % correlate each row of fluxes_allreactions with auc + % model2.c(rxnpos) = epsilon_methylation ; + % [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + % grate_ccle_exp_soft(i,3) = fluxstate_gurobi(rxnpos); %methylation flux + % end + + for rxncount = 1:length(model2.rxns) + if basalflag + %minfluxflag = 1; + %[~, grate, solverobj_ccle] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth + model2.c(rxncount) = epsilon_methylation ; + [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + g_rate(i,rxncount) = fluxstate_gurobi(rxncount); %methylation flux + end + end + else + grate_ccle_exp_soft(i,:) = NaN; + end + end + % Save in local github repo and another local folder. wd is MeCorr + s_var= num2str(name(i_weight)); + s_name= strcat('grate_', s_var, 'FBA_allrxns'); + s_path1= strcat('VariablesSaved/MaxAllRxns_Results/', s_name); + save(s_path1, 'g_rate'); + s_path2= strcat('./../../../MaxAllRxns_local/', s_name); + save(s_path2, 'g_rate'); + + % figure; h = histogram(grate_ccle_exp_soft(:,3),70); + % xlabel('Predicted methylation flux') + % ylabel('Total cell lines') + % title('Distribution of methylation flux among CCLE cell lines','fontweight','normal') + + %fluxesAll_1pFBA= fluxes_allrxns; + % save('VariablesSaved\grate_3FBA_allRxns', 'g_rate'); + % save('VariablesSaved\fluxesAll_1FBA', 'fluxes_allrxns'); + % save('VariablesSaved\grate_1E-6', 'grate_ccle_exp_soft'); + % save('VariablesSaved\fluxstate_1E-6', 'fluxstate_gurobi'); +end +%% Correlation between flux and auc for a reaction +% flux_allrexns: 1031 cell lines by 3777 reactions +% drug_auc_expt: extract auc values for the 8 methyl drugs +% 8 x 3000 rhos +%fluxes_allrxns= g_rate; +rho= NaN(height(hmei_list),length(model2.rxns)); +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); + hmei_auc_dat = drug_auc_expt(ix,:); + + % Only uses the experiments of fluxes_allrxns that also are in hmei_auc_dat + [ix, pos] = ismember(hmei_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix); + v1 = hmei_auc_dat(:,2); v1= table2array(v1); + for nCol=1:size(fluxes_allrxns,2) + v4= fluxes_allrxns(pos,nCol); + [rho(j,nCol), rho_p(j,nCol)]= corr(v4, v1, 'rows', 'complete'); + end +end +disp('rho & p-value calculated') +%% Create Table of Significant Reactions by Correlation value +% Workflow: Change threshold. Change struct field name (e.g. above3) +% accordingly +sigTF= (abs(rho) > 0.3); +% sigTF= (abs(rho_p) < 0.001); +nSigExpt= sum(sigTF); % sum number of signif expts per rxn (sum each column) +sigIndRxn= (nSigExpt >= 1); +sigRxn= model2.rxns(sigIndRxn); disp(sigRxn) +sigRxnFlux= array2table(fluxes_allrxns(:,sigIndRxn)); +VarNames= genvarname(sigRxn); +sigRxnFlux.Properties.VariableNames=VarNames; + +[iDrug,iRxn]= find(sigTF); n= length(iDrug); disp(n) +sigRxnS.above3(1:n,1)= num2cell(rho(sigTF)); +sigRxnS.above3(1:n,2)= num2cell(rho_p(sigTF)); +sigRxnS.above3(1:n,3)= model2.rxns(iRxn); +sigRxnS.above3(1:n,4)= model2.rxnNames(iRxn); +sigRxnS.above3(1:n,5)= model2.subSystems(iRxn); +sigRxnS.above3(1:n,6)= cpd_infoMe(iDrug, 1); +sigRxnS.above3(1:n,7)= cpd_infoMe(iDrug, 3); +sigRxnS.above3(1:n,8)= cpd_infoMe(iDrug, 2); + +sigRxnT= sortrows(sigRxnS.above3, 1); % sort by correlation +sigRxnT= cell2table(sigRxnT); +sigRxnT.Properties.VariableNames{'sigRxnT1'}='Correlation'; +sigRxnT.Properties.VariableNames{'sigRxnT2'}='Pvalue'; +sigRxnT.Properties.VariableNames{'sigRxnT3'}='Rxn'; +sigRxnT.Properties.VariableNames{'sigRxnT4'}='RxnName'; +sigRxnT.Properties.VariableNames{'sigRxnT5'}='Subsystem'; +sigRxnT.Properties.VariableNames{'sigRxnT6'}='Compound'; +sigRxnT.Properties.VariableNames{'sigRxnT7'}='CpdActivity'; +sigRxnT.Properties.VariableNames{'sigRxnT8'}='CpdTarget'; + +% Rename variable and save +% sigRxnT_1pFBA=sigRxnT; +% save('sigRxnT_1pFBA', 'sigRxnT_1pFBA'); +% writetable(sigRxnT,('sigRxnT_4FBA.xlsx')); + +% sigRxnT_maxAllRxns_1FBA=sigRxnT; +% save('sigRxnT_maxAllRxns_1FBA', 'sigRxnT_maxAllRxns_1FBA'); +% writetable(sigRxnT_maxAllRxns_1FBA,('sigRxnT_maxAllRxns_1FBA.xlsx')); + +% sigRxnFlux_1pFBA= sigRxnFlux; +% save('sigRxnFlux_1pFBA', 'sigRxnFlux_1pFBA'); +%% See if fluxes are positive or negative for 6 FBA simulations +% load fluxesAll for all 6 weights (FBA) +clearvars intAll6 model2rxns +intAll6= int.intAll6{1}; model2rxns= cell2table(model2.rxns); +model2rxns.Properties.VariableNames{'Var1'}='Rxn'; +iRxnAll6= find(ismember(intAll6, model2rxns)); +sigRxnFlux6= {fluxesAll_1FBA, fluxesAll_2FBA, fluxesAll_3FBA, fluxesAll_4FBA, ... + fluxesAll_5FBA, fluxesAll_6FBA}; +posneg= zeros(6, height(intAll6)); +for iWeight=1:6 + fluxesi= sigRxnFlux6{iWeight}; + sigRxnFlux= fluxesi(:,iRxnAll6); + for iRxn= 1:height(intAll6) + if sum(sigRxnFlux(:,iRxn)<0) ==0 + posneg(iWeight,iRxn)= 1; + elseif sum(sigRxnFlux(:,iRxn)>0) ==0 + posneg(iWeight,iRxn)= -1; + end % For an element that =1, (rxn, weight) has both pos & neg flux + end +end +disp(posneg) +%% See if flux is pos/neg for each sigRxn of large loop +load('grate_1FBA_allRxns') +load('sigRxnT_maxAllRxns_1FBA') +sigRxn= cell2table(sigRxnT_maxAllRxns_1FBA.Rxn); +model2rxns= cell2table(model2.rxns); +isigRxn= find(ismember(sigRxn, model2rxns)); +sigRxnFlux= g_rate(:,isigRxn); + +posneg= zeros(1,height(sigRxn)); +iWeight=1; +for iRxn=1:height(sigRxn) + if sum(sigRxnFlux(:,iRxn)<0) ==0 + posneg(iWeight,iRxn)= 1; + elseif sum(sigRxnFlux(:,iRxn)>0) ==0 + posneg(iWeight,iRxn)= -1; + end +end +disp(posneg) +%% Calculate correlation between flux and auc. Each reaction maximized. +rho3= NaN(height(hmei_list),length(model2.rxns)); +for j= 1:height(hmei_list) + % Extract the auc for a drug in drug_auc_expt + fx = find(ismember(ctd2compoundidname_name, hmei_list(j,1))); + ix = ismember(drug_auc_expt(:,3), ctd2compoundidname_id(fx,1)); + hmei_auc_dat = drug_auc_expt(ix,:); + % Only uses the experiments of fluxes_allrxns that also are in hmei_auc_dat + [ix, pos] = ismember(hmei_auc_dat(:,1), exptidcelllinemediamatch(:,1)); + v1 = hmei_auc_dat(:,2); + v1= table2array(v1); + + for nCol= 1:size(g_rate,2) + v5= g_rate(pos, nCol); + rho3(j,nCol)= corr(v5, v1, 'rows', 'complete'); + end +end + +sigRho3Ind= false(1,length(model2.rxns)); +for j2= 13:13%length(model2.rxns) + ix= ~isnan(rho3(:,j2)); + sigRho3Ind(1,j2)= (sum(ix)>0); +end +sigRxn3= model2.rxns(sigRho3Ind); +disp(sigRxn3) + +save('VariablesSaved\rho3_1E-1', rho3); +save('VariablesSabed\sigRxn3_1E-1', sigRxn3); +%% predicting sensitivity to hme inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%load('grate_pFBA_1E-2'); load('fluxstate_pFBA_1E-2'); +clearvars r_basal +r_basal(height(hmei_list))= struct('rho',0, 'p',0, 'rlowerbound',0, 'rupperbound',0); +pp_flux= zeros(height(hmei_list), 3); +pp_grate= zeros(height(hmei_list), 3); +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) % 847. ix is a 1D array. Each drug appears multiple times (appears for each medium in drug_auc_expt) + hme_auc_dat = drug_auc_expt(ix,:); + + [ix, pos] = ismember(hme_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix); + v1 = hme_auc_dat(:,2); v1= table2array(v1); + v2 = grate_ccle_exp_soft(pos,:); % growth rate + v3 = fluxstate_gurobi(pos,:); + +% % Scatter plots of methylation (metabolic) flux & calculate correlation +% figure; +% scatter(v1, v3) +% ylabel({'Methylation Flux'},'fontname','helvetica'); %'fontweight','bold') +% xlabel({'Sensitivity (AUC)'; hmei_list.cpd_name{j}},'fontname','helvetica'); +% ylim([0, 0.06]) +% [R, P, RL, RU]= corrcoef(v1, v3); +% r_basal(j).rho=R; r_basal(j).p=P; r_basal(j).rlowerbound=RL; r_basal(j).rupperbound=RU; +% str=sprintf('r= %1.3f',r_basal(j).rho(1,2)); +% T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); +% set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + % Boxplot pairs of methylation flux (Actually graphing metabolic flux, + ...but we hypothesize met flux approximates methyl flux) + ix0 = (v3 < 0.05);sum(ix0); % 0.05 is the default threshold value + ix01 = (v3 > 0.05);sum(ix01); + [hh, pp_flux(j,3)] = ttest2(v1(ix0), v1(ix01)); % + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; + vv = NaN(length(v1), 2); + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); + + figure; + %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); + %hold on; + bh = boxplot(v1, groups,'symbol',''); + set(gca,'xticklabel',{'Low flux','High flux'}) + ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); + xlabel(hmei_list{j,1}); + %ylim([0 20]) + str=sprintf('p= %1.4f',pp_flux(j,3)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + ix0 = (v3 <= prctile(v3, 25)); sum(ix0); + ix01 = (v3 > prctile(v3, 25)); sum(ix01); + [hh, pp_flux(j,1)] = ttest2(v1(ix0), v1(ix01)); + + ix0 = (v3 <= prctile(v3, 50)); sum(ix0); + ix01 = (v3 > prctile(v3, 50)); sum(ix01); + [hh, pp_flux(j,2)] = ttest2(v1(ix0), v1(ix01)); + + % Boxplot pairs of growth rate %%%%%% + ix0 = (v2(:,2) < 0.05);sum(ix0); % 0.05 is the default threshold value + ix01 = (v2(:,2) > 0.05);sum(ix01); + [hh, pp_grate(j,3)] = ttest2(v1(ix0), v1(ix01)); % + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; + vv = NaN(length(v1), 2); + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); + + figure; + %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); + %hold on; + bh = boxplot(v1, groups,'symbol',''); + set(gca,'xticklabel',{'Low growth','High growth'}) + ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); + xlabel(hmei_list{j,1}); + %ylim([0 20]) + str=sprintf('p= %1.4f',pp_grate(j,3)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + ix0 = (v2(:,3) <= prctile(v2(:,3), 25)); sum(ix0); + ix01 = (v2(:,3) > prctile(v2(:,3), 25)); sum(ix01); + [hh, pp_grate(j,1)] = ttest2(v1(ix0), v1(ix01)); + + ix0 = (v2(:,3) <= prctile(v2(:,3), 50)); sum(ix0); + ix01 = (v2(:,3) > prctile(v2(:,3), 50)); sum(ix01); + [hh, pp_grate(j,2)] = ttest2(v1(ix0), v1(ix01)); + +% ix0 = (v2(:,2) <= prctile(v2(:,2), 25)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 25)); sum(ix01); +% [hh, pp_grate(j,1)] = ttest2(v1(ix0), v1(ix01)); +% +% ix0 = (v2(:,2) <= prctile(v2(:,2), 50)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 50)); sum(ix01); +% [hh, pp_grate(j,2)] = ttest2(v1(ix0), v1(ix01)); +end +r_basal= struct2table(r_basal); +r_basal.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +disp(r_basal) +pp_flux= array2table(pp_flux); +pp_flux.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +pp_flux.Properties.VariableNames{'pp_flux3'}= 'pf_hilo'; +pp_flux.Properties.VariableNames{'pp_flux1'}= 'pf_25prctile'; +pp_flux.Properties.VariableNames{'pp_flux2'}= 'pf_50prctile'; +disp(pp_flux) %t-test p-values for each drug +pp_grate= array2table(pp_grate); +pp_grate.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +pp_grate.Properties.VariableNames{'pp_grate3'}= 'pg_hilo'; +pp_grate.Properties.VariableNames{'pp_grate1'}= 'pg_25prctile'; +pp_grate.Properties.VariableNames{'pp_grate2'}= 'pg_50prctile'; +disp(pp_grate) +%save('VariablesSaved\pp_grate_1E-6_pFBA', 'pp_grate'); \ No newline at end of file diff --git a/matlab/MeCorr/MATLAB_CODE_me.m b/matlab/MeCorr/MATLAB_CODE_me.m new file mode 100755 index 0000000..5f8132c --- /dev/null +++ b/matlab/MeCorr/MATLAB_CODE_me.m @@ -0,0 +1,318 @@ + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % this code reproduces all the main text figures from the manuscript (Shen et al) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Bulk methylation model: +% 1) cd ./../scripts. Run make_eGEM to create bulk methylation model +% 2) cd ./../MeCorr. Run section 1 and 3 of this script. Only run Section 2 if you want to +% run the long for-loop +% 3) Parameters of interest to vary: epsilon_methylation, model2 +% 4) Change name of files to which variables are saved +model2 = min_model; +epsilon_methylation = 1E-4; % 1E-2 and 1E-1 are probably accurate (Scott) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % impact of basal metabolic state of CCLE cell lines on sensitivity to demethylase inhibitors +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load supplementary_software_code ctd2celllineidname exptidcelllinemediamatch ctd2celllineidname_id r* +%data from seashore-ludlow study, contains cell line names, growth media +load supplementary_software_code ctd2compoundidname_id drug_auc_expt ctd2compoundidname_name +%data from seashore-ludlow study, contains drug names, drug sensitivity data +load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz +%contains CCLE cell line names, gene expression data (z-transformed) +load supplementary_software_code hcommon_exptdat hcommon1 +%data from seashore-ludlow study, contains drug names, drug sensitivity data for cell lines that were screened against all 8 hme inhibitors +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 +% caompatible with function ismember +ctd2celllineidname= array2table(ctd2celllineidname); +ctd2celllineidname.Properties.VariableNames{'ctd2celllineidname1'}= 'ccl_name'; +ctd2celllineidname_id= array2table(ctd2celllineidname_id); +ctd2celllineidname_id.Properties.VariableNames{'ctd2celllineidname_id'}= 'master_ccl_id'; +ctd2compoundidname_name= array2table(ctd2compoundidname_name); +ctd2compoundidname_name.Properties.VariableNames{'ctd2compoundidname_name'}= 'cpd_name'; +ctd2compoundidname_id= array2table(ctd2compoundidname_id); +ctd2compoundidname_id.Properties.VariableNames{'ctd2compoundidname_id'}= 'cpd_id'; +celllinenames_ccle1= array2table(celllinenames_ccle1); +celllinenames_ccle1.Properties.VariableNames{'celllinenames_ccle1'}= 'ccl_name'; +drug_auc_expt= array2table(drug_auc_expt); +drug_auc_expt.Properties.VariableNames{'drug_auc_expt3'}= 'cpd_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt1'}= 'experiment_id'; +drug_auc_expt.Properties.VariableNames{'drug_auc_expt2'}= 'auc'; +exptidcelllinemediamatch= array2table(exptidcelllinemediamatch); %s4 +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch1'}= 'experiment_id'; +exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch2'}= 'master_ccl_id'; +hcommon1= array2table(hcommon1); +hcommon1.Properties.VariableNames{'hcommon1'}= 'experiment_id'; + +rxnpos = find(ismember(min_model.rxns,'LYSMTF1n')); + +MODE = 1; % reaction (1) or gene list (0) +epsilon = 1E-2; rho = 1; +kappa = 1; % parameters for integrating transcriptomics data. kappa is the strength of down regulation of genes (Chandrasekaran & Price, PNAS, 2010) +minfluxflag = 0; % 0: no Pfba, 1: Pfba +hmetransint = 1; +basalflag = 1; + +% If I don't want to run next for-loop: +%load('VariablesSaved\fluxstate_gurobi'); +%load('VariablesSaved\grate_ccle_exp_soft'); +%% long for-loop (1:1035) +grate_ccle_exp_soft = NaN(height(exptidcelllinemediamatch),2); %contains acetylation flux based on basal metabolic state +%grate_ccle_exp_soft_hdacsign = NaN(height(exptidcelllinemediamatch),8); %contains acetylation flux based on basal metabolic state and impact of each hme inhibitor teatment + +for i = 1:height(exptidcelllinemediamatch) + % match cell line data in CCLE with CTD2 + ii = find(ismember(ctd2celllineidname_id, exptidcelllinemediamatch(i,2))); + iii = find(ismember(celllinenames_ccle1, ctd2celllineidname(ii,1))); + if ~isempty(iii) + iii = iii(1); +% model2 = min_model; +% model2 = metabolicmodel; + + %find up and down-regulated genes in each cell line + ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); + offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % set the glucose uptake based on media + % default glucose is -5 for rpmi + if r1(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi + elseif r2(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% dmem with low glucose + elseif r3(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% Emem.. + elseif r4(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3/2;% mccoy + elseif r5(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*1/2;% mem + elseif r6(i) + model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5*3.15/2;% dmem:f12 + end + + %find reactions from differentially expressed genes + onreactions= findRxnsFromGenes(model2, ongenes); + onreactions= struct2cell(onreactions); + for jj=1:length(onreactions) + onreactions(jj)= onreactions{jj}(1); + end + offreactions= findRxnsFromGenes(model2, offgenes); + offreactions= struct2cell(offreactions); + for jj=1:length(offreactions) + offreactions(jj)= offreactions{jj}(1); + end + % Below 2 lines work for acetyl model, not min methyl model. +% [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); +% [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); + + disp(i) + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % get basal metabolic state based on transcriptome + if basalflag + [~,grate_ccle_exp_soft(i,1), solverobj_ccle(i,1)] =... + constrain_flux_regulation(model2,onreactions,offreactions,... + kappa,rho,epsilon,MODE,[], minfluxflag); %impact on growth + model2.c(rxnpos) = epsilon_methylation ; + [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,... + offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + grate_ccle_exp_soft(i,2) = fluxstate_gurobi(rxnpos); %methylation flux + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % hme inhibitor impact on transcriptome +% if hmetransint +% for kk = 4:-1:1 % match drug order +% ongenesh = hdacexpallgeneids((hdacexpfcs(:,kk) > 1.1) & (hdacexpfcs(:,kk + 4) < 0.01)); +% ongenesh = intersect(ongenesh, model2.genes); +% offgenesh = hdacexpallgeneids((hdacexpfcs(:,kk) < 0.9) & (hdacexpfcs(:,kk + 4) < 0.01)); +% offgenesh = intersect(offgenesh, model2.genes); +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% [~,~,onreactionsh,~] = deleteModelGenes(model2, ongenesh); +% [~,~,offreactionsh,~] = deleteModelGenes(model2, offgenesh); +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % remove conflicting transcripts +% onreactionsh0 = setdiff(onreactionsh, offreactions); +% offreactionsh0 = setdiff(offreactionsh, onreactions); +% onreactions0 = setdiff(onreactions, offreactionsh); +% offreactions0 = setdiff(offreactions, onreactionsh); +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% onreactions1h = [onreactions0; onreactionsh0]; +% offreactions1h = [offreactions0; offreactionsh0]; +% epsilon2 = [repmat(0, size(offreactions1h))]; +% model2.c(rxnpos) = 0; +% % [fluxstate_gurobi,grate_ccle_exp_soft_hdacsign(i,kk)] = constrain_flux_regulation(model2,onreactions1h,offreactions1h,kappa,rho,epsilon,MODE, epsilon2,minfluxflag); % impact on growth +% model2.c(rxnpos) = epsilon_methylation ; +% [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions1h,offreactions1h,kappa,rho,epsilon,MODE,epsilon2,minfluxflag); +% grate_ccle_exp_soft_hdacsign(i,kk + 4) = fluxstate_gurobi(rxnpos); %acetylation flux +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% end +% end + + else + grate_ccle_exp_soft(i,:) = NaN; +% grate_ccle_exp_soft_hdacsign(i,:) = NaN; + end +end + + +figure; h = histogram(grate_ccle_exp_soft(:,2),70); + xlabel('Predicted methylation flux') + ylabel('Total cell lines') +title('Distribution of methylation flux among CCLE cell lines','fontweight','normal') + +save('VariablesSaved\grate_ccle_exp_soft_1E-4', 'grate_ccle_exp_soft'); +save('VariablesSaved\fluxstate_gurobi_1E-4', 'fluxstate_gurobi'); +save('VariablesSaved\solverobj_ccle_1E-4', 'solverobj_ccle'); +%% predicting sensitivity to hme inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +hmei_list = {'BRD-A02303741';'BIX-01294';'methylstat';'QW-BI-011';... + 'UNC0321';'CBB-1007';'UNC0638';'GSK-J4'}; +hmei_list= cell2table(hmei_list); +hmei_list.Properties.VariableNames{'hmei_list'}='cpd_name'; +clearvars r_basal +r_basal(height(hmei_list))= struct('rho',0, 'p',0, 'rlowerbound',0, 'rupperbound',0); +pp_basal= zeros(height(hmei_list), 3); +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) % 847. ix is a 1D array. Each drug appears multiple times (appears for each medium in drug_auc_expt) + hme_auc_dat = drug_auc_expt(ix,:); + + [ix, pos] = ismember(hme_auc_dat(:,1), exptidcelllinemediamatch(:,1)); sum(ix); % + v1 = hme_auc_dat(:,2); v1= table2array(v1); + v2 = grate_ccle_exp_soft(pos,:); % growth rate + v3 = fluxstate_gurobi(pos,:); + + % Scatter plots of methylation (metabolic) flux & calculate correlation + figure; + scatter(v1, v3) + ylabel({'Methylation Flux'},'fontname','helvetica'); %'fontweight','bold') + xlabel({'Sensitivity (AUC)'; hmei_list.cpd_name{j}},'fontname','helvetica'); + ylim([0, 0.06]) + [R, P, RL, RU]= corrcoef(v1, v3); + r_basal(j).rho=R; r_basal(j).p=P; r_basal(j).rlowerbound=RL; r_basal(j).rupperbound=RU; + str=sprintf('r= %1.3f',r_basal(j).rho(1,2)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + % Boxplot pairs of methylation flux (Actually grapphing metabolic flux, + ...but we hypothesize met flux approximates methyl flux) + ix0 = (v3 < 0.05);sum(ix0); % 0.05 is the default threshold value + ix01 = (v3 > 0.05);sum(ix01); + [hh, pp_basal(j,3)] = ttest2(v1(ix0), v1(ix01)); % + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; + vv = NaN(length(v1), 2); + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); + + figure; + %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); + %hold on; + bh = boxplot(v1, groups,'symbol','') + set(gca,'xticklabel',{'Low growth','High growth'}) + ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); + xlabel(hmei_list{j,1}); + ylim([0 20]) + str=sprintf('p= %1.3f',pp_basal(j,3)); + T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str); + set(T, 'fontsize', 14, 'verticalAlignment', 'top', 'horizontalAlignment', 'left'); + + ix0 = (v3 <= prctile(v3, 25)); sum(ix0); + ix01 = (v3 > prctile(v3, 25)); sum(ix01); + [hh, pp_basal(j,1)] = ttest2(v1(ix0), v1(ix01)); + + ix0 = (v3 <= prctile(v3, 50)); sum(ix0); + ix01 = (v3 > prctile(v3, 50)); sum(ix01); + [hh, pp_basal(j,2)] = ttest2(v1(ix0), v1(ix01)); + + % Boxplot pairs of growth rate %%%%%% +% ix0 = (v2(:,2) < 0.05);sum(ix0); % 0.05 is the default threshold value +% ix01 = (v2(:,2) > 0.05);sum(ix01); +% [hh, pp_basal(j,3)] = ttest2(v1(ix0), v1(ix01)); % +% +% groups = NaN(size(ix0)); +% groups(ix0) = 1; groups(ix01) = 2; +% vv = NaN(length(v1), 2); +% vv(1:sum(groups == 1),1) = v1(groups == 1); +% vv(1:sum(groups == 2),2) = v1(groups == 2); +% +% figure; +% %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); +% %hold on; +% bh = boxplot(v1, groups,'symbol','') +% set(gca,'xticklabel',{'Low flux','High flux'}) +% ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); +% xlabel(hmei_list{j,1}); +% ylim([0 20]) +% +% ix0 = (v2(:,2) <= prctile(v2(:,2), 25)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 25)); sum(ix01); +% [hh, pp_basal(j,1)] = ttest2(v1(ix0), v1(ix01)); +% +% ix0 = (v2(:,2) <= prctile(v2(:,2), 50)); sum(ix0); +% ix01 = (v2(:,2) > prctile(v2(:,2), 50)); sum(ix01); +% [hh, pp_basal(j,2)] = ttest2(v1(ix0), v1(ix01)); +end +r_basal= struct2table(r_basal); +r_basal.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +disp(r_basal) +pp_basal= array2table(pp_basal); +pp_basal.Properties.RowNames= {hmei_list.cpd_name{1:8}}; +pp_basal.Properties.VariableNames{'pp_basal3'}= 'pp_highlow'; +pp_basal.Properties.VariableNames{'pp_basal1'}= 'pp_25prctile'; +pp_basal.Properties.VariableNames{'pp_basal2'}= 'pp_50prctile'; +disp(pp_basal) %t-test p-values for each drug + +%% Last section uses grate_ccle_exp_soft_hdacsign, which was not calculated. Do not run this section. +%predicting variation in sensitivity between hme inhibitors based on basal +%metabolic state and drug impact - comparison with drug sensitivity data from seashore-ludlow study +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +g1 = grate_ccle_exp_soft_hdacsign(:, 5:8); +g1 = g1 - repmat(ignoreNaN(g1, @median,2),1,4); % + +[ix pos] = ismember(hcommon1, exptidcelllinemediamatch(:,1)); sum(ix) % match cell lines +v1 = hcommon_exptdat; +v1 = v1 - repmat(median(v1,2),1,4); +v2 = g1(pos,:); +v1 = v1(:); + + +ix0 = (v2(:) <= prctile(v2(:), 25)); sum(ix0) +ix01 = (v2(:) >= prctile(v2(:), 25)); sum(ix01) +[hh pp] = ttest2(v1(ix0), v1(ix01)) % 2e-52.. + +ix0 = (v2(:) <= prctile(v2(:), 50)); sum(ix0) +ix01 = (v2(:) >= prctile(v2(:), 50)); sum(ix01) +[hh pp] = ttest2(v1(ix0), v1(ix01)) % 3e-44 + +ix0 = (v2(:) <= 0.05); sum(ix0) +ix01 = (v2(:) >0.05); sum(ix01) +[hh pp] = ttest2(v1(ix0), v1(ix01)) %2e -62 +mean(v1(ix0)) % -0.01 +mean(v1(ix01)) % -2.61 + + groups = NaN(size(ix0)); + groups(ix0) = 1; groups(ix01) = 2; +vv = NaN(length(v1), 2); + + vv(1:sum(groups == 1),1) = v1(groups == 1); + vv(1:sum(groups == 2),2) = v1(groups == 2); +figure; +bh = boxplot(vv,'symbol','' ,'orientation','horizontal') + set(gca,'yticklabel',{'No','Large'}) + xlabel('Observed differential sensitivity between drugs in a cell line') +title({'Predicted differential sensitivity between drugs' 'vs Observed differential sensitivity (AUC)'}) + + +figure +h1 = histfit(vv(:,1),20,'kernel')%, +set(h1(1),'facecolor','g','facealpha',.15,'edgecolor','none') +set(h1(2),'color','g') +hold on +h2 = histfit(vv(:,2),20,'kernel') +set(h2(1),'facecolor','m','facealpha',.15,'edgecolor','none') +set(h2(2),'color','m') +hh = legend([h1(2),h2(2)],'No difference','Large difference') +title(hh,'Predicted differential acetylation') diff --git a/matlab/MeCorr/MATLAB_CODE_methyl.m b/matlab/MeCorr/MATLAB_CODE_methyl.m index f8e5b1a..cce5560 100755 --- a/matlab/MeCorr/MATLAB_CODE_methyl.m +++ b/matlab/MeCorr/MATLAB_CODE_methyl.m @@ -5,19 +5,18 @@ % originally assigned in 1st module. epsilon_methylation = 1E-2; % or 1E-1 +% Bulk methylation model: +% 1) cd ./../scripts. Run make_eGEM to create bulk methylation model +% model = min; +% 2) cd ./../MeCorr. Run methylVariables.m to create variables for methylation drug data. +% 3) Run the last module of this script, which is split into 3 sections + % Acetylation model: % 1) Load it % load supplementary_software_code acetylation_model %contains metabolic model with nuclear acetylation reaction % model = acetylation_model; % 2) Run methylVariables.m to create variables for methylation drug data. % 3) Run the last module of this script, which is split into 3 sections - -% Bulk methylation model: -% 1) Run make_eGEM to create bulk methylation model -% model = min; -% 2) Run methylVariables.m to create variables for methylation drug data. -% 3) Run the last module of this script, which is split into 3 sections - %% impact of nutrient sources on acetylation - figure 2A load supplementary_software_code acetylation_model %contains metabolic model with nuclear acetylation reaction % load recon1 % contains a methylation rxn, but genes are number ID's, so @@ -25,7 +24,7 @@ load supplementary_software_code labels media_exchange1 mediareactions1 %list of nutrient conditions and uptake rates posgluc = 1385; % glucose uptake reaction in recon1. % Changed the reaction to be found from acetylation to methylation - rxnpos = find(ismember(acetylation_model.rxns,'LYSMTF1n')); % rxnpos = 2451 +rxnpos = find(ismember(acetylation_model.rxns,'LYSMTF1n')); % rxnpos = 2451 objpos = find(acetylation_model.c); %biomass objective minfluxflag = 0; % no PFBA @@ -87,16 +86,14 @@ for i = 2:length(unqgenes) modeltemp = acetylation_model ; - modeltemp = deleteModelGenes(modeltemp,unqgenes(i)); + modeltemp = deleteModelGenes(modeltemp,unqgenes(i)); model3 = modeltemp; - [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 end - - model3.c(rxnpos) = epsilon_acetylation; % epsilon is a weight + model3.c(rxnpos) = epsilon_methylation; % epsilon is a weight [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); if ~isempty(solf.x) && ~isnan(solf.x) acet_screen_rpmi(i,1) = solf.x(rxnpos); % impact on acetylation flux @@ -106,9 +103,9 @@ end acet_screen_rpmi(1,:) = NaN; - - [sx spos] = sort(acet_screen_rpmi(:,1)); - acetgenessorted = unqgenes(spos); + + [sx spos] = sort(acet_screen_rpmi(:,1)); + acetgenessorted = unqgenes(spos); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% impact of culture media components on acetylation - figure 2C @@ -124,7 +121,6 @@ glucflag1 = [ ones(1,6), 0 ,0, 0, 0]; expval1 = [ ones(1,7), 0 , 1, 0]; - [ix, aapos] = ismember( mediareactions1(20:37),acetylation_model.rxns); posgluc = 1385; % glucose uptake reaction in recon1. glnpos = 1386; % glutamine @@ -133,7 +129,7 @@ fattyacidpos = 1445; % linoleic acid model3 = acetylation_model; -model3.c(rxnpos) = epsilon_acetylation ; +model3.c(rxnpos) = epsilon_methylation ; [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); wild_type_acet = solf.x(rxnpos); % default acetylation flux @@ -155,24 +151,18 @@ model3.lb(pyrpos) = -10; end - - model3.c(rxnpos) = epsilon_acetylation; + model3.c(rxnpos) = epsilon_methylation; [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); acet_media_screen_dmem(i,1) = solf.x(rxnpos); % impact on acetylation flux - - end disp('consistency with 10 different experimental conditions - Part I') sum((acet_media_screen_dmem > 0.05*wild_type_acet) == expval') % 9 - - for i = 1:10 model3 = acetylation_model; if ~glucflag1(i) model3.lb(posgluc) = 0; end - switch i case 2 % no ca+ @@ -186,9 +176,7 @@ case 9 % add fatty acid model3.lb(fattyacidpos) = -5; end - - - model3.c(rxnpos) = epsilon_acetylation; + model3.c(rxnpos) = epsilon_methylation; [solf.x,sol11] = constrain_flux_regulation(model3,[],[],0,0,0,[],[],minfluxflag); acet_media_screen_dmem1(i,1) = solf.x(rxnpos); % impact on acetylation flux end @@ -238,11 +226,10 @@ disp(i) - [fluxstate_gurobi,grate_ccle_exp_acetdat(i,1), solverobj_ccle(i,1)] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[], minfluxflag); % impact on growth - - model2.c(rxnpos) = epsilon_acetylation; - [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[],minfluxflag); - grate_ccle_exp_acetdat(i,2) = fluxstate_gurobi(rxnpos); %acetylation flux + [fluxstate_gurobi,grate_ccle_exp_acetdat(i,1), solverobj_ccle(i,1)] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[], minfluxflag); % impact on growth + model2.c(rxnpos) = epsilon_methylation; + [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE ,[],minfluxflag); + grate_ccle_exp_acetdat(i,2) = fluxstate_gurobi(rxnpos); %acetylation flux end end @@ -271,9 +258,8 @@ [solf.x,sol11] = constrain_flux_regulation(model2,[],[],0,0,0,[],[],minfluxflag); disp(i) - growth_basal2(i,1) = solf.x(objpos); % impact on growth - model2.c(rxnpos) = epsilon_acetylation; + model2.c(rxnpos) = epsilon_methylation; [solf.x,sol11] = constrain_flux_regulation(model2,[],[],0,0,0,[],[],minfluxflag); growth_basal2(i,2) = solf.x(rxnpos); %acetylation flux @@ -300,17 +286,19 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% impact of basal metabolic state of CCLE cell lines on sensitivity to demethylase inhibitors %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -load supplementary_software_code ctd2celllineidname exptidcelllinemediamatch ctd2celllineidname_id r* -%data from seashore-ludlow study, contains cell line names, growth media -load supplementary_software_code ctd2compoundidname_id drug_auc_expt ctd2compoundidname_name -%data from seashore-ludlow study, contains drug names, drug sensitivity data +%load supplementary_software_code ctd2celllineidname ctd2celllineidname_id +load supplementary_software_code exptidcelllinemediamatch r* +% data from seashore-ludlow study, contains cell line names, growth media +%load supplementary_software_code ctd2compoundidname_id ctd2compoundidname_name +load supplementary_software_code drug_auc_expt +% data from seashore-ludlow study, contains drug names, drug sensitivity data load supplementary_software_code celllinenames_ccle1 ccleids_met ccle_expression_metz -%contains CCLE cell line names, gene expression data (z-transformed) -load supplementary_software_code hcommon_exptdat hcommon1 -%data from seashore-ludlow study, contains drug names, drug sensitivity data +% contains CCLE cell line names, gene expression data (z-transformed) +%load supplementary_software_code hcommon_exptdat hcommon1 +% data from seashore-ludlow study, contains drug names, drug sensitivity data ...for cell lines that were screened against all 4 hdac inhibitors -load supplementary_software_code hdacexpfcs hdacexpallgeneids -%contains gene expression data after treatment with hdac inhibitors +%load supplementary_software_code hdacexpfcs hdacexpallgeneids +% contains gene expression data after treatment with hdac inhibitors rxnpos = 2451; % loading the above variables changes rxnpos to 3754, which is out of bounds for fluxstate_gurobi % Change data type to be compatible with functions @@ -318,14 +306,15 @@ exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch1'}='index_cpd'; exptidcelllinemediamatch.Properties.VariableNames{'exptidcelllinemediamatch2'}='index_ccl'; celllinenames_ccle1= cell2table(celllinenames_ccle1); -hcommon1= array2table(hcommon1); -hcommon1.Properties.VariableNames{'hcommon1'}= 'index_cpd'; -% Name table variables for clarity +% hcommon1= array2table(hcommon1); +% hcommon1.Properties.VariableNames{'hcommon1'}= 'index_cpd'; +% Name table variables for clarity. draw_auc_expt is not used in code, but +% useful to look at drug_auc_expt_t= array2table(drug_auc_expt); drug_auc_expt_t.Properties.VariableNames{'drug_auc_expt1'}='index_cpd'; drug_auc_expt_t.Properties.VariableNames{'drug_auc_expt2'}='auc'; drug_auc_expt_t.Properties.VariableNames{'drug_auc_expt3'}='index_ccl'; - +%% MODE = 1; % reaction (1) or gene list (0) epsilon = 1E-2; rho = 1; kappa = 1; % parameters for integrating transcriptomics data. kappa is the strength of down regulation of genes (Chandrasekaran & Price, PNAS, 2010) @@ -338,19 +327,20 @@ for i = 1:height(exptidcelllinemediamatch) % match cell line data in CCLE with CTD2 - ii = find(ismember(ctd2celllineidname_id_me, exptidcelllinemediamatch(i,2))); - iii = find(ismember(celllinenames_ccle1, ctd2celllineidname_me(ii,1))); + ii = find(ismember(ctd2clidname_id_me, exptidcelllinemediamatch(i,2))); + iii = find(ismember(celllinenames_ccle1, ctd2clidname_me(ii,1))); if ~isempty(iii) iii = iii(1); - model2 = acetylation_model; - %model2 = min; + + model2 = min; + %model2 = acetylation_model; - %find up and down-regulated genes in each cell line + % find up and down-regulated genes in each cell line ongenes = unique(ccleids_met(ccle_expression_metz(:,iii) > 2)); offgenes = unique(ccleids_met(ccle_expression_metz(:,iii) < -2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % set the glucose uptake based on media + % set the glucose uptake based on media % default glucose is -5 for rpmi if r1(i) model2.lb(find(ismember(model2.rxns, {'EX_glc(e)'}))) = -5;% no change rpmi @@ -367,18 +357,32 @@ end %find reactions from differentially expressed genes - [~,~,onreactions,~] = deleteModelGenes(model2, ongenes); - [~,~,offreactions,~] = deleteModelGenes(model2, offgenes); + onreactions= findRxnsFromGenes(model2, ongenes); + onreactions= struct2cell(onreactions); + for jj=1:length(onreactions) + onreactions(jj)= onreactions{jj}(1); + end + offreactions= findRxnsFromGenes(model2, offgenes); + offreactions= struct2cell(offreactions); + for jj=1:length(offreactions) + offreactions(jj)= offreactions{jj}(1); + end + % Below 2 lines work for acetyl model, not min methyl model. + %[~,~,onreactions,~] = deleteModelGenes(model2, ongenes); + %[~,~,offreactions,~] = deleteModelGenes(model2, offgenes); disp(i) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % get basal metabolic state based on transcriptome - if basalflag - [fluxstate_gurobi,grate_ccle_exp_soft(i,1), solverobj_ccle(i,1)] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[], minfluxflag); % impact on growth - model2.c(rxnpos) = epsilon_methylation ; - [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); - grate_ccle_exp_soft(i,2) = fluxstate_gurobi(rxnpos); %acetylation flux - end + if basalflag + [~,grate_ccle_exp_soft(i,1), solverobj_ccle(i,1)] =... + constrain_flux_regulation(model2,onreactions,offreactions,... + kappa,rho,epsilon,MODE,[], minfluxflag); %impact on growth + model2.c(rxnpos) = epsilon_methylation ; + [fluxstate_gurobi] = constrain_flux_regulation(model2,onreactions,... + offreactions,kappa,rho,epsilon,MODE,[],minfluxflag); + grate_ccle_exp_soft(i,2) = fluxstate_gurobi(rxnpos); %methylation flux + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % hdac inhibitor impact on transcriptome % if hdactransint @@ -417,28 +421,29 @@ figure; h = histogram(grate_ccle_exp_soft(:,2),70); - xlabel('Predicted acetylation flux') - ylabel('Total cell lines') -title('Distribution of acetylation flux among CCLE cell lines','fontweight','normal') +xlabel('Predicted methylation flux') +ylabel('Total cell lines') +title('Distribution of methylation flux among CCLE cell lines','fontweight','normal') -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%predicting sensitivity to hdac inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study +%% I created this section to surpass the long for-loop +% predicting sensitivity to hdac inhibitors based on basal metabolic state - comparison with drug sensitivity data from seashore-ludlow study %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %hdaclist = {'LBH-589','vorinostat','entinostat','belinostat'} hmei_list= {'BRD-A02303741';'BIX-01294';'methylstat';'QW-BI-011';... 'UNC0321';'CBB-1007';'UNC0638';'GSK-J4'}; hmei_list= cell2table(hmei_list); hmei_list.Properties.VariableNames{'hmei_list'}='compound_name'; -%% I created this section to surpass the long for loop + for j = 1:height(hmei_list) -fx = find(ismember(ctd2compoundidname_name_me, hmei_list(j,1))) -ix = ismember(drug_auc_me(:,1), ctd2compoundidname_id_me(fx,1)); +fx = find(ismember(ctd2cpdidname_name_me, hmei_list(j,1))) +ix = ismember(drug_auc_me(:,1), ctd2cpdidname_id_me(fx,1)); sum(ix) %597. ix is 1D logical array & 423 drugs tested -> Each drug ...appears multiple times (appears for each medium in drug_auc_expt) hmei_auc_dat_me = drug_auc_me(ix,:); -[ix pos] = ismember(hmei_auc_dat_me(:,1), exptidcelllinemediamatch(:,1)); sum(ix) % 597 v1 = hmei_auc_dat_me(:,2); v1= table2array(v1); +% Match index_cpd +[ix pos] = ismember(hmei_auc_dat_me(:,3), exptidcelllinemediamatch(:,1)); sum(ix) % 597 v2 = grate_ccle_exp_soft(pos,:); ix0 = (v2(:,2) < 0.05);sum(ix0) % 0.05 is the default threshold value @@ -453,15 +458,20 @@ vv(1:sum(groups == 2),2) = v1(groups == 2); %Col2 is values > 0.05 figure; - %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); - %hold on; - bh = boxplot(v1) -% bh = boxplot(v1, groups,'Symbol','') % groups contains only NaN values. -% No groups found. - set(gca,'xticklabel',{'Low flux','High flux'}) -ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); -xlabel(hmei_list{j,1}); -ylim([0 20]) + scatter(v1, v2(:,2)) + figure; + scatter(KMTi_auc(:,j), grate_ccle_exp_soft(:,2)) + +% figure; +% %clf; UnivarScatter(vv,'Width', 0.3, 'PointSize', 11,'MarkerEdgeColor','w','LineWidth',0.1);%, 'markerfacealpha',0.5); +% %hold on; +% bh = boxplot(v1) +% % bh = boxplot(v1, groups,'Symbol','') % groups contains only NaN values. +% % No groups found. +% set(gca,'xticklabel',{'Low flux','High flux'}) +% ylabel({'Sensitivity (AUC)'},'fontname','helvetica');%,'fontweight','bold'); +% xlabel(hmei_list{j,1}); +% ylim([0 20]) ix0 = (v2(:,2) <= prctile(v2(:,2), 25)); sum(ix0) ix01 = (v2(:,2) > prctile(v2(:,2), 25)); sum(ix01) @@ -480,13 +490,13 @@ g1 = grate_ccle_exp_soft_hdacsign(:, 5:8); g1 = g1 - repmat(ignoreNaN(g1, @median,2),1,4); % -[ix pos] = ismember(hcommon1, exptidcelllinemediamatch(:,1)); sum(ix) % match cell lines +[ix pos] = ismember(ctd2cpdidname_id_me, exptidcelllinemediamatch(:,1)); sum(ix) % match cell lines +%[ix pos] = ismember(hcommon1, exptidcelllinemediamatch(:,1)); sum(ix) % original v1 = KMTi_auc; v1 = v1 - repmat(median(v1,2), 1, height(hmei_list)); v2 = g1(pos,:); v1 = v1(:); - ix0 = (v2(:) <= prctile(v2(:), 25)); sum(ix0) ix01 = (v2(:) >= prctile(v2(:), 25)); sum(ix01) [hh pp] = ttest2(v1(ix0), v1(ix01)) % 2e-52.. diff --git a/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/grate_1FBA_allRxns.mat b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/grate_1FBA_allRxns.mat new file mode 100755 index 0000000..78edd80 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/grate_1FBA_allRxns.mat differ diff --git a/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.mat b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.mat new file mode 100644 index 0000000..34e3c17 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.xlsx b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.xlsx new file mode 100644 index 0000000..703816a Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/MaxAllRxns_Results/sigRxnT_maxAllRxns_1FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-4.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-4.mat new file mode 100755 index 0000000..4882d5e Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-4.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-6.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-6.mat new file mode 100755 index 0000000..c8020c6 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_1E-6.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-1.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-1.mat new file mode 100755 index 0000000..2cf5775 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-1.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-2.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-2.mat new file mode 100755 index 0000000..566cfbc Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-2.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-3.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-3.mat new file mode 100755 index 0000000..25ee5b4 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/fluxstate_pFBA_1E-3.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-4.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-4.mat new file mode 100755 index 0000000..e12036c Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-4.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-6.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-6.mat new file mode 100755 index 0000000..319f5af Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1E-6.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1FBA_allRxns.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1FBA_allRxns.mat new file mode 100755 index 0000000..78edd80 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_1FBA_allRxns.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-1.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-1.mat new file mode 100755 index 0000000..4335e7f Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-1.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-2.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-2.mat new file mode 100755 index 0000000..7a561a3 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-2.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-3.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-3.mat new file mode 100755 index 0000000..2631ebc Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/grate_pFBA_1E-3.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-1_pFBA.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-1_pFBA.mat new file mode 100755 index 0000000..1753b1b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-1_pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-2_pFBA.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-2_pFBA.mat new file mode 100755 index 0000000..7bb1468 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_flux_1E-2_pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-1_pFBA.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-1_pFBA.mat new file mode 100755 index 0000000..9d66b9a Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-1_pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-2_pFBA.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-2_pFBA.mat new file mode 100755 index 0000000..85e9b3d Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-2_pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-4_FBA.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-4_FBA.mat new file mode 100755 index 0000000..359d3f6 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-4_FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-6.mat b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-6.mat new file mode 100755 index 0000000..b549a6b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/VarForPairedBoxplots/pp_grate_1E-6.mat differ diff --git a/matlab/MeCorr/VariablesSaved/analysis.xlsx b/matlab/MeCorr/VariablesSaved/analysis.xlsx new file mode 100755 index 0000000..4687bef Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/analysis.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_1FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_1FBA.mat new file mode 100644 index 0000000..7c2a17f Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_1FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA.mat new file mode 100755 index 0000000..5dd388c Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA_actuallyFBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA_actuallyFBA.mat new file mode 100755 index 0000000..e9ba0e9 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_1pFBA_actuallyFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_2FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_2FBA.mat new file mode 100755 index 0000000..cebaad4 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_2FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_3FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_3FBA.mat new file mode 100755 index 0000000..6ef1463 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_3FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_3pFBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_3pFBA.mat new file mode 100755 index 0000000..1c1f727 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_3pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_4FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_4FBA.mat new file mode 100755 index 0000000..d772410 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_4FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_5FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_5FBA.mat new file mode 100755 index 0000000..1994a56 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_5FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/fluxesAll_6FBA.mat b/matlab/MeCorr/VariablesSaved/fluxesAll_6FBA.mat new file mode 100755 index 0000000..ba5a991 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/fluxesAll_6FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/intersectingRxns.mat b/matlab/MeCorr/VariablesSaved/intersectingRxns.mat new file mode 100755 index 0000000..c915fa6 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/intersectingRxns.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnFlux_1FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnFlux_1FBA.mat new file mode 100755 index 0000000..9a1bb24 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnFlux_1FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnFlux_1pFBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnFlux_1pFBA.mat new file mode 100755 index 0000000..f96e09d Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnFlux_1pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnS_4FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnS_4FBA.mat new file mode 100755 index 0000000..04828b6 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnS_4FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.mat new file mode 100644 index 0000000..48b797b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.xlsx new file mode 100644 index 0000000..2a2a76b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_1FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA.mat new file mode 100755 index 0000000..50b1e13 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA_Correct.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA_Correct.xlsx new file mode 100755 index 0000000..067d9dc Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_1pFBA_Correct.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.mat new file mode 100755 index 0000000..b08f23b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.xlsx new file mode 100755 index 0000000..5ca9853 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_2FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.mat new file mode 100755 index 0000000..ef1ea63 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.xlsx new file mode 100755 index 0000000..a1087e3 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_3FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.mat new file mode 100755 index 0000000..7ba1e6b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.xlsx new file mode 100755 index 0000000..271ea00 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_3pFBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.mat new file mode 100755 index 0000000..e73626b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.xlsx new file mode 100755 index 0000000..dad25f7 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_4FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.mat new file mode 100755 index 0000000..134395e Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.xlsx new file mode 100755 index 0000000..b8d17df Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_5FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.mat new file mode 100755 index 0000000..819d780 Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.xlsx new file mode 100755 index 0000000..1edfc2a Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_6FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.mat b/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.mat new file mode 100755 index 0000000..10dd5fe Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.mat differ diff --git a/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.xlsx b/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.xlsx new file mode 100755 index 0000000..bce3f8b Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/sigRxnT_allRxns_1FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/~$sigRxnT_1pFBA.xlsx b/matlab/MeCorr/VariablesSaved/~$sigRxnT_1pFBA.xlsx new file mode 100644 index 0000000..06872ac Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/~$sigRxnT_1pFBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/~$sigRxnT_2FBA.xlsx b/matlab/MeCorr/VariablesSaved/~$sigRxnT_2FBA.xlsx new file mode 100644 index 0000000..06872ac Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/~$sigRxnT_2FBA.xlsx differ diff --git a/matlab/MeCorr/VariablesSaved/~$sigRxnT_6FBA.xlsx b/matlab/MeCorr/VariablesSaved/~$sigRxnT_6FBA.xlsx new file mode 100644 index 0000000..06872ac Binary files /dev/null and b/matlab/MeCorr/VariablesSaved/~$sigRxnT_6FBA.xlsx differ diff --git a/matlab/MeCorr/constrain_flux_regulation.m b/matlab/MeCorr/constrain_flux_regulation.m index dfbf5ba..5995631 100755 --- a/matlab/MeCorr/constrain_flux_regulation.m +++ b/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)) diff --git a/matlab/MeCorr/methylVariables.m b/matlab/MeCorr/methylVariables.m index ea1ab87..f9b341a 100755 --- a/matlab/MeCorr/methylVariables.m +++ b/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); diff --git a/matlab/models/min.mat b/matlab/models/min.mat new file mode 100755 index 0000000..bf70df5 Binary files /dev/null and b/matlab/models/min.mat differ diff --git a/matlab/models/min2SLC.mat b/matlab/models/min2SLC.mat new file mode 100755 index 0000000..203ec62 Binary files /dev/null and b/matlab/models/min2SLC.mat differ diff --git a/matlab/scripts/6hm_1E-6.fig b/matlab/scripts/6hm_1E-6.fig deleted file mode 100755 index 8de8c8c..0000000 Binary files a/matlab/scripts/6hm_1E-6.fig and /dev/null differ diff --git a/matlab/scripts/6hm_1E-6.jpg b/matlab/scripts/6hm_1E-6.jpg deleted file mode 100755 index 50add7b..0000000 Binary files a/matlab/scripts/6hm_1E-6.jpg and /dev/null differ diff --git a/matlab/scripts/make_eGEM.m b/matlab/scripts/make_eGEM.m index a2341bd..8b43b1e 100644 --- a/matlab/scripts/make_eGEM.m +++ b/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'); \ No newline at end of file +%save('./../models/min.mat', 'min'); +%save('./../models/min2SLC.mat', 'min'); +min_model=min; +clearvars min \ No newline at end of file diff --git a/python/tissue b/python/tissue new file mode 100644 index 0000000..e69de29