本部分选择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处理。
代码语言：javascript复制#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函数就可以
代码语言：javascript复制#批量读取16个10X单细胞转录组数据文件夹
sceList = lapply(folders,function(folder){
CreateSeuratObject(counts = Read10X(folder),
project = folder )
})除此之外的单细胞数据形式有很多种，读取方式也不太一致，后续会介绍其他几种常见数据的读取方式。
三 多样本合并
多样本的单细胞数据合并的时候需要考虑
（1）是否去批次 ?
（2）使用何种方式去批次 ?
这里简单的介绍 不考虑批次直接Merge 以及 使用harmony方法去批次 这2种方式，为了直观的对比一下区别。
1 直接Merge ，不考虑批次效应
代码语言：javascript复制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 表示细胞。
代码语言：javascript复制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官网标准流程 。
代码语言：javascript复制#标准化
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）为方便查看，先根据文献修改样本命名，然后添加组别信息
代码语言：javascript复制#修改样本名 ，因为没有规则，手动修改的
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图
代码语言：javascript复制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去批次分析。后续会介绍其他去批次的方法。
代码语言：javascript复制##############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
去批次后样本分布的比较一致，但是也可能抹掉了 转移组 和 原发组 本身的区别。
关于是否需要去批次以及使用什么方法去批次 ，需要根据数据 以及 组织背景的真实情况来选择。
◆ ◆ ◆ ◆ ◆
精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）
