★ FEATURED ARTICLE 简介SmokingGeneExpression 是一份面向生物信息学与统计分析学习者的 R 语言实战项目资源聚焦用因果推断方法评估吸烟对健康肺组织中基因表达的短期因果效应。项目融合工具变量分析、倾向得分匹配、DAG 等因果推断思路并借助 causalInference、MatchIt 等 R 包控制混杂因素适合具备一定 R 基础、希望深入基因表达差异分析与因果建模的读者参考。资源包共 45 个文件约 3.1MB包含 19 个 png 结果图、13 个 txt 数据说明、8 个 xlsx 表格、2 个 csv、1 个 pdf 文档及 1 个 R 脚本覆盖样本 PCA、基因表达表、差异基因计数与人口统计辅助数据等模块便于按目录结构复现分析流程。目前已有 148 人学习下载。通过该资源可了解从 RNA-seq 预处理、差异表达分析到通路富集与结果验证的完整链路并掌握因果推断在真实健康数据中的落地方式为吸烟相关疾病机制研究提供可借鉴的分析框架。1. 用因果推断拆解吸烟与基因表达的关联一个 R 语言项目能跑出什么吸烟对基因表达的影响很多研究停留在相关性层面——吸烟者某基因表达高不吸烟者低但这不等于吸烟导致了变化。这个项目用因果推断方法专门评估吸烟对健康肺组织中基因表达的短期因果关系。换句话说它要回答的是如果一个人从不吸烟变成吸烟肺组织里的基因表达会发生什么变化项目基于 R 语言构建适合有转录组数据分析基础、想入门因果推断的开发者也适合做生信分析的同学拿来当模板改造成自己的暴露-结局研究。整个流程从数据预处理到因果效应估计再到结果可视化链路完整跑通一遍就能理解因果推断在组学数据里怎么落地。2. 项目架构与数据流从表达矩阵到因果效应估计2.1 目录结构与模块划分拿到源码包后先别急着跑脚本。我一般会花十分钟把目录翻一遍搞清楚数据从哪进、结果从哪出。这个项目的结构比较典型按分析阶段分目录SmokingGeneExpression/ ├── data/ │ ├── raw/ # 原始表达矩阵与样本元数据 │ └── processed/ # 清洗、标准化后的数据 ├── R/ │ ├── 01_preprocess.R # 数据读取、质控、标准化 │ ├── 02_estimate_effect.R # 因果效应估计主逻辑 │ ├── 03_visualize.R # 结果可视化 │ └── utils.R # 公共函数 ├── output/ │ ├── tables/ # 因果效应估计结果表 │ └── figures/ # 火山图、森林图等 ├── config.yml # 参数配置文件 └── run_all.R # 一键运行入口data/raw/放的是原始数据通常是基因表达矩阵行是基因列是样本和样本注释表包含吸烟状态、年龄、性别等协变量。R/下按分析步骤拆成独立脚本每个脚本可以单独跑也可以由run_all.R串起来。config.yml是参数入口比如显著性阈值、协变量列表、输出路径改参数不用翻代码。注意原始数据没有随包提供的话需要自己准备一份表达矩阵和对应的样本信息表格式要求后面会讲。2.2 因果推断的核心思路这个项目要解决的核心问题是吸烟状态暴露对基因表达结局的因果效应。直接做差异表达分析得到的是关联不是因果。项目采用的方法框架是有向无环图DAG指导的协变量调整配合逆概率加权IPW或边际结构模型来估计因果效应。具体来说流程分三步构建因果图根据先验知识确定哪些变量是混杂因素同时影响吸烟和基因表达比如年龄、性别、BMI。这些变量需要调整。估计倾向性得分用逻辑回归建模吸烟概率给定协变量。加权或分层估计效应用倾向性得分做逆概率加权平衡协变量分布后再估计吸烟对每个基因表达的效应。项目里02_estimate_effect.R就是干这件事的。它会对每个基因循环执行上述步骤最终输出每个基因的因果效应估计值、标准误和校正后的 p 值。2.3 数据预处理的关键步骤01_preprocess.R负责把原始数据整理成可分析的格式。核心操作包括# 读取表达矩阵和样本信息 expr_matrix - read.csv(data/raw/expression.csv, row.names 1) sample_info - read.csv(data/raw/sample_metadata.csv) # 确保样本顺序一致 common_samples - intersect(colnames(expr_matrix), sample_info$sample_id) expr_matrix - expr_matrix[, common_samples] sample_info - sample_info[match(common_samples, sample_info$sample_id), ] # 过滤低表达基因 keep_genes - rowSums(expr_matrix 1) 0.5 * ncol(expr_matrix) expr_matrix - expr_matrix[keep_genes, ] # log2 转换如果数据未转换 if (max(expr_matrix, na.rm TRUE) 50) { expr_matrix - log2(expr_matrix 1) } # 标准化分位数归一化 library(preprocessCore) expr_matrix - normalize.quantiles(as.matrix(expr_matrix)) rownames(expr_matrix) - names(keep_genes)[keep_genes] colnames(expr_matrix) - common_samples # 保存处理后的数据 saveRDS(expr_matrix, data/processed/expr_normalized.rds) saveRDS(sample_info, data/processed/sample_info.rds)这段代码做了几件事先对齐样本防止表达矩阵和元数据错位然后过滤掉在超过一半样本里表达量都低于 1 的基因减少多重检验负担接着判断是否需要 log2 转换通常原始 counts 或 FPKM 值跨度大取对数后更接近正态最后用分位数归一化消除样本间的技术偏差。normalize.quantiles来自preprocessCore包是芯片和 RNA-seq 数据常用的标准化方法。参数方面rowSums(expr_matrix 1) 0.5 * ncol(expr_matrix)里的1和0.5可以根据数据质量调整。如果测序深度低阈值可以降到 0.5 甚至 0如果样本量少比例可以放宽到 0.3。标准化方法也可以换成 TMM 或 DESeq2 的vst取决于下游分析假设。2.4 因果效应估计的实现细节02_estimate_effect.R是项目的核心。它读入预处理后的数据对每个基因拟合因果模型。简化后的主循环如下library(tidyverse) library(WeightIt) library(sandwich) expr_matrix - readRDS(data/processed/expr_normalized.rds) sample_info - readRDS(data/processed/sample_info.rds) # 定义协变量 covariates - c(age, sex, bmi, cell_type_fraction) # 创建结果存储 results - data.frame( gene rownames(expr_matrix), effect NA_real_, se NA_real_, p_value NA_real_ ) for (i in seq_len(nrow(expr_matrix))) { gene_expr - expr_matrix[i, ] # 构建分析数据框 df - data.frame( expr gene_expr, smoking sample_info$smoking_status, sample_info[, covariates] ) # 估计倾向性得分并加权 w - weightit( smoking ~ age sex bmi cell_type_fraction, data df, method ps, estimand ATE ) # 加权线性模型估计因果效应 fit - lm(expr ~ smoking, data df, weights w$weights) # 稳健标准误 coef_test - coeftest(fit, vcov vcovHC(fit, type HC0)) results$effect[i] - coef_test[smoking, Estimate] results$se[i] - coef_test[smoking, Std. Error] results$p_value[i] - coef_test[smoking, Pr(|t|)] } # 多重检验校正 results$fdr - p.adjust(results$p_value, method BH) # 保存结果 write.csv(results, output/tables/causal_effects.csv, row.names FALSE)逐段解释weightit函数用来估计倾向性得分并计算权重method ps表示用倾向性得分加权estimand ATE表示估计平均处理效应。加权后的lm拟合的是边际结构模型smoking的系数就是吸烟对基因表达的因果效应估计。coeftest配合vcovHC计算稳健标准误因为加权后残差方差不再恒定。最后用 Benjamini-Hochberg 方法校正 p 值控制错误发现率。这里有个容易翻车的地方weightit默认对协变量中的分类变量做哑变量编码如果sex是因子类型且只有两个水平没问题但如果协变量里有连续变量且分布偏态严重倾向性得分模型可能拟合不好权重极端值多。我一般会检查权重的分布如果某些样本权重超过 10考虑截断或改用匹配方法。参数调整建议estimand可以改成ATT估计处理组的平均处理效应取决于研究问题。协变量列表在covariates向量里改但要注意因果图——不是越多越好调整中介变量会引入偏倚。vcovHC的type参数可以换成HC3小样本下更稳健。2.5 结果可视化与解读03_visualize.R把因果效应估计结果画成火山图和森林图。火山图展示每个基因的效应大小和显著性森林图挑 top 基因展示置信区间。library(ggplot2) library(dplyr) results - read.csv(output/tables/causal_effects.csv) # 火山图 results - results %% mutate(significant ifelse(fdr 0.05 abs(effect) 0.5, Significant, Not significant)) p - ggplot(results, aes(x effect, y -log10(fdr), color significant)) geom_point(alpha 0.6, size 1) scale_color_manual(values c(Significant red, Not significant grey)) labs(x Causal effect (log2 fold change), y -log10(FDR), title Smoking effect on gene expression) theme_minimal() ggsave(output/figures/volcano.pdf, p, width 8, height 6)火山图的横轴是因果效应估计值纵轴是 FDR 的负对数。红点表示 FDR 小于 0.05 且效应绝对值大于 0.5 的基因。这个阈值不是固定的0.5可以根据生物学意义调整——有些研究用0.2捕捉微弱效应有些用1只看强效应。森林图用forestplot包或ggplot2的geom_pointrange实现展示 top 10 基因的效应和 95% 置信区间。提示火山图里显著基因的数量受阈值影响很大报告结果时最好同时给出不同阈值下的基因数避免选择性报告。3. 环境搭建与运行从零跑通全流程3.1 R 环境与依赖包安装项目基于 R 语言建议用 R 4.2 以上版本。依赖包分三类数据处理tidyverse、preprocessCore、因果推断WeightIt、sandwich、可视化ggplot2、forestplot。安装命令# 安装 CRAN 包 install.packages(c(tidyverse, WeightIt, sandwich, ggplot2, forestplot, yaml)) # 安装 Bioconductor 包 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(preprocessCore)WeightIt是因果推断的核心包负责倾向性得分加权。sandwich提供稳健标准误。preprocessCore在 Bioconductor 上用于分位数归一化。如果安装WeightIt时提示依赖CBPS或ebal缺失按提示补装即可。3.2 数据准备与格式要求项目本身不附带数据需要自己准备两份文件data/raw/expression.csv行是基因列是样本第一列是基因名第一行是样本 ID。data/raw/sample_metadata.csv至少包含sample_id、smoking_status0/1 或 never/current、age、sex、bmi等列。表达矩阵的数值可以是 counts、FPKM 或 TPM但要在01_preprocess.R里确认是否需要 log 转换。样本元数据里的smoking_status如果是字符型需要转成因子或数值型weightit才能正确处理。注意样本 ID 在表达矩阵列名和元数据sample_id列里必须完全一致包括大小写和分隔符。我见过因为-和_不一致导致样本全部对不上的情况排查了半天。3.3 一键运行与分步调试run_all.R是入口脚本内容很简单source(R/01_preprocess.R) source(R/02_estimate_effect.R) source(R/03_visualize.R)在 RStudio 里打开项目设置工作目录到项目根目录然后source(run_all.R)即可。如果数据量大比如 2 万个基因02_estimate_effect.R的循环可能跑几十分钟。建议先拿 100 个基因试跑确认流程无误后再全量跑。分步调试时每跑完一个脚本检查data/processed/和output/下有没有生成对应文件。01_preprocess.R跑完应该有expr_normalized.rds和sample_info.rds02_estimate_effect.R跑完应该有causal_effects.csv03_visualize.R跑完应该有volcano.pdf。3.4 参数配置与日志config.yml里可以集中管理参数preprocess: min_expression: 1 min_sample_fraction: 0.5 normalize_method: quantile causal: covariates: - age - sex - bmi - cell_type_fraction estimand: ATE fdr_threshold: 0.05 effect_threshold: 0.5 output: table_dir: output/tables figure_dir: output/figures在 R 脚本里用yaml::read_yaml(config.yml)读取替换硬编码参数。这样改阈值不用动代码也方便记录每次分析的参数组合。日志方面项目没有内置日志系统我一般会在每个脚本开头加sink()把控制台输出重定向到output/logs/下方便回溯。4. 避坑与常见问题跑这个项目容易翻车的几个地方4.1 样本对不上导致表达矩阵错位现象预处理后样本数变少或者因果效应结果全部不显著。原因表达矩阵列名和元数据sample_id格式不一致intersect后只剩少量样本或者顺序错乱导致吸烟状态和表达值张冠李戴。解决在01_preprocess.R里加一行stopifnot(all(colnames(expr_matrix) sample_info$sample_id))对不上直接报错别让它静默通过。另外用match而不是直接赋值确保顺序一致。4.2 倾向性得分模型不收敛现象weightit报错No valid propensity scores或警告Some weights are extreme。原因协变量之间存在完全共线性或者某个协变量在吸烟组和非吸烟组之间完全分离比如某细胞类型只在吸烟组出现。样本量太小也会导致逻辑回归不收敛。解决先检查协变量相关性去掉高度相关的变量。如果某个分类协变量水平太少合并类别。样本量小于 50 时考虑改用匹配方法或贝叶斯模型。权重极端值时在weightit里加trim 0.01截断极端权重。4.3 多重检验校正后没有显著基因现象causal_effects.csv里fdr全部大于 0.05。原因效应本身弱或者检验效能不足。基因表达数据噪声大单个基因的因果效应往往很小多重检验校正后很难存活。解决先看未校正的 p 值分布如果 p 值直方图均匀分布说明没有信号如果富集在 0 附近说明有信号但被校正压掉了。可以放宽 FDR 阈值到 0.1或者用qvalue包估计 pi0 后做校正。另外检查协变量调整是否过度——调整了中介变量会吃掉真实效应。4.4 内存不足导致循环中断现象跑02_estimate_effect.R时 R 崩溃或报cannot allocate vector of size。原因表达矩阵太大循环里反复创建数据框和模型对象内存碎片化。解决在循环里每跑完 1000 个基因调用一次gc()强制垃圾回收。或者把表达矩阵转成稀疏矩阵Matrix包减少内存占用。更彻底的办法是并行化用future包把基因分批并行跑但要注意weightit的随机种子设置。4.5 可视化结果与表格不一致现象火山图里显著基因数量和表格里fdr 0.05的行数对不上。原因火山图里用了额外的效应阈值abs(effect) 0.5而表格里只按 FDR 筛选。或者03_visualize.R读的是旧版结果文件。解决统一筛选条件在config.yml里定义fdr_threshold和effect_threshold表格和图形都用同一套参数。每次跑完02_estimate_effect.R后确认causal_effects.csv的修改时间避免读缓存。5. 二次开发与改造把这个因果推断框架用到自己的数据上5.1 替换暴露与结局变量这个项目的框架不限于吸烟和基因表达。任何二值暴露比如某种治疗、环境暴露和连续结局比如另一个组学指标、临床指标都可以套用。改造步骤在sample_metadata.csv里把smoking_status换成你的暴露变量列名比如treatment_group。在config.yml的covariates列表里换成你的混杂因素。在02_estimate_effect.R里把smoking ~ age sex bmi cell_type_fraction改成treatment_group ~ your_covariates。结局变量从基因表达换成你的指标矩阵行名改成指标名。如果暴露是多分类的比如从不吸烟、曾经吸烟、当前吸烟weightit支持多组倾向性得分但estimand要改成ATE并指定参考组。代码上把smoking ~ ...换成factor(smoking) ~ ...即可。5.2 换因果推断方法项目默认用 IPW 加权线性模型。如果数据里有时间依赖的混杂因素可以换成边际结构模型MSM配合逆概率处理加权IPTW。如果暴露是连续的比如吸烟包年数倾向性得分方法不适用改用广义倾向性得分GPS或工具变量法。换方法时WeightIt包支持method gbm梯度提升、cbps协变量平衡倾向性得分、ebal熵平衡。改一个参数就能切换但每种方法的假设不同选之前想清楚。5.3 并行化加速全基因组因果推断循环慢并行化是刚需。用futurefurrr改造library(future) library(furrr) plan(multisession, workers 4) estimate_effect - function(gene_name, expr_matrix, sample_info, covariates) { gene_expr - expr_matrix[gene_name, ] df - data.frame(expr gene_expr, smoking sample_info$smoking_status, sample_info[, covariates]) w - weightit(smoking ~ ., data df, method ps, estimand ATE) fit - lm(expr ~ smoking, data df, weights w$weights) coef_test - coeftest(fit, vcov vcovHC(fit, type HC0)) data.frame(gene gene_name, effect coef_test[smoking, Estimate], se coef_test[smoking, Std. Error], p_value coef_test[smoking, Pr(|t|)]) } results - future_map_dfr(rownames(expr_matrix), estimate_effect, expr_matrix expr_matrix, sample_info sample_info, covariates covariates, .options furrr_options(seed 123))plan(multisession)在 Windows 上开多进程Linux 和 Mac 可以用multicore。furrr_options(seed 123)保证随机数可复现。4 个 worker 大概能把时间压到原来的三分之一。注意weightit内部有随机成分并行时每个 worker 的种子要独立设置否则结果不可复现。5.4 结果验证与敏感性分析因果推断的结果依赖假设必须做敏感性分析。我一般会跑三组验证E 值衡量未测量混杂需要多强才能推翻结论。用EValue包计算。负对照选一个理论上不受吸烟影响的基因比如管家基因看因果效应是否接近零。不同模型换method gbm或cbps重跑看显著基因是否重叠。library(EValue) # 对 top 基因计算 E 值 top_gene - results %% arrange(fdr) %% slice(1) evalue - evalues.OR(est exp(top_gene$effect), lo exp(top_gene$effect - 1.96 * top_gene$se), hi exp(top_gene$effect 1.96 * top_gene$se), rare FALSE) print(evalue)E 值越大结论越稳健。如果 E 值小于 1.5说明一个中等强度的未测量混杂就能解释掉结果需要谨慎解读。5.5 一个我踩过的坑第一次跑这个项目时我直接把smoking_status当成数值型 0/1 喂给weightit结果它把 0 和 1 当成了连续变量做线性倾向性得分模型权重全乱。后来才发现weightit对暴露变量的类型敏感——二值暴露必须转成因子否则它不报错但结果不可信。从那以后我每次跑因果推断前都强制走一遍str(sample_info)和table(sample_info$smoking_status)确认变量类型和分布没问题再往下跑。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
