jamesbang

V1

2022/09/25阅读：62主题：雁栖湖

🤒 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]]) boxplot(log2(exprSet))``

``exprSet <- normalizeBetweenArrays(exprSet) %>%   log2(.)boxplot(exprSet)``

nice!~ 齐了，接着做吧。

4获取分组数据

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

5整理分组数据

``individuals <- factor(unlist(lapply(pdata\$characteristics_ch1.1,function(x) strsplit(as.character(x),":")[[1]][2])))treatment <- unlist(lapply(pdata\$characteristics_ch1.2,function(x) strsplit(as.character(x),":")[[1]][2]))treatment <- factor(treatment,levels = unique(treatment))``

6非配对处理

6.1 整理分组矩阵

``design_non_paried <- model.matrix(~ 0 + treatment)colnames(design_non_paried) <- c("Control","anti-BTLA")fit1 <- lmFit(exprSet,design_non_paried)fit1 <- eBayes(fit1)``

6.2 差异分析

``allDiff_non_paired <- topTable(fit1,                             adjust = 'BH',                             coef =  "anti-BTLA",                             n = Inf,                             #p.value = 0.05                             )``

7配对处理

7.1 整理分组矩阵

``design_paried <- model.matrix(~ individuals + treatment)fit2 <- lmFit(exprSet,design_paried)fit2 <- eBayes(fit2)``

7.2 差异分析

``allDiff_paired <- topTable(fit2,                             adjust = 'BH',                             coef = "treatment anti-BTLA",                             n = Inf)``

8对比结果

``library(EnhancedVolcano)library(patchwork)p1 <- EnhancedVolcano(allDiff_non_paired,    lab = rownames(allDiff_non_paired),    x = 'logFC',    y =  'P.Value',    title = 'non_paired',    pointSize = 3.0,    labSize = 6.0,    legendPosition = 'right',    pCutoff = 0.05,    FCcutoff = 1)  p2 <- EnhancedVolcano(allDiff_paired,    lab = rownames(allDiff_paired),    x = 'logFC',    y =  'P.Value',    title = 'paired',    pointSize = 3.0,    labSize = 6.0,    legendPosition = 'right',    pCutoff = 0.05,    FCcutoff = 1)p1 + p2``

9小彩蛋

jamesbang
V1

wx🔍: Grassssss 卷起来了