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] TRUEsamples=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 secondstable(sce@meta.data$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(sce@meta.data$seurat_clusters,sce@meta.data$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 30DimPlot(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(sce@meta.data),
cluster =sce@meta.data$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 1p=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(sce@meta.data$seurat_clusters)##
## 0 1 2 3 4 5 6 7 8
## 325 303 275 273 214 145 139 133 80p <- list()
for( i in unique(sce@meta.data$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 <- sce@meta.data$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
## 1plotScoreHeatmap(pred.hesc)

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

sce@meta.data$celltype = "NA"
#先新增列celltype，值均为NA，然后利用下一行代码循环填充
for(i in 1:nrow(celltype)){
sce@meta.data[which(sce@meta.data$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(sce@meta.data),
cluster =clusters,
celltype=sce@meta.data$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 cellsp=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项目数据哦， 大家如果有自己的项目，赶快用起来吧！

写在后面

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

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

赞 (0)

相关推荐

  • 单细胞转录组分析肿瘤异质性

    欢迎来到医科研,这里是白介素2的读书笔记,跟我一起聊临床与科研的故事, 生物医学数据挖掘,R语言,TCGA.GEO, SEER数据挖掘. 肿瘤异质性指的是不同肿瘤细胞可以表现出不同的遗传和表型特征,包 ...

  • 单细胞工具箱|Seurat官网标准流程

    学习单细胞转录组肯定先来一遍Seurat官网的标准流程. 数据来源于Peripheral Blood Mononuclear Cells (PBMC),共2700个单细胞, Illumina Next ...

  • 祖传的单个10x样本的seurat标准代码

    其实单细胞领域进展太快,我那些课程内容关于R包相关的代码基本上过时了,因为R语言本身都经历了一个超级大的变革!考虑到不能把粉丝带歪,我早就全部公开了系列视频课程.还创立了<单细胞天地>这个 ...

  • 祖传的单个10x样本的seurat标准代码（人和鼠需要区别对待）

    昨天发布了:祖传的单个10x样本的seurat标准代码,非常受大家欢迎,我看到有人马上用起来了我的代码,很棒.在10X官网演示数据:https://support.10xgenomics.com/si ...

  • 我的课题只有一个10x样本肿么办？

    前面我们介绍过,如果只有两个10x单细胞转录组样本的数据, 该如何分析,见:两个样品的10x单细胞转录组数据分析策略 ,实际上这个分析策略的文章里面并不是把单细胞转录组数据当做是重点,分析也是很草率, ...

  • Seurat标准流程

    作者 | 单细胞天地小编 刘小泽 课程链接在:http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55 第二单元第7讲:走Seurat标准流程[文 ...

  • 中外合资经营企业参考样本之一项目建议书

    中外合资经营企业参考样本之一项目建议书

  • 三个10X单细胞转录组样本CCA整合

    前面我在单细胞天地分别介绍了如果因为种种原因仅仅是测了一个样本的10X单细胞,或者走经典的2个样本的10X样本该如何分析,并且辅助自己的生物学故事,如下: 我的课题只有一个10x样本肿么办? 两个样品 ...

  • 使用seurat3的merge功能整合8个10X单细胞转录组样本

    本教程演示的数据来源于发表在2017年10月的NC文章:Differentiation dynamics of mammary epithelial cells revealed by single- ...

  • 10X单细胞转录组理论上有3个文件才能被读入R进行seurat分析

    我在单细胞天地教程:表达矩阵逆转为10X的标准输出3个文件,详细介绍过 10X文件的3个标准文件: 比如SRR7722939数据集里面,文件barcodes.tsv 和 genes.tsv,就是表达矩 ...

  • 10X Visium:空间转录组样本制备到数据分析

    男, 一个长大了才会遇到的帅哥, 稳健,潇洒,大方,靠谱. 一段生信缘,一棵技能树. 生信技能树核心成员,单细胞天地特约撰稿人,简书创作者,单细胞数据科学家. 识别复杂生物系统中空间基因表达差异的能力 ...