实验课4 差异基因分析-聚类 R代码整理¶
对应材料：
实验/课件/实验课4-差异基因分析-聚类.pptx回放/生物信息学（协和班）第11周星期4第8,9节_笔记.txt- 老师板书整理
说明¶
- 这份整理以老师板书为主，结合课件和第11周实验课回放做了必要补全。
- 其中“板书直接有的代码主线”我尽量原样保留。
- 个别地方老师板书只写了思路，没有把能运行的细节完全写全，我按课件内容补齐，并单独标出来。
一、结合板书整理后的完整 R 代码¶
library(DESeq2)
library(factoextra)
library(pheatmap)
countData <- read.csv("chose_TCGAcount.csv")
rownames(countData) <- countData$Gene
countData <- countData[, -1]
countData <- countData[rowSums(countData > 0) > 0, ]
meta <- data.frame(
id = colnames(countData),
stringsAsFactors = FALSE
)
meta$type <- sapply(
meta$id,
function(x) strsplit(x, "_")[[1]][4]
)
meta$label <- ifelse(meta$type == "01A", "tumor", "pang")
meta$label <- factor(meta$label)
rownames(meta) <- meta$id
dds <- DESeqDataSetFromMatrix(
countData = countData,
colData = meta,
design = ~ label
)
dds <- DESeq(dds)
res <- results(dds, contrast = c("label", "tumor", "pang"))
res_df <- as.data.frame(res)
diff_gene_up <- subset(
res_df,
padj < 0.05 & log2FoldChange > 1
)
diff_gene_down <- subset(
res_df,
padj < 0.05 & log2FoldChange < -1
)
DEGS <- c(rownames(diff_gene_up), rownames(diff_gene_down))
data_choose <- countData[DEGS, ]
data_choosescale <- scale(t(data_choose))
d <- dist(data_choosescale)
fit1 <- hclust(d, method = "ward.D2")
plot(fit1, hang = -1, cex = 0.3, main = "clustering")
groups <- cutree(fit1, k = 2)
rect.hclust(fit1, k = 2, border = "red")
annotation_row <- data.frame(
label = as.character(meta$label),
stringsAsFactors = FALSE
)
rownames(annotation_row) <- meta$id
pheatmap::pheatmap(
data_choosescale,
annotation_row = annotation_row,
clustering_method = "ward.D2",
show_colnames = FALSE,
show_rownames = FALSE,
main = "Heatmap of DEGs (rows: samples, columns: genes)",
fontsize = 6,
fontsize_row = 6,
fontsize_col = 6
)
二、这份代码对应老师板书的哪几块¶
1. countData 读取与预处理¶
对应第一张板书：
countData <- read.csv("chose_TCGAcount.csv")
rownames(countData) <- countData$Gene
countData <- countData[, -1]
countData <- countData[rowSums(countData > 0) > 0, ]
这里的意思是：
- 读入原始
count matrix - 把
Gene这一列变成行名 - 删除原来的
Gene列 - 去掉完全不表达的低表达基因
这一步和课件最后总流程图是完全对应的：
2. 样本注释 meta¶
对应第一张板书下半部分：
meta <- data.frame(id = colnames(countData))
meta$type <- sapply(meta$id, function(x) strsplit(x, "_")[[1]][4])
meta$label <- ifelse(meta$type == "01A", "tumor", "pang")
meta$label <- factor(meta$label)
这一步就是从样本名里拆出样本类型，再转成分组变量：
01A->tumor- 其他这里课堂里按这套数据记作
pang
也就是课件里的：
Extract sample type (01A / 11A)Assign group label (tumor / pang)Convert to factor
3. DESeq2 差异分析¶
对应第二张板书：
dds <- DESeqDataSetFromMatrix(
countData = countData,
colData = meta,
design = ~ label
)
dds <- DESeq(dds)
res <- results(dds, contrast = c("label", "tumor", "pang"))
这一段的逻辑很固定：
- 用表达矩阵和样本信息构建
dds - 跑
DESeq2 - 提取
tumor vs pang的结果
老师这节实验课板书基本就是照着课件第 14 页那三步在写。
4. 上调/下调差异基因筛选¶
对应第二张板书最下面：
diff_gene_up <- subset(res, padj < 0.05 & log2FoldChange > 1)
diff_gene_down <- subset(res, padj < 0.05 & log2FoldChange < -1)
这对应课件里最关键的筛选原则：
padj < 0.05|log2FoldChange| > 1
并按方向分成：
UP：log2FC > 1DOWN：log2FC < -1
5. 提取 DEGs 并做标准化¶
对应第三张板书前半部分：
DEGS <- c(rownames(diff_gene_up), rownames(diff_gene_down))
data_choose <- countData[DEGS, ]
data_choosescale <- scale(t(data_choose))
这里有两个重要点：
- 先把上调和下调基因合并起来
scale()默认按列标准化，所以老师特意加了一次转置
课件第 17 页专门强调过这一点：
6. 层次聚类¶
对应第三张板书中间部分：
d <- dist(data_choosescale)
fit1 <- hclust(d, method = "ward.D2")
plot(fit1, hang = -1, cex = 0.3, main = "clustering")
groups <- cutree(fit1, k = 2)
rect.hclust(fit1, k = 2, border = "red")
这一段就是：
- 先算样本间距离
dist() - 再用
hclust()做层次聚类 - 再画树状图
- 然后把树切成两类
7. 热图¶
对应第四张板书中间部分：
pheatmap::pheatmap(
data_choosescale,
annotation_row = annotation_row,
clustering_method = "ward.D2",
show_colnames = FALSE,
show_rownames = FALSE,
main = "Heatmap of DEGs (rows: samples, columns: genes)",
fontsize = 6,
fontsize_row = 6,
fontsize_col = 6
)
老师板书里热图参数基本就是照着课件第 21 页写的，重点是：
- 输入是标准化后的
DEGs表达矩阵 annotation_row用来标注样本分组- 聚类方法还是
ward.D2 - 行列名最好别显示，不然太密
三、我按课件和回放补全了哪些地方¶
1. library() 这三句¶
老师板书右上角只写了包名：
DESeq2factoextrapheatmap
我这里补成了：
这样是完整可运行的写法。
2. rownames(meta) <- meta$id¶
老师板书里只写到了：
但真正运行 DESeqDataSetFromMatrix() 时，colData 的行名要和 countData 的列名对上，所以我补了：
这一步是“为了能顺利跑通而补的”，不是板书原句，但逻辑上是必须的。
3. res_df <- as.data.frame(res)¶
老师板书写的是直接：
这在有些环境里也能跑，但为了后面 subset() 和 rownames() 更稳定，我加了：
然后再对 res_df 做筛选。
4. annotation_row¶
老师板书写的是：
我保留了这个思路，只顺手加了：
避免旧版本 R 自动转因子。
四、结合第11周实验课回放，老师口头强调了什么¶
虽然第11周这个回放主体是在讲单细胞实验课，但里面有几句对这份 DESeq2 + 聚类 代码的复习方式很有帮助。
1. 实验课更重要的是“流程”和“图”，不是死记每个命令¶
回放里老师明确说：
- 实验课考试更可能让你画流程图
- 不一定让你完整编码
- 但基础编程思路要会
这对这节差异基因分析实验课的意思就是：
- 你要会写这条主线
count matrix输入
-> 基因名设为行名
-> 去低表达基因
-> 构建样本注释meta
-> DESeq2差异分析
-> 筛选上调/下调DEGs
-> 提取DEGs表达矩阵
-> 转置并标准化
-> dist计算距离
-> hclust聚类
-> 画树状图和热图
2. 写流程图时，不必死记所有函数名¶
回放里老师讲得很直接：
- 写流程图的时候不需要把每个命令都写出来
- 如果命令记得住可以写
- 记不住就写“创建对象”“过滤”“标准化”“聚类”“画图”这些步骤
所以这节课如果考应用题，最稳妥的写法不是硬背整串命令，而是把关键步骤和关键阈值写出来。
3. 关键阈值和关键判断比命令名更重要¶
对这节课来说，最该记的是这些：
padj < 0.05log2FoldChange > 1为上调log2FoldChange < -1为下调ward.D2是层次聚类方法scale()前要注意转置
也就是说，考试时比起死背 DESeqDataSetFromMatrix() 这个名字，更重要的是你知道：
- 输入是什么
- 怎么分组
- 按什么标准筛
DEGs - 为什么要标准化
- 聚类和热图分别在展示什么
五、你复习时最该盯住的代码逻辑¶
如果不要求逐行默写代码，这节课最值得背的是下面这一条主线：
读入count matrix
-> 设置基因为行名
-> 删除Gene列
-> 去掉全零/低表达基因
-> 从样本名提取01A/11A
-> 标注tumor/pang
-> 构建DESeq2对象
-> 提取差异分析结果
-> 按padj和log2FC筛选DEGs
-> 取出DEGs表达矩阵
-> 转置并标准化
-> 计算距离并层次聚类
-> 画树状图和热图
其中最容易单独拿出来考的，是这几个点：
- 为什么
Gene要变成行名 - 为什么
meta$label要转成factor - 为什么
tumor vs pang时log2FoldChange > 0代表在tumor中上调 - 为什么
scale()前要先转置一次 - 为什么树状图和热图都能反映样本聚类关系
六、如果考试要你手写伪代码，可以写成这样¶
输入：count matrix（行是基因，列是样本）
1. 读取count矩阵
2. 将Gene列设置为行名，并删除原Gene列
3. 去除完全不表达或低表达基因
4. 根据样本名提取样本类型
5. 将样本分为tumor和pang，并构建样本注释表
6. 使用DESeq2构建差异分析对象
7. 执行标准化、离散度估计和差异检验
8. 提取tumor与pang的比较结果
9. 按padj < 0.05且|log2FC| > 1筛选差异基因
10. 合并上调和下调基因，提取其表达矩阵
11. 对表达矩阵转置并标准化
12. 计算样本间距离并进行层次聚类
13. 绘制树状图和热图
14. 输出聚类和热图结果
七、一句提醒¶
这节实验课的 R 代码，最重要的不是把每个函数名一字不差背下来，而是你要会解释：
- 输入是什么
- 分组信息是怎么来的
- 差异基因按什么标准筛
- 为什么聚类前要标准化
- 树状图和热图各自在说明什么