前往小程序，Get更优阅读体验！
立即前往
文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
首页
学习
活动
专区
工具
TVP最新优惠活动
返回腾讯云官网
社区首页 >专栏 >seurat标准流程实例之2个10x样本的项目（GSE135927数据集）

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

作者头像
生信技能树
发布于 2020-09-14 16:37:29
5.9K2
发布于 2020-09-14 16:37:29
举报
文章被收录于专栏：生信技能树

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

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

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

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

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

教程目录大纲如下：

  • 1、准备原始分析数据
  • 2、创建Seurat对象
  • 3、过滤质控
  • 4、降维聚类
  • 5、clusters细胞类型注释

1、准备原始分析数据

代码语言：javascript
复制
# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135927
# 两个样本，共6个文件。需要分到两个文件夹里，并且重命名
fs=list.files('./GSE135927_RAW/','^GSM')
fs
代码语言：javascript
复制
# 自行下载GSE135927数据集的GSE135927_RAW压缩包并且解压哦，这样上面的代码就可以运行啦
代码语言：javascript
复制
# 然后获取两个样本信息，因为是批量，所以下面的代码可能不好理解，需要熟练掌握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"))
})
代码语言：javascript
复制
## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
代码语言：javascript
复制
samples=list.files("GSE135927_RAW/")
samples
代码语言：javascript
复制
## [1] "GSM4038043" "GSM4038044"
代码语言：javascript
复制
# 是两个文件夹的名字哦

2、创建Seurat对象

代码语言：javascript
复制
# 循环读取两个文件夹下面的10x的的3个文件
代码语言：javascript
复制
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
代码语言：javascript
复制
## An object of class Seurat
## 33538 features across 2042 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)
代码语言：javascript
复制
table(sce.big$orig.ident)
代码语言：javascript
复制
##
## GSM4038043 GSM4038044
##       1359        683
代码语言：javascript
复制
#save(sce.big,file = 'sce.big.merge.ls_12.Rdata')

3、过滤质控

代码语言：javascript
复制
#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)]
代码语言：javascript
复制
##  [1] "MT-ND1"  "MT-ND2"  "MT-CO1"  "MT-CO2"  "MT-ATP8" "MT-ATP6"
代码语言：javascript
复制
#核糖体蛋白基因
rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]
代码语言：javascript
复制
##   [1] "RPL22"          "RPL11"          "RPS6KA1"        "RPS8"
##   [5] "RPL5"           "RPS27"          "RPS6KC1"        "RPS7"
##   [9] "RPS27A"         "RPL31"          "RPL37A"         "RPL32"
##  [13] "RPL15"          "RPSA"           "RPL14"          "RPL29"
代码语言：javascript
复制
#计算上述两种基因在所有细胞中的比例
代码语言：javascript
复制
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))
代码语言：javascript
复制
VlnPlot(raw_sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
代码语言：javascript
复制
VlnPlot(raw_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
代码语言：javascript
复制
VlnPlot(raw_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)
代码语言：javascript
复制
raw_sce #2042个
代码语言：javascript
复制
## An object of class Seurat
## 33538 features across 2042 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)
代码语言：javascript
复制
#按照三个指标过滤细胞
raw_sce1 <- subset(raw_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
raw_sce1 #1887个
代码语言：javascript
复制
## An object of class Seurat
## 33538 features across 1887 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)

4、降维聚类

代码语言：javascript
复制
sce=raw_sce1
#counts归一化
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",
                     scale.factor = 10000)
GetAssay(sce,assay = "RNA")
代码语言：javascript
复制
## 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
代码语言：javascript
复制
#前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)
代码语言：javascript
复制
#判断最终选取的主成分数，这里我判断18个
ElbowPlot(sce)
代码语言：javascript
复制
sce <- FindNeighbors(sce, dims = 1:18)
sce <- FindClusters(sce, resolution = 0.9)
代码语言：javascript
复制
## 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
代码语言：javascript
复制
table(sce@meta.data$RNA_snn_res.0.9)
代码语言：javascript
复制
##
##   0   1   2   3   4   5   6   7   8
## 325 303 275 273 214 145 139 133  80
代码语言：javascript
复制
# 分成9个cluster

#tSNE可视化
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:18, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
代码语言：javascript
复制
# 在两个样本里9个cluster的分布情况
table(sce@meta.data$seurat_clusters,sce@meta.data$orig.ident)
代码语言：javascript
复制
##
##     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
代码语言：javascript
复制
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')
代码语言：javascript
复制
#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)
代码语言：javascript
复制
##                                   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
代码语言：javascript
复制
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)
代码语言：javascript
复制
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)))
代码语言：javascript
复制
#寻找每个cluster的高变代表基因，并选取前5个，进行可视化
table(sce@meta.data$seurat_clusters)
代码语言：javascript
复制
##
##   0   1   2   3   4   5   6   7   8
## 325 303 275 273 214 145 139 133  80
代码语言：javascript
复制
p <- 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
}
代码语言：javascript
复制
##                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
代码语言：javascript
复制
p[[1]][[2]]
代码语言：javascript
复制
p[[2]][[1]]
代码语言：javascript
复制
#查阅所有的marker基因，可用于人工注释cell type
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25,
                              thresh.use = 0.25)
DT::datatable(sce.markers)
代码语言：javascript
复制
#热图可视化每个cluster的marker基因表达差异
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)
代码语言：javascript
复制
#save(sce,sce.markers,file = 'sce.output.merge.ls_12.Rdata')

5、clusters细胞类型注释

代码语言：javascript
复制
# 这里的代码是针对上一步分的9个cluster注释celltype，并不是jimmy老师那样的针对每个样本独立注释哦，而且呢，因为我电脑网络问题，采取了云盘下载singleR数据库的方式，如果你的网络好，就直接看jimmy老师的教程哈！
代码语言：javascript
复制
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)
代码语言：javascript
复制
##
##  Central memory CD8 T cells Effector memory CD8 T cells
##                           3                           2
##          T regulatory cells                  Th17 cells
##                           2                           1
##                   Th2 cells
##                           1
代码语言：javascript
复制
plotScoreHeatmap(pred.hesc)
代码语言：javascript
复制
#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')
代码语言：javascript
复制
#利用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)
代码语言：javascript
复制
##                                   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
代码语言：javascript
复制
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)
代码语言：javascript
复制
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项目数据哦， 大家如果有自己的项目，赶快用起来吧！

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2020-09-10，如有侵权请联系 cloudcommunity@tencent.com 删除
geo

本文分享自 生信技能树 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ，欢迎热爱写作的你一起参与！

geo
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
LV.
文章
0
获赞
0
目录
  • 1、准备原始分析数据
  • 2、创建Seurat对象
  • 3、过滤质控
  • 4、降维聚类
  • 5、clusters细胞类型注释
相关产品与服务
灰盒安全测试
腾讯知识图谱（Tencent Knowledge Graph，TKG）是一个集成图数据库、图计算引擎和图可视化分析的一站式平台。支持抽取和融合异构数据，支持千亿级节点关系的存储和计算，支持规则匹配、机器学习、图嵌入等图数据挖掘算法，拥有丰富的图数据渲染和展现的可视化方案。
产品介绍产品文档
领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

Copyright © 2013 - 2024 Tencent Cloud.

All Rights Reserved. 腾讯云 版权所有

登录 后参与评论
2