差异基因表达分析是一种常见的生信分析方法，是每个生信人都必须掌握的技术，本文将使用R语言演示如何利用limma包分析TCGA的RNA基因表达矩阵。

首先，准备好所需的数据，如下图所示，基因表达数据为一个包含样品与基因的矩阵。

首先，打开R之后先加载所需的R包。其中，limma是差异基因表达分析的一个常用R包，ggplot2和ggrepel是用来绘图的。

library(limma)

library(ggplot2)

library(ggrepel)

设定好工作目录后，读取基因表达矩阵。因为我在Xena上下载的基因表达矩阵的包含了肿瘤组织样本和癌旁组织样本，因此要区分这两种组织作为分组依据。可以根据TCGA的15位样本编号的最后两位进行区分，01是肿瘤组织样本，11是正常组织样本。

#读取基因表达矩阵

data<-read.table('BLCA_signature_n-t.tsv', row.names = 1, header = T)

#按肿瘤/正常组织分组

Tdata<-subset(data, substr(rownames(data),14,15)=='01')  #统计肿瘤样本数目

Ndata<-subset(data, substr(rownames(data),14,15)=='11')  #统计癌旁组织样本数目

data<-rbind(Tdata, Ndata)  #这一步是为了将上两步筛选出的肿瘤样本和正常样本按顺序排列好，方便后续分析

随后就是对样本进行分组了，我这里按照癌组织和正常组织进行分组，我们需要一个记录了分组信息的列表list，这里的list可以按照自己想要的分组进行设置，但是要和前面的data的排列顺序互相匹配。然后，根据这个list生成一个列名为分组名，只包含0和1的分组矩阵。分组矩阵预览效果如下。

#这里的list可以按照自己想要的分组进行设置，但是要和前面的data的排列顺序互相匹配

list<-c(rep('Tumor', nrow(Tdata)),rep('Normal',nrow(Ndata)))

#生成分组矩阵

design <- model.matrix(~ 0 + list)

colnames(design) <- c('Tumor', 'Normal')

准备工作已经完成了，接下来进行的就是limma的主体部分。注意进行lmFit时的基因表达矩阵的基因名要放到行名，不要搞错了。

在进行makeContrasts的时候，记得改好分组信息，要和前面的分组矩阵保持一致。

#limma

data<-t(data)  #最终矩阵的基因名在行名，记得检查一下不要搞错了

fit <- lmFit(data, design)

contrast.matrix <- makeContrasts(Tumor - Normal, levels = design)  #这里的“Tumor”和“Normal”是我的分组信息，记得和前面的分组保持一致

fit2 <- contrasts.fit(fit, contrast.matrix)

fit2 <- eBayes(fit2)

diff<-topTable(fit2, adjust="BH", sort.by="logFC", n=Inf)

diff矩阵就是我们的最终结果，大致预览效果如下所示。

行名是基因名，logFC（log2 fold change）是两组之间差异表达的倍数，使用log2处理过。AveExpr是基因在所有样本中的平均表达量，t是用于t-test的，可以衡量组间差异显著性，P.value就是P值，adj.P.Val是校正过的P值，这里我用的是“BH”方法进行的校正。B是表示基因表达差异的贝叶斯统计量。这里我们基本上只用到logFC、P.value和adj.P.Val，其它可以不用管。通常我们认为|logFC|>=1，P值<0.05就算是一个差异表达基因，当然，这个具体情况具体分析，不一定按照这个标准筛选。

之后就是做差异基因表达专属的火山图了。这里先把p值转换为负对数形式，再用ggplot就可以画出一幅很基本的图。

#转换p值为-lg坐标

diff$P.Value<-log(diff$P.Value, 0.1)

#ggplot

ggplot(diff, aes(logFC, P.Value))+geom_point(size=2.5)+  #使用点图

 scale_x_continuous(limits = c(-6, 6))+  #调整横坐标，使图形左右对称

 labs(y='-lg(p.value)')+  #修改纵坐标名称

 geom_vline(xintercept = c(1,-1))+geom_hline(yintercept = -log(0.05, 10))  #画辅助线

当然，这张图可以说是非常的丑陋，我们要给它加上标签和颜色。

这一步的统计要在将p值取负对数之前。

#颜色设置

color<-ifelse(diff$P.Value<=0.05 & abs(diff$logFC)>=1,

       ifelse(diff$logFC>0, 'up', 'down'), 'unsignificant')  #ifelse判断是否显著确定分组

#标签名字

label<-ifelse(diff$P.Value<=0.05 & abs(diff$logFC)>=1, rownames(diff), '')  #同理

将这些数据用于作图即可获得比较好看的火山图了。

#ggplot

ggplot(diff, aes(logFC, P.Value, colour=color))+geom_point(size=2.5)+

 scale_x_continuous(limits = c(-6, 6))+

 labs(y='-lg(p.value)')+

 scale_color_manual(values=c('blue','black','red'))+  #设置分组的颜色

 geom_text_repel(aes(logFC, P.Value, label=label))+  #在对应的点上显示基因名

 geom_vline(xintercept = c(1,-1))+geom_hline(yintercept = -log(0.05, 10))+

 theme_bw()  #加个风格，个人习惯这款

以上就是差异基因表达分析和作图的基本方法。

由于up临近期末，因此这次是近期的最后一次更新，希望大家谅解。

参考文献：https://blog.csdn.net/weixin_73043012/article/details/134840337?spm=1001.2014.3001.5501（其实是我本人）