Results from e-QTL
by Brain region
2020-04-29
## eQTL results
get_egene_by_peer = function(tissue){
work_results = paste0("~/pd-omics/katia/Microglia/mic_314s/eQTL_mic/2nd_pass_eQTL/EUR/QTL-mapping-pipeline/results/",tissue,"_eur_expression")
peer0 = read.table(gzfile(paste0(work_results, "/peer0/",tissue,"_eur_expression_peer0.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer0) = peer0$phenotype_id
peer5 = read.table(gzfile(paste0(work_results, "/peer5/",tissue,"_eur_expression_peer5.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer5) = peer5$phenotype_id
peer10 = read.table(gzfile(paste0(work_results, "/peer10/",tissue,"_eur_expression_peer10.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer10) = peer10$phenotype_id
peer15 = read.table(gzfile(paste0(work_results, "/peer15/",tissue,"_eur_expression_peer15.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer15) = peer15$phenotype_id
peer20 = read.table(gzfile(paste0(work_results, "/peer20/",tissue,"_eur_expression_peer20.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer20) = peer20$phenotype_id
egene_by_peer = data.frame(peer = c(0,5,10,15,20), egene = c(sum(peer0$qval < 0.05),
sum(peer5$qval < 0.05),
sum(peer10$qval < 0.05),
sum(peer15$qval < 0.05),
sum(peer20$qval < 0.05)), tissue = tissue)
return(egene_by_peer)
}
MFG = get_egene_by_peer("MFG")
STG = get_egene_by_peer("STG")
THA = get_egene_by_peer("THA")
SVZ = get_egene_by_peer("SVZ")
all_egenes = rbind(MFG, STG, THA, SVZ)
# All e-genes are significative! eGenes by PEER: qval 5%
ggplot2::ggplot(all_egenes, aes(x = peer, y = egene, color=tissue)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(0, 5, 10, 15,20)) +
scale_y_continuous(breaks = c(50, 100, 150, 200, 250,300)) +
labs(x = "PEER factors", y = "Number of e-genes") +
theme_classic()
eGenes table: qval 5%
Number of eGenes by the best PEER factor.
files4eqtl = "~/pd-omics/katia/Microglia/mic_314s/eQTL_mic/files4eQTL/"
key_MFG = read.table(paste0(files4eqtl, "MFG_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_MFG$participant_id))
key_STG = read.table(paste0(files4eqtl, "STG_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_STG$participant_id))
key_THA = read.table(paste0(files4eqtl, "THA_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_THA$participant_id))
key_SVZ = read.table(paste0(files4eqtl, "SVZ_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_SVZ$participant_id))
best_peer_by_sample = all_egenes[all_egenes$peer == 10 ,] #Peer factor with more e-genes
best_peer_by_sample$num_samples = 0
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "MFG"] <- length(unique(key_MFG$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "STG"] <- length(unique(key_STG$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "THA"] <- length(unique(key_THA$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "SVZ"] <- length(unique(key_SVZ$participant_id))
best_peer_by_sample$peer[best_peer_by_sample$tissue == "SVZ"] <- 5
best_peer_by_sample$egene[best_peer_by_sample$tissue == "SVZ"] <- all_egenes$egene[all_egenes$tissue %in% "SVZ" & all_egenes$peer == 5]
createDT(best_peer_by_sample)## eQTL results
get_egene_by_peer = function(tissue){
work_results = paste0("~/pd-omics/katia/Microglia/mic_314s/eQTL_mic/2nd_pass_eQTL/EUR/QTL-mapping-pipeline/results/",tissue,"_eur_expression")
peer0 = read.table(gzfile(paste0(work_results, "/peer0/",tissue,"_eur_expression_peer0.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer0) = peer0$phenotype_id
peer5 = read.table(gzfile(paste0(work_results, "/peer5/",tissue,"_eur_expression_peer5.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer5) = peer5$phenotype_id
peer10 = read.table(gzfile(paste0(work_results, "/peer10/",tissue,"_eur_expression_peer10.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer10) = peer10$phenotype_id
peer15 = read.table(gzfile(paste0(work_results, "/peer15/",tissue,"_eur_expression_peer15.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer15) = peer15$phenotype_id
peer20 = read.table(gzfile(paste0(work_results, "/peer20/",tissue,"_eur_expression_peer20.cis_qtl.txt.gz")), stringsAsFactors=F, header = T)
rownames(peer20) = peer20$phenotype_id
egene_by_peer = data.frame(peer = c(0,5,10,15,20), egene = c(sum(peer0$qval < 0.10),
sum(peer5$qval < 0.10),
sum(peer10$qval < 0.10),
sum(peer15$qval < 0.10),
sum(peer20$qval < 0.10)), tissue = tissue)
return(egene_by_peer)
}
MFG = get_egene_by_peer("MFG")
STG = get_egene_by_peer("STG")
THA = get_egene_by_peer("THA")
SVZ = get_egene_by_peer("SVZ")
all_egenes = rbind(MFG, STG, THA, SVZ)
# All e-genes are significative! eGenes by PEER: qvalue 10%
ggplot2::ggplot(all_egenes, aes(x = peer, y = egene, color=tissue)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(0, 5, 10, 15,20)) +
scale_y_continuous(breaks = c(50, 100, 150, 200, 250,300)) +
labs(x = "PEER factors", y = "Number of e-genes") +
theme_classic()
eGenes table: qvalue 10%
files4eqtl = "~/pd-omics/katia/Microglia/mic_314s/eQTL_mic/files4eQTL/"
key_MFG = read.table(paste0(files4eqtl, "MFG_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_MFG$participant_id))
key_STG = read.table(paste0(files4eqtl, "STG_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_STG$participant_id))
key_THA = read.table(paste0(files4eqtl, "THA_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_THA$participant_id))
key_SVZ = read.table(paste0(files4eqtl, "SVZ_sample_key.txt"), header = T, stringsAsFactors = F, check.names = F)
#length(unique(key_SVZ$participant_id))
best_peer_by_sample = all_egenes[all_egenes$peer == 10 ,] #Peer factor with more e-genes
best_peer_by_sample$num_samples = 0
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "MFG"] <- length(unique(key_MFG$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "STG"] <- length(unique(key_STG$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "THA"] <- length(unique(key_THA$participant_id))
best_peer_by_sample$num_samples[best_peer_by_sample$tissue == "SVZ"] <- length(unique(key_SVZ$participant_id))
best_peer_by_sample$peer[best_peer_by_sample$tissue == "SVZ"] <- 5
best_peer_by_sample$egene[best_peer_by_sample$tissue == "SVZ"] <- all_egenes$egene[all_egenes$tissue %in% "SVZ" & all_egenes$peer == 5]
createDT(best_peer_by_sample)R version 3.6.1 (2019-07-05) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Catalina 10.15.4
Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] stats graphics grDevices utils datasets methods base
other attached packages: [1] ggplot2_3.3.0
loaded via a namespace (and not attached): [1] Rcpp_1.0.4 pillar_1.4.3 compiler_3.6.1 later_1.0.0
[5] tools_3.6.1 digest_0.6.25 jsonlite_1.6.1 evaluate_0.14
[9] lifecycle_0.2.0 tibble_2.1.3 gtable_0.3.0 pkgconfig_2.0.3
[13] rlang_0.4.5 shiny_1.4.0 crosstalk_1.0.0 yaml_2.2.0
[17] xfun_0.11 fastmap_1.0.1 withr_2.1.2 dplyr_0.8.5
[21] stringr_1.4.0 knitr_1.26 htmlwidgets_1.5.1 grid_3.6.1
[25] DT_0.13 tidyselect_0.2.5 glue_1.4.0 R6_2.4.1
[29] rmarkdown_2.0 purrr_0.3.3 farver_2.0.3 magrittr_1.5
[33] promises_1.1.0 scales_1.1.0 htmltools_0.4.0 assertthat_0.2.1 [37] xtable_1.8-4 mime_0.8 colorspace_1.4-1 httpuv_1.5.2
[41] stringi_1.4.6 munsell_0.5.0 crayon_1.3.4