Score

# # 讀取資料
# prs_wt = read.delim("C:/Users/yuting/Desktop/pdc4pm_og/TPMI_GWAS_Stats/PRS_Weight/585.32_prs_weights.tsv", stringsAsFactors = F)
# 
# all156 = read.delim("C:/Users/yuting/Desktop/pdc4pm_og/TPMI_GWAS_Stats/GWAS_Stats/585.32_all.156.txt", stringsAsFactors = F)
# 
# 
# # ====================
# #  GWAS score
# # ====================
# 
# # OR → log(OR)
# all156$ES_GWAS = log(all156$OR)
# 
# all156$ID = paste0(all156$X.CHROM, ":",all156$POS)
# 
# # 轉成risk allele與risk effect size
# all156$EA_GWAS = all156$A1
# 
# idx_es = which(!is.na(all156$ES_GWAS) & all156$ES_GWAS < 0)
# 
# if (length(idx_es) > 0) {
#   all156$ES_GWAS[idx_es] = -all156$ES_GWAS[idx_es]
#   all156$EA_GWAS[idx_es] = all156$OMITTED[idx_es]
# }
# 
# 
# all156$ID_ALT_REF = paste0(
#   all156$X.CHROM, ":",
#   all156$POS, ":",
#   all156$A1, ":",
#   all156$OMITTED
# )
# 
# all156$ID_REF_ALT = paste0(
#   all156$X.CHROM, ":",
#   all156$POS, ":",
#   all156$OMITTED, ":",
#   all156$A1
# )
# 
# 
# # 匯出GWAS score檔
# write.table(all156[, c("ID", "EA_GWAS", "ES_GWAS", "ID_ALT_REF", "ID_REF_ALT")], file = "score_ES_GWAS.txt", sep = "\t", quote = F, row.names = F)
# 
# # ====================
# #  PRS methods
# # ====================
# prs_merged = merge(
#   prs_wt,
#   all156[, c("rsID", "X.CHROM", "POS", "A1", "OMITTED", "P")],
#   by.x = "MarkerID",
#   by.y = "rsID",
#   all.x = T
# )
# 
# 
# # Marker為chr1_189313672_G_A格式者，拆解Chr_Pos_OMITTED_Alt
# library(tidyr)
# library(dplyr)
# 
# idx = !grepl("^rs", prs_merged$MarkerID, ignore.case = TRUE)
# 
# prs_merged[idx, c("X.CHROM", "POS", "OMITTED", "A1")] =
#   separate(
#     data.frame(MarkerID = prs_merged$MarkerID[idx]),
#     col = MarkerID,
#     into = c("X.CHROM", "POS", "OMITTED", "A1"),
#     sep = "_"
#   )
# prs_merged$X.CHROM = sub("^chr", "", prs_merged$X.CHROM)
# 
# prs_merged$ID = paste0(prs_merged$X.CHROM, ":", prs_merged$POS)
# 
# prs_merged$ID_ALT_REF = 
#   paste0(
#     prs_merged$X.CHROM, ":",
#     prs_merged$POS, ":",
#     prs_merged$A1, ":",
#     prs_merged$OMITTED
#   )
# 
# 
# prs_merged$ID_REF_ALT = paste0(
#     prs_merged$X.CHROM, ":",
#     prs_merged$POS, ":",
#     prs_merged$OMITTED, ":",
#     prs_merged$A1
#   )
# 
# 
# # 轉成risk allele與risk effect size
# # 須注意每個TPMI PheWeb的PRS模型，做哪些PRS方法。
# prs_methods = c("LDpred2", "Lassosum2", "PRS.CS", "SBayesR", "MegaPRS")
# 
# for (m in prs_methods) {
#   
#   ea_col = paste0("EA_", m)
#   beta = prs_merged[[m]]
#   
#   prs_merged[[ea_col]] = prs_merged$Effect_Allele
#   
#   # beta < 0 的變異
#   idx = which(
#     !is.na(beta) &
#       beta < 0
#   )
#   
#   if (length(idx) > 0) {
#     
#     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
#     prs_merged[[m]][idx] = -prs_merged[[m]][idx]
#   }
# }
# 
# # 數字不顯示科學記號
# options(scipen = 999)
# 
# # 匯出各PRS score檔
# for (m in prs_methods) {
#   
#   out = prs_merged[, c(
#     "MarkerID",
#     paste0("EA_", m),
#     m,
#     "ID_ALT_REF",
#     "ID_REF_ALT"
#   )]
#   
#   write.table(
#     out,
#     file = paste0("score_", m, ".txt"),
#     sep = "\t",
#     quote = FALSE,
#     row.names = FALSE
#   )
# }