Libraries

library(tidyverse)
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(ggbeeswarm)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(SimComp)
library(readxl)
library(dplyr)
library(pwr)

Dataset

dataset_exp_0_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_0_preliminary_red_yellow")

dataset_exp_0_bm_imprint
## # A tibble: 42 × 11
##    id    batch   run    bm imprint   sex drug    dose   label  colour dose_level
##    <chr> <dbl> <dbl> <dbl>   <dbl> <dbl> <chr>   <chr>  <chr>  <chr>       <dbl>
##  1 id5       1  4076   137     516     1 00_ctrl 0_null 00_re… 0_red           0
##  2 id6       1  2079    54     134     0 00_ctrl 0_null 00_re… 0_red           0
##  3 id7       1  3701   361     578     1 00_ctrl 0_null 00_re… 0_red           0
##  4 id8       1   615   113     589     1 00_ctrl 0_null 00_re… 0_red           0
##  5 id10      4   602   262     584     1 00_ctrl 0_null 00_re… 0_red           0
##  6 id13      4  3073   288     584     0 00_ctrl 0_null 00_re… 0_red           0
##  7 id17      4   511   121     -12     1 00_ctrl 0_null 00_re… 0_red           0
##  8 id19      4  4610   198     437     0 00_ctrl 0_null 00_re… 0_red           0
##  9 id22      5   748   159     591     0 00_ctrl 0_null 00_re… 0_red           0
## 10 id23      5    32   159     447     0 00_ctrl 0_null 00_re… 0_red           0
## # … with 32 more rows
dataset_exp_0_bm_imprint %>% 
  group_by (label) %>% 
  get_summary_stats(bm)
## # A tibble: 4 × 14
##   label    variable     n   min   max median    q1    q3   iqr   mad  mean    sd
##   <chr>    <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 00_red_… bm          10    54   361   159  125    246   121   63.0  185.  93.0
## 2 01_red_… bm          11  -292   564   129  -37    314.  352. 301.   129. 274. 
## 3 02_yell… bm          10  -189   567   444.  97.5  532.  435  162.   325. 284. 
## 4 03_yell… bm          11  -366   504   167   88    343   255  219.   165. 250. 
## # … with 2 more variables: se <dbl>, ci <dbl>
dataset_exp_0_bm_imprint %>% 
  group_by (label) %>% 
  get_summary_stats(imprint)
## # A tibble: 4 × 14
##   label    variable     n   min   max median    q1    q3   iqr   mad  mean    sd
##   <chr>    <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 00_red_… imprint     10   -12   591    547  440. 584    144.  63.8  445.  213.
## 2 01_red_… imprint     11  -311   596    146  -94  378    472  368.   130   292.
## 3 02_yell… imprint     10   -37   522    188   -9  304.   314. 288.   192.  211.
## 4 03_yell… imprint     11  -597   592   -262 -569    0.5  570. 430.  -185.  414.
## # … with 2 more variables: se <dbl>, ci <dbl>

x-y plot in all data

fig <- ggplot(data = dataset_exp_0_bm_imprint, mapping = aes(x=bm, y=imprint, color=label))+
  geom_point(size=3)+
  geom_vline(xintercept=0, linetype = "dotted")+
  geom_hline(yintercept=0, linetype = "dotted")+
  geom_smooth(method = "lm")+
  xlim(-600, 600)+
  ylim(-600, 600)+
  theme_classic()
fig
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 15 rows containing missing values (geom_smooth).

ggsave(plot = fig, filename = "exp_0_imprint_vs_bm_control_vpa.png", dpi = 300, height = 10, width = 13, units = "cm")
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 15 rows containing missing values (geom_smooth).

imprinting of all data, linear fitting

fit_exp_0_imprint_bm <- lm (imprint ~ bm * colour, data = dataset_exp_0_bm_imprint)
summary (fit_exp_0_imprint_bm)
## 
## Call:
## lm(formula = imprint ~ bm * colour, data = dataset_exp_0_bm_imprint)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -641.61 -238.01   23.26  211.82  532.18 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        149.3820    89.1486   1.676   0.1020  
## bm                   0.8369     0.3506   2.387   0.0221 *
## colour1_yellow    -215.2363   130.3292  -1.651   0.1069  
## bm:colour1_yellow   -0.5876     0.4394  -1.337   0.1891  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 322.7 on 38 degrees of freedom
## Multiple R-squared:  0.2805, Adjusted R-squared:  0.2237 
## F-statistic: 4.937 on 3 and 38 DF,  p-value: 0.005418
car::Anova(fit_exp_0_imprint_bm)
## Anova Table (Type II tests)
## 
## Response: imprint
##            Sum Sq Df F value   Pr(>F)   
## bm         499297  1  4.7960 0.034735 * 
## colour    1074371  1 10.3198 0.002681 **
## bm:colour  186161  1  1.7882 0.189100   
## Residuals 3956096 38                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fig <- ggplot(data = dataset_exp_0_bm_imprint, mapping = aes(x=bm, y=imprint, color=label))+
  geom_point(size=3)+
  facet_wrap(~drug)+
  geom_vline(xintercept=0, linetype = "dotted")+
  geom_hline(yintercept=0, linetype = "dotted")+
  geom_smooth(method = "lm")+
  xlim(-400, 600)+
  ylim(-400, 600)+
  theme_classic()
fig
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 15 rows containing missing values (geom_smooth).

ggsave(plot = fig, filename = "exp_2_imprint_vs_bm_all_wrap.png", dpi = 300, height = 10, width = 15, units = "cm")
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 4 rows containing non-finite values (stat_smooth).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 15 rows containing missing values (geom_smooth).

Analysis of BM

BM all data

fig <- ggplot(data = dataset_exp_0_bm_imprint, mapping = aes(x=label, y=bm))+
  geom_hline(yintercept=0, linetype = "dotted")+
  geom_boxplot()+
  ylim(-400, 600)+
  geom_quasirandom(shape=16, size=3, colour="black")+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
fig

ggsave(plot = fig, filename = "exp_0_bm_label_all.png", dpi = 300, height = 10, width = 7.5, units = "cm")

Two-way ANOVA of BM

m <- lm(bm ~ drug * colour, data = dataset_exp_0_bm_imprint)
anova(m)
## Analysis of Variance Table
## 
## Response: bm
##             Df  Sum Sq Mean Sq F value Pr(>F)
## drug         1  121363  121363  2.1140 0.1542
## colour       1   76288   76288  1.3289 0.2562
## drug:colour  1   28105   28105  0.4896 0.4884
## Residuals   38 2181492   57408

Analysis of imprinting

imprinting all data

fig <- ggplot(data = dataset_exp_0_bm_imprint, mapping = aes(x=label, y=imprint))+
  geom_hline(yintercept=0, linetype = "dotted")+
  geom_boxplot()+
  ylim(-400, 600)+
  geom_quasirandom(shape=16, size=3, colour="black")+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
fig
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
## Warning: Removed 4 rows containing missing values (position_quasirandom).

ggsave(plot = fig, filename = "exp_0_imprint_label_all.png", dpi = 300, height = 10, width = 7.5, units = "cm")
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

## Warning: Removed 4 rows containing missing values (position_quasirandom).

Two-way ANOVA of imprinting

m <- lm(imprint ~ drug * colour, data = dataset_exp_0_bm_imprint)
anova(m)
## Analysis of Variance Table
## 
## Response: imprint
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## drug         1 1252323 1252323 14.0827 0.0005838 ***
## colour       1  856572  856572  9.6324 0.0036004 ** 
## drug:colour  1   10032   10032  0.1128 0.7388093    
## Residuals   38 3379199   88926                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] 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] pwr_1.3-0        readxl_1.3.1     SimComp_3.3      rstatix_0.7.0   
##  [5] car_3.0-12       carData_3.0-5    ggbeeswarm_0.6.0 plotly_4.10.0   
##  [9] forcats_0.5.1    stringr_1.4.0    dplyr_1.0.7      purrr_0.3.4     
## [13] readr_2.1.1      tidyr_1.1.4      tibble_3.1.6     ggplot2_3.3.5   
## [17] tidyverse_1.3.1 
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-144         fs_1.5.2             lubridate_1.8.0     
##  [4] httr_1.4.2           tools_3.6.3          backports_1.4.1     
##  [7] utf8_1.2.2           R6_2.5.1             vipor_0.4.5         
## [10] DBI_1.1.1            lazyeval_0.2.2       mgcv_1.8-31         
## [13] colorspace_2.0-2     withr_2.4.3          tidyselect_1.1.1    
## [16] compiler_3.6.3       survPresmooth_1.1-11 mratios_1.4.2       
## [19] cli_3.1.0            rvest_1.0.2          xml2_1.3.3          
## [22] sandwich_3.0-1       labeling_0.4.2       scales_1.1.1        
## [25] mvtnorm_1.1-3        digest_0.6.29        rmarkdown_2.11      
## [28] pkgconfig_2.0.3      htmltools_0.5.2      highr_0.9           
## [31] dbplyr_2.1.1         fastmap_1.1.0        htmlwidgets_1.5.4   
## [34] rlang_0.4.12         rstudioapi_0.13      farver_2.1.0        
## [37] jquerylib_0.1.4      generics_0.1.1       zoo_1.8-9           
## [40] jsonlite_1.7.2       magrittr_2.0.1       Matrix_1.2-18       
## [43] Rcpp_1.0.7           munsell_0.5.0        fansi_0.5.0         
## [46] 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       
## [52] grid_3.6.3           crayon_1.4.2         lattice_0.20-40     
## [55] haven_2.4.3          splines_3.6.3        hms_1.1.1           
## [58] knitr_1.37           pillar_1.6.4         codetools_0.2-16    
## [61] reprex_2.0.1         glue_1.6.0           evaluate_0.14       
## [64] data.table_1.14.2    modelr_0.1.8         vctrs_0.3.8         
## [67] tzdb_0.2.0           cellranger_1.1.0     gtable_0.3.0        
## [70] assertthat_0.2.1     xfun_0.29            broom_0.7.10        
## [73] survival_3.1-8       viridisLite_0.4.0    beeswarm_0.4.0      
## [76] TH.data_1.1-0        ellipsis_0.3.2