第一救拉，同樣的，清除變量+加載之前的數據
> rm(list = ls())
> load(file = "step2output.Rdata")
> group_list
[1] control control control treat treat treat
Levels: control treat
差異分析忱嘹，用limma包來做
需要表達矩陣exp和group_list分組對象映凳，不需要改
{
> library(limma)#加載limma包
#需要自己做好三個數據（表達矩陣，分組矩陣畏鼓，差異比較矩陣）
#總共三個步驟（lmFit,eBayes,topTable）
+ design=model.matrix(~group_list)#把group_list設置成一個model.matrix
#> design
# (Intercept) group_listtreat
#1 1 0
#2 1 0
#3 1 0
#4 1 1
#5 1 1
#6 1 1
#attr(,"assign")
#[1] 0 1
#attr(,"contrasts")
#attr(,"contrasts")$group_list
#[1] "contr.treatment"
+ fit=lmFit(exp,design)#lmFit用于線性擬合酱酬；至少兩個輸入，一個是表達矩陣云矫，一個是分組對象膳沽。
#表達矩陣必須是matrix類數據結構，每一列都是存放一個樣本让禀，每一行是一個探針信息或者是注釋后的基因名
+ fit=eBayes(fit)#根據lmFit的擬合結果進行統(tǒng)計推斷
+ deg=topTable(fit,coef=2,number = Inf)#得出差異比較矩陣
}
#view(deg)
差異表達矩陣deg.png
image.png
image.png
為deg數據框添加幾列
{
#1.加probe_id列：把deg的行名（探針名）賦值給probe_id挑社，并加一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))#mutate增加一列
head(deg)
#2.加symbol列，火山圖要用
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
#按照symbol列去重復#因為存在一個探針映射多個基因的情況
deg <- deg[!duplicated(deg$symbol),]
}
#3.加change列,下面兩行是閾值巡揍，可修改
logFC_t=1
logP_t = 0.01
test1 = deg$P.Value < logP_t
test2 = deg$logFC < -logFC_t
test3 = deg$logFC > logFC_t
{
change = ifelse(test1 & test2 ,
"down" ,
ifelse(test1 & test3 ,
"up",
"stable"))#不滿足down和up時痛阻，就是stable;分成三種
deg <- mutate(deg,change)#加一列有down和up和stable
#4.加ENTREZID列，用于富集分析（symbol轉entrezid腮敌，然后inner_join）
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人類基因ID轉換庫
#其他物種http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
}
save(group_list,deg,file = "step4output.Rdata")
image.png
