KOMPUTE method testing - OF data (v10.1)
Coby Warkentin and Donghyung Lee
2023-06-19
- Load packages
- Prep control phenotype data
- Load Open Field control phenotypes
- Heatmap showing measured phenotypes
- Remove phenotypes with num of obs < 15000
- Remove samples with num of measured phenotypes < 10
- Heapmap of measured phenotypes after filtering
- Reshape the data (long to wide)
- Distribution of each phenotype
- Rank Z transformation
- Principal Variance Component Analysis
- Removing batch effects using ComBat
- PVCA on residuals from ComBat
- Compute correlations between phenotypes
- Prep IMPC summary stat
- Read OF summary stat (IMPCv10.1)
- Duplicates in gene-phenotype pair
- Using Stouffer’s method, merge multiple z-scores of a gene-phenotype pair into a z-score
- Make z-score matrix (long to wide)
- Z-score Distribution
- Estimate genetic correlation matrix between phenotypes using Zscores
- Phenotype Corr VS Genetic Corr btw phenotypes
- Test of the correlation between genetic correlation matrices
- Test KOMPUTE imputation algorithm
Last updated: 2023-06-19
Checks: 7 0
Knit directory: komputeExamples/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0.1). 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(20230110) 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 bb7f5ee. 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: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: code/.DS_Store
Untracked files:
Untracked: analysis/figures.Rmd_old
Unstaged changes:
Deleted: analysis/figures.Rmd
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_OF_v10.1.Rmd)
and HTML (docs/kompute_test_OF_v10.1.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 | 7685a09 | statsleelab | 2023-01-10 | first commit |
| html | 7685a09 | statsleelab | 2023-01-10 | first commit |
Load packages
rm(list=ls())
library(data.table)
library(dplyr)
library(reshape2)
library(ggplot2)
library(tidyr) #spread
library(RColorBrewer)
#library(irlba) # partial PCA
#library(cowplot)
library(circlize)
library(ComplexHeatmap)Prep control phenotype data
Load Open Field control phenotypes
OF.data <- readRDS("data/OF.data.rds")
dim(OF.data)[1] 344844 10Heatmap showing measured phenotypes
This heatmaps show phenotypes measured for each control mouse. Columns represent mice and rows represent phenotypes.
mtest <- table(OF.data$proc_param_name_stable_id, OF.data$biological_sample_id)
mtest <-as.data.frame.matrix(mtest)
dim(mtest)[1] 51 24604if(FALSE){
nmax <-max(mtest)
library(circlize)
col_fun = colorRamp2(c(0, nmax), c("white", "red"))
col_fun(seq(0, nmax))
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)
}Remove phenotypes with num of obs < 15000
mtest <- table(OF.data$proc_param_name, OF.data$biological_sample_id)
dim(mtest)[1] 16 24604#head(mtest[,1:10])
mtest0 <- mtest>0
#head(mtest0[,1:10])
rowSums(mtest0) OF_Center average speed OF_Center distance travelled
22856 23638
OF_Center permanence time OF_Center resting time
24599 17698
OF_Distance travelled - total OF_Latency to center entry
21805 17883
OF_Number of center entries OF_Number of rears - total
17892 12814
OF_Percentage center time OF_Periphery average speed
21730 22857
OF_Periphery distance travelled OF_Periphery permanence time
23639 24600
OF_Periphery resting time OF_Whole arena average speed
17699 24603
OF_Whole arena permanence OF_Whole arena resting time
23821 24593 rmv.pheno.list <- rownames(mtest)[rowSums(mtest0)<15000]
#rmv.pheno.list
dim(OF.data)[1] 344844 10OF.data <- OF.data %>% filter(!(proc_param_name %in% rmv.pheno.list))
dim(OF.data)[1] 332030 10# number of phenotypes left
length(unique(OF.data$proc_param_name))[1] 15Remove samples with num of measured phenotypes < 10
mtest <- table(OF.data$proc_param_name, OF.data$biological_sample_id)
dim(mtest)[1] 15 24604head(mtest[,1:10])
21653 21713 21742 21745 21747 21751 21753 21756
OF_Center average speed 1 1 1 1 1 1 1 1
OF_Center distance travelled 1 1 1 1 1 1 1 1
OF_Center permanence time 1 1 1 1 1 1 1 1
OF_Center resting time 1 1 1 1 1 1 1 1
OF_Distance travelled - total 0 0 0 0 0 0 0 0
OF_Latency to center entry 1 1 1 1 1 1 1 1
21759 21800
OF_Center average speed 1 1
OF_Center distance travelled 1 1
OF_Center permanence time 1 1
OF_Center resting time 1 1
OF_Distance travelled - total 0 0
OF_Latency to center entry 1 1mtest0 <- mtest>0
head(mtest0[,1:10])
21653 21713 21742 21745 21747 21751 21753 21756
OF_Center average speed TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
OF_Center distance travelled TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
OF_Center permanence time TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
OF_Center resting time TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
OF_Distance travelled - total FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
OF_Latency to center entry TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
21759 21800
OF_Center average speed TRUE TRUE
OF_Center distance travelled TRUE TRUE
OF_Center permanence time TRUE TRUE
OF_Center resting time TRUE TRUE
OF_Distance travelled - total FALSE FALSE
OF_Latency to center entry TRUE TRUEsummary(colSums(mtest0)) Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 11.00 15.00 13.41 15.00 15.00 rmv.sample.list <- colnames(mtest)[colSums(mtest0)<10]
length(rmv.sample.list)[1] 1747dim(OF.data)[1] 332030 10OF.data <- OF.data %>% filter(!(biological_sample_id %in% rmv.sample.list))
dim(OF.data)[1] 319816 10# number of observations to use
length(unique(OF.data$biological_sample_id))[1] 22857Heapmap of measured phenotypes after filtering
if(FALSE){
mtest <- table(OF.data$proc_param_name, OF.data$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_OF.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()
}Reshape the data (long to wide)
OF.mat <- OF.data %>%
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(OF.mat) sex phenotyping_center strain_name
21653 female WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
21713 female WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
21742 male WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
21745 male WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
21747 male WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
21751 male WTSI C57BL/6Brd-Tyr<c-Brd> * C57BL/6N
OF_Center average speed OF_Center distance travelled
21653 51.5 4259
21713 40.1 3266
21742 51.0 710
21745 38.0 2580
21747 31.4 3022
21751 14.9 1723
OF_Center permanence time OF_Center resting time
21653 102 20
21713 87 6
21742 17 3
21745 100 33
21747 134 39
21751 240 130
OF_Distance travelled - total OF_Latency to center entry
21653 NA 5.0
21713 NA 6.1
21742 NA 24.0
21745 NA 8.3
21747 NA 6.6
21751 NA 18.7
OF_Number of center entries OF_Percentage center time
21653 193 NA
21713 221 NA
21742 73 NA
21745 165 NA
21747 210 NA
21751 75 NA
OF_Periphery average speed OF_Periphery distance travelled
21653 35.3 12923
21713 31.9 11625
21742 19.1 5769
21745 25.1 7910
21747 26.0 8208
21751 11.5 2150
OF_Periphery permanence time OF_Periphery resting time
21653 498 135
21713 513 152
21742 583 292
21745 500 191
21747 466 154
21751 360 184
OF_Whole arena average speed OF_Whole arena permanence
21653 38.3 600
21713 33.4 600
21742 20.5 600
21745 27.4 600
21747 27.2 600
21751 12.8 600
OF_Whole arena resting time
21653 155
21713 158
21742 295
21745 224
21747 193
21751 314dim(OF.mat)[1] 22857 18summary(colSums(is.na(OF.mat[,-1:-3]))) Min. 1st Qu. Median Mean 3rd Qu. Max.
0 0 1 1677 3882 5159 Distribution of each phenotype
ggplot(melt(OF.mat), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Rank Z transformation
library(RNOmni)
OF.mat.rank <- OF.mat
dim(OF.mat.rank)[1] 22857 18OF.mat.rank <- OF.mat.rank[complete.cases(OF.mat.rank),]
dim(OF.mat.rank)[1] 14863 18dim(OF.mat)[1] 22857 18OF.mat <- OF.mat[complete.cases(OF.mat),]
dim(OF.mat)[1] 14863 18OF.mat.rank <- cbind(OF.mat.rank[,1:3], apply(OF.mat.rank[,-1:-3], 2, RankNorm))
ggplot(melt(OF.mat.rank), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Principal Variance Component Analysis
Here we conducted a PVCA analysis on the phenotype matrix data and measure the proportion of variance explained by each important covariate (sex, phenotyping_center).
source("code/PVCA.R")
meta <- OF.mat.rank[,1:3] ## looking at covariates sex, phenotyping_center, and strain_name
head(meta) sex phenotyping_center strain_name
39638 female MRC Harwell C57BL/6NTac
39639 female HMGU C57BL/6NCrl
39640 female HMGU C57BL/6NTac
39641 male HMGU C57BL/6NCrl
39642 female MRC Harwell C57BL/6NTac
39643 female HMGU C57BL/6NCrldim(meta)[1] 14863 3summary(meta) # variables are still characters sex phenotyping_center strain_name
Length:14863 Length:14863 Length:14863
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:7428 MRC Harwell:4532 C57BL/6N :3655
male :7435 HMGU :3119 C57BL/6NCrl:3510
ICS :2417 C57BL/6NJcl: 459
RBRC :1323 C57BL/6NTac:7239
CCP-IMG :1141
TCP :1093
(Other) :1238 chisq.test(meta[,1],meta[,2])
Pearson's Chi-squared test
data: meta[, 1] and meta[, 2]
X-squared = 3.7637, df = 7, p-value = 0.8066chisq.test(meta[,2],meta[,3])
Pearson's Chi-squared test
data: meta[, 2] and meta[, 3]
X-squared = 29526, df = 21, p-value < 2.2e-16meta<-meta[,-3] # phenotyping_center and strain_name strongly associated and this caused confouding in PVCA analysis so strain_name dropped.
G <- t(OF.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
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
png(file="docs/figure/figures.Rmd/pvca_OF_1.png", width=600, height=350)
pvca.plot
dev.off()quartz_off_screen
2 Removing batch effects using ComBat
We remove the center effect using ComBat.
library(sva)Loading required package: mgcvLoading required package: nlme
Attaching package: 'nlme'The following object is masked from 'package:lme4':
lmListThe following object is masked from 'package:dplyr':
collapseThis is mgcv 1.8-40. For overview type 'help("mgcv-package")'.Loading required package: genefilter
Attaching package: 'genefilter'The following object is masked from 'package:ComplexHeatmap':
dist2Loading required package: BiocParallelcombat_komp = ComBat(dat=G, batch=meta$phenotyping_center, par.prior=TRUE, prior.plots=TRUE, mod=NULL)Found 1 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.Found8batchesAdjusting for0covariate(s) or covariate level(s)Standardizing Data across genesFitting L/S model and finding priorsFinding parametric adjustmentsAdjusting the Data
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
combat_komp[1:5,1:5] 39638 39639 39640 39641
OF_Center average speed 0.59446136 -0.1046728 0.12751142 0.2488662
OF_Center distance travelled 0.49851881 -0.3961273 0.26323301 -0.1860701
OF_Center permanence time 0.05536160 -0.5740468 0.28108074 -0.3840993
OF_Center resting time 0.09818595 -0.9419518 0.58940330 -0.4078182
OF_Distance travelled - total 0.03620577 -0.0705752 0.03808406 0.4042278
39642
OF_Center average speed -0.16244236
OF_Center distance travelled -0.98532415
OF_Center permanence time -0.81014937
OF_Center resting time 0.05024321
OF_Distance travelled - total -0.11281875G[1:5,1:5] # for comparison, combat_komp is same form and same dimensions as G 39638 39639 39640 39641
OF_Center average speed 0.425513501 1.234197021 1.5167530 1.6644355
OF_Center distance travelled 0.129374222 0.997641495 1.5808081 1.1834251
OF_Center permanence time -0.515972599 0.009697299 0.8574740 0.1980120
OF_Center resting time 0.001686461 -1.202492878 0.3483423 -0.6615647
OF_Distance travelled - total -0.470073474 1.298879853 1.4120914 1.7935746
39642
OF_Center average speed -0.44088826
OF_Center distance travelled -1.46426536
OF_Center permanence time -1.37773801
OF_Center resting time -0.04360882
OF_Distance travelled - total -0.61963422PVCA on residuals from ComBat
The center effect should be much lower.
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)
}
# 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
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
png(file="docs/figure/figures.Rmd/pvca_OF_2.png", width=600, height=350)
pvca.plot.nobatch
dev.off()quartz_off_screen
2 Compute correlations between phenotypes
OF.cor.rank <- cor(OF.mat.rank[,-1:-3], use="pairwise.complete.obs") # pearson correlation coefficient
OF.cor <- cor(OF.mat[,-1:-3], use="pairwise.complete.obs", method="spearman") # spearman
OF.cor.combat <- cor(t(combat_komp), use="pairwise.complete.obs")
pheno.list <- rownames(OF.cor)
ht1 = Heatmap(OF.cor, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Spearm. Corr.")
draw(ht1)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
ht2 = Heatmap(OF.cor.rank, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Corr. RankZ")
draw(ht2)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
ht3 = Heatmap(OF.cor.combat, show_column_names = F, row_names_gp = gpar(fontsize = 9), name="Corr. ComBat")
draw(ht3)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Prep IMPC summary stat
Read OF summary stat (IMPCv10.1)
OF.stat <- readRDS("data/OF.stat.rds")
dim(OF.stat)[1] 45326 8table(OF.stat$parameter_name, OF.stat$procedure_name)
OF
Center average speed 3335
Center distance travelled 3591
Center permanence time 3632
Center resting time 1802
Distance travelled - total 3550
Latency to center entry 1818
Number of center entries 1816
Number of rears - total 2772
Percentage center time 3375
Periphery average speed 3335
Periphery distance travelled 3591
Periphery permanence time 3633
Periphery resting time 1803
Whole arena average speed 3637
Whole arena resting time 3636length(unique(OF.stat$marker_symbol)) #3362[1] 3362length(unique(OF.stat$allele_symbol)) #3412[1] 3412length(unique(OF.stat$proc_param_name)) #15 # number of phenotypes in association statistics data set[1] 15length(unique(OF.data$proc_param_name)) #15 # number of phenotypes in final control data[1] 15pheno.list.stat <- unique(OF.stat$proc_param_name)
pheno.list.ctrl <- unique(OF.data$proc_param_name)
sum(pheno.list.stat %in% pheno.list.ctrl)[1] 14sum(pheno.list.ctrl %in% pheno.list.stat)[1] 14## extract common phenotype list
common.pheno.list <- sort(intersect(pheno.list.ctrl, pheno.list.stat))
common.pheno.list [1] "OF_Center average speed" "OF_Center distance travelled"
[3] "OF_Center permanence time" "OF_Center resting time"
[5] "OF_Distance travelled - total" "OF_Latency to center entry"
[7] "OF_Number of center entries" "OF_Percentage center time"
[9] "OF_Periphery average speed" "OF_Periphery distance travelled"
[11] "OF_Periphery permanence time" "OF_Periphery resting time"
[13] "OF_Whole arena average speed" "OF_Whole arena resting time" length(common.pheno.list) # 14 - each data set had one phenotype not present in the other[1] 14## Use summary statistics of common phenotypes
dim(OF.stat)[1] 45326 8OF.stat <- OF.stat %>% filter(proc_param_name %in% common.pheno.list)
dim(OF.stat)[1] 42554 8length(unique(OF.stat$proc_param_name))[1] 14Duplicates in gene-phenotype pair
mtest <- table(OF.stat$proc_param_name, OF.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" "#FFEEE7FF" "#FFDCD0FF" "#FFCBB9FF" "#FFB9A2FF" "#FFA78CFF"
[7] "#FF9576FF" "#FF8161FF" "#FF6D4CFF" "#FF5636FF" "#FF3A1FFF" "#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")`use_raster` is automatically set to TRUE for a matrix with more than
2000 columns You can control `use_raster` argument by explicitly
setting TRUE/FALSE to it.
Set `ht_opt$message = FALSE` to turn off this message.'magick' package is suggested to install to give better rasterization.
Set `ht_opt$message = FALSE` to turn off this message.draw(ht)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Using Stouffer’s method, merge multiple z-scores of a gene-phenotype pair into a z-score
## sum(z-score)/sqrt(# of zscore)
sumz <- function(z){ sum(z)/sqrt(length(z)) }
OF.z = OF.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(OF.z)[1] 35836 3Make z-score matrix (long to wide)
nan2na <- function(df){
out <- data.frame(sapply(df, function(x) ifelse(is.nan(x), NA, x)))
colnames(out) <- colnames(df)
out
}
OF.zmat = dcast(OF.z, marker_symbol ~ proc_param_name, value.var = "zscore",
fun.aggregate = mean) %>% tibble::column_to_rownames(var="marker_symbol")
OF.zmat = nan2na(OF.zmat) #convert nan to na
dim(OF.zmat)[1] 3360 14dim(OF.zmat)[1] 3360 14saveRDS(OF.zmat, file = "data/OF.zmat.rds")
id.mat <- 1*(!is.na(OF.zmat)) # multiply 1 to make this matrix numeric
nrow(as.data.frame(colSums(id.mat)))[1] 14dim(id.mat)[1] 3360 14## heatmap of gene - phenotype (red: tested, white: untested)
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")`use_raster` is automatically set to TRUE for a matrix with more than
2000 columns You can control `use_raster` argument by explicitly
setting TRUE/FALSE to it.
Set `ht_opt$message = FALSE` to turn off this message.'magick' package is suggested to install to give better rasterization.
Set `ht_opt$message = FALSE` to turn off this message.draw(ht)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Z-score Distribution
We plot association Z-score distribution for each phenotype.
ggplot(melt(OF.zmat), aes(x=value)) +
geom_histogram() +
facet_wrap(~variable, scales="free", ncol=5)+
theme(strip.text.x = element_text(size = 6))
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Estimate genetic correlation matrix between phenotypes using Zscores
Here, we estimate the genetic correlations between phenotypes using association Z-score matrix (num of genes:3983, num of phenotypes 19).
OF.zmat <- OF.zmat[,common.pheno.list]
OF.zcor = cor(OF.zmat, use="pairwise.complete.obs")
ht = Heatmap(OF.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)"
)
draw(ht)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
Phenotype Corr VS Genetic Corr btw phenotypes
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.
OF.cor.rank.fig <- OF.cor.rank[common.pheno.list,common.pheno.list]
OF.cor.fig <- OF.cor[common.pheno.list,common.pheno.list]
OF.cor.combat.fig <- OF.cor.combat[common.pheno.list, common.pheno.list]
OF.zcor.fig <- OF.zcor
ht = Heatmap(OF.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)
OF.cor.rank.fig <- OF.cor.rank.fig[pheno.order,pheno.order]
ht1 = Heatmap(OF.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")
OF.cor.fig <- OF.cor.fig[pheno.order,pheno.order]
ht2 = Heatmap(OF.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")
OF.cor.combat.fig <- OF.cor.combat.fig[pheno.order,pheno.order]
ht3 = Heatmap(OF.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")
OF.zcor.fig <- OF.zcor.fig[pheno.order,pheno.order]
ht4 = Heatmap(OF.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)Warning: Heatmap/annotation names are duplicated: CorrWarning: Heatmap/annotation names are duplicated: Corr, CorrWarning: Heatmap/annotation names are duplicated: Corr, Corr, Corr
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
png(file="docs/figure/figures.Rmd/cors_OF.png", width=800, height=250)
draw(ht1+ht2+ht3+ht4)Warning: Heatmap/annotation names are duplicated: CorrWarning: Heatmap/annotation names are duplicated: Corr, CorrWarning: Heatmap/annotation names are duplicated: Corr, Corr, Corrdev.off()quartz_off_screen
2 Test of the correlation between genetic correlation matrices
We use the Mantel’s test for 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(OF.cor.fig)
a2 <- to.upper(OF.cor.rank.fig)
a3 <- to.upper(OF.cor.combat.fig)
a4 <- to.upper(OF.zcor.fig)
plot(a4, a1)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
plot(a4, a2)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
plot(a4, a3)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
mantel.rtest(as.dist(1-OF.cor.fig), as.dist(1-OF.zcor.fig), nrepet = 9999) #nrepet = number of permutationsMonte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.9461638
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
8.253839656 -0.000564017 0.013156434 mantel.rtest(as.dist(1-OF.cor.rank.fig), as.dist(1-OF.zcor.fig), nrepet = 9999)Monte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.9513592
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
8.1916710930 0.0009080079 0.0134621564 mantel.rtest(as.dist(1-OF.cor.combat.fig), as.dist(1-OF.zcor.fig), nrepet = 9999)Monte-Carlo test
Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
Observation: 0.9660052
Based on 9999 replicates
Simulated p-value: 1e-04
Alternative hypothesis: greater
Std.Obs Expectation Variance
8.835872572 -0.001694674 0.011994508 Test KOMPUTE imputation algorithm
Load KOMPUTE package
if(!"kompute" %in% rownames(installed.packages())){
library(devtools)
devtools::install_github("dleelab/kompute")
}
library(kompute)Simulation study - imputed vs measured
We randomly select measured gene-phenotype association z-scores, mask those, impute them using KOMPUTE method. Then we compare the imputed z-scores to the measured ones.
zmat <-t(OF.zmat)
dim(zmat)[1] 14 3360## filter genes with na < 10
zmat0 <- is.na(zmat)
num.na<-colSums(zmat0)
summary(num.na) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 4.000 3.335 5.000 13.000 zmat <- zmat[,num.na<10]
dim(zmat)[1] 14 3244#pheno.cor <- OF.cor.fig
#pheno.cor <- OF.cor.rank.fig
pheno.cor <- OF.cor.combat.fig
#pheno.cor <- OF.zcor.fig
zmat <- zmat[rownames(pheno.cor),,drop=FALSE]
rownames(zmat) [1] "OF_Center distance travelled" "OF_Number of center entries"
[3] "OF_Whole arena average speed" "OF_Periphery average speed"
[5] "OF_Distance travelled - total" "OF_Periphery distance travelled"
[7] "OF_Center average speed" "OF_Percentage center time"
[9] "OF_Center permanence time" "OF_Center resting time"
[11] "OF_Latency to center entry" "OF_Periphery permanence time"
[13] "OF_Whole arena resting time" "OF_Periphery resting time" rownames(pheno.cor) [1] "OF_Center distance travelled" "OF_Number of center entries"
[3] "OF_Whole arena average speed" "OF_Periphery average speed"
[5] "OF_Distance travelled - total" "OF_Periphery distance travelled"
[7] "OF_Center average speed" "OF_Percentage center time"
[9] "OF_Center permanence time" "OF_Center resting time"
[11] "OF_Latency to center entry" "OF_Periphery permanence time"
[13] "OF_Whole arena resting time" "OF_Periphery resting time" colnames(pheno.cor) [1] "OF_Center distance travelled" "OF_Number of center entries"
[3] "OF_Whole arena average speed" "OF_Periphery average speed"
[5] "OF_Distance travelled - total" "OF_Periphery distance travelled"
[7] "OF_Center average speed" "OF_Percentage center time"
[9] "OF_Center permanence time" "OF_Center resting time"
[11] "OF_Latency to center entry" "OF_Periphery permanence time"
[13] "OF_Whole arena resting time" "OF_Periphery resting time" npheno <- nrow(zmat)
## percentage of missing Z-scores in the original data
100*sum(is.na(zmat))/(nrow(zmat)*ncol(zmat)) # 21%[1] 21.92179nimp <- 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)Run KOMPUTE method
kompute.res <- kompute(zmat.imp, pheno.cor, 0.01)
KOMPute running...# of genes: 3244# of phenotypes: 14# of imputed Z-scores: 11956# measured vs imputed
length(org.z)[1] 2000imp.z <- as.matrix(kompute.res$zmat)[mask.i]
imp.info <- as.matrix(kompute.res$infomat)[mask.i]
plot(imp.z, org.z)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
imp <- data.frame(org.z=org.z, imp.z=imp.z, info=imp.info)
dim(imp)[1] 2000 3imp <- imp[complete.cases(imp),]
imp <- subset(imp, info>=0 & info <= 1)
dim(imp)[1] 2000 3cor.val <- round(cor(imp$imp.z, imp$org.z), digits=3)
cor.val[1] 0.938plot(imp$imp.z, imp$org.z)
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
info.cutoff <- 0.8
imp.sub <- subset(imp, info>info.cutoff)
dim(imp.sub)[1] 1573 3summary(imp.sub$imp.z) Min. 1st Qu. Median Mean 3rd Qu. Max.
-8.10593 -0.98431 -0.05220 -0.06551 0.84600 8.89281 summary(imp.sub$info) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.8000 0.8825 0.9284 0.9171 0.9633 0.9878 cor.val <- round(cor(imp.sub$imp.z, imp.sub$org.z), digits=3)
cor.val[1] 0.968g <- ggplot(imp.sub, aes(x=imp.z, y=org.z, col=info)) +
geom_point() +
labs(title=paste0("IMPC Behavior Data (OF), Info>", info.cutoff, ", Cor=",cor.val),
x="Imputed Z-scores", y = "Measured Z-scores", col="Info") +
theme_minimal()
g
| Version | Author | Date |
|---|---|---|
| 7685a09 | statsleelab | 2023-01-10 |
# save plot
png(file="docs/figure/figures.Rmd/sim_results_OF.png", width=600, height=350)
g
dev.off()quartz_off_screen
2
sessionInfo()R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] kompute_0.1.0 ade4_1.7-20 sva_3.44.0
[4] BiocParallel_1.30.3 genefilter_1.78.0 mgcv_1.8-40
[7] nlme_3.1-158 lme4_1.1-31 Matrix_1.5-1
[10] RNOmni_1.0.1 ComplexHeatmap_2.12.1 circlize_0.4.15
[13] RColorBrewer_1.1-3 tidyr_1.2.0 ggplot2_3.4.1
[16] reshape2_1.4.4 dplyr_1.0.9 data.table_1.14.2
[19] workflowr_1.7.0.1
loaded via a namespace (and not attached):
[1] minqa_1.2.5 colorspace_2.1-0 rjson_0.2.21
[4] rprojroot_2.0.3 XVector_0.36.0 GlobalOptions_0.1.2
[7] fs_1.5.2 clue_0.3-62 rstudioapi_0.13
[10] farver_2.1.1 bit64_4.0.5 AnnotationDbi_1.58.0
[13] fansi_1.0.4 codetools_0.2-18 splines_4.2.1
[16] doParallel_1.0.17 cachem_1.0.6 knitr_1.39
[19] jsonlite_1.8.0 nloptr_2.0.3 annotate_1.74.0
[22] cluster_2.1.3 png_0.1-8 compiler_4.2.1
[25] httr_1.4.3 assertthat_0.2.1 fastmap_1.1.0
[28] limma_3.52.4 cli_3.6.0 later_1.3.0
[31] htmltools_0.5.3 tools_4.2.1 gtable_0.3.1
[34] glue_1.6.2 GenomeInfoDbData_1.2.8 Rcpp_1.0.10
[37] Biobase_2.56.0 jquerylib_0.1.4 vctrs_0.5.2
[40] Biostrings_2.64.0 iterators_1.0.14 xfun_0.31
[43] stringr_1.4.0 ps_1.7.1 lifecycle_1.0.3
[46] XML_3.99-0.10 edgeR_3.38.4 getPass_0.2-2
[49] MASS_7.3-58.1 zlibbioc_1.42.0 scales_1.2.1
[52] promises_1.2.0.1 parallel_4.2.1 yaml_2.3.5
[55] memoise_2.0.1 sass_0.4.2 stringi_1.7.8
[58] RSQLite_2.2.15 highr_0.9 S4Vectors_0.34.0
[61] foreach_1.5.2 BiocGenerics_0.42.0 boot_1.3-28
[64] shape_1.4.6 GenomeInfoDb_1.32.3 rlang_1.0.6
[67] pkgconfig_2.0.3 matrixStats_0.62.0 bitops_1.0-7
[70] evaluate_0.16 lattice_0.20-45 purrr_0.3.4
[73] labeling_0.4.2 bit_4.0.4 processx_3.7.0
[76] tidyselect_1.2.0 plyr_1.8.7 magrittr_2.0.3
[79] R6_2.5.1 IRanges_2.30.0 generics_0.1.3
[82] DBI_1.1.3 pillar_1.8.1 whisker_0.4
[85] withr_2.5.0 survival_3.3-1 KEGGREST_1.36.3
[88] RCurl_1.98-1.8 tibble_3.1.8 crayon_1.5.1
[91] utf8_1.2.3 rmarkdown_2.14 GetoptLong_1.0.5
[94] locfit_1.5-9.6 blob_1.2.3 callr_3.7.1
[97] git2r_0.30.1 digest_0.6.29 xtable_1.8-4
[100] httpuv_1.6.5 stats4_4.2.1 munsell_0.5.0
[103] bslib_0.4.0