Last updated: 2021-10-29
Checks: 7 0
Knit directory: komputeExamples/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20211027)
was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 46485bd. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish
or wflow_git_commit
). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Unstaged changes:
Modified: README.md
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/kompute_test_10212021_BC.Rmd
) and HTML (docs/kompute_test_10212021_BC.html
) files. If you’ve configured a remote Git repository (see ?wflow_git_remote
), click on the hyperlinks in the table below to view the files as they were in that past version.
File | Version | Author | Date | Message |
---|---|---|---|---|
Rmd | 46485bd | Coby Warkentin | 2021-10-29 | wflow_publish("analysis/kompute_test_10212021_BC.Rmd") |
Rmd | b7ebb1f | GitHub | 2021-10-29 | Add files via upload |
rm(list=ls())
library(data.table)
library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:data.table':
between, first, last
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
#library(kableExtra)
library(reshape2)
Attaching package: 'reshape2'
The following objects are masked from 'package:data.table':
dcast, melt
library(ggplot2)
library(tidyr) #spread
Attaching package: 'tidyr'
The following object is masked from 'package:reshape2':
smiths
#library(pheatmap)
library(RColorBrewer)
#library(GGally) #ggpairs
library(irlba) # partial PCA
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
#library(gridExtra)
#library(cowplot)
library(circlize)
========================================
circlize version 0.4.13
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/
If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
in R. Bioinformatics 2014.
This message can be suppressed by:
suppressPackageStartupMessages(library(circlize))
========================================
library(ComplexHeatmap)
Loading required package: grid
========================================
ComplexHeatmap version 2.10.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
genomic data. Bioinformatics 2016.
The new InteractiveComplexHeatmap package can directly export static
complex heatmaps into an interactive Shiny app with zero effort. Have a try!
This message can be suppressed by:
suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
#options(max.print = 3000)
# load(file="~/Google Drive Miami/Miami_IMPC/data/v10.1/AllControls_small.Rdata")
load(file="G:/.shortcut-targets-by-id/1SeBOMb4GZ2Gkldxp4QNEnFWHOiAqtRTz/Miami_IMPC/data/v10.1/AllControls_small.Rdata")
dim(allpheno)
[1] 12025349 15
head(allpheno)
biological_sample_id procedure_stable_id
1 242736 IMPC_ACS_003
2 155527 IMPC_GRS_001
3 72957 IMPC_OFD_001
4 183476 IMPC_GRS_001
5 81327 IMPC_OFD_001
6 226658 IMPC_ACS_003
procedure_name parameter_stable_id
1 Acoustic Startle and Pre-pulse Inhibition (PPI) IMPC_ACS_036_001
2 Grip Strength IMPC_GRS_008_001
3 Open Field IMPC_OFD_020_001
4 Grip Strength IMPC_GRS_009_001
5 Open Field IMPC_OFD_007_001
6 Acoustic Startle and Pre-pulse Inhibition (PPI) IMPC_ACS_034_001
parameter_name data_point sex
1 % Pre-pulse inhibition - PPI4 88.0373 female
2 Forelimb grip strength measurement mean 75.3333 female
3 Distance travelled - total 6130.3000 female
4 Forelimb and hindlimb grip strength measurement mean 120.3970 male
5 Whole arena resting time 607.0000 female
6 % Pre-pulse inhibition - PPI2 68.0742 female
age_in_weeks weight phenotyping_center date_of_experiment strain_name
1 9 18.15 JAX 2014-08-18T00:00:00Z C57BL/6NJ
2 9 20.00 BCM 2018-02-20T00:00:00Z C57BL/6N
3 9 21.20 RBRC 2015-08-03T00:00:00Z C57BL/6NTac
4 9 28.40 WTSI 2015-03-02T00:00:00Z C57BL/6N
5 9 22.62 KMPC 2018-01-30T00:00:00Z C57BL/6NTac
6 10 18.54 JAX 2015-12-16T00:00:00Z C57BL/6NJ
developmental_stage_name observation_type data_type
1 Early adult unidimensional FLOAT
2 Early adult unidimensional FLOAT
3 Early adult unidimensional FLOAT
4 Early adult unidimensional FLOAT
5 Early adult unidimensional FLOAT
6 Early adult unidimensional FLOAT
allpheno.org <- allpheno
#allpheno <- allpheno.org
We use BC only.
allpheno = allpheno %>%
filter(procedure_name=="Body Composition (DEXA lean/fat)") %>%
mutate(proc_short_name=recode(procedure_name, "Body Composition (DEXA lean/fat)"="BC")) %>%
#mutate(parameter_name=recode(parameter_name, "Triglyceride"="Triglycerides")) %>%
mutate(proc_param_name=paste0(proc_short_name,"_",parameter_name)) %>%
mutate(proc_param_name_stable_id=paste0(proc_short_name,"_",parameter_name,"_",parameter_stable_id))
## Extract time series data and find out parameter names
ts <- allpheno %>% filter(observation_type=="time_series")
table(ts$proc_param_name)
< table of extent 0 >
# Filter out time series data
allpheno <- allpheno %>% filter(observation_type!="time_series")
table(allpheno$proc_param_name)
BC_BMC/Body weight
24446
BC_Body length
19233
BC_Body weight
25879
BC_Bone Area
24461
BC_Bone Mineral Content (excluding skull)
24918
BC_Bone Mineral Density (excluding skull)
24919
BC_Fat mass
25154
BC_Fat/Body weight
24682
BC_Lean mass
25156
BC_Lean/Body weight
24684
This heatmaps show phenotypes measured for each control mouse. Columns represent mice and rows represent phenotypes.
mtest <- table(allpheno$proc_param_name_stable_id, allpheno$biological_sample_id)
mtest <-as.data.frame.matrix(mtest)
dim(mtest)
[1] 51 25690
if(FALSE){
nmax <-max(mtest)
library(circlize)
col_fun = colorRamp2(c(0, nmax), c("white", "red"))
col_fun(seq(0, nmax))
pdf("~/Google Drive Miami/Miami_IMPC/output/measured_phenotypes_controls_BC.pdf", width = 12, height = 14)
ht = Heatmap(as.matrix(mtest), cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, col = col_fun,
row_names_gp = gpar(fontsize = 8), name="Count")
draw(ht)
dev.off()
}
mtest <- table(allpheno$proc_param_name, allpheno$biological_sample_id)
dim(mtest)
[1] 10 25690
head(mtest[,1:10])
39638 39639 39640 39641 39642 39643
BC_BMC/Body weight 1 1 1 1 1 1
BC_Body length 0 0 0 1 0 1
BC_Body weight 2 1 1 1 1 1
BC_Bone Area 1 1 1 1 1 1
BC_Bone Mineral Content (excluding skull) 2 1 1 1 1 1
BC_Bone Mineral Density (excluding skull) 2 1 1 1 1 1
39650 39651 39652 39657
BC_BMC/Body weight 1 1 1 1
BC_Body length 1 0 0 1
BC_Body weight 1 1 1 1
BC_Bone Area 1 1 1 1
BC_Bone Mineral Content (excluding skull) 1 1 1 1
BC_Bone Mineral Density (excluding skull) 1 1 1 1
mtest0 <- mtest>0
head(mtest0[,1:10])
39638 39639 39640 39641 39642 39643
BC_BMC/Body weight TRUE TRUE TRUE TRUE TRUE TRUE
BC_Body length FALSE FALSE FALSE TRUE FALSE TRUE
BC_Body weight TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Area TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Mineral Content (excluding skull) TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Mineral Density (excluding skull) TRUE TRUE TRUE TRUE TRUE TRUE
39650 39651 39652 39657
BC_BMC/Body weight TRUE TRUE TRUE TRUE
BC_Body length TRUE FALSE FALSE TRUE
BC_Body weight TRUE TRUE TRUE TRUE
BC_Bone Area TRUE TRUE TRUE TRUE
BC_Bone Mineral Content (excluding skull) TRUE TRUE TRUE TRUE
BC_Bone Mineral Density (excluding skull) TRUE TRUE TRUE TRUE
rowSums(mtest0)
BC_BMC/Body weight
24446
BC_Body length
19233
BC_Body weight
25675
BC_Bone Area
24461
BC_Bone Mineral Content (excluding skull)
24714
BC_Bone Mineral Density (excluding skull)
24715
BC_Fat mass
24950
BC_Fat/Body weight
24682
BC_Lean mass
24952
BC_Lean/Body weight
24684
rmv.pheno.list <- rownames(mtest)[rowSums(mtest0)<15000]
rmv.pheno.list
character(0)
dim(allpheno)
[1] 243532 18
allpheno <- allpheno %>% filter(!(proc_param_name %in% rmv.pheno.list))
dim(allpheno)
[1] 243532 18
# number of phenotypes left
length(unique(allpheno$proc_param_name))
[1] 10
mtest <- table(allpheno$proc_param_name, allpheno$biological_sample_id)
dim(mtest)
[1] 10 25690
head(mtest[,1:10])
39638 39639 39640 39641 39642 39643
BC_BMC/Body weight 1 1 1 1 1 1
BC_Body length 0 0 0 1 0 1
BC_Body weight 2 1 1 1 1 1
BC_Bone Area 1 1 1 1 1 1
BC_Bone Mineral Content (excluding skull) 2 1 1 1 1 1
BC_Bone Mineral Density (excluding skull) 2 1 1 1 1 1
39650 39651 39652 39657
BC_BMC/Body weight 1 1 1 1
BC_Body length 1 0 0 1
BC_Body weight 1 1 1 1
BC_Bone Area 1 1 1 1
BC_Bone Mineral Content (excluding skull) 1 1 1 1
BC_Bone Mineral Density (excluding skull) 1 1 1 1
mtest0 <- mtest>0
head(mtest0[,1:10])
39638 39639 39640 39641 39642 39643
BC_BMC/Body weight TRUE TRUE TRUE TRUE TRUE TRUE
BC_Body length FALSE FALSE FALSE TRUE FALSE TRUE
BC_Body weight TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Area TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Mineral Content (excluding skull) TRUE TRUE TRUE TRUE TRUE TRUE
BC_Bone Mineral Density (excluding skull) TRUE TRUE TRUE TRUE TRUE TRUE
39650 39651 39652 39657
BC_BMC/Body weight TRUE TRUE TRUE TRUE
BC_Body length TRUE FALSE FALSE TRUE
BC_Body weight TRUE TRUE TRUE TRUE
BC_Bone Area TRUE TRUE TRUE TRUE
BC_Bone Mineral Content (excluding skull) TRUE TRUE TRUE TRUE
BC_Bone Mineral Density (excluding skull) TRUE TRUE TRUE TRUE
summary(colSums(mtest0))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 9.00 10.00 9.44 10.00 10.00
rmv.sample.list <- colnames(mtest)[colSums(mtest0)<7]
length(rmv.sample.list)
[1] 1829
dim(allpheno)
[1] 243532 18
allpheno <- allpheno %>% filter(!(biological_sample_id %in% rmv.sample.list))
dim(allpheno)
[1] 234305 18
# number of observations to use
length(unique(allpheno$biological_sample_id))
[1] 23861
if(FALSE){
mtest <- table(allpheno$proc_param_name, allpheno$biological_sample_id)
dim(mtest)
mtest <-as.data.frame.matrix(mtest)
nmax <-max(mtest)
library(circlize)
col_fun = colorRamp2(c(0, nmax), c("white", "red"))
col_fun(seq(0, nmax))
pdf("~/Google Drive Miami/Miami_IMPC/output/measured_phenotypes_controls_after_filtering_BC.pdf", width = 10, height = 3)
ht = Heatmap(as.matrix(mtest), cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, col = col_fun,
row_names_gp = gpar(fontsize = 7), name="Count")
draw(ht)
dev.off()
}
ap.mat <- allpheno %>%
dplyr::select(biological_sample_id, proc_param_name, data_point, sex, phenotyping_center, strain_name) %>%
##consider weight or age in weeks
arrange(biological_sample_id) %>%
distinct(biological_sample_id, proc_param_name, .keep_all=TRUE) %>% ## remove duplicates, maybe mean() is better.
spread(proc_param_name, data_point) %>%
tibble::column_to_rownames(var="biological_sample_id")
head(ap.mat)
sex phenotyping_center strain_name BC_BMC/Body weight BC_Body length
39638 female MRC Harwell C57BL/6NTac 0.0567631 NA
39639 female HMGU C57BL/6NCrl 0.0224897 NA
39640 female HMGU C57BL/6NTac 0.0214276 NA
39641 male HMGU C57BL/6NCrl 0.0191929 9.7
39642 female MRC Harwell C57BL/6NTac 0.0242145 NA
39643 female HMGU C57BL/6NCrl 0.0224004 9.4
BC_Body weight BC_Bone Area BC_Bone Mineral Content (excluding skull)
39638 22.12 13.51560 1.25560
39639 23.30 8.11161 0.52401
39640 23.20 6.85683 0.49712
39641 29.70 8.61072 0.57003
39642 25.27 8.42837 0.61190
39643 24.30 8.42616 0.54433
BC_Bone Mineral Density (excluding skull) BC_Fat mass BC_Fat/Body weight
39638 0.0929 6.9362 0.313571
39639 0.0646 6.1000 0.261803
39640 0.0725 3.3000 0.142241
39641 0.0662 7.1000 0.239057
39642 0.0726 3.4382 0.136059
39643 0.0646 6.9000 0.283951
BC_Lean mass BC_Lean/Body weight
39638 11.2569 0.508901
39639 14.7000 0.630901
39640 17.1000 0.737069
39641 21.7000 0.730640
39642 20.5392 0.812790
39643 15.3000 0.629630
dim(ap.mat)
[1] 23861 13
summary(colSums(is.na(ap.mat[,-1:-3])))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 0.00 0.00 532.50 0.75 5323.00
ggplot(melt(ap.mat), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
Using sex, phenotyping_center, strain_name as id variables
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 5325 rows containing non-finite values (stat_bin).
library(RNOmni)
ap.mat.rank <- ap.mat
dim(ap.mat.rank)
[1] 23861 13
ap.mat.rank <- ap.mat.rank[complete.cases(ap.mat.rank),]
dim(ap.mat.rank)
[1] 18538 13
dim(ap.mat)
[1] 23861 13
ap.mat <- ap.mat[complete.cases(ap.mat),]
dim(ap.mat)
[1] 18538 13
#rankZ <- function(x){ qnorm((rank(x,na.last="keep")-0.5)/sum(!is.na(x))) }
#ap.mat.rank <- ap.mat
#dim(ap.mat.rank)
#ap.mat.rank <- ap.mat.rank[complete.cases(ap.mat.rank),]
#dim(ap.mat.rank)
#library(RNOmni)
#ap.mat.rank <- cbind(ap.mat.rank[,1:3], apply(ap.mat.rank[,-1:-3], 2, RankNorm))
ap.mat.rank <- cbind(ap.mat.rank[,1:3], apply(ap.mat.rank[,-1:-3], 2, RankNorm))
ggplot(melt(ap.mat.rank), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
Using sex, phenotyping_center, strain_name as id variables
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Please conduct a PVCA analysis on the phenotype matrix data (op.mat[,-1:-3]). I think you can measure the proportion of variance explained by each important covariate (sex, phenotyping_center, strain_name)
#source("~/Google Drive Miami/Miami_IMPC/reference/PVCA/examples/PVCA.R")
source("G:/.shortcut-targets-by-id/1SeBOMb4GZ2Gkldxp4QNEnFWHOiAqtRTz/Miami_IMPC/reference/PVCA/examples/PVCA.R")
meta <- ap.mat.rank[,1:3] ## looking at covariates sex, phenotyping_center, and strain_name
head(meta)
sex phenotyping_center strain_name
39641 male HMGU C57BL/6NCrl
39643 female HMGU C57BL/6NCrl
39657 male HMGU C57BL/6NCrl
39750 female HMGU C57BL/6NCrl
39763 female HMGU C57BL/6NCrl
39773 female HMGU C57BL/6NCrl
dim(meta)
[1] 18538 3
summary(meta) # variables are still characters
sex phenotyping_center strain_name
Length:18538 Length:18538 Length:18538
Class :character Class :character Class :character
Mode :character Mode :character Mode :character
meta[sapply(meta, is.character)] <- lapply(meta[sapply(meta, is.character)], as.factor)
summary(meta) # now all variables are converted to factors
sex phenotyping_center strain_name
female:9282 WTSI :4573 B6Brd;B6Dnk;B6N-Tyr<c-Brd>: 176
male :9256 JAX :3675 C57BL/6N :9127
ICS :2075 C57BL/6NCrl :2997
BCM :1832 C57BL/6NJ :3675
RBRC :1692 C57BL/6NJcl : 728
UC Davis:1441 C57BL/6NTac :1835
(Other) :3250
chisq.test(meta[,1],meta[,2])
Pearson's Chi-squared test
data: meta[, 1] and meta[, 2]
X-squared = 13.572, df = 10, p-value = 0.1934
chisq.test(meta[,2],meta[,3])
Warning in chisq.test(meta[, 2], meta[, 3]): Chi-squared approximation may be
incorrect
Pearson's Chi-squared test
data: meta[, 2] and meta[, 3]
X-squared = 59770, df = 50, p-value < 2.2e-16
meta<-meta[,-3] # phenotyping_center and strain_name strongly associated and this caused confouding in PVCA analysis so strain_name dropped.
G <- t(ap.mat.rank[,-1:-3]) ## phenotype matrix data
set.seed(09302021)
# Perform PVCA for 10 random samples of size 1000 (more computationally efficient)
pvca.res <- matrix(nrow=10, ncol=3)
for (i in 1:10){
sample <- sample(1:ncol(G), 1000, replace=FALSE)
pvca.res[i,] <- PVCA(G[,sample], meta[sample,], threshold=0.6, inter=FALSE)
}
# Average effect size across samples
pvca.means <- colMeans(pvca.res)
names(pvca.means) <- c(colnames(meta), "resid")
# Plot PVCA
pvca.plot <- PlotPVCA(pvca.means, "PVCA of Phenotype Matrix Data")
pvca.plot
#ggsave(filename = "pvca_plot.png", pvca.plot, width=8, height=6)
If a large proportion of variance is explained by these covariats, we need to remove their effects from the data.
library(sva)
Loading required package: mgcv
Loading required package: nlme
Attaching package: 'nlme'
The following object is masked from 'package:lme4':
lmList
The following object is masked from 'package:dplyr':
collapse
This is mgcv 1.8-36. For overview type 'help("mgcv-package")'.
Loading required package: genefilter
Attaching package: 'genefilter'
The following object is masked from 'package:ComplexHeatmap':
dist2
Loading required package: BiocParallel
combat_komp = ComBat(dat=G, batch=meta$phenotyping_center, par.prior=TRUE, prior.plots=TRUE, mod=NULL)
Found11batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data
combat_komp[,1:3]
39641 39643 39657
BC_BMC/Body weight -0.54465772 0.76242613 -0.8478037
BC_Body length 0.23714642 -0.15635195 0.2371464
BC_Body weight 0.90694779 -0.38659680 0.9844317
BC_Bone Area 0.41881733 0.15104262 1.0646332
BC_Bone Mineral Content (excluding skull) 0.58720884 0.15610494 0.3622582
BC_Bone Mineral Density (excluding skull) 0.37590404 0.04804937 -0.5216066
BC_Fat mass 0.06021668 -0.01420044 0.9679394
BC_Fat/Body weight -0.38016638 0.16405541 0.6313684
BC_Lean mass 1.24681829 -0.49595622 0.1077335
BC_Lean/Body weight 0.72009214 -0.27108379 -0.6575397
G[,1:3] # for comparison, combat_komp is same form and same dimensions as G
39641 39643 39657
BC_BMC/Body weight 0.61537676 1.3908581 0.43552296
BC_Body length 0.04917004 -0.1811585 0.04917004
BC_Body weight 0.72940618 -0.3554534 0.79438976
BC_Bone Area -0.05160698 -0.2731793 0.48277848
BC_Bone Mineral Content (excluding skull) 1.43981495 1.1313384 1.27885140
BC_Bone Mineral Density (excluding skull) 1.37276777 1.2301457 0.98233626
BC_Fat mass 0.70042463 0.6502323 1.31265922
BC_Fat/Body weight 0.54873676 0.9724173 1.33622379
BC_Lean mass 0.64897996 -1.1933776 -0.55519275
BC_Lean/Body weight -0.13426402 -1.1704055 -1.57439333
set.seed(09302021)
# Perform PVCA for 10 samples (more computationally efficient)
pvca.res.nobatch <- matrix(nrow=10, ncol=3)
for (i in 1:10){
sample <- sample(1:ncol(combat_komp), 1000, replace=FALSE)
pvca.res.nobatch[i,] <- PVCA(combat_komp[,sample], meta[sample,], threshold=0.6, inter=FALSE)
}
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
# Average effect size across samples
pvca.means.nobatch <- colMeans(pvca.res.nobatch)
names(pvca.means.nobatch) <- c(colnames(meta), "resid")
# Plot PVCA
pvca.plot.nobatch <- PlotPVCA(pvca.means.nobatch, "PVCA of Phenotype Matrix Data with Reduced Batch Effect")
pvca.plot.nobatch
The batch effect of phenotyping center has been largely removed.
ap.cor.rank <- cor(ap.mat.rank[,-1:-3], use="pairwise.complete.obs") # pearson correlation coefficient
#ap.cor <- cor(ap.mat[,-1:-3], use="pairwise.complete.obs") # pearson correlation coefficient
ap.cor <- cor(ap.mat[,-1:-3], use="pairwise.complete.obs", method="spearman")
ap.cor.combat <- cor(t(combat_komp), use="pairwise.complete.obs")
#ap.cor <- cor(ap.mat[,-1:-3], use="pairwise.complete.obs", method="spearman") # use original phenotype data
#ap.cor <- cor(ap.mat.rank[,-1:-3], use="pairwise.complete.obs", method="spearman") # use rankZ transformed phenotype data
#col <- colorRampPalette(c("steelblue", "white", "darkorange"))(100)
#ap.cor.out <- pheatmap(ap.cor, cluster_rows = T, cluster_cols=T, show_colnames=F, col=col, fontsize = 7)
#col <- colorRampPalette(c("white","darkorange"))(100)
#pheatmap(abs(op.cor), cluster_rows = T, cluster_cols=T, show_colnames=F, col=col)
if(FALSE){
#pdf("~/Google Drive Miami/Miami_IMPC/output/genetic_corr_btw_phenotypes.pdf", width = 11, height = 8)
pdf("~/Google Drive Miami/Miami_IMPC/output/genetic_corr_btw_phenotypes_Pearson_BC.pdf", width = 4.9, height = 2)
ht = Heatmap(ap.cor, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Corr")
draw(ht)
dev.off()
pdf("~/Google Drive Miami/Miami_IMPC/output/genetic_corr_btw_phenotypes_rankZ_BC.pdf", width = 4.9, height = 2)
ht = Heatmap(ap.cor.rank, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Corr")
draw(ht)
dev.off()
pdf("~/Google Drive Miami/Miami_IMPC/output/genetic_corr_btw_phenotypes_ComBat_Adjusted_BC.pdf", width = 4.9, height = 2)
ht = Heatmap(ap.cor.combat, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Corr")
draw(ht)
dev.off()
}
pheno.list <- rownames(ap.cor)
# KOMPv10.1.file = "~/Google Drive Miami/Miami_IMPC/data/v10.1/IMPC_ALL_statistical_results.csv.gz"
KOMPv10.1.file = "G:/.shortcut-targets-by-id/1SeBOMb4GZ2Gkldxp4QNEnFWHOiAqtRTz/Miami_IMPC/data/v10.1/IMPC_ALL_statistical_results.csv.gz"
KOMPv10.1 = fread(KOMPv10.1.file, header=TRUE, sep=",")
KOMPv10.1$parameter_name <- trimws(KOMPv10.1$parameter_name) #remove white spaces
KOMPv10.1$proc_param_name <- paste0(KOMPv10.1$procedure_name,"_",KOMPv10.1$parameter_name)
#head(KOMPv10.1, 10)
#sort(table(KOMPv10.1$procedure_name))
#sort(table(KOMPv10.1$proc_param_name), decreasing = TRUE)[1:100]
#sort(table(KOMPv10.1$procedure_name))
#table(KOMPv10.1$procedure_name, KOMPv10.1$parameter_name)
#table(KOMPv10.1$procedure_name, KOMPv10.1$statistical_method)
table(KOMPv10.1$procedure_name, KOMPv10.1$data_type)
adult-gross-path
Acoustic Startle and Pre-pulse Inhibition (PPI) 0
Acoustic Startle&PPI 0
Allergy (GMC) 0
Anti-nuclear antibody assay 0
Antigen Specific Immunoglobulin Assay 0
Auditory Brain Stem Response 0
Body Composition (DEXA lean/fat) 0
Body Weight 0
Bodyweight (GMC) 0
Bone marrow immunophenotyping 0
Buffy coat peripheral blood leukocyte immunophenotyping 0
Calorimetry 0
Challenge Whole Body Plethysmography 0
Clinical Chemistry 0
Clinical chemistry (GMC) 0
Combined SHIRPA and Dysmorphology 0
Cortical Bone MicroCT 0
DEXA 0
Dexa-scan analysis 0
DSS Histology 0
Dysmorphology 0
Ear epidermis immunophenotyping 0
ECG (Electrocardiogram) (GMC) 0
Echo 0
Electrocardiogram (ECG) 0
Electroconvulsive Threshold Testing 0
Electroretinography 0
ELISA (GMC) 0
ERG (Electroretinogram) (GMC) 0
Eye Morphology 0
Eye size (GMC) 0
FACS (GMC) 0
FACs Analysis 0
Fasted Clinical Chemistry 0
Fear Conditioning 0
Femoral Microradiography 0
Fertility of Homozygous Knock-out Mice 0
Food efficiency (GMC) 0
Grip-Strength 0
Grip Strength 0
Grip Strength (GMC) 0
Gross Morphology Embryo E12.5 0
Gross Morphology Embryo E14.5-E15.5 0
Gross Morphology Embryo E18.5 0
Gross Morphology Embryo E9.5 0
Gross Morphology Placenta E12.5 0
Gross Morphology Placenta E14.5-E15.5 0
Gross Morphology Placenta E18.5 0
Gross Morphology Placenta E9.5 0
Gross Pathology and Tissue Collection 633844
Haematology 0
Haematology (GMC) 0
Haematology test 0
Heart Dissection 0
Heart Weight 0
Heart weight/tibia length 0
Hematology 0
Hole-board Exploration 0
Holeboard (GMC) 0
Hot Plate 0
Immunoglobulin 0
Indirect Calorimetry 0
Indirect ophthalmoscopy 0
Insulin Blood Level 0
Intraperitoneal glucose tolerance test (IPGTT) 0
IPGTT 0
Light-Dark Test 0
Mesenteric Lymph Node Immunophenotyping 0
Modified SHIRPA 0
Nociception Hotplate (GMC) 0
Open-field 0
Open Field 0
Ophthalmoscope 0
Organ Weight 0
pDexa (GMC) 0
Plasma Chemistry 0
Rotarod 0
Rotarod A (GMC) 0
Shirpa (GMC) 0
Simplified IPGTT 0
Sleep Wake 0
Slit Lamp 0
Spleen Immunophenotyping 0
Spontaneous breathing (GMC) 0
Tail Suspension 0
Three-point Bend 0
Trabecular Bone MicroCT 0
Trichuris 0
Urinalysis 0
Vertebra Compression 0
Vertebral Microradiography 0
Viability E12.5 Secondary Screen 0
Viability E14.5-E15.5 Secondary Screen 0
Viability E18.5 Secondary Screen 0
Viability E9.5 Secondary Screen 0
Viability Primary Screen 0
Whole blood peripheral blood leukocyte immunophenotyping 0
X-ray 0
X-Ray 0
X-Ray (GMC) 0
categorical embryo
Acoustic Startle and Pre-pulse Inhibition (PPI) 0 0
Acoustic Startle&PPI 0 0
Allergy (GMC) 0 0
Anti-nuclear antibody assay 984 0
Antigen Specific Immunoglobulin Assay 0 0
Auditory Brain Stem Response 0 0
Body Composition (DEXA lean/fat) 0 0
Body Weight 0 0
Bodyweight (GMC) 0 0
Bone marrow immunophenotyping 0 0
Buffy coat peripheral blood leukocyte immunophenotyping 0 0
Calorimetry 0 0
Challenge Whole Body Plethysmography 0 0
Clinical Chemistry 0 0
Clinical chemistry (GMC) 0 0
Combined SHIRPA and Dysmorphology 156103 0
Cortical Bone MicroCT 0 0
DEXA 0 0
Dexa-scan analysis 0 0
DSS Histology 21 0
Dysmorphology 17289 0
Ear epidermis immunophenotyping 0 0
ECG (Electrocardiogram) (GMC) 0 0
Echo 0 0
Electrocardiogram (ECG) 0 0
Electroconvulsive Threshold Testing 0 0
Electroretinography 0 0
ELISA (GMC) 0 0
ERG (Electroretinogram) (GMC) 0 0
Eye Morphology 85613 0
Eye size (GMC) 0 0
FACS (GMC) 0 0
FACs Analysis 0 0
Fasted Clinical Chemistry 0 0
Fear Conditioning 0 0
Femoral Microradiography 0 0
Fertility of Homozygous Knock-out Mice 0 0
Food efficiency (GMC) 0 0
Grip-Strength 0 0
Grip Strength 0 0
Grip Strength (GMC) 0 0
Gross Morphology Embryo E12.5 0 40122
Gross Morphology Embryo E14.5-E15.5 0 13670
Gross Morphology Embryo E18.5 0 21665
Gross Morphology Embryo E9.5 0 23135
Gross Morphology Placenta E12.5 0 4623
Gross Morphology Placenta E14.5-E15.5 0 2685
Gross Morphology Placenta E18.5 0 2976
Gross Morphology Placenta E9.5 0 3827
Gross Pathology and Tissue Collection 0 0
Haematology 0 0
Haematology (GMC) 0 0
Haematology test 0 0
Heart Dissection 570 0
Heart Weight 0 0
Heart weight/tibia length 103 0
Hematology 0 0
Hole-board Exploration 0 0
Holeboard (GMC) 0 0
Hot Plate 0 0
Immunoglobulin 0 0
Indirect Calorimetry 0 0
Indirect ophthalmoscopy 3040 0
Insulin Blood Level 0 0
Intraperitoneal glucose tolerance test (IPGTT) 0 0
IPGTT 0 0
Light-Dark Test 0 0
Mesenteric Lymph Node Immunophenotyping 0 0
Modified SHIRPA 20876 0
Nociception Hotplate (GMC) 0 0
Open-field 0 0
Open Field 0 0
Ophthalmoscope 8246 0
Organ Weight 0 0
pDexa (GMC) 0 0
Plasma Chemistry 0 0
Rotarod 0 0
Rotarod A (GMC) 0 0
Shirpa (GMC) 0 0
Simplified IPGTT 0 0
Sleep Wake 0 0
Slit Lamp 14988 0
Spleen Immunophenotyping 0 0
Spontaneous breathing (GMC) 0 0
Tail Suspension 1 0
Three-point Bend 0 0
Trabecular Bone MicroCT 0 0
Trichuris 10 0
Urinalysis 0 0
Vertebra Compression 0 0
Vertebral Microradiography 0 0
Viability E12.5 Secondary Screen 0 422
Viability E14.5-E15.5 Secondary Screen 0 210
Viability E18.5 Secondary Screen 0 176
Viability E9.5 Secondary Screen 0 316
Viability Primary Screen 0 0
Whole blood peripheral blood leukocyte immunophenotyping 0 0
X-ray 60705 0
X-Ray 1791 0
X-Ray (GMC) 47 0
line
Acoustic Startle and Pre-pulse Inhibition (PPI) 0
Acoustic Startle&PPI 0
Allergy (GMC) 0
Anti-nuclear antibody assay 0
Antigen Specific Immunoglobulin Assay 0
Auditory Brain Stem Response 0
Body Composition (DEXA lean/fat) 0
Body Weight 0
Bodyweight (GMC) 0
Bone marrow immunophenotyping 0
Buffy coat peripheral blood leukocyte immunophenotyping 0
Calorimetry 0
Challenge Whole Body Plethysmography 0
Clinical Chemistry 0
Clinical chemistry (GMC) 0
Combined SHIRPA and Dysmorphology 0
Cortical Bone MicroCT 0
DEXA 0
Dexa-scan analysis 0
DSS Histology 0
Dysmorphology 0
Ear epidermis immunophenotyping 0
ECG (Electrocardiogram) (GMC) 0
Echo 0
Electrocardiogram (ECG) 0
Electroconvulsive Threshold Testing 0
Electroretinography 0
ELISA (GMC) 0
ERG (Electroretinogram) (GMC) 0
Eye Morphology 0
Eye size (GMC) 0
FACS (GMC) 0
FACs Analysis 0
Fasted Clinical Chemistry 0
Fear Conditioning 0
Femoral Microradiography 0
Fertility of Homozygous Knock-out Mice 7512
Food efficiency (GMC) 0
Grip-Strength 0
Grip Strength 0
Grip Strength (GMC) 0
Gross Morphology Embryo E12.5 0
Gross Morphology Embryo E14.5-E15.5 0
Gross Morphology Embryo E18.5 0
Gross Morphology Embryo E9.5 0
Gross Morphology Placenta E12.5 0
Gross Morphology Placenta E14.5-E15.5 0
Gross Morphology Placenta E18.5 0
Gross Morphology Placenta E9.5 0
Gross Pathology and Tissue Collection 0
Haematology 0
Haematology (GMC) 0
Haematology test 0
Heart Dissection 0
Heart Weight 0
Heart weight/tibia length 0
Hematology 0
Hole-board Exploration 0
Holeboard (GMC) 0
Hot Plate 0
Immunoglobulin 0
Indirect Calorimetry 0
Indirect ophthalmoscopy 0
Insulin Blood Level 0
Intraperitoneal glucose tolerance test (IPGTT) 0
IPGTT 0
Light-Dark Test 0
Mesenteric Lymph Node Immunophenotyping 0
Modified SHIRPA 0
Nociception Hotplate (GMC) 0
Open-field 0
Open Field 0
Ophthalmoscope 0
Organ Weight 0
pDexa (GMC) 0
Plasma Chemistry 0
Rotarod 0
Rotarod A (GMC) 0
Shirpa (GMC) 0
Simplified IPGTT 0
Sleep Wake 0
Slit Lamp 0
Spleen Immunophenotyping 0
Spontaneous breathing (GMC) 0
Tail Suspension 0
Three-point Bend 0
Trabecular Bone MicroCT 0
Trichuris 0
Urinalysis 0
Vertebra Compression 0
Vertebral Microradiography 0
Viability E12.5 Secondary Screen 0
Viability E14.5-E15.5 Secondary Screen 0
Viability E18.5 Secondary Screen 0
Viability E9.5 Secondary Screen 0
Viability Primary Screen 7362
Whole blood peripheral blood leukocyte immunophenotyping 0
X-ray 0
X-Ray 0
X-Ray (GMC) 0
unidimensional
Acoustic Startle and Pre-pulse Inhibition (PPI) 20676
Acoustic Startle&PPI 8356
Allergy (GMC) 27
Anti-nuclear antibody assay 982
Antigen Specific Immunoglobulin Assay 482
Auditory Brain Stem Response 2149
Body Composition (DEXA lean/fat) 42582
Body Weight 7263
Bodyweight (GMC) 27
Bone marrow immunophenotyping 8043
Buffy coat peripheral blood leukocyte immunophenotyping 7560
Calorimetry 731
Challenge Whole Body Plethysmography 46
Clinical Chemistry 99693
Clinical chemistry (GMC) 548
Combined SHIRPA and Dysmorphology 3645
Cortical Bone MicroCT 18
DEXA 5529
Dexa-scan analysis 4516
DSS Histology 133
Dysmorphology 0
Ear epidermis immunophenotyping 3112
ECG (Electrocardiogram) (GMC) 130
Echo 16196
Electrocardiogram (ECG) 32701
Electroconvulsive Threshold Testing 409
Electroretinography 5
ELISA (GMC) 96
ERG (Electroretinogram) (GMC) 148
Eye Morphology 3263
Eye size (GMC) 34
FACS (GMC) 185
FACs Analysis 3644
Fasted Clinical Chemistry 3744
Fear Conditioning 420
Femoral Microradiography 12
Fertility of Homozygous Knock-out Mice 0
Food efficiency (GMC) 135
Grip-Strength 2472
Grip Strength 19766
Grip Strength (GMC) 19
Gross Morphology Embryo E12.5 0
Gross Morphology Embryo E14.5-E15.5 0
Gross Morphology Embryo E18.5 0
Gross Morphology Embryo E9.5 0
Gross Morphology Placenta E12.5 0
Gross Morphology Placenta E14.5-E15.5 0
Gross Morphology Placenta E18.5 0
Gross Morphology Placenta E9.5 0
Gross Pathology and Tissue Collection 0
Haematology 4731
Haematology (GMC) 428
Haematology test 6038
Heart Dissection 1510
Heart Weight 4312
Heart weight/tibia length 776
Hematology 49606
Hole-board Exploration 991
Holeboard (GMC) 513
Hot Plate 2044
Immunoglobulin 929
Indirect Calorimetry 4919
Indirect ophthalmoscopy 0
Insulin Blood Level 2081
Intraperitoneal glucose tolerance test (IPGTT) 12990
IPGTT 408
Light-Dark Test 11077
Mesenteric Lymph Node Immunophenotyping 44829
Modified SHIRPA 1245
Nociception Hotplate (GMC) 91
Open-field 7458
Open Field 50546
Ophthalmoscope 757
Organ Weight 1671
pDexa (GMC) 321
Plasma Chemistry 894
Rotarod 1798
Rotarod A (GMC) 6
Shirpa (GMC) 197
Simplified IPGTT 622
Sleep Wake 5333
Slit Lamp 0
Spleen Immunophenotyping 56372
Spontaneous breathing (GMC) 576
Tail Suspension 534
Three-point Bend 30
Trabecular Bone MicroCT 28
Trichuris 0
Urinalysis 1034
Vertebra Compression 21
Vertebral Microradiography 14
Viability E12.5 Secondary Screen 0
Viability E14.5-E15.5 Secondary Screen 0
Viability E18.5 Secondary Screen 0
Viability E9.5 Secondary Screen 0
Viability Primary Screen 0
Whole blood peripheral blood leukocyte immunophenotyping 0
X-ray 2543
X-Ray 231
X-Ray (GMC) 0
unidimensional-ReferenceRange
Acoustic Startle and Pre-pulse Inhibition (PPI) 0
Acoustic Startle&PPI 0
Allergy (GMC) 0
Anti-nuclear antibody assay 0
Antigen Specific Immunoglobulin Assay 0
Auditory Brain Stem Response 26586
Body Composition (DEXA lean/fat) 0
Body Weight 0
Bodyweight (GMC) 0
Bone marrow immunophenotyping 0
Buffy coat peripheral blood leukocyte immunophenotyping 0
Calorimetry 0
Challenge Whole Body Plethysmography 0
Clinical Chemistry 0
Clinical chemistry (GMC) 0
Combined SHIRPA and Dysmorphology 0
Cortical Bone MicroCT 0
DEXA 0
Dexa-scan analysis 0
DSS Histology 0
Dysmorphology 0
Ear epidermis immunophenotyping 0
ECG (Electrocardiogram) (GMC) 0
Echo 0
Electrocardiogram (ECG) 0
Electroconvulsive Threshold Testing 0
Electroretinography 0
ELISA (GMC) 0
ERG (Electroretinogram) (GMC) 0
Eye Morphology 0
Eye size (GMC) 0
FACS (GMC) 0
FACs Analysis 0
Fasted Clinical Chemistry 0
Fear Conditioning 0
Femoral Microradiography 0
Fertility of Homozygous Knock-out Mice 0
Food efficiency (GMC) 0
Grip-Strength 0
Grip Strength 0
Grip Strength (GMC) 0
Gross Morphology Embryo E12.5 0
Gross Morphology Embryo E14.5-E15.5 0
Gross Morphology Embryo E18.5 0
Gross Morphology Embryo E9.5 0
Gross Morphology Placenta E12.5 0
Gross Morphology Placenta E14.5-E15.5 0
Gross Morphology Placenta E18.5 0
Gross Morphology Placenta E9.5 0
Gross Pathology and Tissue Collection 0
Haematology 0
Haematology (GMC) 0
Haematology test 0
Heart Dissection 0
Heart Weight 0
Heart weight/tibia length 0
Hematology 0
Hole-board Exploration 0
Holeboard (GMC) 0
Hot Plate 0
Immunoglobulin 0
Indirect Calorimetry 0
Indirect ophthalmoscopy 0
Insulin Blood Level 0
Intraperitoneal glucose tolerance test (IPGTT) 0
IPGTT 0
Light-Dark Test 0
Mesenteric Lymph Node Immunophenotyping 0
Modified SHIRPA 0
Nociception Hotplate (GMC) 0
Open-field 0
Open Field 0
Ophthalmoscope 0
Organ Weight 0
pDexa (GMC) 0
Plasma Chemistry 0
Rotarod 0
Rotarod A (GMC) 0
Shirpa (GMC) 0
Simplified IPGTT 0
Sleep Wake 0
Slit Lamp 0
Spleen Immunophenotyping 0
Spontaneous breathing (GMC) 0
Tail Suspension 0
Three-point Bend 0
Trabecular Bone MicroCT 0
Trichuris 0
Urinalysis 0
Vertebra Compression 0
Vertebral Microradiography 0
Viability E12.5 Secondary Screen 0
Viability E14.5-E15.5 Secondary Screen 0
Viability E18.5 Secondary Screen 0
Viability E9.5 Secondary Screen 0
Viability Primary Screen 0
Whole blood peripheral blood leukocyte immunophenotyping 29520
X-ray 10450
X-Ray 414
X-Ray (GMC) 0
#dat <- KOMPv10.1 %>% select(procedure_name=="Gross Pathology and Tissue Collection")
# extract unidimensional data only.
dim(KOMPv10.1)
[1] 1779903 88
KOMPv10.1.ud <- KOMPv10.1 %>% filter(data_type=="unidimensional")
dim(KOMPv10.1.ud)
[1] 580001 88
table(allpheno$procedure_name)
Body Composition (DEXA lean/fat)
234305
#"Auditory Brain Stem Response"
#"Clinical Chemistry"
#"Body Composition (DEXA lean/fat)"
#"Intraperitoneal glucose tolerance test (IPGTT)"
#"Hematology"
# count the number of tests in each phenotype
proc.list <- table(KOMPv10.1.ud$procedure_name)
#proc.list <- proc.list[proc.list>1000]
proc.list
Acoustic Startle and Pre-pulse Inhibition (PPI)
20676
Acoustic Startle&PPI
8356
Allergy (GMC)
27
Anti-nuclear antibody assay
982
Antigen Specific Immunoglobulin Assay
482
Auditory Brain Stem Response
2149
Body Composition (DEXA lean/fat)
42582
Body Weight
7263
Bodyweight (GMC)
27
Bone marrow immunophenotyping
8043
Buffy coat peripheral blood leukocyte immunophenotyping
7560
Calorimetry
731
Challenge Whole Body Plethysmography
46
Clinical Chemistry
99693
Clinical chemistry (GMC)
548
Combined SHIRPA and Dysmorphology
3645
Cortical Bone MicroCT
18
DEXA
5529
Dexa-scan analysis
4516
DSS Histology
133
Ear epidermis immunophenotyping
3112
ECG (Electrocardiogram) (GMC)
130
Echo
16196
Electrocardiogram (ECG)
32701
Electroconvulsive Threshold Testing
409
Electroretinography
5
ELISA (GMC)
96
ERG (Electroretinogram) (GMC)
148
Eye Morphology
3263
Eye size (GMC)
34
FACS (GMC)
185
FACs Analysis
3644
Fasted Clinical Chemistry
3744
Fear Conditioning
420
Femoral Microradiography
12
Food efficiency (GMC)
135
Grip-Strength
2472
Grip Strength
19766
Grip Strength (GMC)
19
Haematology
4731
Haematology (GMC)
428
Haematology test
6038
Heart Dissection
1510
Heart Weight
4312
Heart weight/tibia length
776
Hematology
49606
Hole-board Exploration
991
Holeboard (GMC)
513
Hot Plate
2044
Immunoglobulin
929
Indirect Calorimetry
4919
Insulin Blood Level
2081
Intraperitoneal glucose tolerance test (IPGTT)
12990
IPGTT
408
Light-Dark Test
11077
Mesenteric Lymph Node Immunophenotyping
44829
Modified SHIRPA
1245
Nociception Hotplate (GMC)
91
Open-field
7458
Open Field
50546
Ophthalmoscope
757
Organ Weight
1671
pDexa (GMC)
321
Plasma Chemistry
894
Rotarod
1798
Rotarod A (GMC)
6
Shirpa (GMC)
197
Simplified IPGTT
622
Sleep Wake
5333
Spleen Immunophenotyping
56372
Spontaneous breathing (GMC)
576
Tail Suspension
534
Three-point Bend
30
Trabecular Bone MicroCT
28
Urinalysis
1034
Vertebra Compression
21
Vertebral Microradiography
14
X-ray
2543
X-Ray
231
length(proc.list)
[1] 79
pheno.list <- table(KOMPv10.1.ud$proc_param_name)
pheno.list <- pheno.list[pheno.list>1000] # find list of phenotypes with more than 1000 tests (i.e. 1000 mutants tested)
pheno.list <- names(pheno.list)
pheno.list
[1] "Acoustic Startle and Pre-pulse Inhibition (PPI)_% Pre-pulse inhibition - Global"
[2] "Acoustic Startle and Pre-pulse Inhibition (PPI)_% Pre-pulse inhibition - PPI1"
[3] "Acoustic Startle and Pre-pulse Inhibition (PPI)_% Pre-pulse inhibition - PPI2"
[4] "Acoustic Startle and Pre-pulse Inhibition (PPI)_% Pre-pulse inhibition - PPI3"
[5] "Acoustic Startle and Pre-pulse Inhibition (PPI)_% Pre-pulse inhibition - PPI4"
[6] "Acoustic Startle and Pre-pulse Inhibition (PPI)_Response amplitude - S"
[7] "Body Composition (DEXA lean/fat)_BMC/Body weight"
[8] "Body Composition (DEXA lean/fat)_Body length"
[9] "Body Composition (DEXA lean/fat)_Bone Area"
[10] "Body Composition (DEXA lean/fat)_Bone Mineral Content (excluding skull)"
[11] "Body Composition (DEXA lean/fat)_Bone Mineral Density (excluding skull)"
[12] "Body Composition (DEXA lean/fat)_Fat mass"
[13] "Body Composition (DEXA lean/fat)_Fat/Body weight"
[14] "Body Composition (DEXA lean/fat)_Lean mass"
[15] "Body Composition (DEXA lean/fat)_Lean/Body weight"
[16] "Body Weight_Body Weight"
[17] "Clinical Chemistry_Alanine aminotransferase"
[18] "Clinical Chemistry_Albumin"
[19] "Clinical Chemistry_Alkaline phosphatase"
[20] "Clinical Chemistry_Alpha-amylase"
[21] "Clinical Chemistry_Aspartate aminotransferase"
[22] "Clinical Chemistry_Calcium"
[23] "Clinical Chemistry_Chloride"
[24] "Clinical Chemistry_Creatine kinase"
[25] "Clinical Chemistry_Creatinine"
[26] "Clinical Chemistry_Free fatty acids"
[27] "Clinical Chemistry_Fructosamine"
[28] "Clinical Chemistry_Glucose"
[29] "Clinical Chemistry_Glycerol"
[30] "Clinical Chemistry_HDL-cholesterol"
[31] "Clinical Chemistry_Iron"
[32] "Clinical Chemistry_LDL-cholesterol"
[33] "Clinical Chemistry_Magnesium"
[34] "Clinical Chemistry_Phosphorus"
[35] "Clinical Chemistry_Potassium"
[36] "Clinical Chemistry_Sodium"
[37] "Clinical Chemistry_Total bilirubin"
[38] "Clinical Chemistry_Total cholesterol"
[39] "Clinical Chemistry_Total protein"
[40] "Clinical Chemistry_Triglyceride"
[41] "Clinical Chemistry_Triglycerides"
[42] "Clinical Chemistry_Urea"
[43] "Clinical Chemistry_Urea (Blood Urea Nitrogen - BUN)"
[44] "Combined SHIRPA and Dysmorphology_Locomotor activity"
[45] "Echo_Cardiac Output"
[46] "Echo_Ejection Fraction"
[47] "Echo_Fractional Shortening"
[48] "Echo_HR"
[49] "Echo_LVAWd"
[50] "Echo_LVIDd"
[51] "Echo_LVIDs"
[52] "Echo_LVPWd"
[53] "Echo_LVPWs"
[54] "Echo_Stroke Volume"
[55] "Electrocardiogram (ECG)_CV"
[56] "Electrocardiogram (ECG)_HR"
[57] "Electrocardiogram (ECG)_HRV"
[58] "Electrocardiogram (ECG)_PQ"
[59] "Electrocardiogram (ECG)_PR"
[60] "Electrocardiogram (ECG)_QRS"
[61] "Electrocardiogram (ECG)_QTc"
[62] "Electrocardiogram (ECG)_QTc Dispersion"
[63] "Electrocardiogram (ECG)_rMSSD"
[64] "Electrocardiogram (ECG)_RR"
[65] "Electrocardiogram (ECG)_ST"
[66] "Grip Strength_Forelimb and hindlimb grip strength measurement mean"
[67] "Grip Strength_Forelimb and hindlimb grip strength normalised against body weight"
[68] "Grip Strength_Forelimb grip strength measurement mean"
[69] "Grip Strength_Forelimb grip strength normalised against body weight"
[70] "Heart Weight_Heart weight"
[71] "Hematology_Basophil cell count"
[72] "Hematology_Basophil differential count"
[73] "Hematology_Eosinophil cell count"
[74] "Hematology_Eosinophil differential count"
[75] "Hematology_Hematocrit"
[76] "Hematology_Hemoglobin"
[77] "Hematology_Lymphocyte cell count"
[78] "Hematology_Lymphocyte differential count"
[79] "Hematology_Mean cell hemoglobin concentration"
[80] "Hematology_Mean cell volume"
[81] "Hematology_Mean corpuscular hemoglobin"
[82] "Hematology_Mean platelet volume"
[83] "Hematology_Monocyte cell count"
[84] "Hematology_Monocyte differential count"
[85] "Hematology_Neutrophil cell count"
[86] "Hematology_Neutrophil differential count"
[87] "Hematology_Platelet count"
[88] "Hematology_Red blood cell count"
[89] "Hematology_Red blood cell distribution width"
[90] "Hematology_White blood cell count"
[91] "Hot Plate_Time of first response"
[92] "Indirect Calorimetry_Respiratory Exchange Ratio"
[93] "Insulin Blood Level_Insulin"
[94] "Intraperitoneal glucose tolerance test (IPGTT)_Area under glucose response curve"
[95] "Intraperitoneal glucose tolerance test (IPGTT)_Fasted blood glucose concentration"
[96] "Intraperitoneal glucose tolerance test (IPGTT)_Initial response to glucose challenge"
[97] "Light-Dark Test_Dark side time spent"
[98] "Light-Dark Test_Fecal boli"
[99] "Light-Dark Test_Latency to first transition into dark"
[100] "Light-Dark Test_Light side time spent"
[101] "Light-Dark Test_Percent time in dark"
[102] "Light-Dark Test_Percent time in light"
[103] "Light-Dark Test_Side changes"
[104] "Light-Dark Test_Time mobile dark side"
[105] "Light-Dark Test_Time mobile light side"
[106] "Modified SHIRPA_Locomotor activity"
[107] "Open Field_Center average speed"
[108] "Open Field_Center distance travelled"
[109] "Open Field_Center permanence time"
[110] "Open Field_Center resting time"
[111] "Open Field_Distance travelled - total"
[112] "Open Field_Latency to center entry"
[113] "Open Field_Number of center entries"
[114] "Open Field_Number of rears - total"
[115] "Open Field_Percentage center time"
[116] "Open Field_Periphery average speed"
[117] "Open Field_Periphery distance travelled"
[118] "Open Field_Periphery permanence time"
[119] "Open Field_Periphery resting time"
[120] "Open Field_Whole arena average speed"
[121] "Open Field_Whole arena resting time"
[122] "X-ray_Tibia length"
length(pheno.list) #122
[1] 122
# Use phenotypes with more than 1000 tests (i.e. 1000 mutants tested)
dim(KOMPv10.1.ud)
[1] 580001 88
ap.stat <- KOMPv10.1.ud %>% filter(proc_param_name %in% pheno.list)
dim(ap.stat)
[1] 358088 88
mtest <- table(ap.stat$proc_param_name, ap.stat$marker_symbol)
mtest <-as.data.frame.matrix(mtest)
dim(mtest)
[1] 122 5954
if(FALSE){
nmax <-max(mtest)
library(circlize)
col_fun = colorRamp2(c(0, nmax), c("white", "red"))
col_fun(seq(0, nmax))
pdf("~/Google Drive Miami/Miami_IMPC/output/KMOPv10.1_heatmap_gene_vs_pheno_after_filtering.pdf", width = 10, height = 10)
ht = Heatmap(as.matrix(mtest), cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, col = col_fun,
row_names_gp = gpar(fontsize = 5), name="Count")
draw(ht)
dev.off()
}
table(ap.stat$procedure_name)
Acoustic Startle and Pre-pulse Inhibition (PPI)
20676
Body Composition (DEXA lean/fat)
42582
Body Weight
1242
Clinical Chemistry
95963
Combined SHIRPA and Dysmorphology
3645
Echo
12435
Electrocardiogram (ECG)
32701
Grip Strength
18995
Heart Weight
3657
Hematology
48130
Hot Plate
1329
Indirect Calorimetry
2215
Insulin Blood Level
2081
Intraperitoneal glucose tolerance test (IPGTT)
12990
Light-Dark Test
11077
Modified SHIRPA
1245
Open Field
45326
X-ray
1799
table(allpheno$procedure_name)
Body Composition (DEXA lean/fat)
234305
ap.stat = ap.stat %>%
dplyr::select(phenotyping_center, procedure_name, parameter_name, zygosity, allele_symbol,
genotype_effect_parameter_estimate, genotype_effect_stderr_estimate,
genotype_effect_p_value, phenotyping_center, allele_name, marker_symbol) %>%
filter(procedure_name == "Body Composition (DEXA lean/fat)") %>%
mutate(procedure_name=recode(procedure_name, "Body Composition (DEXA lean/fat)"="BC")) %>%
mutate(z_score = genotype_effect_parameter_estimate/genotype_effect_stderr_estimate,
proc_param_name=paste0(procedure_name,"_",parameter_name),
gene_pheno = paste0(parameter_name, "_", allele_symbol))
table(ap.stat$parameter_name, ap.stat$procedure_name)
BC
BMC/Body weight 4780
Body length 4197
Bone Area 4780
Bone Mineral Content (excluding skull) 4780
Bone Mineral Density (excluding skull) 4780
Fat mass 4816
Fat/Body weight 4815
Lean mass 4817
Lean/Body weight 4817
length(unique(ap.stat$marker_symbol)) #4428
[1] 4428
length(unique(ap.stat$allele_symbol)) #4559
[1] 4559
length(unique(ap.stat$proc_param_name)) #9 # number of phenotypes in association statistics data set
[1] 9
length(unique(allpheno$proc_param_name)) #10 # number of phenotypes in final control data
[1] 10
pheno.list.stat <- unique(ap.stat$proc_param_name)
pheno.list.ctrl <- unique(allpheno$proc_param_name)
sum(pheno.list.stat %in% pheno.list.ctrl)
[1] 9
sum(pheno.list.ctrl %in% pheno.list.stat)
[1] 9
## extract common phenotype list
common.pheno.list <- sort(intersect(pheno.list.ctrl, pheno.list.stat))
common.pheno.list
[1] "BC_BMC/Body weight"
[2] "BC_Body length"
[3] "BC_Bone Area"
[4] "BC_Bone Mineral Content (excluding skull)"
[5] "BC_Bone Mineral Density (excluding skull)"
[6] "BC_Fat mass"
[7] "BC_Fat/Body weight"
[8] "BC_Lean mass"
[9] "BC_Lean/Body weight"
length(common.pheno.list)
[1] 9
# Use summary statistics of common phenotypes
dim(ap.stat)
[1] 42582 13
ap.stat <- ap.stat %>% filter(proc_param_name %in% common.pheno.list)
dim(ap.stat)
[1] 42582 13
length(unique(ap.stat$proc_param_name))
[1] 9
mtest <- table(ap.stat$proc_param_name, ap.stat$marker_symbol)
mtest <-as.data.frame.matrix(mtest)
nmax <-max(mtest)
col_fun = colorRamp2(c(0, nmax), c("white", "red"))
col_fun(seq(0, nmax))
[1] "#FFFFFFFF" "#FFF0EBFF" "#FFE2D7FF" "#FFD3C4FF" "#FFC4B0FF" "#FFB59DFF"
[7] "#FFA68BFF" "#FF9678FF" "#FF8666FF" "#FF7554FF" "#FF6342FF" "#FF4E2FFF"
[13] "#FF351BFF" "#FF0000FF"
ht = Heatmap(as.matrix(mtest), cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, col = col_fun,
row_names_gp = gpar(fontsize = 8), name="Count")
draw(ht)
## sum(z-score)/sqrt(# of zscore)
sumz <- function(z){ sum(z)/sqrt(length(z)) }
ap.z = ap.stat %>%
dplyr::select(marker_symbol, proc_param_name, z_score) %>%
na.omit() %>%
group_by(marker_symbol, proc_param_name) %>%
summarize(zscore = sumz(z_score)) ## combine z-scores
`summarise()` has grouped output by 'marker_symbol'. You can override using the `.groups` argument.
dim(ap.z)
[1] 32174 3
nan2na <- function(df){
out <- data.frame(sapply(df, function(x) ifelse(is.nan(x), NA, x)))
colnames(out) <- colnames(df)
out
}
ap.zmat = dcast(ap.z, marker_symbol ~ proc_param_name, value.var = "zscore",
fun.aggregate = mean) %>% tibble::column_to_rownames(var="marker_symbol")
ap.zmat = nan2na(ap.zmat) #convert nan to na
dim(ap.zmat)
[1] 4412 9
id.mat <- 1*(!is.na(ap.zmat)) # multiply 1 to make this matrix numeric
nrow(as.data.frame(colSums(id.mat)))
[1] 9
dim(id.mat)
[1] 4412 9
# heatmap of gene - phenotype (red: tested, white: untested)
if(FALSE){
pdf("~/Google Drive Miami/Miami_IMPC/output/missing_tests_after_filtering_BC.pdf", width = 6, height = 2.7)
ht = Heatmap(t(id.mat),
cluster_rows = T, clustering_distance_rows ="binary",
cluster_columns = T, clustering_distance_columns = "binary",
show_row_dend = F, show_column_dend = F, # do not show dendrogram
show_column_names = F, col = c("white","red"),
row_names_gp = gpar(fontsize = 10), name="Missing")
draw(ht)
dev.off()
}
We plot a association Z-score distribution for each phenotype.
ggplot(melt(ap.zmat), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
No id variables; using all as measure variables
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 7534 rows containing non-finite values (stat_bin).
Here, we estimate the genetic correlations between phenotypes using association Z-score matrix (num of genes:5479, num of phenotypes 14).
ap.zmat <- ap.zmat[,common.pheno.list]
ap.zcor = cor(ap.zmat, use="pairwise.complete.obs", method="spearman")
#col <- colorRampPalette(c("steelblue", "white", "darkorange"))(100)
#pheatmap(op.zcor, cluster_rows = T, cluster_cols=T, show_colnames=F, col=col)
#op.cor.order <- op.cor.out$tree_row[["order"]]
#op.zcor.org <- op.zcor # this will be used in correlation matrix test
#op.zcor <- op.zcor[op.cor.order,]
#op.zcor <- op.zcor[,op.cor.order]
#pheatmap(ap.zcor, cluster_rows = F, cluster_cols=F, show_colnames=F, col=col)
ht = Heatmap(ap.zcor, cluster_rows = T, cluster_columns = T, show_column_names = F, #col = col_fun,
row_names_gp = gpar(fontsize = 10),
#name="Genetic corr (Z-score)"
name="Genetic Corr (Zscore)"
)
draw(ht)
#pheno.order <- row_order(ht)
#ap.zcor <- ap.zcor[pheno.order,pheno.order]
We compare a correlation matrix obtained using control mice phenotype data v.s. a genetic correlation matrix estimated using association Z-scores. As you can see, both correlation heatmaps have similar correlation pattern.
ap.cor.rank.fig <- ap.cor.rank[common.pheno.list,common.pheno.list]
ap.cor.fig <- ap.cor[common.pheno.list,common.pheno.list]
ap.cor.combat.fig <- ap.cor.combat[common.pheno.list, common.pheno.list]
ap.zcor.fig <- ap.zcor
ht = Heatmap(ap.cor.rank.fig, cluster_rows = TRUE, cluster_columns = TRUE, show_column_names = F, #col = col_fun,
show_row_dend = F, show_column_dend = F, # do not show dendrogram
row_names_gp = gpar(fontsize = 8), column_title="Phenotype Corr (RankZ, Pearson)", column_title_gp = gpar(fontsize = 8),
name="Corr")
pheno.order <- row_order(ht)
Warning: The heatmap has not been initialized. You might have different results
if you repeatedly execute this function, e.g. when row_km/column_km was
set. It is more suggested to do as `ht = draw(ht); row_order(ht)`.
draw(ht)
if(FALSE){
pdf("~/Google Drive Miami/Miami_IMPC/output/comp_pheno_corr_gene_corr_combat_BC.pdf", width = 9, height = 2)
ap.cor.rank.fig <- ap.cor.rank.fig[pheno.order,pheno.order]
ht1 = Heatmap(ap.cor.rank.fig, cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, #col = col_fun,
show_row_dend = F, show_column_dend = F, # do not show dendrogram
row_names_gp = gpar(fontsize = 8), column_title="Phenotype Corr (RankZ, Pearson)", column_title_gp = gpar(fontsize = 8),
name="Corr")
ap.cor.fig <- ap.cor.fig[pheno.order,pheno.order]
ht2 = Heatmap(ap.cor.fig, cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, #col = col_fun,
row_names_gp = gpar(fontsize = 8), column_title="Phenotype Corr (Spearman)", column_title_gp = gpar(fontsize = 8),
name="Corr")
ap.cor.combat.fig <- ap.cor.combat.fig[pheno.order,pheno.order]
ht3 = Heatmap(ap.cor.combat.fig, cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, #col = col_fun,
row_names_gp = gpar(fontsize = 8), column_title="Phenotype Corr (Combat, Pearson)", column_title_gp = gpar(fontsize = 8),
name="Corr")
ap.zcor.fig <- ap.zcor.fig[pheno.order,pheno.order]
ht4 = Heatmap(ap.zcor.fig, cluster_rows = FALSE, cluster_columns = FALSE, show_column_names = F, #col = col_fun,
row_names_gp = gpar(fontsize = 8), column_title="Genetic Corr (Pearson)", column_title_gp = gpar(fontsize = 8),
name="Corr"
)
draw(ht1+ht2+ht3+ht4)
dev.off()
}
It looks like Jenrich (1970) test is too conservative here. Instead, we use Mantel test testing the correlation between two distance matrices.
####################
# Use Mantel test
# https://stats.idre.ucla.edu/r/faq/how-can-i-perform-a-mantel-test-in-r/
# install.packages("ade4")
library(ade4)
to.upper<-function(X) X[upper.tri(X,diag=FALSE)]
a1 <- to.upper(ap.cor.fig)
a2 <- to.upper(ap.cor.rank.fig)
a3 <- to.upper(ap.cor.combat.fig)
a4 <- to.upper(ap.zcor.fig)
plot(a4, a1)
plot(a4, a2)
plot(a4, a3)
mantel.rtest(as.dist(1-ap.cor.fig), as.dist(1-ap.zcor.fig), nrepet = 9999) #nrepet = number of permutations
Monte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.7725431
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
4.421228997 -0.001000792 0.030611450
mantel.rtest(as.dist(1-ap.cor.rank.fig), as.dist(1-ap.zcor.fig), nrepet = 9999)
Monte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.7667761
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
4.2448574362 0.0001120909 0.0326199968
mantel.rtest(as.dist(1-ap.cor.combat.fig), as.dist(1-ap.zcor.fig), nrepet = 9999)
Monte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.9087556
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
4.8503970246 0.0005979198 0.0350564523
Impute z-scores of untested gene-pheno pair using phenotype correlation matrix
if(FALSE){
library(devtools)
devtools::install_github("dleelab/kompute")
}
library(kompute)
We randomly select measured gene-phenotype association z-scores, mask those, impute them using kompute algorithm. Then we compare the imputed z-scores to the measured ones.
zmat <-t(ap.zmat)
dim(zmat)
[1] 9 4412
#filter genes with na < 20
zmat0 <- is.na(zmat)
num.na<-colSums(zmat0)
summary(num.na)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 1.000 1.708 3.000 8.000
zmat <- zmat[,num.na<10]
dim(zmat)
[1] 9 4412
#pheno.cor <- ap.cor.fig
#pheno.cor <- ap.cor.rank.fig
pheno.cor <- ap.cor.combat.fig
#pheno.cor <- ap.zcor.fig
zmat <- zmat[rownames(pheno.cor),,drop=FALSE]
rownames(zmat)
[1] "BC_BMC/Body weight"
[2] "BC_Body length"
[3] "BC_Bone Area"
[4] "BC_Bone Mineral Content (excluding skull)"
[5] "BC_Bone Mineral Density (excluding skull)"
[6] "BC_Fat mass"
[7] "BC_Fat/Body weight"
[8] "BC_Lean mass"
[9] "BC_Lean/Body weight"
rownames(pheno.cor)
[1] "BC_BMC/Body weight"
[2] "BC_Body length"
[3] "BC_Bone Area"
[4] "BC_Bone Mineral Content (excluding skull)"
[5] "BC_Bone Mineral Density (excluding skull)"
[6] "BC_Fat mass"
[7] "BC_Fat/Body weight"
[8] "BC_Lean mass"
[9] "BC_Lean/Body weight"
colnames(pheno.cor)
[1] "BC_BMC/Body weight"
[2] "BC_Body length"
[3] "BC_Bone Area"
[4] "BC_Bone Mineral Content (excluding skull)"
[5] "BC_Bone Mineral Density (excluding skull)"
[6] "BC_Fat mass"
[7] "BC_Fat/Body weight"
[8] "BC_Lean mass"
[9] "BC_Lean/Body weight"
npheno <- nrow(zmat)
# percentage of missing Z-scores in the original data
100*sum(is.na(zmat))/(nrow(zmat)*ncol(zmat)) # 43%
[1] 18.97351
nimp <- 2000 # # of missing/imputed Z-scores
set.seed(1111)
# find index of all measured zscores
all.i <- 1:(nrow(zmat)*ncol(zmat))
measured <- as.vector(!is.na(as.matrix(zmat)))
measured.i <- all.i[measured]
# mask 2000 measured z-scores
mask.i <- sort(sample(measured.i, nimp))
org.z = as.matrix(zmat)[mask.i]
zvec <- as.vector(as.matrix(zmat))
zvec[mask.i] <- NA
zmat.imp <- matrix(zvec, nrow=npheno)
rownames(zmat.imp) <- rownames(zmat)
kompute.res <- kompute(zmat.imp, pheno.cor, 0.01)
KOMPute running...
# of genes: 4412
# of phenotypes: 9
# of imputed Z-scores: 9534
# measured vs imputed
length(org.z)
[1] 2000
imp.z <- as.matrix(kompute.res$zmat)[mask.i]
imp.info <- as.matrix(kompute.res$infomat)[mask.i]
plot(imp.z, org.z)
imp <- data.frame(org.z=org.z, imp.z=imp.z, info=imp.info)
dim(imp)
[1] 2000 3
imp <- imp[complete.cases(imp),]
imp <- subset(imp, info>=0 & info <= 1)
dim(imp)
[1] 1995 3
cor.val <- round(cor(imp$imp.z, imp$org.z), digits=3)
cor.val
[1] 0.856
plot(imp$imp.z, imp$org.z)
info.cutoff <- 0.9
imp.sub <- subset(imp, info>info.cutoff)
dim(imp.sub)
[1] 582 3
summary(imp.sub$imp.z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.6750 -1.3848 -0.3271 -0.4539 0.6020 6.1920
summary(imp.sub$info)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.9001 0.9327 0.9588 0.9494 0.9691 0.9715
cor.val <- round(cor(imp.sub$imp.z, imp.sub$org.z), digits=3)
cor.val
[1] 0.946
g <- ggplot(imp.sub, aes(x=imp.z, y=org.z, col=info)) +
geom_point() +
labs(title=paste0("IMPC Behavior Data, Info>", info.cutoff, ", Cor=",cor.val),
x="Imputed Z-scores", y = "Measured Z-scores", col="Info") +
theme_minimal()
g
#filename <- "~/Google Drive Miami/Miami_IMPC/output/realdata_measured_vs_imputed_info_BC.pdf"
#ggsave(filename, plot=g, height=4, width=5)
sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] kompute_0.1.0 ade4_1.7-18 sva_3.42.0
[4] BiocParallel_1.28.0 genefilter_1.76.0 mgcv_1.8-36
[7] nlme_3.1-152 lme4_1.1-27.1 RNOmni_1.0.0
[10] ComplexHeatmap_2.10.0 circlize_0.4.13 irlba_2.3.3
[13] Matrix_1.3-4 RColorBrewer_1.1-2 tidyr_1.1.4
[16] ggplot2_3.3.5 reshape2_1.4.4 dplyr_1.0.7
[19] data.table_1.14.2 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] bitops_1.0-7 matrixStats_0.61.0 fs_1.5.0
[4] bit64_4.0.5 doParallel_1.0.16 httr_1.4.2
[7] GenomeInfoDb_1.30.0 rprojroot_2.0.2 tools_4.1.1
[10] utf8_1.2.2 R6_2.5.1 DBI_1.1.1
[13] BiocGenerics_0.40.0 colorspace_2.0-2 GetoptLong_1.0.5
[16] withr_2.4.2 tidyselect_1.1.1 bit_4.0.4
[19] compiler_4.1.1 git2r_0.28.0 Biobase_2.54.0
[22] labeling_0.4.2 scales_1.1.1 stringr_1.4.0
[25] digest_0.6.28 minqa_1.2.4 R.utils_2.11.0
[28] rmarkdown_2.11 XVector_0.34.0 pkgconfig_2.0.3
[31] htmltools_0.5.2 limma_3.50.0 fastmap_1.1.0
[34] highr_0.9 rlang_0.4.12 GlobalOptions_0.1.2
[37] RSQLite_2.2.8 shape_1.4.6 jquerylib_0.1.4
[40] farver_2.1.0 generics_0.1.1 R.oo_1.24.0
[43] RCurl_1.98-1.5 magrittr_2.0.1 GenomeInfoDbData_1.2.7
[46] Rcpp_1.0.7 munsell_0.5.0 S4Vectors_0.32.0
[49] fansi_0.5.0 R.methodsS3_1.8.1 lifecycle_1.0.1
[52] edgeR_3.36.0 stringi_1.7.5 whisker_0.4
[55] yaml_2.2.1 zlibbioc_1.40.0 MASS_7.3-54
[58] plyr_1.8.6 blob_1.2.2 parallel_4.1.1
[61] promises_1.2.0.1 crayon_1.4.1 lattice_0.20-44
[64] Biostrings_2.62.0 splines_4.1.1 annotate_1.72.0
[67] KEGGREST_1.34.0 locfit_1.5-9.4 knitr_1.36
[70] pillar_1.6.4 boot_1.3-28 rjson_0.2.20
[73] codetools_0.2-18 stats4_4.1.1 XML_3.99-0.8
[76] glue_1.4.2 evaluate_0.14 png_0.1-7
[79] vctrs_0.3.8 nloptr_1.2.2.2 httpuv_1.6.3
[82] foreach_1.5.1 gtable_0.3.0 purrr_0.3.4
[85] clue_0.3-60 assertthat_0.2.1 cachem_1.0.6
[88] xfun_0.27 xtable_1.8-4 later_1.3.0
[91] survival_3.2-11 tibble_3.1.5 iterators_1.0.13
[94] memoise_2.0.0 AnnotationDbi_1.56.0 IRanges_2.28.0
[97] cluster_2.1.2 ellipsis_0.3.2