生信面试 R 语言统计题精选¶
一句话概述:生信面试中 R 语言题集中在数据可视化(ggplot2)、差异分析(DESeq2/edgeR)、统计检验(t检验/ANOVA/非参)和生存分析四个方向,掌握 tidyverse 生态就能应对大部分考题。
核心知识点速查表¶
| 概念 | 白话解释 |
|---|---|
| tidyverse | R语言数据科学工具集(dplyr+ggplot2+tidyr等) |
| ggplot2 | R的王牌绑图包,图层叠加语法 |
| dplyr | 数据操作五大动词:filter/select/mutate/arrange/summarise |
| DESeq2 | RNA-seq差异表达分析的标准工具 |
| Bioconductor | R的生物信息专用包管理平台 |
| 正态检验 | 判断数据是否符合正态分布(决定用什么检验) |
| 参数检验 | 假设数据符合某种分布(如t检验假设正态) |
| 非参检验 | 不假设分布(如Wilcoxon检验) |
一、R 基础操作(快速回顾)¶
# ========== R基础操作 ==========
# 数据类型
x <- 5 # 赋值(推荐用<-)
name <- "TP53" # 字符串
flag <- TRUE # 逻辑值
vec <- c(1, 2, 3, 4, 5) # 向量(一维数组)
mat <- matrix(1:12, nrow=3, ncol=4) # 矩阵(二维数组)
df <- data.frame( # 数据框(最常用!类似Excel表)
gene = c("TP53", "BRCA1", "EGFR"),
logFC = c(2.5, -1.8, 3.2),
padj = c(0.001, 0.03, 0.0001)
)
lst <- list(name="John", age=25) # 列表(可以混合类型)
# 向量操作
length(vec) # 长度
vec[1] # 第1个元素(R从1开始!不是0!)
vec[c(1,3)] # 第1和第3个元素
vec[vec > 3] # 筛选>3的元素
# 数据框操作
df$gene # 提取gene列
df[1, ] # 第1行
df[, 2] # 第2列
df[df$padj < 0.05, ] # 筛选padj<0.05的行
nrow(df) # 行数
ncol(df) # 列数
str(df) # 查看结构
summary(df) # 统计摘要
# 读写文件
df <- read.csv("data.csv") # 读CSV
df <- read.delim("data.tsv") # 读TSV(tab分隔)
write.csv(df, "output.csv", row.names=FALSE) # 写CSV
二、数据处理(dplyr,高频考点)¶
# ========== dplyr数据处理五大动词 ==========
library(dplyr) # 加载dplyr
# 示例数据
expr <- data.frame(
gene = paste0("Gene", 1:100), # 100个基因
logFC = rnorm(100, mean=0, sd=2), # 随机log2FC
padj = runif(100, 0, 1), # 随机padj
baseMean = rpois(100, lambda=500) # 随机表达量
)
# filter — 行筛选(选行)
sig_genes <- expr %>% # %>% 管道符(把结果传给下一步)
filter(abs(logFC) > 1 & padj < 0.05) # 显著差异基因
nrow(sig_genes) # 多少个显著基因
# select — 列选择(选列)
expr %>% select(gene, logFC, padj) # 只保留3列
expr %>% select(-baseMean) # 去掉baseMean列
# mutate — 新增/修改列
expr <- expr %>%
mutate(
direction = ifelse(logFC > 0, "Up", "Down"), # 新增方向列
neg_log_padj = -log10(padj) # 计算-log10(padj)
)
# arrange — 排序
expr %>% arrange(padj) # 按padj升序
expr %>% arrange(desc(abs(logFC))) # 按|logFC|降序
# summarise + group_by — 分组统计
expr %>%
filter(padj < 0.05) %>% # 先筛选显著的
group_by(direction) %>% # 按方向分组
summarise(
count = n(), # 每组计数
mean_logFC = mean(logFC), # 每组平均logFC
median_padj = median(padj) # 每组中位padj
)
# 面试经典:合并两个表格
gene_info <- data.frame(
gene = paste0("Gene", 1:50),
pathway = sample(c("apoptosis","cell_cycle","metabolism"), 50, replace=TRUE)
)
merged <- expr %>%
inner_join(gene_info, by="gene") # 按gene列合并(只保留交集)
# left_join 保留左表所有行
# right_join 保留右表所有行
# full_join 保留所有行
三、数据可视化(ggplot2,必考)¶
火山图(Volcano Plot)¶
# ========== 火山图(差异表达可视化经典) ==========
library(ggplot2) # 加载ggplot2
# 准备数据
expr <- expr %>%
mutate(
significance = case_when( # 根据条件分类
padj < 0.05 & logFC > 1 ~ "Up", # 显著上调
padj < 0.05 & logFC < -1 ~ "Down", # 显著下调
TRUE ~ "NS" # 不显著
)
)
# 画火山图
ggplot(expr, aes(x=logFC, y=-log10(padj), # X轴logFC,Y轴-log10(padj)
color=significance)) + # 按显著性着色
geom_point(alpha=0.6, size=1.5) + # 散点图
scale_color_manual(values=c( # 自定义颜色
"Up"="red", "Down"="blue", "NS"="grey"
)) +
geom_hline(yintercept=-log10(0.05), # padj=0.05水平线
linetype="dashed", color="grey50") +
geom_vline(xintercept=c(-1, 1), # logFC=±1垂直线
linetype="dashed", color="grey50") +
labs(title="Volcano Plot", # 标题
x="log2(Fold Change)", # X轴标签
y="-log10(adjusted p-value)") + # Y轴标签
theme_classic() # 经典主题
ggsave("volcano_plot.pdf", width=8, height=6) # 保存为PDF
热图(Heatmap)¶
# ========== 热图(表达谱可视化) ==========
library(pheatmap) # 加载pheatmap
# 模拟表达矩阵
set.seed(42) # 设置随机种子(可重现)
expr_matrix <- matrix(rnorm(200), # 随机数据
nrow=20, ncol=10)
rownames(expr_matrix) <- paste0("Gene", 1:20) # 行名=基因
colnames(expr_matrix) <- paste0("Sample", 1:10) # 列名=样本
# 注释信息
annotation_col <- data.frame(
Group = rep(c("Control", "Treatment"), each=5), # 分组
row.names = colnames(expr_matrix)
)
# 画热图
pheatmap(expr_matrix,
scale = "row", # 按行标准化(Z-score)
clustering_method = "ward.D2", # 聚类方法
annotation_col = annotation_col, # 列注释
color = colorRampPalette(c("blue","white","red"))(100), # 颜色
fontsize_row = 8, # 行字体大小
fontsize_col = 10, # 列字体大小
main = "Gene Expression Heatmap") # 标题
箱线图 + 统计检验¶
# ========== 箱线图 + 统计检验 ==========
library(ggplot2)
library(ggpubr) # 加载ggpubr(添加统计检验)
# 模拟数据
data <- data.frame(
group = rep(c("Control", "Treatment"), each=30), # 两组各30个
expression = c(rnorm(30, mean=5, sd=1), # 对照组
rnorm(30, mean=7, sd=1.5)) # 处理组
)
# 箱线图 + t检验p值
ggplot(data, aes(x=group, y=expression, fill=group)) +
geom_boxplot(alpha=0.7) + # 箱线图
geom_jitter(width=0.2, alpha=0.3) + # 散点抖动
stat_compare_means(method="t.test", # 添加t检验
label="p.format") + # 显示p值
labs(title="Gene Expression Comparison",
y="Expression (TPM)") +
theme_classic() +
theme(legend.position="none") # 隐藏图例
ggsave("boxplot_comparison.pdf", width=5, height=5)
四、统计检验(核心考点)¶
# ========== 统计检验完整流程 ==========
# === 第 1 步:正态性检验(决定用参数还是非参) ===
# Shapiro-Wilk检验(p>0.05 = 正态分布)
shapiro.test(data$expression[data$group == "Control"]) # 检验对照组
shapiro.test(data$expression[data$group == "Treatment"]) # 检验处理组
# === 第 2 步:方差齐性检验 ===
# Levene检验(p>0.05 = 方差齐)
car::leveneTest(expression ~ group, data=data)
# === 第 3 步:选择检验方法 ===
# 情况1:正态+两组 → t检验
t.test(expression ~ group, data=data) # 独立样本t检验
t.test(expression ~ group, data=data, var.equal=TRUE) # 方差齐时用等方差t检验
t.test(before, after, paired=TRUE) # 配对t检验
# 情况2:不正态+两组 → Wilcoxon检验(非参数)
wilcox.test(expression ~ group, data=data) # Mann-Whitney U检验
# 情况3:正态+三组及以上 → ANOVA
# 模拟三组数据
data3 <- data.frame(
group = rep(c("A", "B", "C"), each=20),
value = c(rnorm(20,5,1), rnorm(20,7,1), rnorm(20,6,1))
)
aov_result <- aov(value ~ group, data=data3) # 单因素方差分析
summary(aov_result) # 查看结果
TukeyHSD(aov_result) # 事后两两比较
# 情况4:不正态+三组及以上 → Kruskal-Wallis检验
kruskal.test(value ~ group, data=data3) # Kruskal-Wallis检验
pairwise.wilcox.test(data3$value, data3$group, # 事后两两比较
p.adjust.method="BH") # BH多重检验校正
# === 多重检验校正 ===
pvalues <- c(0.01, 0.04, 0.03, 0.001, 0.06) # 原始p值
p.adjust(pvalues, method="BH") # BH校正(FDR)
p.adjust(pvalues, method="bonferroni") # Bonferroni校正(更保守)
# 面试要点:
# BH(Benjamini-Hochberg)= 控制FDR(假发现率),生信最常用
# Bonferroni = 控制FWER(家族错误率),最保守
# 为什么要多重检验校正?因为同时做很多检验,偶然发现显著的概率增加
五、差异表达分析(DESeq2)¶
# ========== DESeq2差异表达分析流程 ==========
# BiocManager::install("DESeq2") # 安装
library(DESeq2) # 加载
# 准备输入:原始counts矩阵 + 样本信息
# counts矩阵:行=基因,列=样本,值=raw read counts
count_matrix <- as.matrix(read.csv("counts.csv", row.names=1))
# 样本信息
coldata <- data.frame(
condition = factor(c(rep("Control",3), # 3个对照
rep("Treatment",3)), # 3个处理
levels=c("Control","Treatment")), # Control为参考
row.names = colnames(count_matrix) # 行名=样本名
)
# 创建DESeq2对象
dds <- DESeqDataSetFromMatrix(
countData = count_matrix, # counts矩阵
colData = coldata, # 样本信息
design = ~ condition # 设计公式
)
# 预过滤低表达基因(加速运算)
keep <- rowSums(counts(dds)) >= 10 # 至少总计10个reads
dds <- dds[keep, ] # 过滤
# 运行差异分析
dds <- DESeq(dds) # 核心函数(做了3步)
# 1. estimateSizeFactors — 标准化
# 2. estimateDispersions — 估计离散度
# 3. nbinomWaldTest — 负二项分布Wald检验
# 提取结果
res <- results(dds,
contrast=c("condition","Treatment","Control"), # 比较组
alpha=0.05) # 显著性阈值
res <- res[order(res$padj), ] # 按padj排序
summary(res) # 结果摘要
# 保存结果
write.csv(as.data.frame(res), "DESeq2_results.csv")
# 可视化:MA图
plotMA(res, ylim=c(-5,5), # MA图
main="DESeq2 MA Plot")
# 可视化:PCA图
vsd <- vst(dds, blind=FALSE) # 方差稳定化转换
plotPCA(vsd, intgroup="condition") + # PCA图
theme_classic()
六、生存分析¶
# ========== Kaplan-Meier生存分析 ==========
library(survival) # 生存分析核心包
library(survminer) # 漂亮的生存曲线
# 模拟临床数据
clinical <- data.frame(
patient = paste0("P", 1:100),
time = rexp(100, rate=0.02), # 生存时间(指数分布)
status = sample(0:1, 100, # 1=事件(死亡), 0=删失
replace=TRUE, prob=c(0.3,0.7)),
gene_group = sample(c("High","Low"), 100, # 基因高/低表达组
replace=TRUE)
)
# 拟合生存模型
surv_obj <- Surv(time=clinical$time, # 生存时间
event=clinical$status) # 事件状态
fit <- survfit(surv_obj ~ gene_group, # 按基因分组
data=clinical)
print(fit) # 打印中位生存时间
# 画Kaplan-Meier生存曲线
ggsurvplot(fit,
data=clinical,
pval=TRUE, # 显示log-rank p值
risk.table=TRUE, # 风险表
conf.int=TRUE, # 置信区间
palette=c("red","blue"), # 颜色
xlab="Time (months)", # X轴
ylab="Survival Probability",# Y轴
title="Kaplan-Meier Survival Curve",
ggtheme=theme_classic()) # 主题
# Cox回归(多变量)
cox_model <- coxph(surv_obj ~ gene_group + age + stage,
data=clinical) # 多因素Cox回归
summary(cox_model) # 查看HR、p值、CI
七、常见报错与解决¶
| 问题 | 原因 | 解决方法 |
|---|---|---|
| object not found | 变量名拼写错误或未赋值 | 检查变量名,ls() 列出所有变量 |
| package not found | 包未安装 | install.packages("pkg") 或 BiocManager |
| non-numeric argument | 对非数值数据做数学运算 | as.numeric() 转换,检查NA |
| Error in plot | 数据中有NA或Inf | na.omit() 或 is.finite() 过滤 |
| factor levels错误 | 分组变量顺序不对 | factor(x, levels=c("Control","Treatment")) |
| 中文乱码 | 编码问题 | Sys.setlocale("LC_ALL","UTF-8") |
八、面试高频问题¶
Q1:R 中 data.frame 和 matrix 的区别?
data.frame 可以混合不同类型(数值+字符),matrix 只能存同一类型。data.frame 有列名可以用
$访问,适合存储样本信息等异质数据。matrix 运算速度更快,适合存储表达矩阵等同质数据。
Q2:什么时候用 t 检验,什么时候用 Wilcoxon?
数据正态分布 + 方差齐 → t检验(参数检验,统计效力更强)。数据不正态或样本量很小 → Wilcoxon(非参数检验,不依赖分布假设)。先做 Shapiro-Wilk 正态性检验再决定。
Q3:DESeq2 和 edgeR 的区别?
都用负二项分布建模counts数据。DESeq2用中位数比值法(median of ratios)标准化,Wald检验;edgeR用TMM标准化,精确检验或QLF检验。样本量<3时edgeR更稳定,大样本两者结果接近。一般推荐DESeq2(更保守、假阳性更低)。
Q4:p值和FDR(q值)的区别?
p值是单次检验犯错的概率。当同时检验上万个基因时,即使p<0.05也会有大量假阳性(5%×20000=1000个假阳性)。FDR(False Discovery Rate)控制的是"你说显著的里面有多少比例是假的",padj<0.05意味着预期假阳性不超过5%。生信中必须用FDR校正后的padj。
九、速查表¶
# === 生信 R 语言速查 ===
# 数据操作 (dplyr)
df %>% filter(padj < 0.05) # 筛选行
df %>% select(gene, logFC) # 选择列
df %>% mutate(new_col = logFC * 2) # 新增列
df %>% arrange(padj) # 排序
df %>% group_by(group) %>% summarise(n=n()) # 分组统计
# 可视化 (ggplot2)
ggplot(df, aes(x, y, color=group)) + # 基本语法
geom_point() + # 散点
geom_boxplot() + # 箱线图
geom_bar(stat="identity") + # 柱状图
theme_classic() # 主题
# 统计检验选择
# 2组正态 → t.test()
# 2组非正态 → wilcox.test()
# 3+组正态 → aov() + TukeyHSD()
# 3+组非正态 → kruskal.test()
# 多重校正 → p.adjust(method="BH")
# DESeq2
dds <- DESeqDataSetFromMatrix(counts, coldata, ~condition)
dds <- DESeq(dds)
res <- results(dds, alpha=0.05)
# 生存分析
fit <- survfit(Surv(time, status) ~ group, data)
ggsurvplot(fit, pval=TRUE, risk.table=TRUE)
参考资料:R for Data Science (Wickham 2023)、DESeq2 vignette、ggplot2: Elegant Graphics (Wickham 2024)、Bioconductor