Load required packages
suppressPackageStartupMessages(library(package = "knitr")) suppressPackageStartupMessages(library(package = "Biobase")) suppressPackageStartupMessages(library(package = "biomaRt")) suppressPackageStartupMessages(library(package = "ggbeeswarm")) suppressPackageStartupMessages(library(package = "tidyverse"))
Set session options
dataCacheDir <- here::here("data_cache", "2021_03_08") opts_chunk$set(tidy = FALSE, fig.path = "Figure/") options(stringsAsFactors = FALSE, width = 80)
Read virtual study from Evan
geFile <- list.files(path = file.path(workDir, "Input"), pattern = "noResponse_norm_young.rds", full.names = TRUE) esetYngNoResp <- readRDS(file = geFile) geFile <- list.files(path = file.path(workDir, "Input"), pattern = "noResponse_norm_older.rds", full.names = TRUE) esetOldNoResp <- readRDS(file = geFile)
Combine young and elderly ExpressionSet
featuresCommon <- intersect(featureNames(esetYngNoResp), featureNames(esetOldNoResp)) eset <- ExpressionSet(assayData = cbind(exprs(esetYngNoResp)[featuresCommon, ], exprs(esetOldNoResp)[featuresCommon, ]), phenoData = AnnotatedDataFrame(data = rbind(pData(esetYngNoResp), pData(esetOldNoResp))))
Average chrY genes intensities
human <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl") # create gene 2 chromosome data.frame gene2chr <- getBM(attributes = c("hgnc_symbol", "chromosome_name"), filters = "hgnc_symbol", values = featuresCommon, mart = human) chrYgenes <- filter(gene2chr, chromosome_name %in% "Y") %>% .$hgnc_symbol %>% as.vector()
Plot chrY average as function of gender
plotDF <- colMeans(exprs(eset)[chrYgenes, ], na.rm = TRUE) %>% data.frame(chry = .) %>% rownames_to_column() %>% merge(y = pData(eset), by.x = "rowname", by.y = "uid") %>% mutate(matrix = ifelse(test = study_accession %in% "SDY1328" & age_reported < 50, yes = paste0(matrix, "_young"), no = paste0(matrix, "_old"))) # exclude studies with only one gender or missing gender annotation singleGender <- plotDF %>% group_by(study_accession, matrix) %>% summarize(n.gender = length(unique(setdiff(gender, c("Not Specified", "Unknown"))))) %>% filter(n.gender <= 1) %>% mutate(matrix = as.vector(matrix)) singleGender %>% select(study_accession, n.gender) %>% distinct() %>% arrange(n.gender, study_accession) %>% print() # find outlier (1.5 x IGR from Q1 and Q3) outlierDF <- plotDF %>% filter(gender != "Not Specified") %>% group_by(study_accession, matrix, gender) %>% mutate(up = chry >= quantile(chry, probs = 0.75) + 1.5 * IQR(chry), dn = chry <= quantile(chry, probs = 0.25) - 1.5 * IQR(chry)) quartileDF <- plotDF %>% filter(gender != "Not Specified") %>% group_by(study_accession, matrix, gender) %>% summarize(q1 = quantile(chry, probs = 0.25) - 1.5 * IQR(chry), q3 = quantile(chry, probs = 0.75) + 1.5 * IQR(chry)) %>% mutate(gender = c("Female" = "Male", "Male" = "Female")[gender]) # flag outlier samples (possible swap) flagDF <- merge(x = outlierDF, y = quartileDF, by = c("study_accession", "matrix", "gender")) %>% mutate(flag = FALSE, flag = ifelse(test = gender %in% "Female" & up & chry >= q1, yes = TRUE, no = flag), flag = ifelse(test = gender %in% "Male" & dn & chry <= q3, yes = TRUE, no = flag)) plotDF <- merge(x = plotDF, y = select(flagDF, rowname, flag), by = "rowname", all.x = TRUE) %>% # if gender not specified, set flag as false mutate(flag = ifelse(test = is.na(flag), yes = FALSE, no = flag)) ggplot(data = filter(plotDF, !(matrix %in% singleGender$matrix)), mapping = aes(x = gender, y = chry)) + geom_boxplot(outlier.color = "transparent", fill = "grey") + geom_jitter(height = 0, width = 0.25, mapping = aes(color = flag)) + labs(y = "Average probe intensities (chrY)") + scale_color_manual(values = c("FALSE" = "black", "TRUE" = "red")) + facet_wrap(facets = ~study_accession+matrix, scale = "free") + theme_bw() + theme(plot.title = element_text(hjust = 0.5), legend.pos = "none", strip.text = element_text(size = 5))
Print possible possible mislabelled samples
filter(plotDF, flag) %>% select(rowname, chry, matrix, study_accession) %>% rename(uid = rowname) %>% arrange(study_accession, uid) %>% print()
Plot jitter plot for four example studies
plotTemp <- filter(plotDF, (study_accession %in% "SDY1260" & matrix %in% "SDY1260_PBMC_MCV4_Geo_old") | (study_accession %in% "SDY400" & matrix %in% "SDY400_PBMC_Young_Geo_old") | (study_accession %in% "SDY984" & matrix %in% "SDY984_PBMC_Elderly_Geo_old") | study_accession %in% "SDY1291") ggplot(data = plotTemp, mapping = aes(x = gender, y = chry)) + geom_boxplot(outlier.color = "transparent", fill = "grey") + geom_beeswarm(mapping = aes(color = flag), cex = 2.5, size = 0.8) + labs(y = "Average probe intensities (chrY)", x = "Sex") + scale_color_manual(values = c("FALSE" = "black", "TRUE" = "red")) + facet_wrap(facets = ~study_accession+matrix, scale = "free") + theme_bw() + theme(plot.title = element_text(hjust = 0.5), legend.pos = "none", strip.text = element_text(size = 6))
Print session info
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.