PRS 概述
PRS 计算的步骤
获取GWAS数据：首先，从GWAS研究中获取SNP与目标性状的效应大小和p值。 数据的质量控制（QC）：确保基础GWAS数据和目标数据的质量，例如基因型的缺失率和杂合度。 LD剪枝：通过剪枝减少SNP之间的连锁不平衡（LD），以确保PRS计算的SNP是独立的。 PRS计算：通过将效应大小（如beta值）与个体的基因型相乘并求和，计算每个个体的多基因风险评分。 关联分析和结果解释：通过回归分析或其他统计方法，评估PRS与目标性状的关联性，并解释结果。
PRS 计算中的关键点
SNP选择：应选择与目标性状最相关的SNP，这通常基于GWAS中的p值和效应大小。 效应修正：使用效应修正方法（如贝叶斯方法或LASSO）来减少随机误差。 LD问题：需要通过剪枝来处理SNP之间的LD问题，确保每个SNP独立贡献效应。
PRS 计算的挑战
样本偏差：基础数据和目标数据之间的种群差异可能会影响PRS的计算结果，因此应尽量匹配相似种群。 遗传与环境的交互作用：PRS不仅反映遗传因素，还可能受到环境因素的影响，因此在分析时需要考虑这些交互效应。 数据质量：确保GWAS数据和目标样本数据的高质量是PRS计算的关键。
# 加载所需的R包library(bigsnpr) # 处理大型SNP数据library(data.table) # 快速处理大规模数据表格library(ggplot2) # 可视化工具library(dplyr) # 数据操作library(MASS) # 执行回归分析# 第一步：加载GWAS汇总统计数据gwas_data <- fread("gwas_summary_stats.csv") # 从文件加载GWAS数据# 第二步：加载参考基因组数据ref_data <- snp_attach("data.rds") # 加载参考基因组# 第三步：数据质量控制# 过滤掉低频SNP和高p值SNP，确保数据质量qc_gwas <- gwas_data %>%filter(MAF > 0.01 & p < 1e-5) # 过滤小等位基因频率<1%和p值>1e-5的SNP# 第四步：LD剪枝以确保SNP独立性clumped_snp <- snp_clumping(G = ref_data$genotypes, # 基因型矩阵infos.chr = ref_data$map$chromosome, # 染色体信息infos.pos = ref_data$map$physical.pos, # SNP的物理位置thr.r2 = 0.1, # LD相关性阈值size = 500 # LD窗口大小设置为500kb)# 第五步：计算多基因风险评分（PRS）prs <- snp_PRS(G = ref_data$genotypes, # 使用参考数据的基因型betas.keep = qc_gwas$beta, # GWAS数据中的效应大小ind.keep = clumped_snp, # 经过剪枝保留的SNPlpS.keep = -log10(qc_gwas$p) # 使用p值进行SNP选择)# 第六步：加载目标样本数据并将PRS分数添加进去target_data <- fread("target_data.csv") # 加载目标样本数据target_data$PRS <- prs # 将PRS分数加入目标数据# 第七步：回归分析，评估PRS与目标性状的关联性# 对连续性状进行线性回归model <- lm(phenotype ~ PRS + age + sex + PC1 + PC2, data = target_data)summary(model) # 显示回归分析结果# 第八步：评估PRS对二分类结果的预测能力（例如疾病状态）# 使用逻辑回归进行二分类性状的分析logit_model <- glm(disease_status ~ PRS + age + sex + PC1 + PC2,family = binomial, data = target_data)summary(logit_model) # 显示逻辑回归分析结果# 第九步：PRS分布的可视化ggplot(target_data, aes(x = PRS)) +geom_histogram(binwidth = 0.05, fill = "blue", alpha = 0.7) +labs(title = "PRS Distribution", x = "Polygenic Risk Score", y = "Frequency") +theme_minimal()# 第十步：PRS与性状的关系可视化ggplot(target_data, aes(x = PRS, y = phenotype)) +geom_point(alpha = 0.5) +geom_smooth(method = "lm", se = FALSE, color = "red") +labs(title = "PRS vs Phenotype", x = "Polygenic Risk Score", y = "Phenotype") +theme_minimal()# 第十一步：进行模型预测和评估# 对模型的预测进行评估，如伪R方(pseudo-R2)用于二分类回归模型null_model <- glm(disease_status ~ age + sex + PC1 + PC2, family = binomial, data = target_data)pseudo_r2 <- 1 - (logLik(logit_model) / logLik(null_model)) # 计算伪R方print(paste("Pseudo-R2: ", pseudo_r2)) # 打印伪R方# 第十二步：绘制ROC曲线，评估PRS对二分类结果的预测性能library(pROC)roc_curve <- roc(target_data$disease_status, predict(logit_model, type = "response"))plot(roc_curve, main = "ROC Curve for PRS Prediction")
