本文详细介绍了使用limma、Glimma和edgeR进行RNA-seq数据分析的过程，包括读取featureCounts表达矩阵、预处理、构建DGEList对象、标准化、差异表达计算（如logFC、P值和热图展示）以及实验设计。
目录
5.5.1 pheatmap显示差异最大的50个基因（上调25个，下调25个）的热图
#安装 R 包
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("RNAseq123")
library(limma)
library(Glimma)
library(edgeR)
library(org.Mm.eg.db)
#设置工作目录,xxx 为表达矩阵数据目录
#setwd("xxx")
主要参考了以下文章
https://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html#%E6%91%98%E8%A6%811.读取数据
1.1读取featureCounts表达矩阵
确保已经生成并导入数据count3_sorted
#读取数据
#读取featureCounts表达矩阵
count3 <- read.table("counts_matrix3.csv" ,
header = TRUE,
row.names = 1)
#去掉1-5行
count3 <- count3[,-c(1:5)]
#将行名改成grpdata中sample列列名
colnames(count3) <- grpdata$sample
#按行名从小到大顺序排列
count3_sorted <- count3[order(rownames(count3)), ]
确保含有分组数据grpdata
1.2设置实验组和对照组
grpdata$condition <- factor(grpdata$condition, levels = c("Etoh_8h", "ABA_8h"))2.表达值的预处理
2.1过滤
#过滤，选择那些在至少2个样本中计数大于3的基因
x <- count3_sorted[rowSums(count3_sorted >=3) >=2,]
#rowSums(count3_sorted >= 3 每个基因，计算所有计数大于或等于 3 的样本的总和
#rowSums(count3_sorted >= 3) >= 2 创建一个逻辑向量，如果上一步计算的总和大于或等于2，则逻辑值为TRUE
#count3_sorted[]将 count3_sorted 数据框中的行筛选出来，只保留那些满足前面逻辑向量为 TRUE 的行。2.2构建DGEList对象
#构建DGEList对象
#确保已经生成并导入数据count3_sorted
#确保含有分组数据grpdata
dge <- DGEList(counts = x,
genes=row.names(x),
group=as.factor(grpdata$condition))
#备份，防止后续分析失败
dge2 <- dge3.DGEList对象探索
3.1 barplot展示DGEList每个样本reads计数总和
barplot(colSums(dge$counts),las=2,col = rep(rainbow(2),each=2))#colSums(dge$counts)：colSums 函数计算 dge$counts 矩阵中每一列的总和
#las=2 表示标签将平行于条形图的轴，即标签将水平显示，
#col = rep(rainbow(4), each=2)：col 参数指定条形的颜色。rainbow(2) 生成一个包含 2 种颜色的彩虹调色板。rep 函数将这个颜色向量复制扩展，使得每个颜色被重复两次（因为 each=2）。
4.计算标准化因子
标准化方法：“TMM" "TMMwsp" "RLE" "upperquartile" "none"，下面以TMM为例
4.1使用TMM方法计算样本的归一化因子
dge <- calcNormFactors(dge,method = "TMM")4.2plotMDS可视化样本之间的相似性和差异性
#可视化样本之间的相似性和差异性(按样本名称显示结果)
plotMDS(dge,
+ labels = grpdata$sample,
+ col =rainbow(2))
#可视化样本之间的相似性和差异性(按分组显示结果)
plotMDS(dge.tmm,
+ labels = grpdata$condition,
+ col =rainbow(2))
5.差异表达计算
5.1生成实验矩阵
design.matrix <- model.matrix(~ group, data = dge.$samples)生成设计矩阵如图，design.matrix 现在包含了与 group 变量相关的设计矩阵，每一行对应一个样本，每一列对应 group 变量的一个水平（例如，对照组和处理组）。对照组元素为0，实验组元素为1。这个设计矩阵可以用于多种统计模型，包括线性模型（lm）、广义线性模型（glm）以及 DESeq2 中的差异表达分析。
5.2估计离散系数（差异表达分析关键步骤）
执行 estimateDisp 函数后，dge 对象将被更新，包含了每个基因的方差估计值。这些方差估计值将用于后续的差异表达分析。
dge <- estimateDisp(dge,design = design.matrix) 5.2.1 BCV查看离散系数
plotBCV(dge, cex = 0.8)
参数 cex = 0.8 用于调整图中点的大小。
- 如果某些基因的 BCV 远大于 1，这可能表明这些基因的表达量在样本间的变异性较高。
- 果某些基因的 BCV 接近 0，这可能表明这些基因的表达量在样本间几乎没有变异性。
- 图中可能会出现一些远离趋势线的点，这些点可能代表异常值或异常表达的基因。异常值可能需要进一步调查，因为它们可能代表了重要的生物学现象，或者是数据处理中的错误。
5.2.2plotMeanVar查看离散系数
plotMeanVar(dge, show.raw=TRUE,
show.tagwise=TRUE,
show.binned=TRUE)
- x轴通常表示基因表达量的均值（或某种形式的表达量估计），而y轴表示方差（或方差的估计）。
- 图中可能会出现一些远离趋势线的点，这些点可能代表异常值或异常表达的基因
5.3 exacTest检验
5.3.1 提取logFC，PValue
#对deg精确检验，确定基因表达量在不同条件下是否显著差异，得到精确p值和foldchange
dge_res <- exactTest(dge)
5.3.2提取显著差异表达基因
#按logFC排序
dge_results <- as.data.frame(topTags(dge_res,
n=nrow(dge$counts),
sort.by = "logFC"))
#n=nrow(dge$counts) 表示提取与 dge$counts 矩阵中的基因数量相同的显著基因数量
#sort.by "logFC" 表示根据对数变化倍数（log fold change）对结果进行排序
5.3 最大似然法检测（方法2）
# glmFit 函数来拟合一个负二项分布的线性模型
fit <- glmFit(y = dge, design = design.matrix)
# glmLRT 函数执行似然比检验
lrt <- glmLRT(fit, coef=2)
#按logFC排序
lrt_results <- as.data.frame(topTags(lrt,
n=nrow(dge$counts),
sort.by = "logFC"))
5.4筛选差异表达基因
#将结果以表格形式导出
write.csv(x = dge_results,file = "dge_results.csv")
#筛选出主要差异表达基因，logFC>=2,并且q值<0.05
logFC_indices2 <- which(abs(dge_results$logFC) >= 1.5)
FDR_indices2 <- which(dge_results$FDR<0.05)
sig3_indices <- intersect(logFC_indices2,FDR_indices2)
sig3 <- dge_results[sig3_indices,]5.5差异表达基因热图
5.5.1 pheatmap显示总体差异基因热图（基因分为4簇）
#使用 cpm 函数计算每个基因的每百万转录本数
cpm <- cpm(dge$counts, log=FALSE)
#拷贝cpm
cpm2 <- cpm
#绘图
pheatmap(cpm,
scale = "row",
kmeans_k= 4, #先将基因聚为4类
clustering_method = "ward.D",
cellwidth = 40, #指定宽度
cellheight = 40, #指定高度
main = "title", #指定标题
frontsize_row = 8, #指定行标签（基因）大小
frontsize_col = 12, # 指定列标签（样本）大小
angle_col = 45, #指定行角度
labels_col = paste(grpdata$condition," : ",grpdata$sample,sep = ""),
colorRampPalette(c("blue", "white", "red"))(50),
)
5.5.2 pheatmap显示差异最大的50个基因（上调25个，下调25个）的热图
#使用 cpm 函数计算每个基因的每百万转录本数
cpm <- cpm(dge$counts, log=FALSE)
#拷贝cpm
cpm2 <- cpm
#对sig3按logFC降序排序
sig3 <- sig3[order(sig3$logFC, decreasing = TRUE),]
#创建geneid向量，提取sig3的行名储存其中
geneid <- rownames(sig3)
#创建一个新的基因ID向量 sig50_geneid，它包含 geneid向量中的前25个基因ID和最后25个基因ID。
sig50_geneid <- c(geneid[1:25],geneid[(length(geneid)-24):length(geneid)])
#将sig50_geneid转化为factor类型，并作为行索引选择对应的行（基因名称）
cpm2 <- cpm2[as.factor(sig50_geneid),]
#绘图
pheatmap(cpm2,
scale = "row",
#kmeans_k= 4, #先将基因聚为4类,和下一项不兼容，完整绘图前面加#
clustering_method = "ward.D", # 指定聚类方法"wadr.D"，可换
cellwidth = 8, #指定宽度
cellheight = 40, #指定高度
main = "title", #指定标题
frontsize_row = 6, #指定行标签（基因）大小
frontsize_col = 12, # 指定列标签（样本）大小
angle_col = 45, #指定行角度
labels_col = paste(grpdata$condition," : ",grpdata$sample,sep = ""),
colorRampPalette(c("blue", "white", "red"))(50),
)
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
