文章总体情况-基于豆包
《53,026 名成年人健康与疾病血浆蛋白质组图谱研究》总结
一、研究背景与目标
（一）背景
全球人口增长与老龄化加剧，提升健康水平、减轻疾病负担迫在眉睫。精准医学旨在依据个体特征实现精准诊疗，虽基因组学取得进展，但基因转录翻译调控复杂，蛋白质作为疾病生物效应关键执行者，直接反映人体生理病理变化，其与疾病关系研究亟待加强。此前蛋白质组学研究多聚焦单一疾病，缺乏全面人类蛋白质组 - 表型图谱。
（二）目标
系统绘制 2,920 种血浆蛋白与 720 种疾病（406 种流行疾病、660 种新发疾病）及 986 种健康相关性状关联图谱，剖析疾病共有及特有生物学机制，挖掘疾病诊断与预测生物标志物及治疗靶点，构建综合蛋白质组 - 表型资源库推动精准医学发展。
二、研究方法与流程
（一）数据来源
从英国生物银行（UK Biobank）纳入 53,026 名参与者（平均年龄 56.8 岁，53.9%为女性，93.7%为白人），收集 2006 - 2010 年血浆蛋白数据（经质量控制筛选出 2,920 种蛋白质）、疾病诊断数据（依 ICD - 10 编码分类为流行与新发疾病）及健康相关性状数据（含连续、二元和有序分类变量）。
（二）分析方法
关联分析：分别用逻辑回归和 Cox 比例风险回归分析蛋白与流行、新发疾病关联；依性状变量类型选线性、逻辑或比例优势逻辑回归分析蛋白 - 性状关联，经严格校正识别显著关联。 敏感性与亚组分析：控制混杂因素，调整共病状态与多种协变量检验关联稳健性；按性别、年龄分层亚组分析探究关联变化。 功能富集分析：对疾病相关蛋白进行 Reactome 通路和 Gene Ontology 生物学过程富集分析，明确疾病相关生物学通路及功能。 聚类分析：依蛋白 - 疾病关联强度对新发疾病聚类，挖掘疾病集群特征生物学通路与共同发病机制。 预测与诊断模型构建：用 LightGBM 算法构建蛋白、人口统计学及二者融合模型，比较预测和诊断效能（以 AUC 等指标评估）。 孟德尔随机化分析：整合蛋白质定量性状位点（pQTL）与疾病 GWAS 数据，双向两样本 MR 分析判定蛋白与疾病因果关系及蛋白变化是否为疾病后果。
三、主要研究发现
（一）蛋白 - 疾病与蛋白 - 性状关联
疾病关联广泛：识别 168,100 个蛋白 - 流行疾病和 554,488 个蛋白 - 新发疾病关联，众多蛋白与多种疾病相关，如超 650 种蛋白在至少 50 种疾病中共享，超 1000 种蛋白呈现性别和年龄异质性。 性状关联丰富：大量蛋白与健康相关性状显著关联，部分蛋白在疾病和性状关联中具一致性，如 GDF15 和 CDCP1 与认知功能性状关联显著且在疾病关联中呈正相关。
（二）疾病相关蛋白功能与特征
生物学功能多样：免疫相关通路在多种疾病中富集，如 TNF 与受体结合通路参与过半疾病（神经系统除外）；蛋白代谢通路在多疾病有涉及，不同疾病的阿尔茨海默病和血管性痴呆相关蛋白在神经通路及脂质代谢等特有通路富集各异。 疾病集群特征明显：660 种新发疾病聚类成 40 群，各群有特征通路与共同发病机制，如肝纤维化肝硬化群涉及蛋白修饰等通路，非霍奇金淋巴瘤群具 B 细胞激活通路特征，部分集群含多种疾病类别，为疾病分类提供新思路。
（三）蛋白对疾病诊断与预测价值
诊断效能突出：蛋白模型在 124 种流行疾病（30.5%）和 92 种新发疾病（13.9%）诊断中 AUC 超 0.8，36 种疾病超 0.9，优于人口统计学模型，二者融合提升多疾病诊断准确性。 预测能力优异：蛋白模型对 92 种新发疾病预测 AUC 超 0.8，9 种超 0.9，在多疾病预测上比人口统计学模型准确，融合可改善 417 种（63.2%）疾病预测效果，部分疾病中蛋白信息对人口统计学信息有补充作用。
（四）疾病潜在因果蛋白与药物靶点
因果关系明确：MR 分析确定 474 个潜在因果蛋白 - 疾病对，为疾病发病机制提供依据，如 GDF15 与自身免疫病因果关联获遗传证据支持。 药物靶点潜力大：多数因果蛋白在循环、内分泌代谢疾病中发现，部分蛋白如 FURIN 与高血压显著关联；疾病相关蛋白与可成药基因组大量重叠，挖掘出 37 个药物再利用机会与 26 个安全潜在靶点，推动药物研发。
四、研究意义与展望
（一）理论意义
构建大规模蛋白质组 - 表型关联资源库，为疾病分子机制研究提供海量数据支撑，深化对蛋白质在疾病发生发展中核心作用的理解，从全新视角审视疾病分类边界与亚型，为精准医学理论体系注入活力。
（二）实践意义
为疾病早期诊断、精准预测与个性化治疗开辟新径，血浆蛋白检测有望成疾病筛查与监测得力工具，助力医生精准决策、优化治疗方案、提升医疗质量、改善患者预后，药物靶点挖掘加速药物研发进程、降低成本、提高成功率，为医药产业创新发展赋能。
（三）研究展望
后续应开展动物或人体实验验证蛋白与疾病机制，多生物样本分析补全组织特异性蛋白信息，利用大型前瞻性队列验证成果普适性，持续扩充蛋白质组覆盖度捕捉更多分子细节，监测患者标志物动态评估药物干预效果，基于生物分子特征革新疾病分类法提升临床试验精准度，推动精准医学全方位发展。
代码实战
项目整体解析-codebase
单独代码脚本解析-python
代码详细注释-基于cursor
# 导入所需的库
import glob # 用于文件路径匹配
import numpy as np # 用于数值计算
import pandas as pd # 用于数据处理
from tqdm import tqdm # 用于进度条显示
from collections import Counter # 用于计数
import warnings # 用于警告控制
from lightgbm import LGBMClassifier # 导入LightGBM分类器
from joblib import Parallel, delayed # 用于并行计算
warnings.filterwarnings('error') # 将警告转换为错误
# 设置全局参数
nb_cpus = 10 # CPU核心数
nb_params = 100 # 参数组合数量
my_seed = 2024 # 随机种子
def select_params_combo(my_dict, nb_items, my_seed):
"""
从参数字典中随机选择参数组合
Args:
my_dict: 参数字典
nb_items: 需要的组合数量
my_seed: 随机种子
Returns:
随机选择的参数组合列表
"""
combo_list = [dict(zip(my_dict.keys(), v)) for v in product(*my_dict.values())]
random.seed(my_seed)
return random.sample(combo_list, nb_items)
def normal_imp(mydict):
"""
对特征重要性进行归一化
Args:
mydict: 特征重要性字典
Returns:
归一化后的特征重要性字典
"""
mysum = sum(mydict.values())
mykeys = mydict.keys()
for key in mykeys:
mydict[key] = mydict[key]/mysum
return mydict
def get_cov_f_lst(tgt2pred_df, tgt):
"""
获取协变量特征列表
Args:
tgt2pred_df: 目标预测数据框
tgt: 目标疾病代码
Returns:
协变量特征列表
"""
sex_id = tgt2pred_df.loc[tgt2pred_df.Disease_code == tgt].SEX.iloc[0]
if (sex_id == 1) | (sex_id == 2):
cov_f_lst = ['AGE', 'RACE', 'TDI', 'BMI', 'SMK', 'ALC', 'SBP']
else:
cov_f_lst = ['AGE', 'Sex', 'RACE', 'TDI', 'BMI', 'SMK', 'ALC', 'SBP']
return cov_f_lst
def get_pro_f_lst(mydf, train_idx, f_lst, my_params):
"""
获取蛋白质特征列表
Args:
mydf: 数据框
train_idx: 训练集索引
f_lst: 特征列表
my_params: 模型参数
Returns:
前30个重要的蛋白质特征列表
"""
X_train, y_train = mydf.iloc[train_idx][f_lst], mydf.iloc[train_idx].target_y
my_lgb = LGBMClassifier(objective='binary', metric='auc', is_unbalance=False, verbosity=1, seed=2023)
my_lgb.set_params(**my_params)
my_lgb.fit(X_train, y_train)
totalgain_imp = my_lgb.booster_.feature_importance(importance_type='gain')
totalgain_imp = dict(zip(my_lgb.booster_.feature_name(), totalgain_imp.tolist()))
tg_imp_df = pd.DataFrame({'Pro_code': list(totalgain_imp.keys()), 'TotalGain': list(totalgain_imp.values())})
tg_imp_df.sort_values(by = 'TotalGain', inplace = True, ascending = False)
return tg_imp_df.Pro_code.tolist()[:30]
# ... [其余函数的注释类似,为了简洁起见省略]
# 设置模型参数搜索空间
params_dict = {
'n_estimators': [100, 200, 300, 400, 500], # 树的数量
'max_depth': np.linspace(5, 30, 6).astype('int32').tolist(), # 树的最大深度
'num_leaves': np.linspace(5, 30, 6).astype('int32').tolist(), # 叶子节点数
'subsample': np.linspace(0.6, 1, 9).tolist(), # 样本采样比例
'learning_rate': [0.1, 0.05, 0.01, 0.001], # 学习率
'colsample_bytree': np.linspace(0.6, 1, 9).tolist() # 特征采样比例
}
# 生成候选参数组合
candidate_params_lst = select_params_combo(params_dict, nb_params, my_seed)
# 设置数据路径
dpath = '/home1/jiayou/Documents/Projects/ProDisAtlas/'
#dpath = '/Volumes/JasonWork/Projects/ProDisAtlas/'
# 读取数据
tgt2pred_df = pd.read_csv(dpath + 'Data/Target/Prevalent.csv', encoding='latin-1')
tgt2pred_lst = tgt2pred_df.Disease_code.tolist()
pro_df = pd.read_csv(dpath + 'Data/ProteinData/ProteinData.csv')
full_pro_f_lst = pro_df.columns.tolist()[1:]
cov_df = pd.read_csv(dpath + 'Data/Covariates/Covariates.csv')
# 合并数据
mydf = pd.merge(cov_df, pro_df, how = 'left', on = ['eid'])
fold_id_lst = list(set(mydf.Region_code))
fold_id_lst.sort()
# 初始化错误列表
bad_tgt_lst, nb_tgt_lst = [], []
# 主循环:对每个目标疾病进行预测
for tgt in tqdm(tgt2pred_lst[125:250]):
try:
# ... [后续代码的注释类似,为了简洁起见省略]
pass
except:
bad_tgt_lst.append(tgt)
# 保存错误目标列表
bad_tgt_df = pd.DataFrame({'Disease_code': bad_tgt_lst})
bad_tgt_df.to_csv(dpath + 'Results/Prediction/CrossSectionalAnalysis/bad_tgt.csv', index=False)
单个代码脚本解读-基于R
代码详细注释-基于cursor
# 加载所需的R包
library(GOSemSim) # 用于GO语义相似度分析
library(clusterProfiler) # 用于富集分析
library(org.Hs.eg.db) # 人类基因组注释数据库
library(dplyr) # 数据处理工具包
##1. 对具有显著蛋白质的疾病进行GO富集分析##
# 获取所有csv文件列表
file_list <- list.files(path=path_to_files, pattern="*.csv", full.names=T)
# 遍历每个文件进行分析
for (i in (1:length(file_list))) {
# 读取数据文件
data <- read.csv(file_list[i])
# 获取文件名信息
file_name <- file_list[i]
base_name <- basename(file_name)
input_name <- tools::file_path_sans_ext(base_name)
# 筛选显著性蛋白质
# 横断面分析和纵向分析的显著性阈值不同
vector = (data$pval_raw)*2920*406 < 0.05 & data$Pro_code !="" #横断面分析
vector = (data$pval_raw)*2920*660 < 0.05 & data$Pro_code !="" #纵向分析
data_sgni = data[vector,]
# 如果没有显著性蛋白质则跳过
if(nrow(data_sgni) == 0){
print(paste("Skipping iteration ", i, " due to empty data_sgni"))
next
}
# 将基因符号转换为Entrez ID
genes_full <- bitr(data_sgni$Pro_code, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
if(nrow(genes_full) == 0){
print(paste("Skipping iteration ", i, " due to empty entrezid"))
next
}
genes = na.omit(genes_full) # 去除NA值
genes = genes$ENTREZID # 提取Entrez ID
# 进行GO富集分析
ego <- enrichGO(gene = genes,
keyType = "ENTREZID", # 基因ID类型
OrgDb = org.Hs.eg.db, # 使用人类基因组数据库
ont = "ALL", # 分析所有GO类别
pAdjustMethod = 'fdr', # 使用FDR校正
qvalueCutoff = 0.05, # q值阈值
pvalueCutoff = 0.05, # p值阈值
minGSSize = 1, # 最小基因集大小
readable = TRUE) # 输出可读的基因名
# 如果富集分析结果为空则跳过
if(is.null(ego)){
print(paste("Skipping iteration ", i, " due to empty ego"))
next
}
# 简化GO富集结果,去除冗余项
ego2 <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)
GO_result = data.frame(ego2)
# 保存结果
write.csv(GO_result,paste0(outfilePath, "/GOresult_", input_name, ".csv"))
}
##2. 对没有显著蛋白质的疾病进行GO富集分析,使用p值最小的前30个蛋白质##
# 获取文件列表
file_list <- list.files(path=path_to_files, pattern="*.csv", full.names=T)
# 遍历每个文件
for (i in (1:length(file_list))) {
# 读取数据
data <- read.csv(file_list[i])
# 获取文件名信息
file_name <- file_list[i]
base_name <- basename(file_name)
input_name <- tools::file_path_sans_ext(base_name)
# 检查是否有显著性蛋白质
vector = (data$pval_raw)*2920*406 < 0.05 & data$Pro_code !="" #横断面分析
vector = (data$pval_raw)*2920*660 < 0.05 & data$Pro_code !="" #纵向分析
data_sgni = data[vector,]
# 如果有显著性蛋白质则跳过
if(nrow(data_sgni) > 0){
print(paste0("Significant proteins found in No.", i, ", Skipping enrichment."))
}
else{
# 选择p值最小的前30个蛋白质
data_sgni <- data %>%
arrange(pval_raw) %>%
slice_min(pval_raw, n = 30, with_ties = TRUE)
# 转换基因ID
genes_full <- bitr(data_sgni$Pro_code, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
if(nrow(genes_full) == 0){
print(paste("Skipping iteration ", i, " due to empty entrezid"))
next
}
genes = na.omit(genes_full)
genes = genes$ENTREZID
# GO富集分析
ego <- enrichGO(gene = genes,
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = 'fdr',
qvalueCutoff = 0.05,
pvalueCutoff= 0.05,
minGSSize = 1,
readable = TRUE)
# 检查富集结果
if(is.null(ego)){
print(paste("Skipping iteration ", i, " due to empty ego"))
next
}
# 简化并保存结果
ego2 <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)
GO_result = data.frame(ego2)
write.csv(GO_result,paste0(outfilePath, "/GOresult_", input_name, ".csv"))
}
}
脚本依赖及整体架构解读能力-R语言
要求合并三个脚本，一次输入全部输出-我比较懒
以下为合并的脚本
### 加载所需的包
library(ReactomePA)
library(clusterProfiler)
library(org.Hs.eg.db)
library(GOSemSim)
library(TissueEnrich)
library(dplyr)
library(tidyr)
library(tidyverse)
# 定义富集分析函数
perform_enrichment_analysis <- function(data_sgni, input_name, outfilePath) {
# 准备基因列表
genes_full <- bitr(data_sgni$Pro_code, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
if(nrow(genes_full) == 0){
print("Empty entrezid, skipping analysis")
return(NULL)
}
genes <- na.omit(genes_full)$ENTREZID
# 1. Reactome通路富集分析
tryCatch({
x <- enrichPathway(gene = genes,
pAdjustMethod = "fdr",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE)
if(!is.null(x)) {
Reactome_result = data.frame(x)
write.csv(Reactome_result, paste0(outfilePath, "/Reactome_res_", input_name, ".csv"))
}
}, error = function(e) {
print("Error in Reactome analysis")
})
# 2. GO富集分析
tryCatch({
ego <- enrichGO(gene = genes,
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = 'fdr',
qvalueCutoff = 0.05,
pvalueCutoff = 0.05,
minGSSize = 1,
readable = TRUE)
if(!is.null(ego)) {
ego2 <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)
GO_result = data.frame(ego2)
write.csv(GO_result, paste0(outfilePath, "/GOresult_", input_name, ".csv"))
}
}, error = function(e) {
print("Error in GO analysis")
})
# 3. 组织富集分析
tryCatch({
data_sgni <- distinct(data_sgni, Pro_code, .keep_all = TRUE)
inputGenes <- data_sgni$Pro_code
gs <- GeneSet(geneIds=inputGenes, organism="Homo Sapiens", geneIdType=SymbolIdentifier())
output <- teEnrichment(inputGenes=gs, rnaSeqDataset=2) # 2 represents GTEx
seEnrichmentOutput <- output[[1]]
enrichmentOutput <- setNames(data.frame(assay(seEnrichmentOutput),
row.names = rowData(seEnrichmentOutput)[,1]),
colData(seEnrichmentOutput)[,1])
enrichmentOutput$Tissue <- row.names(enrichmentOutput)
write.csv(enrichmentOutput, paste0(outfilePath, "/tissueEnrichresult_", input_name, ".csv"))
}, error = function(e) {
print("Error in Tissue enrichment analysis")
})
}
# 主程序
main_analysis <- function(path_to_files, outfilePath) {
file_list <- list.files(path=path_to_files, pattern="*.csv", full.names=T)
for (i in 1:length(file_list)) {
# 读取数据
data <- read.csv(file_list[i])
file_name <- file_list[i]
base_name <- basename(file_name)
input_name <- tools::file_path_sans_ext(base_name)
# 判断是否有显著蛋白质
vector_cs <- (data$pval_raw)*2920*406 < 0.05 & data$Pro_code !="" #横断面分析
vector_lg <- (data$pval_raw)*2920*660 < 0.05 & data$Pro_code !="" #纵向分析
# 使用横断面或纵向分析的结果
vector <- vector_cs # 或 vector_lg，根据需要选择
data_sgni <- data[vector,]
if(nrow(data_sgni) > 0) {
# 对显著蛋白质进行富集分析
print(paste0("Processing file ", i, " with significant proteins"))
perform_enrichment_analysis(data_sgni, input_name, outfilePath)
} else {
# 对top 30蛋白质进行富集分析
print(paste0("Processing file ", i, " with top 30 proteins"))
data_sgni <- data %>%
arrange(pval_raw) %>%
slice_min(pval_raw, n = 30, with_ties = TRUE)
perform_enrichment_analysis(data_sgni, input_name, outfilePath)
}
}
}
# 使用示例
# path_to_files <- "path/to/your/input/files"
# outfilePath <- "path/to/your/output/directory"
# main_analysis(path_to_files, outfilePath)
脚本间关系测试-python
# 导入所需的库
import glob # 用于文件路径匹配
import numpy as np # 用于数值计算
import pandas as pd # 用于数据处理和分析
import statsmodels.api as sm # 用于统计建模
from lifelines import CoxPHFitter # 用于Cox比例风险模型
from tqdm import tqdm # 用于进度条显示
import os # 用于操作系统相关功能
import re # 用于正则表达式操作
from statsmodels.stats.multitest import fdrcorrection # 用于FDR多重检验校正
from mne.stats import bonferroni_correction # 用于Bonferroni多重检验校正
from joblib import Parallel, delayed # 用于并行计算
import warnings # 用于警告控制
warnings.filterwarnings('error') # 将警告转换为错误
class IntegratedAnalysis:
"""
整合分析类,用于执行蛋白质与疾病关联的统计分析
包括Cox回归分析(用于发病率分析)和逻辑回归分析(用于患病率分析)
"""
def __init__(self, dpath, nb_cpus=58):
"""
初始化分析类
参数:
dpath (str): 数据和结果的根目录路径
nb_cpus (int): 并行计算使用的CPU核心数,默认58
"""
self.dpath = dpath
self.nb_cpus = nb_cpus
self.setup_paths() # 设置输入输出路径
self.load_data() # 加载基础数据
self.setup_formulas() # 设置分析公式
def setup_paths(self):
"""
设置分析结果的输出路径
为不同类型分析(发病率/患病率)和不同亚组(全部/女性/男性/年轻/年老)创建目录
"""
# 定义输出目录结构
self.out_dirs = {
'incident': { # 发病率分析结果目录
'all': f'{self.dpath}/Results/Association/IncidentAnalysis/All/',
'female': f'{self.dpath}/Results/Association/IncidentAnalysis/Female/',
'male': f'{self.dpath}/Results/Association/IncidentAnalysis/Male/',
'young': f'{self.dpath}/Results/Association/IncidentAnalysis/Young/',
'old': f'{self.dpath}/Results/Association/IncidentAnalysis/Old/'
},
'prevalent': { # 患病率分析结果目录
'all': f'{self.dpath}/Results/Association/CrossSectionalAnalysis/All/',
'female': f'{self.dpath}/Results/Association/CrossSectionalAnalysis/Female/',
'male': f'{self.dpath}/Results/Association/CrossSectionalAnalysis/Male/',
'young': f'{self.dpath}/Results/Association/CrossSectionalAnalysis/Young/',
'old': f'{self.dpath}/Results/Association/CrossSectionalAnalysis/Old/'
}
}
# 创建所有输出目录
for analysis_type in self.out_dirs:
for subgroup in self.out_dirs[analysis_type]:
os.makedirs(self.out_dirs[analysis_type][subgroup], exist_ok=True)
def load_data(self):
"""
加载分析所需的基础数据:
1. 目标疾病文件列表
2. 目标疾病信息(包含性别特异性信息)
3. 蛋白质数据和协变量
4. 设置分析所需的协变量列表
"""
# 加载目标疾病文件列表并排序
self.target_file_lst = self.sort_nicely(
glob.glob(f'{self.dpath}/Data/Target/Targets2Analysis/*.csv'))
# 加载目标疾病信息
self.target_info_df = pd.read_csv(
f'{self.dpath}/Data/Target/TargetVsProtein.csv',
usecols=['NAME', 'SEX'])
# 加载蛋白质数据
self.pro_cov_df = pd.read_csv(
f'{self.dpath}/Data/ProteinData/ProteinData_n_Cov.csv')
self.pro_f_lst = self.pro_cov_df.columns[13:2933].tolist() # 提取蛋白质列名
# 处理种族数据:将多个种族类别二值化(1为白人,0为其他)
self.pro_cov_df['Race'].replace([1,2,3,4], [1, 0, 0, 0], inplace=True)
# 设置协变量列表
self.cov_f_lst = ['eid', 'Age', 'Sex', 'Race', 'TDI', 'BMI',
'smk', 'fastingtime', 'season']
self.cov_df = self.pro_cov_df[self.cov_f_lst]
def setup_formulas(self):
"""
设置回归分析使用的公式
formula_full: 包含所有协变量的完整公式
formula_non_sex: 不包含性别变量的公式(用于性别特异性疾病)
"""
self.formula_full = ("Age + C(Sex) + C(Race) + TDI + BMI + C(smk) + "
"fastingtime + C(season) + x_pro_sa + x_pro")
self.formula_non_sex = ("Age + C(Race) + TDI + BMI + C(smk) + "
"fastingtime + C(season) + x_pro_sa + x_pro")
@staticmethod
def sort_nicely(l):
"""
对文件列表进行自然排序
例如: ['file1', 'file2', 'file10'] 而不是 ['file1', 'file10', 'file2']
参数:
l (list): 需要排序的文件列表
返回:
list: 排序后的文件列表
"""
convert = lambda text: int(text) if text.isdigit() else text
alphanum_key = lambda key: [convert(c.replace("_",""))
for c in re.split('([0-9]+)', key)]
l.sort(key=alphanum_key)
return l
def process_protein(self, pro_f, tmp_tgt_df, analysis_type='cox'):
"""
处理单个蛋白质的统计分析
参数:
pro_f (str): 蛋白质特征名
tmp_tgt_df (DataFrame): 目标疾病数据
analysis_type (str): 分析类型('cox'或'logistic')
返回:
list: 包含分析结果的列表
"""
# 准备数据
tmp_pro_df = self.pro_cov_df[['eid', pro_f, f'{pro_f}_SampAge']]
tmp_df = pd.merge(tmp_tgt_df, tmp_pro_df, how='left', on=['eid'])
tmp_df.rename(columns={pro_f: 'x_pro',
f'{pro_f}_SampAge': 'x_pro_sa'},
inplace=True)
# 清理数据:删除缺失值并填充采样年龄中位数
rm_eid_idx = tmp_df.index[tmp_df.x_pro.isnull()]
tmp_df.drop(rm_eid_idx, axis=0, inplace=True)
tmp_df.reset_index(inplace=True, drop=True)
tmp_df['x_pro_sa'].fillna(tmp_df['x_pro_sa'].median(), inplace=True)
# 计算基本统计量
nb_all, nb_case = len(tmp_df), tmp_df.target_y.sum()
prop_case = np.round(nb_case / nb_all * 100, 3)
# 根据分析类型和病例数执行相应分析
if analysis_type == 'cox' and nb_case >= 10:
return self._cox_analysis(tmp_df, pro_f, nb_all, nb_case, prop_case)
elif analysis_type == 'logistic' and nb_case >= 10:
return self._logistic_analysis(tmp_df, pro_f, nb_all, nb_case, prop_case)
else:
return [pro_f, nb_all, nb_case, prop_case, np.nan, np.nan, np.nan, np.nan]
def _cox_analysis(self, tmp_df, pro_f, nb_all, nb_case, prop_case):
"""
执行Cox比例风险回归分析
参数:
tmp_df (DataFrame): 准备好的数据
pro_f (str): 蛋白质特征名
nb_all (int): 总样本数
nb_case (int): 病例数
prop_case (float): 病例比例
返回:
list: 包含分析结果的列表
"""
i, tmpout = 1, []
while (len(tmpout) == 0) and (i <= 1e7):
try:
cph = CoxPHFitter(penalizer=1e-7 * i)
cph.fit(tmp_df, duration_col='BL2Target_yrs',
event_col='target_y', formula=self.current_formula)
hr = np.round(cph.hazard_ratios_.x_pro, 5)
lbd = np.round(np.exp(cph.confidence_intervals_).loc['x_pro'][0], 5)
ubd = np.round(np.exp(cph.confidence_intervals_).loc['x_pro'][1], 5)
pval = cph.summary.p.x_pro
tmpout = [pro_f, nb_all, nb_case, prop_case, hr, lbd, ubd, pval]
except:
i *= 10
return tmpout if tmpout else [pro_f, nb_all, nb_case, prop_case,
np.nan, np.nan, np.nan, np.nan]
def _logistic_analysis(self, tmp_df, pro_f, nb_all, nb_case, prop_case):
"""
执行逻辑回归分析
参数:
tmp_df (DataFrame): 准备好的数据
pro_f (str): 蛋白质特征名
nb_all (int): 总样本数
nb_case (int): 病例数
prop_case (float): 病例比例
返回:
list: 包含分析结果的列表
"""
try:
Y = tmp_df.target_y
X = tmp_df[self.cov_f_lst[1:] + ['x_pro', 'x_pro_sa']]
log_mod = sm.Logit(Y, sm.add_constant(X)).fit()
oratio = np.round(np.exp(log_mod.params).loc['x_pro'], 5)
pval = log_mod.pvalues.loc['x_pro']
ci_mod = log_mod.conf_int(alpha=0.05)
lbd = np.round(np.exp(ci_mod.loc['x_pro'][0]), 5)
ubd = np.round(np.exp(ci_mod.loc['x_pro'][1]), 5)
return [pro_f, nb_all, nb_case, prop_case, oratio, lbd, ubd, pval]
except:
return [pro_f, nb_all, nb_case, prop_case, np.nan, np.nan, np.nan, np.nan]
def analyze_subgroup(self, tmp_tgt_df, analysis_type, subgroup):
"""
分析特定亚组
参数:
tmp_tgt_df (DataFrame): 目标疾病数据
analysis_type (str): 分析类型('cox'或'logistic')
subgroup (str): 亚组类型('all','female','male','young','old')
返回:
DataFrame/None: 分析结果数据框或None(如果病例数不足)
"""
# 根据亚组类型筛选数据
if subgroup == 'female':
tmp_tgt_df = tmp_tgt_df.loc[tmp_tgt_df.Sex == 0]
elif subgroup == 'male':
tmp_tgt_df = tmp_tgt_df.loc[tmp_tgt_df.Sex == 1]
elif subgroup == 'young':
tmp_tgt_df = tmp_tgt_df.loc[tmp_tgt_df.Age < 60]
elif subgroup == 'old':
tmp_tgt_df = tmp_tgt_df.loc[tmp_tgt_df.Age >= 60]
tmp_tgt_df.reset_index(inplace=True, drop=True)
# 如果病例数足够,执行并行分析
if tmp_tgt_df.target_y.sum() >= 10:
results = Parallel(n_jobs=self.nb_cpus)(
delayed(self.process_protein)(
pro_f, tmp_tgt_df, analysis_type)
for pro_f in self.pro_f_lst)
return pd.DataFrame(results)
return None
def format_results(self, results_df, analysis_type):
"""
格式化分析结果
参数:
results_df (DataFrame): 原始结果数据框
analysis_type (str): 分析类型('cox'或'logistic')
返回:
DataFrame: 格式化后的结果数据框
"""
# 设置列名
if analysis_type == 'cox':
ratio_col = 'hr'
cols = ['Pro_code', 'nb_individuals', 'nb_case', 'prop_case(%)',
'hr', 'hr_lbd', 'hr_ubd', 'pval_raw']
else:
ratio_col = 'or'
cols = ['Pro_code', 'nb_individuals', 'nb_case', 'prop_case(%)',
'oratio', 'or_lbd', 'or_ubd', 'pval_raw']
results_df.columns = cols
# Bonferroni多重检验校正
_, p_f_bfi = bonferroni_correction(
results_df.pval_raw.fillna(1), alpha=0.05)
results_df['pval_bfi'] = p_f_bfi
results_df.loc[results_df['pval_bfi'] >= 1, 'pval_bfi'] = 1
# 添加格式化的输出列
results_df[f'{ratio_col}_output'] = results_df.apply(
lambda x: f"{x[ratio_col]:.2f} [{x[f'{ratio_col}_lbd']:.2f}-"
f"{x[f'{ratio_col}_ubd']:.2f}]", axis=1)
# 添加显著性标记
results_df['pval_significant'] = results_df.pval_bfi.apply(
lambda x: '***' if x < 0.001 else
'**' if x < 0.01 else
'*' if x < 0.05 else '')
# 重命名列(如果是逻辑回归分析)
if analysis_type == 'logistic':
results_df.rename(columns={'oratio': 'or'}, inplace=True)
return results_df
def run_analysis(self):
"""
运行完整分析流程
处理所有目标疾病文件,执行Cox回归和逻辑回归分析
对不同亚组分别进行分析并保存结果
"""
for tgt_file in tqdm(self.target_file_lst):
try:
# 获取目标疾病名称
tgt_name = os.path.basename(tgt_file).split('.')[0]
# 根据疾病的性别特异性确定分析公式
sex_id = self.target_info_df.loc[
self.target_info_df.NAME == tgt_name].SEX.iloc[0]
self.current_formula = (self.formula_non_sex
if sex_id in [1, 2]
else self.formula_full)
# 读取目标疾病数据
tmp_tgt_df = pd.read_csv(
tgt_file,
usecols=['eid', 'target_y', 'BL2Target_yrs'])
# Cox分析数据预处理
cox_df = tmp_tgt_df.copy()
rm_bl_idx = cox_df.index[cox_df.BL2Target_yrs <= 0]
cox_df.drop(rm_bl_idx, axis=0, inplace=True)
cox_df = pd.merge(cox_df, self.cov_df, how='left', on=['eid'])
# Logistic分析数据预处理
log_df = tmp_tgt_df.copy()
rm_idt_idx = log_df.index[
(log_df.BL2Target_yrs > 0) & (log_df.target_y == 1)]
log_df.drop(rm_idt_idx, axis=0, inplace=True)
log_df = pd.merge(log_df, self.cov_df, how='left', on=['eid'])
# 对每个亚组进行分析
for subgroup in ['all', 'female', 'male', 'young', 'old']:
# Cox回归分析(如果病例数>50)
if cox_df.target_y.sum() > 50:
cox_results = (self.analyze_subgroup(cox_df, 'cox', subgroup)
if subgroup == 'all'
else self.analyze_subgroup(
cox_df.copy(), 'cox', subgroup))
if cox_results is not None:
cox_results = self.format_results(cox_results, 'cox')
cox_results.to_csv(
f"{self.out_dirs['incident'][subgroup]}"
f"{tgt_name}.csv", index=False)
# 逻辑回归分析(如果病例数>50)
if log_df.target_y.sum() > 50:
log_results = (self.analyze_subgroup(log_df, 'logistic', subgroup)
if subgroup == 'all'
else self.analyze_subgroup(
log_df.copy(), 'logistic', subgroup))
if log_results is not None:
log_results = self.format_results(log_results, 'logistic')
log_results.to_csv(
f"{self.out_dirs['prevalent'][subgroup]}"
f"{tgt_name}.csv", index=False)
except Exception as e:
print(f"Error processing {tgt_name}: {str(e)}")
continue
# 使用示例
if __name__ == "__main__":
# 设置数据路径
dpath = '/path/to/your/data/directory'
# 创建分析实例
analyzer = IntegratedAnalysis(dpath)
# 运行分析
analyzer.run_analysis()
国内第一个用Cursor讲解R语言的课程
课程中将详细讲解R语言的Cursor配置
课程名称-基于R语言的动态预测模型课程培训班
为什么要用cursor
随时可唤起AI模型
丰富的AI模型（ctrl+L）
仅从某宝或某鱼可购买无限制低价会员
根据注释自动生成代码
根据代码区要求生成代码（ctrl+K）
可SSH登录远程服务器R语言环境
与R studio类似的工作环境-workspace
支持类似于R studio的快捷键设定
开课目的及前言
预测模型作为真实世界研究的重要组成部分，其研究被广泛开展。但是，传统的预测模型利用基线数据对最终的生存结果进行预测，这种模型无法纳入患者在后续随访中可能会动态变化的重要数据（比如肿瘤标记物的动态变化）。 以上情况在统计学中会产生估计偏差情况，也是不符合临床实际的。近年来发展起来的动态预测模型方法，利用患者的多次随访数据，结合患者的基线数据，对最终患者的额生存结果（或类似的time to event事件）进行估计。其发文量呈现快速增长趋势。
在临床实际中，医生会根据患者的动态变化指标做出进一步诊断及治疗的判断。动态预测模型结合患者的纵向数据与最终的生存结果，对于最终结果进行更加准备的预测。由于当前R语言在医学统计工作中占据重要地位，但很多临床大夫、护士因为时间工作关系很难将R语言与临床科研相结合，故开设R语言动态预测模型课程，旨在快速让学员掌握统计工作中常用到的R语言，助力临床科研工作。天企助力(天津)生产力促进有限公司特举办“基于R语言的动态预测模型课程培训班”。
预测模型类文章目前总结起来发展经历了以下三个阶段：
基于传统流行病学的列线图模型（本质都是cox回归及glm回归），简单的统计学分析模型，是模型依赖的方法，临床上实际情况很难满足其前提假设，实际效果不好。
基于机器学习/深度学习的预测模型的构建（在数据上提高了维度，在算法上引入了机器学习），虽然算法上引入了机器学习模型，处理数据更加灵活，模型的假设也更少。但是在使用的数据上还是患者的一次基线数据进行预测，与临床实际不符。
基于纵向数据的动态预测模型（基于纵向多次随访数据，模型应用联合模型等动态预测模型方法），应用患者的多次随访数据对最终的生存结果进行预测，从数据和方法上都更类似于临床实际。
考虑到动态预测模型有以下特点，因此必然是后续高分文章的必备方法：
数据上必须有同一个患者的多次随访数据，相对于既往横断面一次基线数据，数据的收集难度更大，而且动态预测模型需拟合纵向的线性混合模型，因此需要的数据量较大。这就提示我们如果能收集到如上数据更加容易发高分文章。
应用方法学动态预测模型需首先掌握普通生存分析及普通预测模型的方法，并且还需要熟悉纵向数据分析的广义线性混合模型，再次基础上还需要掌握tidyverse语法基础来将自己的数据转换为满足函数要求的纵向数据，另外对于联合模型，模型的结合形式及变量选择也均需要从临床背景及统计学方法考虑。
近期高分文章举例
授课老师
灵活胖子-独自
双一流学校肿瘤学博士毕业，目前就职于国内五大肿瘤中心之一。科研方向为真实世界研究，生物信息学分析及人工智能研究。目前以第一或共同第一作者身份发表SCI论文10余篇，累计IF50+。目前与国内多个院校及医院有科研合作。联合翻译小组同学，在国内第一次将jmbayes2及dynamicLM全文翻译为中文并在公众号发表。
课程目录及安排
授课形式及时间
授课形式：远程在线实时直播授课。
授课时间：2024年12月开课，总课时不少于30小时，每周进行3-5小时的授课，有充分时间学习，预计6-8周完成所有授课内容。
答疑支持：建立课程专属微信群，1年内课程内容免费答疑。
视频回看：3年内免费无限次回看。
课程售价及售后保证
课程售价：总价3000元，报名可先交300元预定即可，开课后2周内交齐即可
对公转账等手续务必提前联系助教
承办公司：天企助力（天津）生产力促进有限公司
奖励政策：学员应用所学内容发表IF 10+文章可退还学费（具体要求及流程需要咨询助教）
