Adding 5 individual figure scripts - #44
Conversation
| library(ggpubr) | ||
| library(patchwork) | ||
|
|
||
| syn <- synapseLogin() |
There was a problem hiding this comment.
@CamiloPosso where is this function defined? My running fails on this
|
Just added |
sgosline
left a comment
There was a problem hiding this comment.
i had some issues. I'm wondering if it's my setup so i'm including my sessionInfo() below:
R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.4.1
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] survminer_0.5.2 survival_3.8-6 GSVA_2.4.8 msigdbr_26.1.0 glmnet_4.1-10
[6] Matrix_1.7-4 data.table_1.18.2.1 patchwork_1.3.2 ggpubr_0.6.3 ggplot2_4.0.2
[11] MSnSet.utils_0.2.0 MSnbase_2.36.0 ProtGenerics_1.42.0 S4Vectors_0.48.0 mzR_2.44.0
[16] Rcpp_1.1.1 Biobase_2.70.0 BiocGenerics_0.56.0 generics_0.1.4 dplyr_1.2.0
[21] reticulate_1.45.0
loaded via a namespace (and not attached):
[1] fs_1.6.7 matrixStats_1.5.0 bitops_1.0-9 httr_1.4.8
[5] RColorBrewer_1.1-3 doParallel_1.0.17 dynamicTreeCut_1.63-1 tools_4.5.2
[9] backports_1.5.0 R6_2.6.1 HDF5Array_1.38.0 lazyeval_0.2.2
[13] mgcv_1.9-4 rhdf5filters_1.22.0 GetoptLong_1.1.0 withr_3.0.2
[17] gridExtra_2.3 preprocessCore_1.72.0 WGCNA_1.74 cli_3.6.5
[21] labeling_0.4.3 S7_0.2.1 randomForest_4.7-1.2 genefilter_1.92.0
[25] askpass_1.2.1 yulab.utils_0.2.4 foreign_0.8-91 DOSE_4.4.0
[29] R.utils_2.13.0 MetaboCoreUtils_1.18.1 invgamma_1.2 limma_3.66.0
[33] rstudioapi_0.18.0 impute_1.84.0 RSQLite_2.4.6 shape_1.4.6.1
[37] gtools_3.9.5 car_3.1-5 zip_2.3.3 GO.db_3.22.0
[41] MALDIquant_1.22.3 abind_1.4-8 R.methodsS3_1.8.2 lifecycle_1.0.5
[45] yaml_2.3.12 edgeR_4.8.2 carData_3.0-6 SummarizedExperiment_1.40.0
[49] rhdf5_2.54.1 gplots_3.3.0 qvalue_2.42.0 SparseArray_1.10.9
[53] BiocFileCache_3.0.0 grid_4.5.2 blob_1.3.0 crayon_1.5.3
[57] PSMatch_1.14.0 lattice_0.22-9 beachmat_2.26.0 cowplot_1.2.0
[61] annotate_1.88.0 KEGGREST_1.50.0 magick_2.9.1 pillar_1.11.1
[65] knitr_1.51 ComplexHeatmap_2.26.1 fgsea_1.36.2 GenomicRanges_1.62.1
[69] rjson_0.2.23 boot_1.3-32 codetools_0.2-20 fastmatch_1.1-8
[73] glue_1.8.0 outliers_0.15 memuse_4.2-3 pcaMethods_2.2.0
[77] MultiAssayExperiment_1.36.1 vctrs_0.7.1 png_0.1-8 Rdpack_2.6.6
[81] gtable_0.3.6 assertthat_0.2.1 cachem_1.1.0 xfun_0.56
[85] openxlsx_4.2.8.1 rbibutils_2.4.1 S4Arrays_1.10.1 Seqinfo_1.0.0
[89] reformulas_0.4.4 ncdf4_1.24 SingleCellExperiment_1.32.0 iterators_1.0.14
[93] statmod_1.5.1 ROCR_1.0-12 nlme_3.1-168 bit64_4.6.0-1
[97] filelock_1.0.3 R.cache_0.17.0 irlba_2.3.7 affyio_1.80.0
[101] KernSmooth_2.23-26 otel_0.2.0 rpart_4.1.24 colorspace_2.1-2
[105] DBI_1.3.0 Hmisc_5.2-5 nnet_7.3-20 tidyselect_1.2.1
[109] bit_4.6.0 compiler_4.5.2 curl_7.0.0 graph_1.88.1
[113] httr2_1.2.2 htmlTable_2.4.3 h5mread_1.2.1 DelayedArray_0.36.0
[117] checkmate_2.3.4 scales_1.4.0 caTools_1.18.3 affy_1.88.0
[121] varSelRF_0.7-9 rappdirs_0.3.4 SpatialExperiment_1.20.0 stringr_1.6.0
[125] MSnID_1.44.2 digest_0.6.39 minqa_1.2.8 rmarkdown_2.30
[129] XVector_0.50.0 htmltools_0.5.9 pkgconfig_2.0.3 base64enc_0.1-6
[133] lme4_2.0-1 umap_0.2.10.0 sparseMatrixStats_1.22.0 MatrixGenerics_1.22.0
[137] dbplyr_2.5.2 fastmap_1.2.0 rlang_1.1.7 GlobalOptions_0.1.3
[141] htmlwidgets_1.6.4 farver_2.1.2 jsonlite_2.0.0 BiocParallel_1.44.0
[145] mzID_1.48.0 GOSemSim_2.36.0 R.oo_1.27.1 BiocSingular_1.26.0
[149] magrittr_2.0.4 Formula_1.2-5 Rhdf5lib_1.32.0 babelgene_22.9
[153] MsCoreUtils_1.22.1 vsn_3.78.1 stringi_1.8.7 pROC_1.19.0.1
[157] MASS_7.3-65 AnnotationHub_4.0.0 plyr_1.8.9 FField_0.1.0
[161] parallel_4.5.2 ggrepel_0.9.7 Biostrings_2.78.0 splines_4.5.2
[165] Spectra_1.20.1 circlize_0.4.17 locfit_1.5-9.12 Boruta_9.0.0
[169] igraph_2.2.2 fastcluster_1.3.0 QFeatures_1.20.0 ggsignif_0.6.4
[173] pkgload_1.5.0 ScaledMatrix_1.18.0 reshape2_1.4.5 Heatplus_3.18.0
[177] BiocVersion_3.22.0 XML_3.99-0.22 evaluate_1.0.5 BiocManager_1.30.27
[181] nloptr_2.2.1 foreach_1.5.2 tidyr_1.3.2 openssl_2.3.5
[185] purrr_1.2.1 clue_0.3-67 BiocBaseUtils_1.12.0 rsvd_1.0.5
[189] broom_1.0.12 xtable_1.8-8 AnnotationFilter_1.34.0 RSpectra_0.16-2
[193] rstatix_0.7.3 tibble_3.3.1 memoise_2.0.1 AnnotationDbi_1.72.0
[197] IRanges_2.44.0 cluster_2.1.8.2 sva_3.58.0 GSEABase_1.72.0
| df.u$subcohort = pData(m_exp28)[, "subcohort"] | ||
| df.v = xx$data.v | ||
|
|
||
| p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], |
There was a problem hiding this comment.
This fails for me - it feels like df.u is trying to subset a plot object, since xx is actually a plot?
| left_join(., p_data, by = "sample_name") %>% filter(sample_name %in% sampleNames(m_exp28)) | ||
|
|
||
| plots = list() | ||
| for (cell_type in unique(maturation_df$X)){ |
There was a problem hiding this comment.
Since you are doing very similar plots, it might be useful to make a function to compute protein scores and weights across the cohort. A For loop is fine, but if you end up using it later, a function outside would be optimal.
| yy = plot_pca(m_, phenotype = "Treatment", output_type = "full") | ||
| df.u = yy$data.u | ||
| df.u$Treatment = pData(m_)[, "Treatment"] | ||
| df.u = df.u %>% |
There was a problem hiding this comment.
df.u is not defined here, I"m guessing this is from Figure 0 code. However, it fails for the same reason - yy is a graphics object, not a PCA object.
| ```{r} | ||
|
|
||
|
|
||
| diffexp = diffexp_helper(m_exp28_[complete_features, ], contrast_var = "refractory", |
There was a problem hiding this comment.
m_exp28_ is not defined. I changed this to m_exp28 which is defined above, but then it can't find the refractory column of the pData.
| return(res_) | ||
| } | ||
|
|
||
| hazards_res = hazards_(msnset[features_50perc, ]) |
There was a problem hiding this comment.
This line fails, I think it's the subsetting of msnset:
> msnset[features_50perc, ]
Error in orig[[nm]][i, , ..., drop = drop] : subscript out of bounds
| source("../scripts/predict_helper.R") | ||
| source("../scripts/s2n_model.R") | ||
|
|
||
| m_exp28 = readRDS("msnset_exp28_global.RDS") |
There was a problem hiding this comment.
cannot find this file - maybe download from synapse?
|
|
||
| complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] | ||
| m_210$new = paste(m_210$InitialAMLDiagnosis, m_210$PostChemotherapy) | ||
| beataml_clinical = openxlsx::read.xlsx("../../PNNL_clinical_summary_12_08_2021_updated_OS_4patients_02_28_2022.xlsx") |
|
I don't think it was on your end, just some minor typos/bugs and a function (modification of plot_pca) I forgot to include in the commits. All the scripts should be working now (famous last words). |
| PCA of SSGSEA | ||
|
|
||
| ```{r} | ||
| treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") |
There was a problem hiding this comment.
I think it might be helpful early on to agree on nomenclature for this sample, perhaps calling then 'Untreated' and 'Ven + HMA treated' throughout would be easier.
| Diffexp of SSGSEA | ||
|
|
||
| ```{r} | ||
| diffexp_results = diffexp_helper(m_exp28, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) |
There was a problem hiding this comment.
Same naming conventions go here - "Refractory", "Short response", "Long response" I think is what we agreed on.
| Showing the signature in the exvivo data (training data) | ||
|
|
||
| ```{r} | ||
| # model <- logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = TRUE, |
There was a problem hiding this comment.
I had to uncomment this to run the next line.
|
|
||
|
|
||
| ```{r} | ||
|
|
There was a problem hiding this comment.
i'd try to comment each block to describe what you're doing, and also prose/plain text above.
| p | ||
|
|
||
| model_features = results_ttest %>% filter(pval_cutoff == 0.001) %>% pull(feature) %>% unique() | ||
| plot_pca(m_exp28[model_features, ], phenotype = "subcohort") + ggtitle("PCA of 21 model features") + |
| # feat2 = all_features_df_ # using the ven+hma model | ||
| common = intersect(feat1$feature, feat2$feature) | ||
|
|
||
| feat1 = left_join(feat1, diffexp_results_venhma %>% select(logFC, feature, t_test_pval), by = "feature") |
There was a problem hiding this comment.
this fails because feat1 is not defined. if i define it above, then it failsb ecause the diffexp_results_venhma is not defined.
| complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] | ||
| m_210$new = paste(m_210$InitialAMLDiagnosis, m_210$PostChemotherapy) | ||
| beataml_clinical = openxlsx::read.xlsx(syn$get("syn74822940")$path) | ||
| ``` |
There was a problem hiding this comment.
is the above section in every markdown? if so would it be helpful to have a 'load data' script?
| refractory_markers = model$full_model$features %>% filter(group_ == "DV") %>% pull(feature) | ||
| longresponse_markers = model$full_model$features %>% filter(group_ != "DV") %>% pull(feature) | ||
|
|
||
| xx = predict_helper_(m_md_anderson, refractory_markers, longresponse_markers) %>% |
There was a problem hiding this comment.
is m_md_anderson defined? this fails for me
This script uses only the data from synapse, and generate all the panels for the exvivo DecVen and Ven signatures. It runs LOOCV to find the number of markers for each signature, and then applies the signature in the invivo data, showing boxplots of the results.
Added the stats table to synapse
Figure 0 - Introductory panels for the two datasets, clinical + exvivo
Figure 1 - Applying existing signatures in clinical dataset
Figure 2 - Deriving exvivo signature of resistance and applying to clinical data
Figure 3 - Deriving clinical signature of resistance and applying to exvivo data
Figure 4 - Markers from Cox Hazard with survival time and also with remission time
Figure 5 - Applying signatures to BeatAML and correlating with drug AUC.
scripts - Necessary helper scripts, loading data, extracting signatures, applying them, etc.