前面的教程：把火山圖做成煙花圖 ，我們批評了一個學員做差異分析時候的錯誤的思路，就是把每個分組裡面的單個樣品簡單粗暴的複製粘貼3次充當生物學重複，會導致計算得到的p值如此的詭異。
exprSet=cbind(
exprSet[,1:2],
exprSet[,1:2],
exprSet[,1:2],
exprSet[,1:2]
)
colnames(exprSet)=1:8那麼， 假如你的兩個分組真的就是都有且僅有一個樣品，你該如何進行差異分析呢？
當然也是有辦法的，大家在前面的教程：把火山圖做成煙花圖 也看到了大量的留言給出了解決方案！我們就摘選其中最簡單的edgeR包來分享一下。
首先仍然是從查看表達量矩陣和分組信息這個時候需要參考前面的教程：把火山圖做成煙花圖 ，自己製作 Step01-airwayData.Rdata 這個文件，它裡面存儲了表達量矩陣和分組信息：
rm(list = ls())
options(stringsAsFactors = F)
# 讀取基因表達矩陣信息
lname <- load(file = "./data/Step01-airwayData.Rdata")
lname
# 查看分組信息和表達矩陣數據
exprSet <- filter_count
dim(exprSet)
exprSet[1:4,1:4]
group_list
table(group_list)
# 加載包
library(DESeq2)
exprSet = exprSet[,1:2]
group_list=group_list[1:2]可以看到的是，我們簡單粗暴的選取了這個表達量矩陣的前兩個樣品而已。這樣的話，我們的兩個分組就各自都有且僅有一個樣品啦，這樣的無生物學重複的實驗設計其實也可以勉強做差異分析。
使用 edgeR包 設置 bcv參數來做無重複的差異分析代碼很簡單， 如下所示：
#加載包
library(edgeR)
#設置分組信息
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.1#設置bcv為0.1
et <- exactTest(exprSet, dispersion=bcv^2)
DEG_edgeR=as.data.frame(topTags(et, n = nrow(exprSet$counts)))
head(DEG_edgeR)
write.csv(DEG_edgeR, 'single-sample-deg-result.csv', quote = FALSE)前面的教程：把火山圖做成煙花圖 裡面我們展示了 DEseq2和limma兩個包，這個時候又出來了 edgeR包 ，其實就大滿貫啦。轉錄組差異分析最常用的3個包，就是它們了。
仍然是繪製火山圖代碼如下所示：
#篩選上下調，設定閾值
fc_cutoff <- 1.5
pvalue <- 0.05
DEG_edgeR$regulated <- "normal"
loc_up <- intersect(which(DEG_edgeR$logFC>log2(fc_cutoff)),
which(DEG_edgeR$PValue<pvalue))
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
which(DEG_edgeR$PValue<pvalue))
DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"
head(DEG_edgeR)
library(AnnoProbe)
ag=annoGene(rownames(DEG_edgeR),
ID_type = 'ENSEMBL',species = 'human'
)
head(ag)
DEG_edgeR$ENSEMBL=rownames(DEG_edgeR)
deg_anno=merge(ag,DEG_edgeR,by='ENSEMBL')
head(deg_anno)
library(EnhancedVolcano)
colnames(deg_anno)
EnhancedVolcano(deg_anno,
lab =deg_anno$SYMBOL,
x = 'logFC',
y = 'FDR')出圖如下所示：
正常火山圖可以看到這個時候的火山圖其實是非常正常的！
可以跟原始的差異分析做對比首先查看兩次差異分析：
reulsts_edgeR = deg_anno
colnames(reulsts_edgeR)
load( file = "./data/Step03-DESeq2_nrDEG.Rdata")
reulsts_DEseq2 = deg_anno
colnames(reulsts_DEseq2)
df=merge(reulsts_edgeR[,c(1,2,7,10,11)],
reulsts_DEseq2[,c(1,8,12,13)],
by='ENSEMBL')
table(df$regulated.x,df$regulated.y)如下所示：
> colnames(reulsts_edgeR)
[1] "ENSEMBL" "SYMBOL" "biotypes" "chr" "start" "end"
[7] "logFC" "logCPM" "PValue" "FDR" "regulated"
> colnames(reulsts_DEseq2)
[1] "ENSEMBL" "SYMBOL" "biotypes" "chr"
[5] "start" "end" "baseMean" "log2FoldChange"
[9] "lfcSE" "stat" "pvalue" "padj"
[13] "regulated"需要把兩個差異分析結果矩陣合併，然後對比！
down normal up
down 0 355 670
normal 424 13643 860
up 719 789 0可以看到， 上下調居然是反過來了，哈哈哈，可能是前面的代碼沒有設置因子。
我們可以散點圖更加直觀查看，代碼如下所示：
library(ggplot2)
library(ggpubr)
library(ggstatsplot)
p <- ggplot(df, aes(x=df$logFC,
df$log2FoldChange))+
geom_point()+
labs(x = "logFC",
y = "log2FoldChange")+
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept = c( -1.5,0,1.5),lty=4,col="#ec183b",lwd=0.8) +
geom_hline(yintercept = c(-1.5,0,1.5 ),lty=4,col="#18a5ec",lwd=0.8) +
#xlim(-5,5)+
#ylim(-3,3)+
theme_ggstatsplot()+
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank())
p效果如下：
散點圖更加直觀查看兩個差異分析結果的對比還有個韋恩圖，也可以查看兩次差異分析結果的區別：
代碼如下所示：
library(data.table)
library(ggplot2)
library(ggstatsplot)
library(ggvenn)
a <- list( 'edgeR_up'=split(df$ENSEMBL,df$regulated.x)$up,
'edgeR_down'=split(df$ENSEMBL,df$regulated.x)$down,
'DEseq2_up'=split(df$ENSEMBL,df$regulated.y)$up,
'DEseq2_down'=split(df$ENSEMBL,df$regulated.y)$down
)
p1 <- ggvenn(a,
fill_alpha = 0.9,
fill_color = c('#7fc97f','#beaed4','#fdc086','#ffff99'),
set_name_size = 5,
stroke_size = 0.5)
p1出圖如下：
韋恩圖查看兩次 差異分析的結果的區別因為我們前面兩次 差異分析的因子沒有設置，所以這個上下調是反過來的，但是仍然是可以看到，兩次差異分析的結果的一致性比較好！
學徒作業運行我的代碼，把因子設置好，重新進行散點圖，韋恩圖的繪製！