Libraries
library(tidyverse)
library(ggplot2)
library(plotrix)
library(car)
library(rstatix)
library(SimComp)
library(readxl)
library(dplyr)
dataset of the control data
d_BM_exp_2_3_bm_imprint <- read_excel(path = "/mnt/c/Users/Toshiya Matsushima/OneDrive/R projects/BM_project/data analysis using Rstudio/BM_project_dataset.xlsx", sheet = "exp_2_3_bm_imprint")
d_BM_exp_2_3_bm_imprint_00_ctrl <- d_BM_exp_2_3_bm_imprint %>% filter(drug=="00_ctrl") %>%
arrange(batch) %>%
mutate(id_num=seq(1:nrow(.)))
d_BM_exp_2_3_bm_imprint_00_ctrl
mean(d_BM_exp_2_3_bm_imprint_00_ctrl$bm)
[1] 136.3089
std.error(d_BM_exp_2_3_bm_imprint_00_ctrl$bm)
[1] 16.13828
mean(d_BM_exp_2_3_bm_imprint_00_ctrl$imprint)
[1] 419.5203
std.error(d_BM_exp_2_3_bm_imprint_00_ctrl$imprint)
[1] 20.25447
BM score
Bootstrap computations
d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- NULL
for (i in 1:10000){
set.seed(i)
d_BM_exp_2_3_bm_imprint_00_ctrl_test <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% sample_n(size = 10, replace = F) %>% arrange(id_num)
bm_mean_test_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_test[,4], 2, mean)
d_BM_exp_2_3_bm_imprint_00_ctrl_train <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% setdiff(d_BM_exp_2_3_bm_imprint_00_ctrl_test)
bm_mean_train_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_train[,4], 2, mean)
temp_diff <- bm_mean_test_data - bm_mean_train_data
d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp <- tibble(n_bootstrap = i, mean_diff = temp_diff)
d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>% bind_rows(d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp)
}
d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap
90%, 95%, and 99% confidence intervals (one-sided testing is assumed)
quantile(d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap$mean_diff, probs = c(0.001, 0.01, 0.05, 0.10))
0.1% 1% 5% 10%
-181.13551 -142.16726 -97.65323 -77.08611
density distribution curve of BM calculated from the control dataset
fig <- d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>%
ggplot(mapping = aes(x = mean_diff), size =3)+
geom_vline(xintercept=0, linetype = "solid", color ="blue", size = 1)+
geom_vline(xintercept=-97.65323, linetype = "solid", color ="dark red", size = 1)+
geom_vline(xintercept=-142.1673, linetype = "solid", color ="red", size = 1)+
geom_hline(yintercept=0, linetype = "solid", color = "grey")+
geom_density()+
theme_classic(base_size=15)
fig
ggsave(plot = fig, filename = "exp_2_BM_distribution.png", dpi = 300, height = 10, width = 15, units = "cm")
imprint score
Bootstrap computations
d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- NULL
for (i in 1:10000){
set.seed(i)
d_BM_exp_2_3_bm_imprint_00_ctrl_test <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% sample_n(size = 10, replace = F) %>% arrange(id_num)
imprint_mean_test_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_test [,5], 2, mean)
d_BM_exp_2_3_bm_imprint_00_ctrl_train <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% setdiff(d_BM_exp_2_3_bm_imprint_00_ctrl_test)
imprint_mean_train_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_train [,5], 2, mean)
temp_diff <- imprint_mean_test_data - imprint_mean_train_data
d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp <- tibble(n_bootstrap = i, mean_diff = temp_diff)
d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>% bind_rows(d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp)
}
d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap
99.9%, 99%, 95%, and 90% confidence intervals (one-sided testing is assumed)
quantile(d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap$mean_diff, probs = c(0.001, 0.01, 0.05, 0.10))
0.1% 1% 5% 10%
-289.45321 -201.61150 -135.32757 -99.51062
density distribution curve of BM calculated from the control dataset
fig <- d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>%
ggplot(mapping = aes(x = mean_diff), size =3)+
geom_vline(xintercept=0, linetype = "solid", color ="blue", size = 1)+
geom_vline(xintercept=-135.32757, linetype = "solid", color ="dark red", size = 1)+
geom_vline(xintercept=-201.61150, linetype = "solid", color ="red", size = 1)+
geom_hline(yintercept=0, linetype = "solid", color = "grey")+
geom_density()+
theme_classic(base_size=15)
fig
ggsave(plot = fig, filename = "exp_2_IMPRINTING_distribution.png", dpi = 300, height = 10, width = 15, units = "cm")
SessionInfo
sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] plotly_4.10.0 ggsci_2.9 readxl_1.3.1 SimComp_3.3 rstatix_0.7.0 car_3.0-12 carData_3.0-5 plotrix_3.8-2 ggbeeswarm_0.6.0 forcats_0.5.1 stringr_1.4.0
[12] dplyr_1.0.7 purrr_0.3.4 readr_2.1.1 tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] nlme_3.1-144 fs_1.5.2 lubridate_1.8.0 httr_1.4.2 tools_3.6.3 backports_1.4.1 utf8_1.2.2 R6_2.5.1
[9] vipor_0.4.5 lazyeval_0.2.2 DBI_1.1.1 mgcv_1.8-31 colorspace_2.0-2 withr_2.4.3 tidyselect_1.1.1 compiler_3.6.3
[17] survPresmooth_1.1-11 mratios_1.4.2 cli_3.1.0 rvest_1.0.2 xml2_1.3.3 sandwich_3.0-1 labeling_0.4.2 scales_1.1.1
[25] mvtnorm_1.1-3 digest_0.6.29 rmarkdown_2.11 pkgconfig_2.0.3 htmltools_0.5.2 dbplyr_2.1.1 fastmap_1.1.0 htmlwidgets_1.5.4
[33] rlang_0.4.12 rstudioapi_0.13 jquerylib_0.1.4 generics_0.1.1 farver_2.1.0 zoo_1.8-9 jsonlite_1.7.2 gtools_3.9.2
[41] magrittr_2.0.1 Matrix_1.2-18 Rcpp_1.0.7 munsell_0.5.0 fansi_0.5.0 abind_1.4-5 lifecycle_1.0.1 stringi_1.7.6
[49] multcomp_1.4-18 yaml_2.2.1 MASS_7.3-51.5 grid_3.6.3 crayon_1.4.2 lattice_0.20-40 haven_2.4.3 splines_3.6.3
[57] hms_1.1.1 knitr_1.37 pillar_1.6.4 codetools_0.2-16 reprex_2.0.1 glue_1.6.0 drc_3.0-1 evaluate_0.14
[65] data.table_1.14.2 modelr_0.1.8 vctrs_0.3.8 tzdb_0.2.0 cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1 xfun_0.29
[73] broom_0.7.10 viridisLite_0.4.0 survival_3.1-8 beeswarm_0.4.0 TH.data_1.1-0 ellipsis_0.3.2
LS0tCnRpdGxlOiAiQk1fcHJvamVjdF9leHBfMl9ib290c3RyYXAiCm91dHB1dDogaHRtbF9ub3RlYm9vawphdXRob3I6ICJUb3NoaXlhIE1hdHN1c2hpbWEiCmRhdGU6ICJgciBmb3JtYXQoU3lzLnRpbWUoKSlgIgotLS0KCiMjIExpYnJhcmllcwpgYGB7cn0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShwbG90cml4KQpsaWJyYXJ5KGNhcikKbGlicmFyeShyc3RhdGl4KQpsaWJyYXJ5KFNpbUNvbXApCmxpYnJhcnkocmVhZHhsKQpsaWJyYXJ5KGRwbHlyKQpgYGAKCgojIGRhdGFzZXQgb2YgdGhlIGNvbnRyb2wgZGF0YQpgYGB7cn0KZF9CTV9leHBfMl8zX2JtX2ltcHJpbnQgPC0gcmVhZF9leGNlbChwYXRoID0gIi9tbnQvYy9Vc2Vycy9Ub3NoaXlhIE1hdHN1c2hpbWEvT25lRHJpdmUvUiBwcm9qZWN0cy9CTV9wcm9qZWN0L2RhdGEgYW5hbHlzaXMgdXNpbmcgUnN0dWRpby9CTV9wcm9qZWN0X2RhdGFzZXQueGxzeCIsIHNoZWV0ID0gImV4cF8yXzNfYm1faW1wcmludCIpCmRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwgPC0gZF9CTV9leHBfMl8zX2JtX2ltcHJpbnQgJT4lIGZpbHRlcihkcnVnPT0iMDBfY3RybCIpICU+JQogIGFycmFuZ2UoYmF0Y2gpICU+JQogIG11dGF0ZShpZF9udW09c2VxKDE6bnJvdyguKSkpCmRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwKCm1lYW4oZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybCRibSkKc3RkLmVycm9yKGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwkYm0pCgptZWFuKGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwkaW1wcmludCkKc3RkLmVycm9yKGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwkaW1wcmludCkKYGBgCgojIEJNIHNjb3JlCiMjIEJvb3RzdHJhcCBjb21wdXRhdGlvbnMKYGBge3J9CmRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfQk1fYm9vdHN0cmFwIDwtIE5VTEwKZm9yIChpIGluIDE6MTAwMDApewogICAgc2V0LnNlZWQoaSkKICAKICAgIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfdGVzdCA8LSBkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsICU+JSBzYW1wbGVfbihzaXplID0gMTAsIHJlcGxhY2UgPSBGKSAlPiUgYXJyYW5nZShpZF9udW0pCiAgICBibV9tZWFuX3Rlc3RfZGF0YSA8LSBhcHBseShkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX3Rlc3RbLDRdLCAyLCBtZWFuKQogICAgCiAgICBkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX3RyYWluIDwtIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwgJT4lIHNldGRpZmYoZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF90ZXN0KQogICAgYm1fbWVhbl90cmFpbl9kYXRhIDwtIGFwcGx5KGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfdHJhaW5bLDRdLCAyLCBtZWFuKQogICAgCiAgICB0ZW1wX2RpZmYgPC0gYm1fbWVhbl90ZXN0X2RhdGEgLSBibV9tZWFuX3RyYWluX2RhdGEKICAgIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfQk1fYm9vdHN0cmFwX3RlbXAgPC0gdGliYmxlKG5fYm9vdHN0cmFwID0gaSwgbWVhbl9kaWZmID0gdGVtcF9kaWZmKQogICAgZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAgPC0gZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAgJT4lIGJpbmRfcm93cyhkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX0JNX2Jvb3RzdHJhcF90ZW1wKQp9CgpkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX0JNX2Jvb3RzdHJhcApgYGAKCiMjIDkwJSwgOTUlLCBhbmQgOTklIGNvbmZpZGVuY2UgaW50ZXJ2YWxzIChvbmUtc2lkZWQgdGVzdGluZyBpcyBhc3N1bWVkKQpgYGB7cn0KcXVhbnRpbGUoZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAkbWVhbl9kaWZmLCBwcm9icyA9IGMoMC4wMDEsIDAuMDEsIDAuMDUsIDAuMTApKQpgYGAKCgojIyBkZW5zaXR5IGRpc3RyaWJ1dGlvbiBjdXJ2ZSBvZiBCTSBjYWxjdWxhdGVkIGZyb20gdGhlIGNvbnRyb2wgZGF0YXNldApgYGB7cn0KZmlnIDwtIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfQk1fYm9vdHN0cmFwICU+JQogICAgZ2dwbG90KG1hcHBpbmcgPSBhZXMoeCA9IG1lYW5fZGlmZiksIHNpemUgPTMpKwogICAgICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9MCwgbGluZXR5cGUgPSAic29saWQiLCBjb2xvciA9ImJsdWUiLCBzaXplID0gMSkrCiAgICAgIGdlb21fdmxpbmUoeGludGVyY2VwdD0tOTcuNjUzMjMsIGxpbmV0eXBlID0gInNvbGlkIiwgY29sb3IgPSJkYXJrIHJlZCIsIHNpemUgPSAxKSsKICAgICAgZ2VvbV92bGluZSh4aW50ZXJjZXB0PS0xNDIuMTY3MywgbGluZXR5cGUgPSAic29saWQiLCBjb2xvciA9InJlZCIsIHNpemUgPSAxKSsKICAgIGdlb21faGxpbmUoeWludGVyY2VwdD0wLCBsaW5ldHlwZSA9ICJzb2xpZCIsIGNvbG9yID0gImdyZXkiKSsKICAgIGdlb21fZGVuc2l0eSgpKwogICAgdGhlbWVfY2xhc3NpYyhiYXNlX3NpemU9MTUpCmZpZwpnZ3NhdmUocGxvdCA9IGZpZywgZmlsZW5hbWUgPSAiZXhwXzJfQk1fZGlzdHJpYnV0aW9uLnBuZyIsIGRwaSA9IDMwMCwgaGVpZ2h0ID0gMTAsIHdpZHRoID0gMTUsIHVuaXRzID0gImNtIikKYGBgCgoKIyBpbXByaW50IHNjb3JlCiMjIEJvb3RzdHJhcCBjb21wdXRhdGlvbnMKYGBge3J9CmRfSU1QUklOVF9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAgPC0gTlVMTApmb3IgKGkgaW4gMToxMDAwMCl7CiAgICBzZXQuc2VlZChpKQogIAogICAgZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF90ZXN0IDwtIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmwgJT4lIHNhbXBsZV9uKHNpemUgPSAxMCwgcmVwbGFjZSA9IEYpICU+JSBhcnJhbmdlKGlkX251bSkKICAgIGltcHJpbnRfbWVhbl90ZXN0X2RhdGEgPC0gYXBwbHkoZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF90ZXN0IFssNV0sIDIsIG1lYW4pCiAgICAKICAgIGRfQk1fZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfdHJhaW4gPC0gZF9CTV9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybCAlPiUgc2V0ZGlmZihkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX3Rlc3QpCiAgICBpbXByaW50X21lYW5fdHJhaW5fZGF0YSA8LSBhcHBseShkX0JNX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX3RyYWluIFssNV0sIDIsIG1lYW4pCiAgICAKICAgIHRlbXBfZGlmZiA8LSBpbXByaW50X21lYW5fdGVzdF9kYXRhIC0gaW1wcmludF9tZWFuX3RyYWluX2RhdGEKICAgIGRfSU1QUklOVF9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXBfdGVtcCA8LSB0aWJibGUobl9ib290c3RyYXAgPSBpLCBtZWFuX2RpZmYgPSB0ZW1wX2RpZmYpCiAgICBkX0lNUFJJTlRfZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfQk1fYm9vdHN0cmFwIDwtIGRfSU1QUklOVF9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAgJT4lIGJpbmRfcm93cyhkX0lNUFJJTlRfZXhwXzJfM19ibV9pbXByaW50XzAwX2N0cmxfQk1fYm9vdHN0cmFwX3RlbXApCn0KCmRfSU1QUklOVF9leHBfMl8zX2JtX2ltcHJpbnRfMDBfY3RybF9CTV9ib290c3RyYXAKYGBgCgoKIyMgOTkuOSUsIDk5JSwgOTUlLCBhbmQgOTAlIGNvbmZpZGVuY2UgaW50ZXJ2YWxzIChvbmUtc2lkZWQgdGVzdGluZyBpcyBhc3N1bWVkKQpgYGB7cn0KcXVhbnRpbGUoZF9JTVBSSU5UX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX0JNX2Jvb3RzdHJhcCRtZWFuX2RpZmYsIHByb2JzID0gYygwLjAwMSwgMC4wMSwgMC4wNSwgMC4xMCkpCmBgYAoKCiMjIGRlbnNpdHkgZGlzdHJpYnV0aW9uIGN1cnZlIG9mIEJNIGNhbGN1bGF0ZWQgZnJvbSB0aGUgY29udHJvbCBkYXRhc2V0CmBgYHtyfQpmaWcgPC0gZF9JTVBSSU5UX2V4cF8yXzNfYm1faW1wcmludF8wMF9jdHJsX0JNX2Jvb3RzdHJhcCAlPiUKICAgIGdncGxvdChtYXBwaW5nID0gYWVzKHggPSBtZWFuX2RpZmYpLCBzaXplID0zKSsKICAgICAgZ2VvbV92bGluZSh4aW50ZXJjZXB0PTAsIGxpbmV0eXBlID0gInNvbGlkIiwgY29sb3IgPSJibHVlIiwgc2l6ZSA9IDEpKwogICAgICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9LTEzNS4zMjc1NywgbGluZXR5cGUgPSAic29saWQiLCBjb2xvciA9ImRhcmsgcmVkIiwgc2l6ZSA9IDEpKwogICAgICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9LTIwMS42MTE1MCwgbGluZXR5cGUgPSAic29saWQiLCBjb2xvciA9InJlZCIsIHNpemUgPSAxKSsKICAgIGdlb21faGxpbmUoeWludGVyY2VwdD0wLCBsaW5ldHlwZSA9ICJzb2xpZCIsIGNvbG9yID0gImdyZXkiKSsKICAgIGdlb21fZGVuc2l0eSgpKwogICAgdGhlbWVfY2xhc3NpYyhiYXNlX3NpemU9MTUpCmZpZwpnZ3NhdmUocGxvdCA9IGZpZywgZmlsZW5hbWUgPSAiZXhwXzJfSU1QUklOVElOR19kaXN0cmlidXRpb24ucG5nIiwgZHBpID0gMzAwLCBoZWlnaHQgPSAxMCwgd2lkdGggPSAxNSwgdW5pdHMgPSAiY20iKQpgYGAKCgoKIyBTZXNzaW9uSW5mbwpgYGB7cn0Kc2Vzc2lvbkluZm8oKQpgYGA=