This script performs non-response analyses comparing those included in tissue-dependent subsets to those excluded.
Load packages
library(tidyverse)
library(DNAmArray)
library(lme4)
library(lmerTest)Load data
load("../GOTO_Data/GOTO_targets-filtered.Rdata")
pheno_data <- read_csv('../GOTO_Data/Pheno/GOTO_PhenoData.csv')## Rows: 327 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): visit_date, date_of_birth, age, sex, bmi
## dbl (13): IOP2_ID, timepoint, length, weight, wc, perc_body_fat, fasting_ins...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Pheno variables to numeric
pheno_data$age <- as.numeric(pheno_data$age)
pheno_data$sex <- ifelse(pheno_data$sex == 'male', 0, 1)
pheno_data$bmi <- as.numeric(pheno_data$bmi)Entire GOTO population
before <- pheno_data %>%
filter(timepoint == 0) %>%
arrange(IOP2_ID)
after <- pheno_data %>%
filter(timepoint == 1, IOP2_ID %in% before$IOP2_ID) %>%
arrange(IOP2_ID)
all_ids <- unique(after$IOP2_ID)Muscle
targets_muscle <- targets %>% filter(tissue == 'muscle')
muscle_ids <- as.numeric(as.character(unique(targets_muscle$IOP2_ID)))Adipose
targets_fat <- targets %>% filter(tissue == 'fat')
fat_ids <- as.numeric(as.character(unique(targets_fat$IOP2_ID)))Blood
targets_blood <- targets %>% filter(tissue == 'fasted blood')
blood_ids <- as.numeric(as.character(unique(targets_blood$IOP2_ID)))Overlap
overlap_ids <- muscle_ids[muscle_ids %in% fat_ids]
overlap_ids <- overlap_ids[overlap_ids %in% blood_ids]Traits
traits <- colnames(pheno_data)[9:18]
traits## [1] "bmi" "wc" "perc_body_fat" "fasting_insulin"
## [5] "leptin" "adiponectin" "interleukin_6" "hdl_size_fasted"
## [9] "hdl_c" "systolic_bp"for(i in traits){
change_in <- pheno_data %>%
pivot_wider(id_cols='IOP2_ID',
names_from='timepoint',
values_from=i) %>%
mutate(delta = `1` - `0`)
included <- (change_in %>% filter(IOP2_ID %in% overlap_ids))$delta
excluded <- (change_in %>% filter(!IOP2_ID %in% overlap_ids))$delta
fit <- t.test(included, excluded)
out <- data.frame(trait = i,
delta = fit$estimate[1] - fit$estimate[2],
se = fit$stderr,
p = fit$p.value)
if(i == traits[1]){
res <- out
} else {
res <- rbind(res, out)
}
row.names(res) <- NULL
}Adjust p-values
res$padj <- p.adjust(res$p, method='fdr')
res_full <- resfor(i in traits){
change_in <- pheno_data %>%
pivot_wider(id_cols='IOP2_ID',
names_from='timepoint',
values_from=i) %>%
mutate(delta = `1` - `0`)
included <- (change_in %>% filter(IOP2_ID %in% muscle_ids))$delta
excluded <- (change_in %>% filter(!IOP2_ID %in% muscle_ids))$delta
fit <- t.test(included, excluded)
out <- data.frame(delta = fit$estimate[1] - fit$estimate[2],
se = fit$stderr,
p = fit$p.value)
if(i == traits[1]){
res <- out
} else {
res <- rbind(res, out)
}
row.names(res) <- NULL
}Adjust p-values
res$padj <- p.adjust(res$p, method='fdr')
res_full <- cbind(res_full, res)for(i in traits){
change_in <- pheno_data %>%
pivot_wider(id_cols='IOP2_ID',
names_from='timepoint',
values_from=i) %>%
mutate(delta = `1` - `0`)
included <- (change_in %>% filter(IOP2_ID %in% fat_ids))$delta
excluded <- (change_in %>% filter(!IOP2_ID %in% fat_ids))$delta
fit <- t.test(included, excluded)
out <- data.frame(delta = fit$estimate[1] - fit$estimate[2],
se = fit$stderr,
p = fit$p.value)
if(i == traits[1]){
res <- out
} else {
res <- rbind(res, out)
}
row.names(res) <- NULL
}Adjust p-values
res$padj <- p.adjust(res$p, method='fdr')
res_full <- cbind(res_full, res)for(i in traits){
change_in <- pheno_data %>%
pivot_wider(id_cols='IOP2_ID',
names_from='timepoint',
values_from=i) %>%
mutate(delta = `1` - `0`)
included <- (change_in %>% filter(IOP2_ID %in% blood_ids))$delta
excluded <- (change_in %>% filter(!IOP2_ID %in% blood_ids))$delta
fit <- t.test(included, excluded)
out <- data.frame(delta = fit$estimate[1] - fit$estimate[2],
se = fit$stderr,
p = fit$p.value)
if(i == traits[1]){
res <- out
} else {
res <- rbind(res, out)
}
row.names(res) <- NULL
}Adjust p-values
res$padj <- p.adjust(res$p, method='fdr')
res_full <- cbind(res_full, res)write_csv(res_full, file='../GOTO_Data/Tables/ST02.csv')sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.10 (Green Obsidian)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] lmerTest_3.1-3
## [2] lme4_1.1-30
## [3] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [4] snpStats_1.44.0
## [5] survival_3.5-5
## [6] ggrepel_0.9.1
## [7] ggfortify_0.4.14
## [8] irlba_2.3.5.1
## [9] Matrix_1.5-4.1
## [10] omicsPrint_1.14.0
## [11] MASS_7.3-60
## [12] DNAmArray_2.0.0
## [13] pls_2.8-2
## [14] FDb.InfiniumMethylation.hg19_2.2.0
## [15] org.Hs.eg.db_3.14.0
## [16] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [17] GenomicFeatures_1.46.5
## [18] AnnotationDbi_1.56.2
## [19] IlluminaHumanMethylationEPICmanifest_0.3.0
## [20] minfi_1.40.0
## [21] bumphunter_1.36.0
## [22] locfit_1.5-9.8
## [23] iterators_1.0.14
## [24] foreach_1.5.2
## [25] Biostrings_2.62.0
## [26] XVector_0.34.0
## [27] SummarizedExperiment_1.24.0
## [28] Biobase_2.58.0
## [29] MatrixGenerics_1.10.0
## [30] matrixStats_1.0.0
## [31] GenomicRanges_1.46.1
## [32] GenomeInfoDb_1.34.9
## [33] IRanges_2.32.0
## [34] S4Vectors_0.36.2
## [35] BiocGenerics_0.44.0
## [36] BiocParallel_1.32.6
## [37] MethylAid_1.28.0
## [38] forcats_0.5.2
## [39] stringr_1.5.0
## [40] dplyr_1.1.3
## [41] purrr_0.3.4
## [42] readr_2.1.2
## [43] tidyr_1.2.1
## [44] tibble_3.2.1
## [45] ggplot2_3.4.3
## [46] tidyverse_1.3.2
## [47] rmarkdown_2.16
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.3 tidyselect_1.2.0
## [3] RSQLite_2.2.17 grid_4.2.2
## [5] munsell_0.5.0 codetools_0.2-19
## [7] preprocessCore_1.60.2 withr_2.5.0
## [9] colorspace_2.1-0 filelock_1.0.2
## [11] highr_0.10 knitr_1.43
## [13] rstudioapi_0.14 labeling_0.4.2
## [15] GenomeInfoDbData_1.2.9 farver_2.1.1
## [17] bit64_4.0.5 rhdf5_2.42.1
## [19] vctrs_0.6.3 generics_0.1.3
## [21] xfun_0.39 timechange_0.2.0
## [23] BiocFileCache_2.2.1 R6_2.5.1
## [25] illuminaio_0.40.0 bitops_1.0-7
## [27] rhdf5filters_1.10.1 cachem_1.0.8
## [29] reshape_0.8.9 DelayedArray_0.24.0
## [31] assertthat_0.2.1 vroom_1.5.7
## [33] promises_1.2.0.1 BiocIO_1.8.0
## [35] scales_1.2.1 googlesheets4_1.0.1
## [37] gtable_0.3.3 rlang_1.1.1
## [39] genefilter_1.76.0 splines_4.2.2
## [41] rtracklayer_1.54.0 gargle_1.5.0
## [43] GEOquery_2.62.2 htm2txt_2.2.2
## [45] hexbin_1.28.3 broom_1.0.1
## [47] yaml_2.3.7 reshape2_1.4.4
## [49] RaggedExperiment_1.18.0 modelr_0.1.9
## [51] backports_1.4.1 httpuv_1.6.11
## [53] tools_4.2.2 gridBase_0.4-7
## [55] nor1mix_1.3-0 ellipsis_0.3.2
## [57] jquerylib_0.1.4 RColorBrewer_1.1-3
## [59] siggenes_1.68.0 MultiAssayExperiment_1.20.0
## [61] Rcpp_1.0.10 plyr_1.8.8
## [63] sparseMatrixStats_1.10.0 progress_1.2.2
## [65] zlibbioc_1.44.0 RCurl_1.98-1.12
## [67] prettyunits_1.1.1 openssl_2.0.6
## [69] haven_2.5.1 fs_1.6.2
## [71] magrittr_2.0.3 data.table_1.14.8
## [73] reprex_2.0.2 googledrive_2.0.0
## [75] hms_1.1.2 mime_0.12
## [77] evaluate_0.21 xtable_1.8-4
## [79] XML_3.99-0.14 mclust_6.0.0
## [81] readxl_1.4.1 gridExtra_2.3
## [83] compiler_4.2.2 biomaRt_2.50.3
## [85] crayon_1.5.2 minqa_1.2.5
## [87] htmltools_0.5.5 later_1.3.1
## [89] tzdb_0.4.0 lubridate_1.9.2
## [91] DBI_1.1.3 dbplyr_2.2.1
## [93] rappdirs_0.3.3 boot_1.3-28.1
## [95] cli_3.6.1 quadprog_1.5-8
## [97] pkgconfig_2.0.3 GenomicAlignments_1.30.0
## [99] numDeriv_2016.8-1.1 xml2_1.3.4
## [101] annotate_1.72.0 bslib_0.5.0
## [103] rngtools_1.5.2 multtest_2.50.0
## [105] beanplot_1.3.1 rvest_1.0.3
## [107] doRNG_1.8.6 scrime_1.3.5
## [109] digest_0.6.31 base64_2.0.1
## [111] cellranger_1.1.0 DelayedMatrixStats_1.16.0
## [113] restfulr_0.0.15 curl_5.0.1
## [115] shiny_1.7.2 Rsamtools_2.10.0
## [117] nloptr_2.0.3 rjson_0.2.21
## [119] lifecycle_1.0.3 nlme_3.1-162
## [121] jsonlite_1.8.5 Rhdf5lib_1.20.0
## [123] askpass_1.1 limma_3.54.2
## [125] fansi_1.0.4 pillar_1.9.0
## [127] lattice_0.21-8 KEGGREST_1.34.0
## [129] fastmap_1.1.1 httr_1.4.6
## [131] glue_1.6.2 png_0.1-8
## [133] bit_4.0.5 stringi_1.7.12
## [135] sass_0.4.6 HDF5Array_1.22.1
## [137] blob_1.2.4 memoise_2.0.1Clear
rm(list=ls())