seurat标准流程实例之2个10x样本的项目（GSE135927数据集）

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《上海中医药大学研究生》的分享

前面jimmy老师分享了两个祖传的单细胞转录组数据分析代码，非常给力，是标准流程：

其中有一个环节是需要比较seurat分群以及singleR的分群，这样就可以合理的命名啦。在jimmy老师的督促下，我使用老师的代码处理了GSE135927数据集，直接套用了jimmy老师的标准代码，希望对所有的初学者有帮助！

首先进入GEO可以看到是两个10X的样本：

教程目录大纲如下：

• 1、准备原始分析数据

• 2、创建Seurat对象

• 3、过滤质控

• 4、降维聚类

• 5、clusters细胞类型注释

1、准备原始分析数据

# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135927
# 两个样本，共6个文件。需要分到两个文件夹里，并且重命名
fs=list.files('./GSE135927_RAW/','^GSM')
fs

# 自行下载GSE135927数据集的GSE135927_RAW压缩包并且解压哦，这样上面的代码就可以运行啦

# 然后获取两个样本信息，因为是批量，所以下面的代码可能不好理解，需要熟练掌握R语言哦

samples=str_split(fs,'_',simplify = T)[,1]

lapply(unique(samples),function(x){
  y=fs[grepl(x,fs)]
  folder=paste0("GSE135927_RAW/", str_split(y[1],'_',simplify = T)[,1])
  dir.create(folder,recursive = T)
  #为每个样本创建子文件夹
  file.rename(paste0("GSE135927_RAW/",y[1]),file.path(folder,"barcodes.tsv.gz"))
  #重命名文件，并移动到相应的子文件夹里
  file.rename(paste0("GSE135927_RAW/",y[2]),file.path(folder,"features.tsv.gz"))
  file.rename(paste0("GSE135927_RAW/",y[3]),file.path(folder,"matrix.mtx.gz"))
})

## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE

samples=list.files("GSE135927_RAW/")
samples

## [1] "GSM4038043" "GSM4038044"# 是两个文件夹的名字哦

2、创建Seurat对象

# 循环读取两个文件夹下面的10x的的3个文件

sceList = lapply(samples,function(pro){
  folder=file.path("GSE135927_RAW/",pro)
  CreateSeuratObject(counts = Read10X(folder),
                     project = pro )
})
sce.big <- merge(sceList[[1]],
                 y = c(sceList[[2]]),
                 add.cell.ids = samples,
                 project = "ls_12")
sce.big

## An object of class Seurat
## 33538 features across 2042 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)

table(sce.big$orig.ident)

##
## GSM4038043 GSM4038044
##       1359        683

#save(sce.big,file = 'sce.big.merge.ls_12.Rdata')

3、过滤质控

#raw_sce=get(load(file = 'sce.big.merge.ls_12.Rdata'))
#线粒体基因
raw_sce=sce.big
rownames(raw_sce)[grepl('^mt-',rownames(raw_sce),ignore.case = T)]

## [1] "MT-ND1" "MT-ND2" "MT-CO1" "MT-CO2" "MT-ATP8" "MT-ATP6"

#核糖体蛋白基因
rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]

##   [1] "RPL22"          "RPL11"          "RPS6KA1"        "RPS8"
##   [5] "RPL5"           "RPS27"          "RPS6KC1"        "RPS7"
##   [9] "RPS27A"         "RPL31"          "RPL37A"         "RPL32"
##  [13] "RPL15"          "RPSA"           "RPL14"          "RPL29"

#计算上述两种基因在所有细胞中的比例

raw_sce[["percent.mt"]] <- PercentageFeatureSet(raw_sce, pattern = "^MT-")

rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce))]
C<-GetAssayData(object = raw_sce, slot = "counts")
percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
raw_sce <- AddMetaData(raw_sce, percent.ribo, col.name = "percent.ribo")

plot1 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

VlnPlot(raw_sce, features = c("percent.ribo", "percent.mt"), ncol = 2)VlnPlot(raw_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)VlnPlot(raw_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)raw_sce #2042个

## An object of class Seurat
## 33538 features across 2042 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)

#按照三个指标过滤细胞
raw_sce1 <- subset(raw_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
raw_sce1 #1887个

## An object of class Seurat
## 33538 features across 1887 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)

4、降维聚类

sce=raw_sce1
#counts归一化
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",
                     scale.factor = 10000)
GetAssay(sce,assay = "RNA")

## Assay data with 33538 features for 1887 cells
## First 10 features:
##  MIR1302-2HG, FAM138A, OR4F5, AL627309.1, AL627309.3, AL627309.2,
## AL627309.4, AL732372.1, OR4F29, AC114498.1

#前2000个高变feature RNA
sce <- FindVariableFeatures(sce,
                            selection.method = "vst", nfeatures = 2000)
# 中心化，为下一步PCA做准备
sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
#看看前12个主成分
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)

#判断最终选取的主成分数，这里我判断18个
ElbowPlot(sce)

sce <- FindNeighbors(sce, dims = 1:18)
sce <- FindClusters(sce, resolution = 0.9)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1887
## Number of edges: 67622
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8140
## Number of communities: 9
## Elapsed time: 0 seconds

table([e-mail]$RNA_snn_res.0.9)

##
##   0   1   2   3   4   5   6   7   8
## 325 303 275 273 214 145 139 133  80

# 分成9个cluster

#tSNE可视化
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:18, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)

# 在两个样本里9个cluster的分布情况
table([e-mail]$seurat_clusters,[e-mail]$orig.ident)

##
##     GSM4038043 GSM4038044
##   0        108        217
##   1         75        228
##   2        237         38
##   3        271          2
##   4        201         13
##   5         92         53
##   6        137          2
##   7         74         59
##   8         50         30

DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')

#save(dat,file='merge_for_tSNE.pos.Rdata')
#load(file='merge_for_tSNE.pos.Rdata)

#再尝试用ggplot2可视化上图
phe=data.frame(cell=rownames([e-mail]),
               cluster =[e-mail]$seurat_clusters)
tsne_pos=Embeddings(sce,'tsne')
dat=cbind(tsne_pos,phe)
head(dat)

##                                   tSNE_1     tSNE_2
## GSM4038043_AAACCTGAGGGCTCTC-1 -13.779792   1.304764
## GSM4038043_AAACCTGAGTGAACGC-1  20.375308   5.068255
## GSM4038043_AAACCTGGTATGAATG-1 -25.223356  -3.429902
## GSM4038043_AAACCTGTCTGCGACG-1   6.194084  14.146905
## GSM4038043_AAACGGGAGGTGCAAC-1   5.695203 -28.687411
## GSM4038043_AAACGGGTCGCTAGCG-1   5.739440 -10.725417
##                                                        cell cluster
## GSM4038043_AAACCTGAGGGCTCTC-1 GSM4038043_AAACCTGAGGGCTCTC-1       0
## GSM4038043_AAACCTGAGTGAACGC-1 GSM4038043_AAACCTGAGTGAACGC-1       2
## GSM4038043_AAACCTGGTATGAATG-1 GSM4038043_AAACCTGGTATGAATG-1       0
## GSM4038043_AAACCTGTCTGCGACG-1 GSM4038043_AAACCTGTCTGCGACG-1       3
## GSM4038043_AAACGGGAGGTGCAAC-1 GSM4038043_AAACGGGAGGTGCAAC-1       5
## GSM4038043_AAACGGGTCGCTAGCG-1 GSM4038043_AAACGGGTCGCTAGCG-1       1

p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype =2,show.legend = F)+coord_fixed()
print(p)

theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black"))
p+theme+guides(colour = guide_legend(override.aes = list(size=5)))

#寻找每个cluster的高变代表基因，并选取前5个，进行可视化
table([e-mail]$seurat_clusters)

##
##   0   1   2   3   4   5   6   7   8
## 325 303 275 273 214 145 139 133  80

p <- list()
for( i in unique([e-mail]$seurat_clusters) ){
  markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25)
  print(x = head(markers_df))
  markers_genes =  rownames(head(x = markers_df, n = 4))
  p1 <- VlnPlot(object = sce, features =markers_genes,log =T,ncol = 2)
  p[[i]][[1]] <- p1
  p2 <- FeaturePlot(object = sce, features=markers_genes,ncol = 2)
  p[[i]][[2]] <- p2
}

##                p_val avg_logFC pct.1 pct.2     p_val_adj
## NRP2   1.051573e-129 1.1153915 0.585 0.059 3.526767e-125
## RPL32  1.930788e-125 0.7562735 1.000 0.976 6.475476e-121
## RPS3A  2.407233e-125 0.7527394 1.000 0.985 8.073378e-121
## EEF1A1 1.509065e-122 0.7042242 1.000 0.997 5.061102e-118
## RPL18A 6.076286e-119 0.7034375 1.000 0.983 2.037865e-114
## RPS14  8.819300e-115 0.5999948 1.000 0.981 2.957817e-110

p[[1]][[2]]p[[2]][[1]]

#查阅所有的marker基因，可用于人工注释cell type
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25,
                              thresh.use = 0.25)
DT::datatable(sce.markers)

#热图可视化每个cluster的marker基因表达差异
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)

#save(sce,sce.markers,file = 'sce.output.merge.ls_12.Rdata')

5、clusters细胞类型注释

# 这里的代码是针对上一步分的9个cluster注释celltype，并不是jimmy老师那样的针对每个样本独立注释哦，而且呢，因为我电脑网络问题，采取了云盘下载singleR数据库的方式，如果你的网络好，就直接看jimmy老师的教程哈！

refdata <- get(load("ref_Monaco_114s.RData"))
sce_for_SingleR <- GetAssayData(sce, slot="data")
clusters <- [e-mail]$seurat_clusters
pred.hesc <- SingleR(test = sce_for_SingleR, ref = refdata,
                     labels = refdata$label.fine,
                     #因为样本主要为免疫细胞（而不是全部细胞），因此设置为label.fine
                     method = "cluster", clusters = clusters,
                     #这里我们为上一步分的9个cluster注释celltype
                     assay.type.test = "logcounts", assay.type.ref = "logcounts")
table(pred.hesc$labels)

##
##  Central memory CD8 T cells Effector memory CD8 T cells
##                           3                           2
##          T regulatory cells                  Th17 cells
##                           2                           1
##                   Th2 cells
##                           1

plotScoreHeatmap(pred.hesc)

#tSNE可视化
celltype = data.frame(ClusterID=rownames(pred.hesc), celltype=pred.hesc$labels, stringsAsFactors = F)
#如下为sce对象注释细胞cluster鉴定结果。

[e-mail]$celltype = "NA"
#先新增列celltype，值均为NA，然后利用下一行代码循环填充
for(i in 1:nrow(celltype)){
  [e-mail][which([e-mail]$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}

DimPlot(sce, group.by="celltype", label=T, label.size=5, reduction='tsne')

#利用ggplot2再可视化上图
phe=data.frame(cell=rownames([e-mail]),
               cluster =clusters,
               celltype=[e-mail]$celltype)
tsne_pos=Embeddings(sce,'tsne')
dat=cbind(tsne_pos,phe)
head(dat)

##                                   tSNE_1     tSNE_2
## GSM4038043_AAACCTGAGGGCTCTC-1 -13.779792   1.304764
## GSM4038043_AAACCTGAGTGAACGC-1  20.375308   5.068255
## GSM4038043_AAACCTGGTATGAATG-1 -25.223356  -3.429902
## GSM4038043_AAACCTGTCTGCGACG-1   6.194084  14.146905
## GSM4038043_AAACGGGAGGTGCAAC-1   5.695203 -28.687411
## GSM4038043_AAACGGGTCGCTAGCG-1   5.739440 -10.725417
##                                                        cell cluster
## GSM4038043_AAACCTGAGGGCTCTC-1 GSM4038043_AAACCTGAGGGCTCTC-1       0
## GSM4038043_AAACCTGAGTGAACGC-1 GSM4038043_AAACCTGAGTGAACGC-1       2
## GSM4038043_AAACCTGGTATGAATG-1 GSM4038043_AAACCTGGTATGAATG-1       0
## GSM4038043_AAACCTGTCTGCGACG-1 GSM4038043_AAACCTGTCTGCGACG-1       3
## GSM4038043_AAACGGGAGGTGCAAC-1 GSM4038043_AAACGGGAGGTGCAAC-1       5
## GSM4038043_AAACGGGTCGCTAGCG-1 GSM4038043_AAACGGGTCGCTAGCG-1       1
##                                                  celltype
## GSM4038043_AAACCTGAGGGCTCTC-1                   Th2 cells
## GSM4038043_AAACCTGAGTGAACGC-1 Effector memory CD8 T cells
## GSM4038043_AAACCTGGTATGAATG-1                   Th2 cells
## GSM4038043_AAACCTGTCTGCGACG-1          T regulatory cells
## GSM4038043_AAACGGGAGGTGCAAC-1  Central memory CD8 T cells
## GSM4038043_AAACGGGTCGCTAGCG-1  Central memory CD8 T cells

p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=celltype))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=celltype,color=celltype),
                 geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
print(p)

theme= theme(panel.grid =element_blank()) +   ## 删去网格
  theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
  theme(axis.line = element_line(size=1, colour = "black"))
p+theme+guides(colour = guide_legend(override.aes = list(size=5)))

理论上jimmy老师的代码是可以适用于绝大部分10x项目数据哦， 大家如果有自己的项目，赶快用起来吧！

写在后面

这些代码大家都可以测试一下，也欢迎加入我们的其它分门别类的学习交流群，如下：

至少在关键时刻，有人可以告诉你该如何搜索。