介紹

常用的差異基因分析軟件主要有DESeq2汁果、edgeR以及Limma驾中。其中，DESeq2適合有重復(fù)的樣本(官方推薦4個(gè)以上)分衫，edgeR可以實(shí)現(xiàn)單個(gè)樣本的差異基因分析蔫仙。但兩者需要輸入的均為原始的read_counts矩陣，并需要gene length信息丐箩，因此只能在同一套參考基因組下進(jìn)行比較摇邦。而Limma是其中唯一支持tpm矩陣進(jìn)行差異基因計(jì)算的，因此Limma可以完成跨物種的差異基因篩選屎勘。

軟件安裝

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

我這里是報(bào)錯(cuò)的施籍，提示版本不合適。在這里我進(jìn)行了手動(dòng)安裝概漱。

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

提示缺少依賴的包"statmod"丑慎，缺啥裝啥

library(statmod)

statmod比較順利，接下來(lái)再次安裝時(shí)提示缺少'make'瓤摧。
這就比較麻煩竿裂，因?yàn)闆](méi)有make這個(gè)包，它時(shí)Rtools下的工具照弥。需要安裝Rtools腻异。
這里可以參考一篇大神的方法Rtools安裝方法
正常在線安裝

# 檢查有沒(méi)有'make'命令
Sys.which("make")
                              make 
""
# 表示命令不存在
# 安裝Rtools，需先安裝前兩個(gè)包
library(installr)
library(stringr)
install.Rtools()

正常安裝還是一如既往失敗了这揣。
那么從官網(wǎng)下載包Rtools,直接下載Rtools一路點(diǎn)擊下一步安裝(安裝到任意位置悔常，除了R本身的文件夾以外)影斑，就可以了。這個(gè)Rtools和R版本要適配机打，如果不適配需要?jiǎng)h除后重新安裝矫户。
安裝好之后，重新打開(kāi)R残邀，從新檢測(cè)make是否存在皆辽。

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

重新通過(guò)本地安裝limma后成功。

Limma的使用

首先是讀入數(shù)據(jù)芥挣，將數(shù)據(jù)整理成便于后續(xù)分析的格式

# 設(shè)置路徑
getwd()
setwd("C:/Users/1/Documents/R/LNN/TSE_PS/")
# 讀入數(shù)據(jù)
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")
# 篩選得到對(duì)應(yīng)樣本的tpm驱闷，以及對(duì)應(yīng)的樣本信息表
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)

展示一下處理好的用以分析的原始數(shù)據(jù)

> 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

設(shè)置分組信息

#### limma計(jì)算差異基因
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轉(zhuǎn)換，并將log轉(zhuǎn)換后無(wú)窮值轉(zhuǎn)換為0
data_tpm2 <- data_tpm1[which(rowSums(data_tpm1)!=0),] 
 # 這里可以log粘我，也可以使用自帶的voom函數(shù)進(jìn)行歸一化
#expr_data = log2(data_tpm2)
#expr_data[expr_data == -Inf] = 0
expr_data <- voom(data_tpm2,design,plot = F)
# 設(shè)置分組信息
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中只管設(shè)置分組即可鼓蜒，countrast.matrix中設(shè)置“treat vs control”，即如果調(diào)換位置征字，改為PS-TSE都弹，就成了PSvsTSE，這點(diǎn)較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

差異基因計(jì)算

這個(gè)DESeq2等流程類似畅厢，只要前面的分組和tpm設(shè)置正確就沒(méi)問(wèn)題

#
fit <- lmFit(expr_data,design) #非線性最小二乘法
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)#用經(jīng)驗(yàn)貝葉斯調(diào)整t-test中方差的部分
DEG <- topTable(fit2, coef = 1,n = Inf,sort.by="logFC")
DEG <- na.omit(DEG)
#

數(shù)據(jù)整理

得到的DEG結(jié)果和其它兩個(gè)軟件類似，結(jié)果如下

> 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

接下來(lái)需要按照自己的需求整理表格氮昧。無(wú)外乎刪掉結(jié)果中不需要的列框杜，增加上調(diào)、下調(diào)標(biāo)識(shí)的列袖肥，聯(lián)合表達(dá)量矩陣咪辱，聯(lián)合注釋信息等等。

# 刪除不需要的列椎组，修改剩余列油狂，名
DEGs_data <- DEG%>%
  mutate(t = NULL,AveExpr=NULL,B=NULL)%>%
  dplyr::rename(P.adj = adj.P.Val)%>%
# 以P.adj < 0.05為標(biāo)準(zhǔn)，可調(diào)寸癌。增加上調(diào)专筷、下調(diào)標(biāo)識(shí)列
  mutate(Direction = if_else(P.adj > 0.05, "NS",
                             if_else(logFC > 1,"UP",
                                     if_else(logFC < -1, "DOWN","NS"))))%>%
# 列名統(tǒng)一為Gene_ID
  rownames_to_column(var = "Gene_ID")%>%
# 加入表達(dá)量信息，列名統(tǒng)一為Gene_ID
  left_join(data_tpm1%>%rownames_to_column(var = "Gene_ID"))%>%
# 加入注釋文件
  left_join(TSE_KEGG_annotation)
# 查看差異基因數(shù)目
DEGs_stat <- DEGs_data%>%
group_by(Direction)%>%
  summarise(gene_number = n())
DEGs_stat
# 寫(xiě)出差異基因集
write.csv(DEGs_data,"Limma_DEGs_TSE_vs_PS_Female.csv")
# End

Limma差異基因篩選就結(jié)束了蒸苇。
上述時(shí)Female組的整個(gè)流程磷蛹，接下來(lái)時(shí)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)#用經(jīng)驗(yàn)貝葉斯調(diào)整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")