这20个题目其实是从用芯片数据做生信分析的初始部分流程拆解而成，里面涉及到的技术内容可能已经过时，比如现在已经用测序数据代替了芯片数据，但生信分析背后的逻辑其实还是类似。
- 对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ， 然后画所有样本的 这些图
- 参考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
- 理解ggplot2的绘图语法，数据和图形元素的映射关系
ggplot2的绘图基本套路：数据+映射aes+几何对象geom_+boxplot、violin、his、density等+图表排列facet+统计变换stat_summary
# 将原始表达矩阵转换为长矩阵
library(reshape2)
exprSet_L=melt(exprSet)
head(exprSet_L)
colnames(exprSet_L)=c('probe','sample','value')
# 根据矩阵转换特点，原矩阵每个列转换成了nrow(exprSet)行，所以每列的group重复赋值nrow(exprSet)次
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
# probe sample value group
#1 MAPK3 CLL11.CEL 5.743132 progres.
#2 TIE1 CLL11.CEL 2.285143 progres.
#3 CYP2C19 CLL11.CEL 3.309294 progres.
#4 CXCR5 CLL11.CEL 1.085264 progres.
#5 CXCR5 CLL11.CEL 7.544884 progres.
#6 DUSP1 CLL11.CEL 5.083793 progres.
# 经过处理，experSet为原始宽矩阵，experSet_L为转换后的长矩阵，一行一个value# 使用ggplot2
library(ggplot2)
# boxplot，数据映射+几何对象
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
# violinplot，数据映射+几何对象
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
# 直方图，数据映射+几何对象+多面板图
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
# 密度分布图，数据映射+几何对象+多面板图
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
# 密度分布图，数据映射+几何对象
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print(p)
# 组合；箱式图+统计点图+网格白色主题+主题设置（文字加粗，转角，空标题）
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)# 拓展， 使用gpairs
library(gpairs)
gpairs(exprSet)
# 聚类直方图
out.dist=dist(t(exprSet),method='euclidean')
out.hclust=hclust(out.dist,method='complete')
plot(out.hclust)
# PCA降唯
pc <- prcomp(t(exprSet),scale=TRUE)
pcx=data.frame(pc$x)
pcr=cbind(samples=rownames(pcx),group_list, pcx)
p=ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)
print(p)- 理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。
- 注意：这个题目出的并不合规，请仔细看。
# 平均值
# 在exprSet按行方向做mean计算（按基因计算），按此排序（默认升序），取最后50个值（即top max 50）。下同。
g_mean <- tail(sort(apply(exprSet,1,mean)),50)
# 中位值
g_median <- tail(sort(apply(exprSet,1,median)),50)
# 最大值
g_max <- tail(sort(apply(exprSet,1,max)),50)
# 最小值
g_min <- tail(sort(apply(exprSet,1,min)),50)
# standard deviation标准差
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
# variance方差
g_var <- tail(sort(apply(exprSet,1,var)),50)
# median absolute deviation，绝对中位差，与中位数偏差值的中位数，描述离散趋势的度量值
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
names(g_mad)
- 根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。
library(pheatmap)
# choose_gene=names(sort(apply(exprSet, 1, mad),decreasing = T)[1:50])
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
# scale()函数用于去中心化和标准化.对每个探针的表达量进行去中心化和标准化
# scale是按column做标准化处理
choose_matrix=scale(choose_matrix)
# heatmap(choose_matrix)
pheatmap(choose_matrix)# 网例美化热图
annotation_col = data.frame(group_list)
rownames(annotation_col)=colnames(exprSet)
colnames(annotation_col)[1]<-"group"
ann_colors = list(group = c(case = "#1B9E77",control = "#D95F02"))
pheatmap(choose_matrix, #pheatmap包简单作图中心化和标准化的数据scale = "row",#参数对行进行归一化color = colorRampPalette(colors = c('#11427C','white','#C31E1F'))(100),#设置热图的颜色cluster_rows = F,cluster_cols = T,#聚类，有些期刊你会发现行和列都进行聚类，或者是只有行或列进行聚类，通过这两个参数控制annotation_col = annotation_col,gaps_col = 8,#对行列进行划分为不同的块, 只有不聚类的情况下才能进行分组，8 代表你想在哪列进行分块show_colnames = T,show_rownames = T,# annotation_legend = FALSE #命令去掉注释图例
)- 取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。
该图由三个子图组成： - 上方：表示交集大小的柱状图 - 下左：表示集合大小的条形图 - 下右：表示集合之间的交叠矩阵，矩阵的列表示每种交集组合，对应于柱状图的横坐标；矩阵的行表示集合，对应于条形图的纵坐标 - 通过这样一张图，可以展示多个集合之间的交叠关系，且很容易从图中看出集合之间的交集信息 - https://zhuanlan.zhihu.com/p/370210775
## UpSetR
# https://cran.r-project.org/web/packages/UpSetR/README.html
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,g_mean=ifelse(g_all %in% names(g_mean) ,1,0),g_median=ifelse(g_all %in% names(g_median) ,1,0),g_max=ifelse(g_all %in% names(g_max) ,1,0),g_min=ifelse(g_all %in% names(g_min) ,1,0),g_sd=ifelse(g_all %in% names(g_sd) ,1,0),g_var=ifelse(g_all %in% names(g_var) ,1,0),g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
- 在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
# [1] "progres." "stable" "progres." "progres." "progres." "progres." "stable" "stable" "progres." "stable" "progres." "stable" "progres." "stable"
# [15] "stable" "progres." "progres." "progres." "progres." "progres." "progres." "stable"
dim(exprSet)
# [1] 8776 22
exprSet[1:5,1:5]
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# MAPK3 5.743132 6.219412 5.523328 5.340477 5.229904
# TIE1 2.285143 2.291229 2.287986 2.295313 2.662170
# CYP2C19 3.309294 3.318466 3.354423 3.327130 3.365113
# CXCR5 7.544884 7.671801 7.474025 7.152482 6.902932
# DUSP1 5.083793 7.610593 7.631311 6.518594 5.059087- 对所有样本的表达矩阵进行聚类并且绘图，然后添加样本的临床表型数据信息(更改样本名)
## hclust
colnames(exprSet)=paste(group_list,1:22,sep='')
# Define nodePar
nodePar <- list(lab.cex = 0.6, # 树上标签文字的缩放比例（0.6倍）pch = c(NA, 19), # 节点符号：第一个 NA 对应分支节点不画点；19 是叶节点画实心圆cex = 0.7, # 节点符号的大小（0.7倍）col = "blue")# 节点符号和标签文字的颜色
# 聚类函数要求行为样本列表，列为基因列表
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10)) #调整绘图边距，四个方向的边距行数（单位是文本行高度），右边距设大一些为了横向绘图时有足够空间显示叶节点标签。
plot(as.dendrogram(hc), # 把 hclust 对象转成 dendrogram，才能接受 nodeParnodePar = nodePar, # 上面定义的节点/标签参数horiz = TRUE # 横向绘制：根在左，叶在右
)
- 对所有样本的表达矩阵进行PCA分析并且绘图，同样要添加表型信息。
-
ggfortify包已经很长时间没有更新维护了，其调用ggplot2包的函数aes_string()已被弃用，原代码未更新。
-
factoextra包也是很久没维护，其fviz_pca_ind函数设置图例标签时，遇到了它无法识别的linetype美学属性，原代码也是一时半会儿不会更新。
-
这个问题在回溯学习的时候经常会遇到，咨询DS老师后，最好的解决办法还是手搓实现，使用系统stats原生计算PCA，和用ggplot2绘图，完全自主可控，也能及时升级代码。核心步骤包括：
- 提取PCA结果：从 prcomp对象中获取样本在主成分上的坐标（得分）。
- 计算方差贡献：了解每个主成分的重要性。
- 使用 ggplot2绘图：将坐标数据转换为数据框，并用 ggplot进行可视化。
-
参考来源： https://www.cnblogs.com/ljbguanli/p/18958694 和 https://algo.ac.cn/archives/shen-ru-qian-chu-zhu-cheng-fen-fen-xi-pca-cong-yuan-li-dao-tui-dao-de-wan-quan-zhi-nan
- PCA是最最重要的降维思想，在研究中不起眼，在工业中可是百试不爽，后续可根据上述参考内容整理PCA单篇。
# PCA方法1：ggfortify
# 两种安装方式都可以
# install.packages("ggfortify")
# BiocManager::install('ggfortify')
library(ggfortify)
# 这个作业要求所有样本降维，所以重新取值
library(CLL)
data(sCLLex)
exprSet <- exprs(sCLLex)
# PCA函数要求行为样本列表，列为基因列表
df=as.data.frame(t(exprSet))
pData = pData(sCLLex)
group_list = as.character(pData[,2])
df$group=group_list
# autoplot uses ggplot2 to draw a particular plot for an object of a particular class in a single command.
autoplot(prcomp(df[, 1:(ncol(df)-1)]), # 只用数值型列做 PCA，不包含最后一列 groupdata = df, # 指定映射来源colour = 'group' # 按 group 列给样本点着色
)
# ggfortify包调用已经被弃用的ggplot的aes_string()函数，一时半会儿还不会更新。还是手搓上价。# PCA方法2
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
df=as.data.frame(t(exprSet))
dat.pca <- PCA(df, graph = FALSE)#这个矩阵是不含有分组信息，重新赋值df用于计算dat.pca
fviz_pca_ind(dat.pca,geom.ind = "point", # show points only (nbut not "text")habillage = as.factor(group_list), # 使用habillage而不是col.ind# palette = c("#00AFBB", "#E7B800"),addEllipses = TRUE, # Concentration ellipseslegend.title = "Groups"
)
# factoextra包的fviz_pca_ind函数在尝试设置图例标签时，遇到了它无法识别的美学属性（这里是 linetype），也是一时半会儿不会更新。得手搓。# PCA方法3:手搓绘图
# 这个作业要求所有样本降维，所以重新取值
exprSet <- exprs(sCLLex)
# PCA函数要求行为样本列表，列为基因列表
df=as.data.frame(t(exprSet))
# 1. 执行PCA，提取主成分得分，转换为数据框
pca_result <- prcomp(df)
pca_scores <- as.data.frame(pca_result$x)
# 将分组信息添加到数据框中
pData = pData(sCLLex)
pca_scores$Group <- as.character(pData[,2])
# 2. 计算主成分的方差贡献比例
percent_var <- pca_result$sdev^2 / sum(pca_result$sdev^2) * 100
# 格式化显示，保留两位小数
percent_var <- round(percent_var, 2)
# 3. 使用ggplot2绘制PCA图
ggplot(pca_scores, aes(x = PC1, y = PC2, color = Group)) +geom_point(size = 3) + # 绘制点stat_ellipse(level = 0.95) + # 添加95%置信椭圆labs( # 设置标签，包含方差贡献x = paste0("PC1 (", percent_var[1], "%)"),y = paste0("PC2 (", percent_var[2], "%)"),color = "Group",# title = "PCA Plot of Iris Dataset") +theme_bw() + # 使用黑白主题coord_fixed(ratio = 1) # 保持X轴和Y轴比例一致- 根据表达矩阵及样本分组信息进行批量T检验，得到检验结果表格
- p值
本数据集总共样本数偏少，22个样本分成2组。所以要比较两组样本之间的是否存在显著差异可以用t检验计算p-value及p-adjust。基本步骤如下。
假设H0：两组样本（progress组和stable组）间无显著差异。每个基因（每行）在两组样本有不同的表达量，对其进行t-test检验，得到p-value（p-adj）。p-value（p-adj）即原假设（或更极端情况）成立的概率，如果小于0.05（显著性水平），则怀疑原假设的正确性，认为原假设H0不成立，即该基因在两个组的表达差异明显。
另外，从工程上说，p.adj值很小，一般取对数来表达差异趋势，可以得到p.adj<0.05 <=> -log10p.adj>1.3 <=> 统计学显著。
- logFC值
logFC（log fold change）即log之后的差异倍数，计算方式：log_2(FC)=log_2(实验组表达量/对照组表达量) ，从计算值的内涵来说|log_2(FC)|>1 <=> 表达水平变化大于2倍
不过，本例中有个疑问，如何判断数据集是否已经过log变换。
- 显著性判断
只有当p.adj < 显著性阈值（一般0.05）且abs(log2FC) > 显著性阈值（一般1）,认为显著。
# t.test检验的组别需要时factor变量
group_list=as.factor(group_list)
# 1. p值计算
pvals = apply(exprSet, 1, function(x){t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")
# 2. 计算logFC值。
dat = exprSet
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
# 3. 将计算得到的p.adj和log2FC形成差异数据表
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
# avg_1 avg_2 log2FC pvals p.adj
# 36129_at 7.875615 8.791753 0.9161377 1.629755e-05 0.2057566
# 37676_at 6.622749 7.965007 1.3422581 4.058944e-05 0.2436177
# 33791_at 7.616197 5.786041 -1.8301554 6.965416e-05 0.2436177
# 39967_at 4.456446 2.152471 -2.3039752 8.993339e-05 0.2436177
# 34594_at 5.988866 7.058738 1.0698718 9.648226e-05 0.2436177
# 32198_at 4.157971 3.407405 -0.7505660 2.454557e-04 0.3516678
# 优化后的批量T检验函数
batch_t_test <- function(expr_matrix, group_vector, control_group = NULL, p_adjust_method = "BH",log2FC_threshold = 1, pval_threshold = 0.05) {# 参数检查与数据预处理if (!is.factor(group_vector)) {group_vector <- as.factor(group_vector)}group_levels <- levels(group_vector)if (length(group_levels) != 2) {stop("分组变量必须包含且仅包含两个水平")}# 确定对照组和处理组if (is.null(control_group)) {control_group <- group_levels[1]treatment_group <- group_levels[2]} else {if (!control_group %in% group_levels) {stop("指定的对照组不存在于分组变量中")}treatment_group <- group_levels[group_levels != control_group]}# 获取组别索引control_idx <- which(group_vector == control_group)treatment_idx <- which(group_vector == treatment_group)# 检查样本量是否足够if (length(control_idx) < 2 | length(treatment_idx) < 2) {stop("每组至少需要2个样本才能进行T检验")}# 提取组别数据control_data <- expr_matrix[, control_idx, drop = FALSE]treatment_data <- expr_matrix[, treatment_idx, drop = FALSE]# 计算基本统计量avg_control <- rowMeans(control_data, na.rm = TRUE)avg_treatment <- rowMeans(treatment_data, na.rm = TRUE)log2FC <- avg_treatment - avg_control# 批量进行T检验（带错误处理），法1# pvals <- sapply(1:nrow(expr_matrix), function(i) {# tryCatch({# control_vals <- as.numeric(control_data[i, ])# treatment_vals <- as.numeric(treatment_data[i, ])# # # 检查方差是否为零（避免T检验错误）# if (sd(control_vals) == 0 & sd(treatment_vals) == 0) {# return(1) # 两组方差都为零，返回p=1# }# # t.test(treatment_vals, control_vals, var.equal = FALSE)$p.value# }, error = function(e) {# warning(paste("基因", rownames(expr_matrix)[i], "T检验出错:", e$message))# return(NA) # 出错时返回NA# })# 批量进行T检验（带错误处理），法2pvals <- apply(expr_matrix, 1, function(x){t.test(as.numeric(x)~group_vector)$p.value})# P值校正p.adj <- p.adjust(pvals, method = p_adjust_method)# 生成完整的差异分析结果表格DEG_results <- data.frame(gene_id = rownames(expr_matrix),avg_control = avg_control,avg_treatment = avg_treatment,log2FC = log2FC,pvalue = pvals,p.adj = p.adj,significant = p.adj < pval_threshold & abs(log2FC) > log2FC_threshold,regulation = ifelse(p.adj < pval_threshold & abs(log2FC) > log2FC_threshold,ifelse(log2FC > 0, "up", "down"), "not significant"),stringsAsFactors = FALSE)# 按P值排序DEG_results <- DEG_results[order(DEG_results$pvalue), ]return(DEG_results)
}# 使用示例
# DEG_results <- batch_t_test(exprSet, group_list, control_group = "stable")
# head(DEG_results)# 结果统计摘要函数
summary_DEG <- function(DEG_df, pval_threshold = 0.05, log2FC_threshold = 1) {sig_up <- sum(DEG_df$p.adj < pval_threshold & DEG_df$log2FC > log2FC_threshold, na.rm = TRUE)sig_down <- sum(DEG_df$p.adj < pval_threshold & DEG_df$log2FC < -log2FC_threshold, na.rm = TRUE)total_genes <- nrow(DEG_df)cat("差异表达基因统计:\n")cat("总基因数:", total_genes, "\n")cat("上调基因数:", sig_up, "\n")cat("下调基因数:", sig_down, "\n")cat("显著基因数:", sig_up + sig_down, "\n")cat("显著比例:", round((sig_up + sig_down) / total_genes * 100, 2), "%\n")
}# 使用示例
summary_DEG(DEG_results)
- 使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图。)。
-
sCLLex是RNA-seq芯片数据，按主流主要使用limma包进行基因差异分析。基因差异分析主要是比较各基因在样本组间的表现差别，涉及指标主要包括log2FC、p-value、adjusted p-value等。
-
limma包是专门针对芯片数据开发的差异分析工具包，虽然现在已经是测序数据的天下，但分析思想一脉相承，且有limma-voom包用于转录组测序数据的差异分析。
-
limma包需要的输入文件：
- 表达矩阵 （exprSet）：芯片数据的ExpressionSet对象可以通过exprs()获得 ，而常规的转录组数据（如count数据）可以通过read.csv()等导入
- 分组矩阵 (design) ：就是将表达矩阵的列（各个样本）分成几组（如case - control，或者时间序列day0, day1, day2 ...），可通过model.matrix()得到
- 比较矩阵（contrast）：定义期望进行的组间比较，可通过makeContrasts(比较关系，levels=design)得到，比较关系定义了谁比较谁。
- 注意：
- 减法的意义：在 group2 - group1中，正值的logFC表示基因在group2中的表达量高于group1（即上调），负值则表示表达量低于group1（下调），logFC的正负号关乎对基因是“上调”还是“下调”的正确解读 。
- 组名的准确性：makeContrasts函数中使用的组名（如 group1）必须与 colnames(design)中设置的名称完全一致，包括大小写。
- 多组比较：当同时进行多个比较时，topTable函数默认会输出一个合并的F检验p值。若要获取特定两两比较的结果，需要通过 coef参数（=2或3）来指定 。
-
limma包分析的主要流程：
- lmFit()：线性拟合模型构建，需要两个数据：表达矩阵exprSet和分组矩阵design ，得到的结果再和比较矩阵contrast一起导入contrasts.fit()函数
- eBayes()：利用上一步contrasts.fit()的结果
- topTable()：利用上一步eBayes()的结果，并最终导出差异分析结果
-
参考来源：https://www.jianshu.com/p/6d18e39e9ba3 和 https://zhuanlan.zhihu.com/p/677665271
- 虽然现在很少用芯片数据了，但结合上述两篇博文可以整理一个limma包的使用集锦单篇。
# 通过实际案例比较分组矩阵有无截距的实际区别
suppressMessages(library(limma))
# ~0的作用为去掉截距（全1的列）
model.matrix( ~ factor(group_list) )
# 输出为：
# (Intercept) factor(group_list)stable
# 1 1 0
# 2 1 1
# 3 1 0
# 4 1 0
# 5 1 0
# 6 1 0
# 7 1 1
# 8 1 1
# 9 1 0
# 10 1 1
# 11 1 0
# 12 1 1
# 13 1 0
# 14 1 1
# 15 1 1
# 16 1 0
# 17 1 0
# 18 1 0
# 19 1 0
# 20 1 0
# 21 1 0
# 22 1 1
# attr(,"assign")
# [1] 0 1
# attr(,"contrasts")
# attr(,"contrasts")$`factor(group_list)`
# [1] "contr.treatment"
model.matrix( ~0+ factor(group_list) )
# 输出为
# factor(group_list)progres. factor(group_list)stable
# 1 1 0
# 2 0 1
# 3 1 0
# 4 1 0
# 5 1 0
# 6 1 0
# 7 0 1
# 8 0 1
# 9 1 0
# 10 0 1
# 11 1 0
# 12 0 1
# 13 1 0
# 14 0 1
# 15 0 1
# 16 1 0
# 17 1 0
# 18 1 0
# 19 1 0
# 20 1 0
# 21 1 0
# 22 0 1
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`factor(group_list)`
# [1] "contr.treatment"
# DEG by limma
suppressMessages(library(limma))
# step0. limma分析需要以下矩阵：表达矩阵exprSet，分组矩阵design，比较矩阵contrast.matrix
# 定义design分组矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
# 定义比较矩阵contrast.matrix，本例中要把progres.组跟stable组进行差异分析比较
# 本例方法可根据实际组名自动设置比较关系
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)##step1，非线形最小二乘法构建线性拟合模型，入参包括表达矩阵+分组矩阵
fit <- lmFit(exprSet,design)
##step2，t-test，入参为setp1结果+比较矩阵
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!用经验贝叶斯调整t-test中的方差
##eBayes() with trend=TRUE
##step3
tempOutput <- topTable(fit2, coef=1, n=Inf)
nrDEG <- na.omit(tempOutput)# write.csv(nrDEG,"limma_notrend.results.csv",quote = F)
head(nrDEG)
# logFC AveExpr t P.Value adj.P.Val B
# 39400_at -1.0284628 5.620700 -5.835799 8.340576e-06 0.03344118 3.233915
# 36131_at 0.9888221 9.954273 5.771526 9.667514e-06 0.03344118 3.116707
# 33791_at 1.8301554 6.950685 5.736161 1.048765e-05 0.03344118 3.051940
# 1303_at -1.3835699 4.463438 -5.731733 1.059523e-05 0.03344118 3.043816
# 36122_at 0.7801404 7.259612 5.141064 4.205709e-05 0.10619415 1.934581
# 36939_at 2.5471980 6.915045 5.038301 5.362353e-05 0.11283285 1.736846
# 使用limma进行差异分析后，可以得到一个DEG矩阵，包括各基因的logFC、p.adj等## 基因差异分析结果的图形化展示方式：火山图
DEG=nrDEG
# 取logFC绝对值的平均值和logFC标准差两倍之和为阈值，高出阈值定义为差异表达显著
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
# {if (pval<0.5&|logFC|>cutoff)--> [if(logFC > cutoff)--> up else down] else not changed}
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
head(DEG)
# ggplot2绘图
library(ggplot2)
# 绘图标题
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=logFC, y=-log10(P.Value), color=change)) + # 数据映射+美学aesgeom_point(alpha=0.4, size=1.75) + # 几何图形==点图theme_set(theme_bw(base_size=20))+ # bw==网格白色主题xlab("log2 fold change") + ylab("-log10 p-value") + # 横纵坐标值ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+ # 图标题+主题title字体设置scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)- 对T检验结果的P值和limma包差异分析的P值画散点图，看看哪些基因相差很大。
## 绘制散点图，观察两种方法得到的logFC、P值、log10(FC)值
# nrDEG是limma包的分析结果，DEG_t.test是t检验的结果
head(nrDEG)
head(DEG_t.test)
# 使两个结果的基因相同，都取nrDEG的，nrDEG有omit过所以应该少一点
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1]) ## 比较logFC：两种方法值正负相反，可猜测实验-对照对比关系是相反的
# 两种方式中，组别对比设置方式不同，后期需要探索函数细节来检验
plot(DEG_t.test[,4],nrDEG[,4]) # 比较pvals：趋势相同，差异尚可接受
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))# 这几个基因在所有样本中的表达量都相近
# DLEU1，淋巴细胞白血病缺失基因，deleted in lymphocytic leukemia
exprSet['GAPDH',]
exprSet['ACTB',]
exprSet['DLEU1',]## 用箱线图观察两组样本的DLEU1基因表达差异
library(ggplot2)
library(ggpubr)
my_comparisons <- list(c("stable", "progres.")
)
# 取DLEU1基因相关数据保存dat：group为progres或stable，基因名称ID，表达量值
dat=data.frame(group=group_list,sampleID= names(exprSet['DLEU1',]),values= as.numeric(exprSet['DLEU1',]))
# 数据=DLEU1表达量，按group分组（值+颜色）
ggboxplot(dat, x = "group", y = "values",color = "group",add = "jitter"
)+
# 添加显著性水平stat_compare_means(comparisons = my_comparisons, method = "t.test")## 用热图heatmap观察基因表达差异（分析）
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
# scale函数默认是按列标准化的
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
有编程基础的话，入门这个领域还是可以上手的，关键的是生物学、细胞学背景知识以及数学基础知识，再就是结合当下流程的机器学习、深度学习和大模型去拔高概念。细胞注释、PCA、差异分析等都可以做下一步的学习内容。