image.png
差异表达分析内容：
• 基因表达量的标准化方法及可视化
➢ counts，RPKM，FPKM，TPM
➢ PCA图、热图等
• 差异表达分析及可视化
➢ limma-voom，edgeR，DESeq2
➢ 差异基因的热图和火山图
• 三个软件包的差异分析结果比较及筛选
➢ logFC含义
➢ 相关性图
1.常用的基因表达的标准化方法
• 现在常用的基因定量方法包括：RPM, RPKM, FPKM, TPM。
• 这些表达量的主要区别是：通过不同的标准化方法为转录本丰度提供一个
数值表示，以便于后续差异分析。
• 标准化的主要目的是去除测序数据的技术偏差：测序深度和基因长度。
• 测序深度：同一条件下，测序深度越深，基因表达的read读数越多。
• 基因长度：同一条件下，不同的基因长度产生不对等的read读数，基
因越长，该基因的read读数越高。
https://mp.weixin.qq.com/s/KSMzgKBlgF2qIadME5nWhw
Counts值:对给定的基因组参考区域，计算比对上的read数，又称为raw count（RC）。计数结果的差异的影响因素：落在参考区域上下限的read是否需要被统计，按照什么样的标准进行统计。
(1)RPM (Reads per million mapped reads)
image.png
(2)RPKM/FPKM
image.png
(3)TPM (Transcript per million)
image.png
三种表达量的比较:
image.png
当我们拿到一个counts之后，首先要对他进行PCA,热图这样的一个初步探索，看看数据的表达量，差异情况等。
这是一个比较粗糙的PCA图
image.png
相关性图
相关性热图
加上分组信息
以上图用到的脚本顺序
以上就是对数据合理性的检查，下面进行差异分析：
2.差异表达分析及可视化
差异分析的目标:
image.png
差异表达分析的统计学模型:在测序出现之前一般都是用芯片
基因芯片是连续型数据：limma来处理
测序数据是离散型数据
学习统计学的网站：https://seeing-theory.brown.edu/
泊松分布：
对于泊松分布而言，其均值和方差是相等的，但是我们的
数据确不符合这样的规律。
紫色实线是泊松分布的拟合结果。
橙色实线是负二项分布的拟合结果。（DESeq2）
橙色虚线是edgeR软件的拟合结果。
image.png
表达差异分析的工具:
image.png
(1)Limma差异分析
image.png
image.png
limma数据准备：转录组测序数据→表达矩阵→分组矩阵
### 读取基因表达矩阵
rm(list = ls())
options(stringsAsFactors = F)
load(file = "data/Step01-airway2exprSet.Rdata")
dim(exprSet)
table(group_list)
library(limma)
library(edgeR)
#### 第一步，创建设计矩阵和对比
# 假设数据符合正态分布，构建线性模型
# 0代表x线性模型的截距为0
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
design
#### 第二步，进行差异表达分析
# 将表达矩阵转换为edgeR的DGEList对象
dge <- DGEList(counts=exprSet)
# 进行标准化和表达量计算
dge <- calcNormFactors(dge)
# 在进行对数转换前会给所有基因的计数加上3
# 以避免对零取对数，且可减小低表达基因之间的差异
logCPM <- cpm(dge, log=TRUE, prior.count=3)
logCPM[1:4,1:4]
# 设置需要进行对比的分组，需要修改
comp='trt-untrt'
f <- 'data/Step03-limma_fit2.Rdata'
if(!file.exists(f)){
v <- voom(dge,design,plot=TRUE, normalize="quantile")
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(contrasts=c(comp),levels = design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
save(fit2,file = "data/Step03-limma_fit2.Rdata")
}
#### 第三步，提取过滤差异分析结果
load(file = "data/Step03-limma_fit2.Rdata")
tmp <- topTable(fit2, coef=comp, n=Inf)
DEG_limma_voom <- na.omit(tmp)
head(DEG_limma_voom)
差异分析结果
# 取表达差异倍数和p值两列
nrDEG <- DEG_limma_voom[,c(1,4)]
colnames(nrDEG) <- c('log2FoldChange','pvalue')
save(nrDEG, DEG_limma_voom, file = "data/Step03-limma_voom_nrDEG.Rdata")
image.png
（2）edgeR差异分析
image.png
image.png
image.png
rm(list = ls())
options(stringsAsFactors = F)
# 读取基因表达矩阵信息,需要修改
load(file = "data/Step01-airway2exprSet.Rdata")
# 查看分组信息和表达矩阵数据
dim(exprSet)
table(group_list)
suppressMessages(library(edgeR))
#### 第一步，构建edgeR的DGEList对象，并过滤
DEG <- DGEList(counts=exprSet,group=factor(group_list))
# 保留在至少在两个样本中cpm值大约1的基因
keep <- rowSums(cpm(DEG)>1) >= 2
table(keep)
DEG <- DEG[keep, , keep.lib.sizes=FALSE]
DEG$samples$lib.size <- colSums(DEG$counts)
# 归一化基因表达分布
DEG <- calcNormFactors(DEG)
# 增加一列$norm.factors
DEG$samples
# 假设数据符合正态分布，构建线性模型
# 0代表x线性模型的截距为0，需要修改
design <- model.matrix(~0+factor(group_list))
rownames(design) <- colnames(DEG)
colnames(design) <- levels(factor(group_list))
#### 第二步，差异表达分析
f <- 'data/Step03-edgeR_lrt.Rdata'
if(!file.exists(f)){
# 计算线性模型的参数
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)
# 拟合线性模型
fit <- glmFit(DEG, design)
# 进行差异分析
# 1,-1意味着前比后
lrt <- glmLRT(fit, contrast=c(1,-1))
# 参考链接：https://www.biostars.org/p/110861/
save(lrt, file = "data/Step03-edgeR_lrt.Rdata")
}
#### 第三步，提取过滤差异分析结果
load("data/Step03-edgeR_lrt.Rdata")
edgeR_DEG <- topTags(lrt, n=nrow(DEG))
edgeR_DEG <- as.data.frame(edgeR_DEG)
head(edgeR_DEG)
image.png
# 取表达差异倍数和p值两列
nrDEG <- edgeR_DEG[,c(1,5)]
colnames(nrDEG) <- c('log2FoldChange','pvalue')
save(nrDEG, edgeR_DEG, file = "data/Step03-edgeR_nrDEG.Rdata")
image.png
（3）DESeq2差异分析
image.png
image.png
image.png
rm(list = ls())
options(stringsAsFactors = F)
# 读取基因表达矩阵信息
load(file = "data/Step01-airway2exprSet.Rdata")
# 查看分组信息和表达矩阵数据
dim(exprSet)
table(group_list)
suppressMessages(library(DESeq2))
#### 第一步，构建DESeq2的DESeq对象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,
design = ~ group_list)
#### 第二步，进行差异表达分析
f <- 'data/Step03-DESeq2_dds2.Rdata'
if(!file.exists(f)){
dds2 <- DESeq(dds)
# 保存差异表达分析结果
save(dds2, file = "data/Step03-DESeq2_dds2.Rdata")
}
#### 第二步，提取差异分析结果
load(file = "data/Step03-DESeq2_dds2.Rdata")
# 提取差异分析结果，trt组对untrt组的差异分析结果
res <- results(dds2,contrast=c("group_list","trt","untrt"))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG <- as.data.frame(resOrdered)
# 去除差异分析结果中包含NA值的行
DESeq2_DEG = na.omit(DEG)
image.png
## 取出包含logFC和P值的两列
nrDEG=DESeq2_DEG[,c(2,6)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
save(nrDEG, DESeq2_DEG, file = "data/Step03-DESeq2_nrDEG.Rdata")
image.png
三大差异分析包的比较:
• limma，edgeR，DESeq2三大包基本是做转录组差异分析的金标准，
大多数转录组的文章都是用这三个R包进行差异分析。
• edgeR差异分析速度快，得到的基因数目比较多，假阳性高（实际不
差异结果差异）。DESeq2差异分析速度慢，得到的基因数目比较少，
假阴性高（实际差异结果不差异）。
• 需要注意的是制作分组信息的因子向量是，因子水平的前后顺序，在
R的很多模型中，默认将因子向量的第一个水平看作对照组。
https://mp.weixin.qq.com/s/SxoP_Br6a8rtCB9Zorc4aA
下面来做一个比较：
image.png
rm(list = ls())
options(stringsAsFactors = F)
# 读取3个软件的差异分析结果
load(file = "data/Step03-limma_voom_nrDEG.Rdata")
limma_DEG <- nrDEG
load(file = "data/Step03-DESeq2_nrDEG.Rdata")
DESeq2_DEG <- nrDEG
load(file = "data/Step03-edgeR_nrDEG.Rdata")
edgeR_DEG <- nrDEG
# 提取所有差异表达的基因名
geneLists <- unique(c(rownames(limma_DEG),rownames(DESeq2_DEG),rownames(edgeR_DEG)))
# 将三个差异分析结果合并成数据框
DEGLists <- data.frame(limma=limma_DEG[geneLists,1],
DESeq2=DESeq2_DEG[geneLists,1],
edgeR=edgeR_DEG[geneLists,1])
library(corrplot)
# 去除NA值
DEGLists <- na.omit(DEGLists)
# 计算相关性
DEGLists_cor <- cor(DEGLists)
DEGLists_cor
# limma DESeq2 edgeR
#limma 1.0000000 0.9820845 0.9821562
#DESeq2 0.9820845 1.0000000 0.9999761
#edgeR 0.9821562 0.9999761 1.0000000
corrplot(DEGLists_cor)
相关性分析：可以看到基本上是没有差异的
logFC(log2foldchange)值的含义：
➢ 简单来讲正的logFC值表示上调的对数倍数，负的logFC值
表示下调的对数倍数。
• logFC= log(expr1)-log(expr2)
• 例如：log2FC = log2(处理组的表达量) -log2(对照组的
表达量)
• 如果一个基因在对照组中的表达量为16，在处理组的表
达量为4，log2FC为-2。
➢ limma/voom、edgeR和DESeq2的logFC值的区别
• 算术平均数和几何平均数那个更适合计算logFC值。
• limma使用对数化的表达矩阵作为输入，所以将零假设
指定在平均log值（对数几何平均数）。
• edgeR和DESeq2使用原始的count矩阵作为输入，所以
将原假设指定在count的平均值上（算术平均数）。
logFC阈值的取舍:
➢ 通过logFC筛选差异基因的结果的标准有两个：
logFC的阈值和矫正后的p值（多重检验会导致假
阳性偏高）。
• 对于p值，业内标准0.05和0.01，不在赘述。
• 一般是在1.2到2之间，筛选到差异基因的数目
在500-1000左右为宜
• 根据logFC的统计指标确定阈值的方法是计算
logFC绝对值的平均数与2倍标准差之和。
（正态分布）
3.差异基因的热图和火山图：
(1)火山图
rm(list = ls())
options(stringsAsFactors = F)
library(pheatmap)
# 加载原始表达矩阵
load(file = "data/Step01-airway2exprSet.Rdata")
# 读取3个软件的差异分析结果
load(file = "data/Step03-limma_voom_nrDEG.Rdata")
limma_DEG <- nrDEG
load(file = "data/Step03-DESeq2_nrDEG.Rdata")
DESeq2_DEG <- nrDEG
load(file = "data/Step03-edgeR_nrDEG.Rdata")
edgeR_DEG <- nrDEG
#### 用DEG数据来作火山图
library(ggplot2)
# 根据需要修改DEG的值
DEG <- limma_DEG #这里用的是limma的DEG，可以换成另外两个
colnames(DEG)
# 使用基础函数plot绘图
plot(DEG$log2FoldChange,-log2(DEG$pvalue))
基础火山图：粗暴简陋
# 确定差异表达倍数
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )
# 取前两位小数
logFC_cutoff <- round(logFC_cutoff, 2)
# 确定上下调表达基因
DEG$change = as.factor(ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'STABLE'))
table(DEG$change)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
g <- ggplot(data=DEG,
aes(x=log2FoldChange,
y=-log10(pvalue),color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20))) +
xlab("log2 fold change") +
ylab("-log10 p-value") +
ggtitle( this_tile ) +
theme(plot.title = element_text(size=15,hjust = 0.5)) +
scale_colour_manual(values = c('blue','black','red'))## 与分组一一对应
print(g)
image.png
### 自定义火山图
## 第一种方法适合于少数基因的展示
library(ggrepel)
DEG$symbol <- rownames(DEG)
DEG$label <- ifelse(DEG$pvalue < 0.00001 & abs(DEG$log2FoldChange) >= 3,DEG$symbol,"")
g1 <- g + geom_text_repel(data = DEG, aes(x = DEG$log2FoldChange,
y = -log10(DEG$pvalue),
label = label),
size = 3,box.padding = unit(0.5, "lines"),
point.padding = unit(0.8, "lines"),
segment.color = "black",
show.legend = FALSE)
print(g1)
image.png
## 第二种方法适合于大量的基因名展示
library(dplyr)
for_label <- DEG %>%
filter(abs(log2FoldChange) > 4 & -log10(pvalue) > -log10(0.0001))
g2 <- g +
geom_point(size = 3, shape = 1, data = for_label) +
ggrepel::geom_label_repel(
aes(label = symbol),
data = for_label,
color="black"
)
print(g2)
ggsave(g2,filename = 'figures/Step04-DEG_volcano_senior2.png')
image.png
(2)热图
rm(list = ls())
options(stringsAsFactors = F)
library(pheatmap)
# 加载原始表达矩阵
load(file = "data/Step01-airway2exprSet.Rdata")
# 读取3个软件的差异分析结果
load(file = "data/Step03-limma_voom_nrDEG.Rdata")
limma_DEG <- nrDEG
load(file = "data/Step03-DESeq2_nrDEG.Rdata")
DESeq2_DEG <- nrDEG
load(file = "data/Step03-edgeR_nrDEG.Rdata")
edgeR_DEG <- nrDEG
dat <- exprSet
dat[1:4,1:4]
table(group_list)
# 提取差异倍数
# 根据需要修改DEG的值
DEG <- limma_DEG
FC <- DEG$log2FoldChange
names(FC) <- rownames(DEG)
class(FC)
FC
# 排序差异倍数，提取前100和后100的基因名
DEG_200 <- c(names(head(sort(FC),100)),names(tail(sort(FC),100)))
head(DEG_200)
# 提取基因的归一化
dat <- t(scale(t(dat[DEG_200,])))
dat[1:4,1:4]
group <- data.frame(group=group_list)
rownames(group)=colnames(dat)
pheatmap(dat,show_colnames =F,show_rownames = F, cluster_cols = T,
annotation_col=group,
#filename = "figures/Step04-heatmap_DEG.png",
# 增加color
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
image.png
# 大于2的值赋值为2
dat[dat > 2] <- 2
# 低于-2的值赋值为-2
dat[dat < -2] <- -2
pheatmap(dat,show_colnames =F,show_rownames = F, cluster_cols = T,
annotation_col=group,
#filename = "figures/Step04-heatmap_DEG_abs_upto2.png",
# 增加color
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
把极值去掉
