介绍

常用的差异基因分析软件主要有DESeq2、edgeR以及Limma。其中，DESeq2适合有重复的样本(官方推荐4个以上)，edgeR可以实现单个样本的差异基因分析。但两者需要输入的均为原始的read_counts矩阵，并需要gene length信息，因此只能在同一套参考基因组下进行比较。而Limma是其中唯一支持tpm矩阵进行差异基因计算的，因此Limma可以完成跨物种的差异基因筛选。

软件安装

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("limma")

我这里是报错的，提示版本不合适。在这里我进行了手动安装。

install.packages("~/R/LNN/TSE_PS/limma_3.58.1.tar.gz", repos = NULL, type = "source")

提示缺少依赖的包"statmod"，缺啥装啥

library(statmod)

statmod比较顺利，接下来再次安装时提示缺少'make'。
这就比较麻烦，因为没有make这个包，它时Rtools下的工具。需要安装Rtools。
这里可以参考一篇大神的方法Rtools安装方法
正常在线安装

# 检查有没有'make'命令
Sys.which("make")
                              make 
""
# 表示命令不存在
# 安装Rtools，需先安装前两个包
library(installr)
library(stringr)
install.Rtools()

正常安装还是一如既往失败了。
那么从官网下载包Rtools,直接下载Rtools一路点击下一步安装(安装到任意位置，除了R本身的文件夹以外)，就可以了。这个Rtools和R版本要适配，如果不适配需要删除后重新安装。
安装好之后，重新打开R，从新检测make是否存在。

> Sys.which("make")
                              make 
"D:\\rtools43\\usr\\bin\\make.exe"
# 表示make命令存在

重新通过本地安装limma后成功。

Limma的使用

首先是读入数据，将数据整理成便于后续分析的格式

# 设置路径
getwd()
setwd("C:/Users/1/Documents/R/LNN/TSE_PS/")
# 读入数据
library(tidyverse)
#
tpm <- read.csv("~/R/LNN/TSE_PS/PS_TSE_tpm.csv")
# 建立样样本信息表
list <- tpm%>%colnames()
write.csv(list,"list.csv")
sample_list <- read.csv("~/R/LNN/TSE_PS/sample_list.csv")
# 筛选得到自己需要的样本信息
target_list <- sample_list%>%filter(tissue == "Inner_ear")
# 筛选得到对应样本的tpm，以及对应的样本信息表
data_tpm <- tpm%>%select(TSE,target_list$sample)
data_list <- target_list%>%select(sample,species,gender)
# 读入注释信息
TSE_KEGG_annotation <- read.delim("~/R/LNN/TSE_PS/TSE_KEGG_annotation.txt", header=FALSE)%>%
  dplyr::rename(Gene_ID = V1,Gene_name = V2)

展示一下处理好的用以分析的原始数据

> head(data_tpm)
           TSE      PS_FI1      PS_FI2      PS_FI3     PS_MI1     PS_MI2      PS_MI3     TSE_FI1
1         KIF6   0.0000000   0.4112962   0.1405429   1.209058   2.457624   1.0145210   6.7368448
2         GOT2  76.6149964 157.7452471 193.3648708 186.903353 396.991380 154.7511390  34.3573896
3 LOC117870645 197.8900227 181.9492006 173.3878021 171.557738 151.897520 190.8632881   1.0875280
4 LOC117870646   0.2811523   0.5090459   0.9132099   2.525185  55.511191   0.6726615   0.2330417
5       DNMT3B   1.1336946   2.9204753   3.6900862   3.968118   3.280471   3.6138520   2.8685177
6       MAPRE1  26.8055168 120.9675223 116.8502857 132.910083 150.910039 113.1194543 135.0258033
     TSE_FI2    TSE_FI3    TSE_MI1     TSE_MI2    TSE_MI3
1   8.313084  11.926323   9.292628   8.8822036  16.877889
2  79.602145  74.369763  53.897337  75.5693887  85.030344
3  11.588775   5.631727   4.972733   6.7368331   6.483457
4   0.000000   0.000000   0.000000   0.7022953   3.704832
5   3.693706   4.203426   5.258760   2.7606219   3.781259
6 177.052411 191.303851 189.669810 112.6413581 108.236785
> head(data_list)
  sample species gender
1 PS_FI1      PS Female
2 PS_FI2      PS Female
3 PS_FI3      PS Female
4 PS_MI1      PS   Male
5 PS_MI2      PS   Male
6 PS_MI3      PS   Male

设置分组信息

#### limma计算差异基因
library(limma)
#### Female
# 筛选样本信息表
expr_list <- data_list%>%
  filter(gender == "Female")
# 筛选tpm样本，并以gene id为行名
data_tpm1 <- data_tpm%>%
  column_to_rownames(var = "TSE")%>%
  select(expr_list$sample)
# 去除0值，log转换，并将log转换后无穷值转换为0
data_tpm2 <- data_tpm1[which(rowSums(data_tpm1)!=0),] 
 # 这里可以log，也可以使用自带的voom函数进行归一化
#expr_data = log2(data_tpm2)
#expr_data[expr_data == -Inf] = 0
expr_data <- voom(data_tpm2,design,plot = F)
# 设置分组信息
group <- data_list%>%
  filter(gender == "Female")%>%
  column_to_rownames(var = "sample")%>%
  select(species)
#coldata <- data.frame(group = factor(rep(c("PS","TSE"), each = 3)))
design <- model.matrix(~0+ factor(group$species))
colnames(design) <- levels(factor(group$species))
rownames(design) <- colnames(expr_data)
#
contrast.matrix <- makeContrasts(TSE-PS,levels = design)

最终得到的分组信息如下：
这里design中只管设置分组即可，countrast.matrix中设置“treat vs control”，即如果调换位置，改为PS-TSE，就成了PSvsTSE，这点较DESeq2中更加人性化。

> head(design)
        PS TSE
PS_FI1   1   0
PS_FI2   1   0
PS_FI3   1   0
TSE_FI1  0   1
TSE_FI2  0   1
TSE_FI3  0   1
> head(contrast.matrix)
      Contrasts
Levels TSE - PS
   PS        -1
   TSE        1

差异基因计算

这个DESeq2等流程类似，只要前面的分组和tpm设置正确就没问题

#
fit <- lmFit(expr_data,design) #非线性最小二乘法
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)#用经验贝叶斯调整t-test中方差的部分
DEG <- topTable(fit2, coef = 1,n = Inf,sort.by="logFC")
DEG <- na.omit(DEG)
#

数据整理

得到的DEG结果和其它两个软件类似，结果如下

> head(DEG)
                 logFC  AveExpr         t      P.Value    adj.P.Val
LOC117871706 -12.23683 6.118416 -32.77278 2.265359e-08 0.0002418326
LOC117876280 -11.97748 2.130500 -29.89722 4.063730e-08 0.0002418326
RPS29         10.91706 5.458532  24.47780 1.447263e-07 0.0003445064
PDIA2        -10.71380 4.612360 -15.83034 2.259883e-06 0.0008560078
PPDPFL       -10.62249 7.106056 -24.51989 1.431585e-07 0.0003445064
ANXA10       -10.32320 4.555014 -17.05096 1.418560e-06 0.0008560078
                    B
LOC117871706 7.295384
LOC117876280 7.114594
RPS29        6.632612
PDIA2        5.117458
PPDPFL       6.637300
ANXA10       5.421359

接下来需要按照自己的需求整理表格。无外乎删掉结果中不需要的列，增加上调、下调标识的列，联合表达量矩阵，联合注释信息等等。

# 删除不需要的列，修改剩余列，名
DEGs_data <- DEG%>%
  mutate(t = NULL,AveExpr=NULL,B=NULL)%>%
  dplyr::rename(P.adj = adj.P.Val)%>%
# 以P.adj < 0.05为标准，可调。增加上调、下调标识列
  mutate(Direction = if_else(P.adj > 0.05, "NS",
                             if_else(logFC > 1,"UP",
                                     if_else(logFC < -1, "DOWN","NS"))))%>%
# 列名统一为Gene_ID
  rownames_to_column(var = "Gene_ID")%>%
# 加入表达量信息，列名统一为Gene_ID
  left_join(data_tpm1%>%rownames_to_column(var = "Gene_ID"))%>%
# 加入注释文件
  left_join(TSE_KEGG_annotation)
# 查看差异基因数目
DEGs_stat <- DEGs_data%>%
group_by(Direction)%>%
  summarise(gene_number = n())
DEGs_stat
# 写出差异基因集
write.csv(DEGs_data,"Limma_DEGs_TSE_vs_PS_Female.csv")
# End

Limma差异基因筛选就结束了。
上述时Female组的整个流程，接下来时Male组的完整脚本。

#### Male
#### Male
expr_list <- data_list%>%
  filter(gender == "Male")
data_tpm1 <- data_tpm%>%
  column_to_rownames(var = "TSE")%>%
  select(expr_list$sample)
data_tpm2 <- data_tpm1[which(rowSums(data_tpm1)!=0),] 
#expr_data = log2(data_tpm2)
#expr_data[expr_data == -Inf] = 0
expr_data <- voom(data_tpm2,design,plot = F)
#
group <- data_list%>%
  filter(gender == "Male")%>%
  column_to_rownames(var = "sample")%>%
  select(species)
#coldata <- data.frame(group = factor(rep(c("PS","TSE"), each = 3)))
design <- model.matrix(~0+ factor(group$species))
colnames(design) <- levels(factor(group$species))
rownames(design) <- colnames(expr_data)
#
contrast.matrix <- makeContrasts(TSE-PS,levels = design)
#
fit <- lmFit(expr_data,design) #非线性最小二乘法
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)#用经验贝叶斯调整t-test中方差的部分
DEG <- topTable(fit2, coef = 1,n = Inf,sort.by="logFC")
DEG <- na.omit(DEG)
#
DEGs_data <- DEG%>%
  mutate(t = NULL,AveExpr=NULL,B=NULL)%>%
  dplyr::rename(P.adj = adj.P.Val)%>%
  mutate(Direction = if_else(P.adj > 0.05, "NS",
                             if_else(logFC > 1,"UP",
                                     if_else(logFC < -1, "DOWN","NS"))))%>%
  rownames_to_column(var = "Gene_ID")%>%
  left_join(data_tpm1%>%rownames_to_column(var = "Gene_ID"))%>%
  left_join(TSE_KEGG_annotation)
DEGs_stat <- DEGs_data%>%
  group_by(Direction)%>%
  summarise(gene_number = n())
DEGs_stat
write.csv(DEGs_data,"Limma_DEGs_TSE_vs_PS_Male.csv")