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

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

R语言单细胞求差异基因 r语言差异基因筛选_科技

首先，打开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')

R语言单细胞求差异基因 r语言差异基因筛选_R语言单细胞求差异基因_02

准备工作已经完成了，接下来进行的就是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矩阵就是我们的最终结果，大致预览效果如下所示。

R语言单细胞求差异基因 r语言差异基因筛选_数据_03

行名是基因名，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))    #画辅助线

R语言单细胞求差异基因 r语言差异基因筛选_科技_04

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

这一步的统计要在将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()    #加个风格，个人习惯这款

R语言单细胞求差异基因 r语言差异基因筛选_科技_05

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

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