scRNA分析|多样本merge 和 harmony去批次

生信补给站

发布于 2022-11-11 18:56:40

4.1K0

发布于 2022-11-11 18:56:40

文章被收录于专栏：生信补给站

本部分选择2020年发表在Genome Med 中的单细胞文章（Single-cell transcriptome analysis of tumor and stromal compartments of pancreatic ductal adenocarcinoma primary tumors and metastatic lesions ），其中含有pancreatic ductal adenocarcinoma的10个原发以及 6个转移样本数据。

后续会依据此数据集进行一些单细胞常见分析以及图表绘制的分享。

一数据下载，整理

在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154778路径下下载GSE154778_RAW.tar文件。

在R的工作路径下新建 data文件夹，然后将GSE154778_RAW.tar文件复制到data中，解压打开。

发现为标准的cellranger输出的三个文件，因此只需要将每个样本的三个文件整理到一个文件夹中，然后更改名字即可。

可以在linux中使用命令行处理，可以在window下手动处理，也可以使用R处理。

#scRNA
library(Seurat) 
library(Matrix)
library(stringr)

#更改目录至 解压后数据所在的路径
setwd("./data")

#查看当前有哪些文件
fs=list.files('./','^GSM')
fs
#提取样本名
samples=str_split(fs,'_',simplify = T)[,2] #按照下划线分割，然后四二个元素作为样本名
unique(samples)

#建立样本文件夹，然后更改对应的三个文件的名字
lapply(unique(samples),function(x){
  y=fs[grepl(x,fs)]
  folder=paste(str_split(y[1],'_',simplify = T)[,2],collapse = '')
  dir.create(folder,recursive = T)
  file.rename(y[1],file.path(folder,"barcodes.tsv.gz"))
  file.rename(y[2],file.path(folder,"features.tsv.gz")) #注意版本
  file.rename(y[3],file.path(folder,"matrix.mtx.gz"))
})

folders=list.files('./')
folders

这里需要注意根据seurat版本号确定是 features.tsv.gz 还是 genes.tsv.gz ？是压缩文件即可还是需要再解压一下。

二批量数据读取

输入文件为标准的cellranger 三个文件，那么使用Read10X函数就可以

#批量读取16个10X单细胞转录组数据文件夹
sceList = lapply(folders,function(folder){ 
  CreateSeuratObject(counts = Read10X(folder), 
                     project = folder )
})

除此之外的单细胞数据形式有很多种，读取方式也不太一致，后续会介绍其他几种常见数据的读取方式。

三多样本合并

多样本的单细胞数据合并的时候需要考虑

（1）是否去批次 ?

（2）使用何种方式去批次 ?

这里简单的介绍不考虑批次直接Merge 以及使用harmony方法去批次这2种方式，为了直观的对比一下区别。

1 直接Merge ，不考虑批次效应

sce.all <- merge(sceList[[1]], 
                 y = c(sceList[[2]],sceList[[3]],sceList[[4]],sceList[[5]],
                       sceList[[6]],sceList[[7]],sceList[[8]],sceList[[9]],sceList[[10]],
                       sceList[[11]],sceList[[12]],sceList[[13]],sceList[[14]],sceList[[15]],
                       sceList[[16]]), 
                 add.cell.ids = folders, #添加样本名
                 project = "scRNA")
sce.all #查看

An object of class Seurat  
51911 features across 17086 samples within 1 assay  
Active assay: RNA (51911 features, 0 variable features)

可以看到合并后的sce.all 是一个seurat对象，然后有51911个基因，17086个细胞。这里的sample 表示细胞。

head(sce.all@meta.data)
table(sce.all@meta.data$orig.ident)

sce.all@meta.data 是数据框，里面存储着每个cell的基本信息，随着后续分析可以添加很多维度的信息，来完成文献中常见的图表绘制。

1.2 过滤

同其他组学数据一样，单细胞数据也需要进行过滤，主要是根据每个细胞的线粒体比例、检测到的基因数、检测的UMI总数这三个方面来对细胞进行简单的过滤，没有固定的标准。

前面head(sce.all@meta.data) 可以看到已经得到了ncount 和 nfeature ，此外seurat 提供了PercentageFeatureSet函数可以很方便的计算某一类特征基因的比例，线粒体基因均为MT- 开头， pattern = "^MT-" 即可，注意后面的横线不要漏掉。

#线粒体
sce.all[["percent.mt"]] <- PercentageFeatureSet(sce.all, pattern = "^MT-")
#核糖体 
sce.all[["percent.rp"]] <- PercentageFeatureSet(sce.all, pattern = "^RP")
#指定基因集
HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ") # 红细胞基因
sce.all[["percent.HB"]]<-PercentageFeatureSet(sce.all,features=HB.genes)

head(sce.all@meta.data)

orig.ident nCount_RNA nFeature_RNA percent.mt percent.rp percent.HB
K16733_AAACATACTCGTTT-1     K16733       2464          965 12.6623377   13.35227          0
K16733_AAACCGTGGGTAGG-1     K16733        689          336  2.1770682   29.02758          0
K16733_AAAGCAGAACGTTG-1     K16733       7145         1919  2.4492652   36.72498          0
K16733_AAAGCAGACTGAGT-1     K16733       1655          621  2.1752266   35.52870          0
K16733_AAAGGCCTGCTCCT-1     K16733      14272         2771  1.5414798   38.81026          0
K16733_AAATACTGTGGATC-1     K16733      13832         2541  0.3181029   40.73164          0

过滤参数阈值可以根据数据的情况和状态进行设置，使用VlnPlot展示，

sce <- subset(sce.all, subset = nFeature_RNA > 500 & 
                nFeature_RNA < 5000 & 
                percent.mt < 30)
VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

1.3 Seurat 标准流程

为了查看样本直接Merge的umap分布情况，直接走一遍seurat的标准流程，详细的参数以及函数细节见单细胞工具箱|Seurat官网标准流程。

#标准化
sce <- NormalizeData(sce)
#高变基因
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(sce)
#归一化
sce <- ScaleData(sce, features = all.genes)
#降维聚类
sce <- RunPCA(sce, npcs = 50)
sce <- FindNeighbors(sce, dims = 1:30)
sce <- FindClusters(sce, resolution = 0.5)
sce <- RunUMAP(sce, dims = 1:30)
sce <- RunTSNE(sce, dims = 1:30)

1.4 绘制umap 图

（1）为方便查看，先根据文献修改样本命名，然后添加组别信息

#修改样本名 ，因为没有规则，手动修改的
sce@meta.data$sample[sce@meta.data$orig.ident == "K16733"] <- "P01"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00006"] <- "P02"
sce@meta.data$sample[sce@meta.data$orig.ident == "T2"] <- "P03"
sce@meta.data$sample[sce@meta.data$orig.ident == "T3"] <- "P04"
sce@meta.data$sample[sce@meta.data$orig.ident == "T4"] <- "P05"
sce@meta.data$sample[sce@meta.data$orig.ident == "T5"] <- "P06"
sce@meta.data$sample[sce@meta.data$orig.ident == "T6"] <- "P07"
sce@meta.data$sample[sce@meta.data$orig.ident == "T8"] <- "P08"
sce@meta.data$sample[sce@meta.data$orig.ident == "T9"] <- "P09"
sce@meta.data$sample[sce@meta.data$orig.ident == "T10"] <- "P10"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00008"] <- "MET01"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00013"] <- "MET02"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00014"] <- "MET03"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00016"] <- "MET04"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00019"] <- "MET05"
sce@meta.data$sample[sce@meta.data$orig.ident == "Y00027"] <- "MET06"
#添加组别sce@meta.data$group <- ifelse( grepl("MET",sce@meta.data$sample ) ,"MET" ,"PT" )

可以使用修改idents的方式就行修改。

（2）绘制umap ，分样本umap ，分组别umap图

p1 <- DimPlot(sce, reduction = "umap", pt.size=0.5, label = F,repel = TRUE)
p2 <- DimPlot(sce, reduction = "umap",group.by = "sample", pt.size=0.5, label = F,repel = TRUE)
p3 <- DimPlot(sce, reduction = "umap",group.by = "group", pt.size=0.5, label = F,repel = TRUE)

p1 + p2 +p3

注：label为是否在umap中添加标签，repel = T 则可以减少文字重叠 。

可以看到转移组和原发组的样本分布区别很大，就需要考虑是否去批次。

原文中是转移组和原发组分别单独分析，没有去批次。

2 harmony 去批次

去批次的方法有很多，本文对比展示使用harmony去批次分析。后续会介绍其他去批次的方法。

##############harmony 去批次################
library(harmony)
#按照样本去批次，也可以选择根据group
sce2 = sce %>% RunHarmony("sample", plot_convergence = TRUE)
sce2
#标准流程，参数不变
sce2 <- sce2 %>% 
  RunUMAP(reduction = "harmony", dims = 1:30) %>% 
  FindNeighbors(reduction = "harmony", dims = 1:30) %>% 
  FindClusters(resolution = 0.5) %>% 
  identity()

p11 <- DimPlot(sce2, reduction = "umap", pt.size=0.5, label = F,repel = TRUE)
p22 <- DimPlot(sce2, reduction = "umap",group.by = "sample", pt.size=0.5, label = F,repel = TRUE)
p33 <- DimPlot(sce2, reduction = "umap",group.by = "group", pt.size=0.5, label = F,repel = TRUE)
p11 + p22 +p33