# 類別變數 (多PRSmix.方法)
# prs_wt <- read.delim(
# "C:/Users/yuting/Desktop/TPMI_GWAS_Stats/PRS_Weight/250.2_prs_weights.tsv",
# stringsAsFactors = FALSE)
# dim(prs_wt) # 1022838 7
#
# all156 <- read.delim(
# "C:/Users/yuting/Desktop/TPMI_GWAS_Stats//GWAS_Stats/250.2_all.156.txt",
# stringsAsFactors = FALSE)
# dim(all156) # 1739048 17
# all156$ES = log(all156$OR)
# all156_subset <- all156[, c("rsID", "X.CHROM", "POS", "A1", "OMITTED", "ES","P")]
#
# prs_merged <- merge(
# prs_wt,
# all156_subset,
# by.x = "MarkerID",
# by.y = "rsID",
# all.x = TRUE
# )
# dim(prs_merged) # 1022838 11
# prs_merged$ID <- paste0(prs_merged$X.CHROM, ":", prs_merged$POS)
# prs_merged <- prs_merged[!is.na(prs_merged$X.CHROM), ]
#
# prs_methods <- c("LDpred2", "Lassosum2", "PRS.CS", "SBayesR", "MegaPRS","PRSmix.")
#
# # ======================================
# # 3️⃣ 生成各方法 effect allele，並翻轉 beta
# # ======================================
# for (m in prs_methods) {
#
# ea_col <- paste0("EA_", m)
# beta <- prs_merged[[m]]
#
# # 預設 effect allele = Effect_Allele
# prs_merged[[ea_col]] <- prs_merged$Effect_Allele
#
# # beta < 0 的 SNP
# idx <- which(!is.na(beta) & beta < 0)
#
# if (length(idx) > 0) {
# # 翻轉 effect allele：選擇 A1 / OMITTED 中不等於 Effect_Allele 的
# alt <- ifelse(
# prs_merged$A1[idx] != prs_merged$Effect_Allele[idx],
# prs_merged$A1[idx],
# prs_merged$OMITTED[idx]
# )
# prs_merged[[ea_col]][idx] <- alt
#
# # 翻正 beta
# prs_merged[[m]][idx] <- -prs_merged[[m]][idx]
# }
# }
#
# # ======================================
# # 4️⃣ ES 處理（只看 A1 / OMITTED）
# # ======================================
# prs_merged$GWAS_ES <- prs_merged$A1
# idx_es <- which(!is.na(prs_merged$ES) & prs_merged$ES < 0)
#
# prs_merged$ES[idx_es] <- -prs_merged$ES[idx_es]
# prs_merged$GWAS_ES[idx_es] <- prs_merged$OMITTED[idx_es]
#
# # ======================================
# # 5️⃣ 安全檢查
# # ======================================
# stopifnot(
# all(prs_merged$A1 != prs_merged$OMITTED) # 確保兩個 allele 不同
# )
#
# # ======================================
# # 6️⃣ 輸出各方法 score 檔（可直接用 PLINK2）
# # ======================================
# for (m in prs_methods) {
# out <- prs_merged[, c("ID", paste0("EA_", m), m)]
# colnames(out) <- c("ID", "A1", "BETA")
# write.table(
# out,
# file = paste0("score_", m, ".txt"),
# sep = "\t",
# quote = FALSE,
# row.names = FALSE
# )
# }
#
# names(prs_merged)
# prs_merged = prs_merged[,c(15,22,13,14,16,3,17,4,18,5,19,6,20,7,21,8,1,9,10)]
# names(prs_merged)
#
#
# write.table(
# prs_merged,
# file = "C:/Users/yuting/Desktop/TPMI_GWAS_Stats/250.2_prs_weights_risk_PRS.txt",
# sep = "\t",
# quote = FALSE,
# row.names = FALSE
# )
#
#
#
#
# # 數值變數
# prs_wt <- read.delim(
# "C:/Users/yuting/Desktop/TPMI_GWAS_Stats/PRS_Weight/BUN_prs_weights.tsv",
# stringsAsFactors = FALSE)
# dim(prs_wt) # 1022838 7
#
# all156 <- read.delim(
# "C:/Users/yuting/Desktop/TPMI_GWAS_Stats//GWAS_Stats/BUN_all.156.txt",
# stringsAsFactors = FALSE)
# dim(all156) # 1739048 17
# # all156$ES = log(all156$OR)
# all156_subset <- all156[, c("rsID", "X.CHROM", "POS", "A1", "OMITTED", "BETA","P")]
#
# prs_merged <- merge(
# prs_wt,
# all156_subset,
# by.x = "MarkerID",
# by.y = "rsID",
# all.x = TRUE
# )
# dim(prs_merged) # 1022838 11
# prs_merged$ID <- paste0(prs_merged$X.CHROM, ":", prs_merged$POS)
# prs_merged <- prs_merged[!is.na(prs_merged$X.CHROM), ]
#
# prs_methods <- c("LDpred2", "Lassosum2", "PRS.CS", "SBayesR", "MegaPRS","PRSmix.")
#
# # ======================================
# # 3️⃣ 生成各方法 effect allele，並翻轉 beta
# # ======================================
# for (m in prs_methods) {
#
# ea_col <- paste0("EA_", m)
# beta <- prs_merged[[m]]
#
# # 預設 effect allele = Effect_Allele
# prs_merged[[ea_col]] <- prs_merged$Effect_Allele
#
# # beta < 0 的 SNP
# idx <- which(!is.na(beta) & beta < 0)
#
# if (length(idx) > 0) {
# # 翻轉 effect allele：選擇 A1 / OMITTED 中不等於 Effect_Allele 的
# alt <- ifelse(
# prs_merged$A1[idx] != prs_merged$Effect_Allele[idx],
# prs_merged$A1[idx],
# prs_merged$OMITTED[idx]
# )
# prs_merged[[ea_col]][idx] <- alt
#
# # 翻正 beta
# prs_merged[[m]][idx] <- -prs_merged[[m]][idx]
# }
# }
#
# # ======================================
# # 4️⃣ ES 處理（只看 A1 / OMITTED）
# # ======================================
# prs_merged$GWAS_ES <- prs_merged$A1
# idx_es <- which(!is.na(prs_merged$BETA) & prs_merged$BETA < 0)
#
# prs_merged$BETA[idx_es] <- -prs_merged$BETA[idx_es]
# prs_merged$GWAS_ES[idx_es] <- prs_merged$OMITTED[idx_es]
#
# # ======================================
# # 5️⃣ 安全檢查
# # ======================================
# stopifnot(
# all(prs_merged$A1 != prs_merged$OMITTED) # 確保兩個 allele 不同
# )
#
# # ======================================
# # 6️⃣ 輸出各方法 score 檔（可直接用 PLINK2）
# # ======================================
# for (m in prs_methods) {
# out <- prs_merged[, c("ID", paste0("EA_", m), m)]
# colnames(out) <- c("ID", "A1", "BETA")
# write.table(
# out,
# file = paste0("score_", m, ".txt"),
# sep = "\t",
# quote = FALSE,
# row.names = FALSE
# )
# }
#
# names(prs_merged)
# prs_merged = prs_merged[,c(15,22,13,14,16,3,17,4,18,5,19,6,20,7,21,8,1,9,10)]
#
#
# write.table(
# prs_merged,
# file = "C:/Users/yuting/Desktop/TPMI_GWAS_Stats/BUN_prs_weights_risk_PRS.txt",
# sep = "\t",
# quote = FALSE,
# row.names = FALSE
# )
############################
## Load libraries
############################
library(dplyr)
## Warning: 套件 'dplyr' 是用 R 版本 4.4.3 來建造的
##
## 載入套件：'dplyr'
## 下列物件被遮斷自 'package:stats':
##
## filter, lag
## 下列物件被遮斷自 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
## Warning: 套件 'ggplot2' 是用 R 版本 4.4.3 來建造的
############################
## Input files (each PRS method)
############################
files <- c(
Paper = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_Paper.chrALL.score.Merge.txt",
MegaPRS = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_MegaPRS.chrALL.score.Merge.txt",
Lassosum2 = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_Lassosum2.chrALL.score.Merge.txt",
LDpred2 = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_LDpred2.chrALL.score.Merge.txt",
SBayesR = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_SBayesR.chrALL.score.Merge.txt",
PRS_CS = "C:/Users/yuting/Desktop/temp/hg38_risk_PRS_PRS_CS.chrALL.score.Merge.txt"
)
############################
## Read & merge all methods
############################
df_all <- bind_rows(
lapply(names(files), function(m) {
df <- read.table(files[m], header = TRUE)
df$Method <- m
df
})
)
############################
## Z-score standardization
## (done separately per method)
############################
df_all <- df_all %>%
group_by(Method) %>%
mutate(Score_z = as.numeric(scale(PRS_SUM))) %>%
ungroup()
############################
## Prepare legend labels
############################
n_df <- df_all %>%
distinct(V1, Group) %>%
count(Group)
labels <- setNames(
paste0(n_df$Group, " (n = ", n_df$n, ")"),
n_df$Group
)
############################
## Force Paper on top & custom label
############################
# 先設定 Method 順序
df_all$Method <- factor(df_all$Method,
levels = c("Paper", "MegaPRS", "Lassosum2", "LDpred2", "SBayesR", "PRS_CS"))
# 自訂 facet 標籤
method_labels <- c(
"Paper" = "Ho, Weang-Kee, et al. Nature communications 11.1 (2020): 3833",
"MegaPRS" = "MegaPRS (TPMI PheWeb 174.1)",
"Lassosum2" = "Lassosum2 (TPMI PheWeb 174.1)",
"LDpred2" = "LDpred2 (TPMI PheWeb 174.1)",
"SBayesR" = "SBayesR (TPMI PheWeb 174.1)",
"PRS_CS" = "PRS-CS (TPMI PheWeb 174.1)"
)
############################
## Plot: all methods together (Paper on top with custom title)
############################
ggplot(df_all, aes(x = Score_z, fill = Group)) +
geom_density(alpha = 0.5, color = NA) +
facet_wrap(~ Method, ncol = 1, labeller = labeller(Method = method_labels)) + # 用自訂標籤
theme_classic() +
labs(
x = "Standardized PRS (Z-score)",
y = "Density",
title = "Standardized PRS Distribution",
fill = "Group"
) +
scale_fill_manual(
values = c("Case" = "#E64B35", "Control" = "#4DBBD5"),
labels = labels
) +
theme(
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
strip.text = element_text(size = 14, face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 13),
legend.background = element_blank(),
legend.key = element_blank(),
panel.spacing = unit(1, "lines") # panel 間距，直排時很好用
)
