jamesbang

V1

2022/09/26阅读：41主题：雁栖湖

# 🤒 limma | 分层样本的差异分析这样搞！（二）

## 2用到的包

``rm(list = ls())library(tidyverse)library(limma)library(GEOquery)``

## 3示例数据

3个样本中对`T细胞``B细胞`分别进行了转录组分析。

``GSE194314 <- getGEO('GSE194314', destdir=".",getGPL = F)exprSet <- exprs(GSE194314[[1]]) exprSet <- normalizeBetweenArrays(exprSet) %>%   log2(.)``

## 4获取分组数据

``pdata <- pData(GSE194314[[1]])``

## 5整理分组数据

``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)``

## 6Multi-level的处理与理解

### 6.1 整理分组矩阵

``measures <- factor(paste(treatment, cell_type,sep="."))design <- model.matrix(~ 0 + measures)colnames(design) <- levels(measures)``

### 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                )``

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                )``

3️⃣ 再比较一下`Control`组内`T细胞``B细胞`之间的差异表达基因。

``res_Bcell_vs_Tcell_For_Control <- topTable(fit2,                 coef = "Bcell_vs_Tcell_For_Control",                n = Inf,                #p.value = 0.05                )``

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                )``

## 7可视化

``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)``

##### jamesbang
V1

wx🔍: Grassssss 卷起来了