本部分选择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个转移样本数据。
后续会依据此数据集进行一些单细胞常见分析以及图表绘制的分享。
一 数据下载，整理
在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展示，
过滤参数阈值可以根据数据的情况和状态进行设置，使用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
去批次后样本分布的比较一致，但是也可能抹掉了 转移组 和 原发组 本身的区别。
关于是否需要去批次以及使用什么方法去批次 ，需要根据数据 以及 组织背景的真实情况来选择。
