![alt](https://img-blog.csdnimg.cn/img_convert/acb5aec682680c759b307a8990932668.png)
1写在前面
上期介绍了用limma
包做配对样本的差异分析。
本期介绍一下Multi-level
如何处理吧。 🥳
应用场景:
Control
和Diseased
的T细胞
和B细胞
分层对比。
2用到的包
rm(list = ls())
library(tidyverse)
library(limma)
library(GEOquery)
3示例数据
这里我们还是利用上期介绍的GEO
数据库上的dataset
。😘
在3个样本中对T细胞
和B细胞
分别进行了转录组分析。
每个样本的细胞都分为Control
或anti-BTLA
组。
GSE194314 <- getGEO('GSE194314', destdir=".",getGPL = F)
exprSet <- exprs(GSE194314[[1]])
exprSet <- normalizeBetweenArrays(exprSet) %>%
log2(.)
![alt](https://img-blog.csdnimg.cn/img_convert/e33363b85aed8ea47a73ef7597c93ebd.png)
4获取分组数据
pdata <- pData(GSE194314[[1]])
![alt](https://img-blog.csdnimg.cn/img_convert/a4c2c8d435b6f3f9268a98c937e0a7d7.png)
5整理分组数据
这里我们提取出分组数据后转为factor
。
individuals <- factor(unlist(lapply(pdata$characteristics_ch1.1,function(x) strsplit(as.character(x),": ")[[1]][2])))
cell_type <- factor(unlist(lapply(pdata$characteristics_ch1,function(x) strsplit(as.character(x),": ")[[1]][2])))
cell_type <- gsub("[[:punct:]]","", cell_type)
cell_type <- gsub("\\s","_", cell_type)
treatment <- unlist(lapply(pdata$characteristics_ch1.2,function(x) strsplit(as.character(x),": ")[[1]][2]))
treatment <- factor(treatment,levels = unique(treatment))
treatment<- gsub("-", "_", treatment)
![alt](https://img-blog.csdnimg.cn/img_convert/ae168282ea6f9d12b73d90e3904f6889.png)
6Multi-level的处理与理解
6.1 整理分组矩阵
如果我们的目的是比较两种细胞类型
间的差异,可以在样本内部进行,因为每个样本都会产生两种细胞类型
的值。 🧐
如果我们想将Control
与anti-BTLA
进行比较,这种比较就是在不同样本之间进行了。
这里大家可以理解为,需要进行组内和组间比较,处理样本时需要用到random effect
,在limma
包中需要调用duplicateCorrelation
函数进行处理。😲
measures <- factor(paste(treatment, cell_type,sep="."))
design <- model.matrix(~ 0 + measures)
colnames(design) <- levels(measures)
![alt](https://img-blog.csdnimg.cn/img_convert/8286e2daf8eba546f71381063ea18374.png)
6.2 相关性计算
然后,我们在同一样本上计算不同组学间的相关性。
corfit <- duplicateCorrelation(exprSet, design, block = individuals)
corfit$consensus
6.3 模型拟合
fit <- lmFit(exprSet, design, block = individuals, correlation = corfit$consensus)
6.4 设置比较组
这里我们设置一下比较组,用-
分隔。👐
cm <- makeContrasts(antiBLTA_vs_Control_For_Bcell =
anti_BTLA.Purified_B_cell-Control.Purified_B_cell,
antiBLTA_vs_Control_For_Tcell =
anti_BTLA.Purified_CD4_T_cell-Control.Purified_CD4_T_cell,
Bcell_vs_Tcell_For_Control =
Control.Purified_B_cell-Control.Purified_CD4_T_cell,
Bcell_vs_Tcell_For_antiBLTA =
anti_BTLA.Purified_B_cell-anti_BTLA.Purified_CD4_T_cell,
levels = design )
6.5 差异分析
1️⃣ 现在我们比较一下在B细胞
中anti-BLTA
组和Control
组之间的差异表达基因。
fit2 <- contrasts.fit(fit, cm)
fit2 <- eBayes(fit2)
res_antiBLTA_vs_Control_For_Bcell <- topTable(fit2,
coef = "antiBLTA_vs_Control_For_Bcell",
n = Inf,
#p.value = 0.05
)
![alt](https://img-blog.csdnimg.cn/img_convert/73ad41452dc8f949d22c0b4a7a2e4fa2.png)
2️⃣ 再比较一下在T细胞
中anti-BLTA
组和Control
组之间的差异表达基因。
res_antiBLTA_vs_Control_For_Tcell <- topTable(fit2,
coef = "antiBLTA_vs_Control_For_Tcell",
n = Inf,
#p.value = 0.05
)
![alt](https://img-blog.csdnimg.cn/img_convert/5103f46e2fb6424c44a8b38a77617890.png)
3️⃣ 再比较一下Control
组内T细胞
和B细胞
之间的差异表达基因。
res_Bcell_vs_Tcell_For_Control <- topTable(fit2,
coef = "Bcell_vs_Tcell_For_Control",
n = Inf,
#p.value = 0.05
)
![alt](https://img-blog.csdnimg.cn/img_convert/5c183d85b90345d2874a10ab7f1db21f.png)
4️⃣ 再比较一下anti-BLTA
组内T细胞
和B细胞
之间的差异表达基因。
res_Bcell_vs_Tcell_For_antiBLTA <- topTable(fit2,
coef = "Bcell_vs_Tcell_For_antiBLTA",
n = Inf,
#p.value = 0.05
)
![alt](https://img-blog.csdnimg.cn/img_convert/b9ed249010dd74d0806aaf1cd043272c.png)
7可视化
我们直接用火山图
可视化对比一下吧。🤒 这里我们把阈值
设置为|logFC|>1
,pvalue<0.05
。
library(EnhancedVolcano)
library(patchwork)
p1 <- EnhancedVolcano(res_antiBLTA_vs_Control_For_Bcell,
lab = rownames(res_antiBLTA_vs_Control_For_Bcell),
x = 'logFC',
y = 'P.Value',
title = 'B cell (anti-BTLA vs Control)',
pointSize = 3.0,
labSize = 6.0,
legendPosition = 'right',
pCutoff = 0.05,
FCcutoff = 1)
p2 <- EnhancedVolcano(res_antiBLTA_vs_Control_For_Tcell,
lab = rownames(res_antiBLTA_vs_Control_For_Tcell),
x = 'logFC',
y = 'P.Value',
title = 'T cell (anti-BTLA vs Control)',
pointSize = 3.0,
labSize = 6.0,
legendPosition = 'right',
pCutoff = 0.05,
FCcutoff = 1)
p3 <- EnhancedVolcano(res_Bcell_vs_Tcell_For_Control,
lab = rownames(res_Bcell_vs_Tcell_For_Control),
x = 'logFC',
y = 'P.Value',
title = 'B cell vs T cell (Control)',
pointSize = 3.0,
labSize = 6.0,
legendPosition = 'right',
pCutoff = 0.05,
FCcutoff = 1)
p4 <- EnhancedVolcano(res_Bcell_vs_Tcell_For_antiBLTA,
lab = rownames(res_Bcell_vs_Tcell_For_antiBLTA),
x = 'logFC',
y = 'P.Value',
title = 'B cell vs T cell (anti-BLTA)',
pointSize = 3.0,
labSize = 6.0,
legendPosition = 'right',
pCutoff = 0.05,
FCcutoff = 1)
wrap_plots(p1,p2,p3,p4)
![alt](https://img-blog.csdnimg.cn/img_convert/786305806246a8f8c6b791b396943814.png)
![](https://img-blog.csdnimg.cn/img_convert/b358e8b901eba9c526053c83eea685f3.png)
点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
![alt](https://img-blog.csdnimg.cn/img_convert/8dd44e7fc35e58afcc45606727289ddf.png)
本文由 mdnice 多平台发布