批次效应去不去除呢?这是个问题
考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求，我开通了一个专栏《100个单细胞转录组数据降维聚类分群图表复现》，也亲自示范了几个，不过自己带娃，读博，时间精力有限，所以把剩余的90多个任务安排了学徒，实习生，学员。真的是太棒了，群策群力！
今天介绍的文章是2021年3月发表在Nature Communications的《Single-cell transcriptional changes associated with drug tolerance and response to combination therapies in cancer》，链接是：https://www.nature.com/articles/s41467-021-21884-z
本研究的目的在单细胞分辨率下观察非小细胞肺癌细胞系体内药物耐受状态，寻找预后和治疗分子靶标。本次文献图表复现的任务是对PC9细胞系药物处理前后的单细胞测序数据分析，发现作者并没有进行整合。
数据在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134836
作者提供了标准的10x文件，如下所示：
这两个样品各自都是一个压缩包，解压后就是标准的10x文件。
分别是：PC9 细胞系和2 µM erlotinib 处理后的PC9细胞系;
下面我们开始复现：
在运行环境下创建一个GSE134836_RAW文件夹，里面包含下面两个文件夹：
step1：读取表达矩阵文件，然后构建Seurat对象
构建Seurat对象后保存为 scesce.all_raw.rds 文件，供后续使用
###### step1:导入数据 ######
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
samples=list.files('GSE134836_RAW/')
samples
sceList = lapply(samples,function(pro){
#pro=samples[1]
folder=file.path('GSE134836_RAW',pro)
print(pro)
print(folder)
print(list.files(folder))
sce=CreateSeuratObject(counts = Read10X(folder),
project = pro )
return(sce)
})
names(sceList)
samples
names(sceList) = samples
sceList
sce <- merge(sceList[[1]], y= sceList[ -1 ] ,
add.cell.ids=samples)
#计算线粒体基因比例
sce=PercentageFeatureSet(sce, "^MT-", col.name = "percent_mito")
#计算核糖体基因比例
sce=PercentageFeatureSet(sce, "^RP[SL]", col.name = "percent_ribo")
#计算红血细胞基因比例
sce=PercentageFeatureSet(sce, "^HB[^(P)]", col.name = "percent_hb")
#根据上述指标，过滤低质量细胞/基因
#过滤指标1:最少表达基因数的细胞&最少表达细胞数的基因
selected_c <- WhichCells(sce, expression = nFeature_RNA > 300)
selected_f <- rownames(sce)[Matrix::rowSums(sce@assays$RNA@counts > 0 ) > 3]
sce.all.filt <- subset(sce, features = selected_f, cells = selected_c)
#过滤指标2:线粒体/核糖体基因比例(根据上面的violin图)
selected_mito <- WhichCells(sce.all.filt, expression = percent_mito < 20)
selected_ribo <- WhichCells(sce.all.filt, expression = percent_ribo > 3)
selected_hb <- WhichCells(sce.all.filt, expression = percent_hb < 0.1)
sce.all.filt <- subset(sce.all.filt, cells = selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_ribo)
sce.all.filt <- subset(sce.all.filt, cells = selected_hb)
save(sce.all.filt, file = "sce.all_filt.Rdata")
step2：降维聚类分群
##### step2 降维聚类分群#####
#读入数据
load("sce.all_filt.Rdata")
#数据预处理
sce.all=sce.all.filt
sce.all = NormalizeData(sce.all)
sce.all <- FindVariableFeatures(sce.all, nfeatures = 3000)
sce.all
# 默认后续使用 HVGs 进行 特征寻找
sce.all = ScaleData(sce.all
#,vars.to.regress = c("nFeature_RNA", "percent_mito")
) #加上vars.to.regress参数后运行so slow!!!
# 降维处理
sce.all <- RunPCA(sce.all, npcs = 30)
sce.all <- RunTSNE(sce.all, dims = 1:30)
sce.all <- RunUMAP(sce.all, dims = 1:30)
# 聚类分析
sce.all <- FindNeighbors(sce.all, dims = 1:30)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1,1.2,1.5,1.6,2)) {
sce.all <- FindClusters(sce.all, graph.name = "RNA_snn", resolution = res, algorithm = 1)
}
p1 <- clustree(sce.all@meta.data, prefix = "RNA_snn_res.")
p1
p2<- DimPlot(sce.all, reduction = "umap", group.by = "RNA_snn_res.1.6", label = T,label.box = T)+NoAxes()
p3<-DimPlot(sce.all, reduction = "umap", group.by = "orig.ident")+NoAxes()
p3+p2
可以看到，如果两个数据集纯粹的使用降维聚类分群，会出现两个样品完全分离的情况，这个也是作者文章里面出现的图表。
但是考虑到 PC9 细胞系和2 µM erlotinib 处理后的PC9细胞系，这两个样品的差异仅仅是药物处理，不太可能导致两个样品里面完全就没有同样的细胞亚群。
所以选择使用数据整合思路
step3: 数据整合
这里采用最常见的CCA整合，代码如下：
###### step3:合并2个数据集######
#拆分为 个seurat子对象
sce.all.list <- SplitObject(sce.all.filt, split.by = "orig.ident")
#数据整合之前要对每个样本的seurat对象进行数据标准化和选择高变基因
for (i in 1:length(sce.all.list)) {
print(i)
sce.all.list[[i]] <- NormalizeData(sce.all.list[[i]], verbose = FALSE)
sce.all.list[[i]] <- FindVariableFeatures(sce.all.list[[i]], selection.method = "vst",
nfeatures = 2000, verbose = FALSE)
}
###以VariableFeatures为基础寻找锚点，运行时间较长
alldata.anchors <- FindIntegrationAnchors(object.list = sce.all.list, dims = 1:30,
reduction = "cca")
##利用锚点整合数据，运行时间较长
sce.all.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:30, new.assay.name = "CCA")
sce.all.int=ScaleData(sce.all.int)
sce.all.int=RunPCA(sce.all.int, npcs = 30)
sce.all.int=RunTSNE(sce.all.int, dims = 1:30)
sce.all.int=RunUMAP(sce.all.int, dims = 1:30)
p4.compare=plot_grid(ncol = 3,
DimPlot(sce.all, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"),
DimPlot(sce.all, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"),
DimPlot(sce.all, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"),
DimPlot(sce.all.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
DimPlot(sce.all.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")
)
p4.compare
整合前后的数据对比：
可以看到，整合前，两个样品泾渭分明的界限，但是整合后，两个样品的单细胞就混合在一起了。这个时候其实把批次效应和生物学差异都给抹去了，非常的尴尬！
然后仍然是常规的降维聚类分群：
# 聚类分析
sce.all.int=FindNeighbors(sce.all.int, dims = 1:30, k.param = 60, prune.SNN = 1/15)
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1,1.2,1.5,1.6,2)) {
sce.all.int=FindClusters(sce.all.int, graph.name = "CCA_snn", resolution = res, algorithm = 1)
}
apply(sce.all.int@meta.data[,grep("CCA_snn_res",colnames(sce.all.int@meta.data))],2,table)
p5_tree <- clustree(sce.all.int@meta.data, prefix = "CCA_snn_res.")
p5_tree
p6<- DimPlot(sce.all.int, reduction = "umap", group.by = "CCA_snn_res.1.6", label = T,label.box = T)+NoAxes()
p7<-DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident")+NoAxes()
p7+p6
可以看到，这次的降维聚类分群，两个样品（ PC9 细胞系和2 µM erlotinib 处理后的PC9细胞系）,有共同细胞亚群，也有各自独有的细胞亚群。
这个时候，其实CCA整合会把批次效应和生物学差异都给抹去了，非常的尴尬！这个问题其实目前是无解的，我们的公众号生信技能树之前也写了很多关于批次效应的推文，例如
大家可以当做是背景知识读一下！
尤其是单细胞领域 ：
也就是说急于纠正一切看起来像是批次效应的事情，有可能会掩盖了有趣的生物学，也许背后藏着潜在的生物学意义。
写在文末
咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！见：急！计划招募100个单细胞爱好者，免费学全套单细胞降维聚类分群和生物学亚群注释