查看: 5694|回复: 3

单细胞转录组学习笔记之Seurat 3.0(一)

[复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
402
注册时间
2018.4.19
在线时间
896 小时

推广达人宣传达人


发表于 2021.7.5 11:26:04 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
在之前的文章中,已经为大家分享了几个R语言的教程,今天再为大家分享R语言的seurat包的学习笔记。

如果你对10X单细胞转录组感兴趣的话,点 这里 报名参加基迪奥6月份的10X单细胞转录组生信培训班,届时会给大家介绍更多相关内容。
1数据导入

本文的范例数据为seurat官网的pbmc-3k数据,文末有下载链接。当然也可以直接使用 基迪奥10X转录组结题报告中的表达量文件,如下图。




[Shell] 纯文本查看 复制代码
#指定数据所在目录;
data_dir <-\"C:/Users/MHY/Desktop/filtered_gene_bc_matrices/hg19\" 

#载入seurat包;
library(Seurat) 

#读入pbmc数据;
pbmc.data <- Read10X(data.dir =data_dir) 


#查看稀疏矩阵的维度,即基因数和细胞数;
dim(pbmc.data)
#[1]32738  2700 


#预览稀疏矩阵(1~10行,1~6列),. 表示0;
pbmc.data[1:10,1:6] 




2创建Seurat对象与数据过滤

在使用CreateSeuratObject()创建对象的同时,过滤数据质量差的细胞。保留在>=3个细胞中表达的基因;保留能检测到>=200个基因的细胞。
pbmc <- CreateSeuratObject(counts =pbmc.data, project = \"pbmc2700\", min.cells = 3, min.features = 200)

计算每个细胞的线粒体基因转录本数的百分比(%),使用[[ ]] 操作符存放到metadata中;


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

过滤细胞:保留gene数大于200小于2500的细胞;目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据;而保留线粒体基因的转录本数低于5%的细胞,为了过滤掉死细胞等低质量的细胞数据。

pbmc <- subset(pbmc, subset =nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
对过滤后的 QC metrics进行可视化(绘制散点图);

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


表达量数据标准化,LogNormalize的算法:
A = log( 1 + ( UMIA ÷ UMITotal ) × 10000 )
pbmc <- NormalizeData(pbmc,normalization.method = \"LogNormalize\", scale.factor = 10000) 鉴定
表达高变基因(2000个),用于下游分析,如PCA;
pbmc <- FindVariableFeatures(pbmc,selection.method = \"vst\", nfeatures = 2000)
提取表达量变变化最高的10个基因;
top10 <- head(VariableFeatures(pbmc),10)
top10
[1]\"PPBP\"   \"LYZ\"    \"S100A9\" \"IGLL5\"  \"GNLY\"   \"FTL\"    \"PF4\"  
[8]\"FTH1\"   \"GNG11\"  \"S100A8\"

plot3 <- VariableFeaturePlot(pbmc)+NoLegend()
plot4 <- LabelPoints(plot = plot3,points = top10, repel = TRUE,xnudge=0,ynudge=0)
plot4



3PCA分析

PCA分析数据准备,使用ScaleData()进行数据归一化;默认只是标准化高变基因(2000个),速度更快,不影响PCA和分群,但影响热图的绘制。
pbmc <- ScaleData(pbmc,vars.to.regress =\"percent.mt\")

而对所有基因进行标准化的方法如下:
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features =all.genes,vars.to.regress = \"percent.mt\")

线性降维(PCA),默认用高变基因集,但也可通过features参数自己指定;
pbmc <- RunPCA(pbmc, features =VariableFeatures(object = pbmc))

检查PCA分群结果, 这里只展示前12个PC,每个PC只显示3个基因;
print(pbmc[[\"pca\"]], dims = 1:12,nfeatures = 3)


绘制pca散点图;
DimPlot(pbmc, reduction = \"pca\")+NoLegend()



画前2个主成分的热图;
DimHeatmap(pbmc, dims = 1:2, cells = 500,balanced = TRUE)




4确定数据集的分群个数

方法1:Jackstraw置换检验算法;重复取样(原数据的1%),重跑PCA,鉴定p-value较小的PC;计算‘null distribution’(即零假设成立时)时的基因scores。
pbmc <- JackStraw(pbmc, num.replicate =100)
pbmc <- ScoreJackStraw(pbmc, dims =1:20)JackStrawPlot(pbmc, dims = 1:15)




方法2:肘部图(碎石图),基于每个主成分对方差解释率的排名。ElbowPlot(pbmc)



分群个数这里选择10,建议尝试选择多个主成分个数做下游分析,对整体影响不大;在选择此参数时,建议选择偏高的数字(为了获取更多的稀有分群,“宁滥勿缺”);有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。

5非线性降维(UMAP/tSNE)

基于PCA空间中的欧氏距离计算nearest neighbor graph,优化任意两个细胞间的距离权重(输入上一步得到的PC维数)。

pbmc <- FindNeighbors(pbmc, dims = 1:10)

接着优化模型,resolution参数决定下游聚类分析得到的分群数,对于3K左右的细胞,设为0.4-1.2 能得到较好的结果(官方说明);如果数据量增大,该参数也应该适当增大。

pbmc <- FindClusters(pbmc, resolution =0.5)

使用Idents()函数可查看不同细胞的分群;
head(Idents(pbmc), 8)



Seurat提供了几种非线性降维的方法进行数据可视化(在低维空间把相似的细胞聚在一起),比如UMAP和t-SNE,运行UMAP需要先安装\'umap-learn\'包,这里不做介绍。
pbmc <- RunTSNE(pbmc, dims = 1:10)
用DimPlot()函数绘制散点图,reduction = \"tsne\",指定绘制类型;如果不指定,默认先从搜索 umap,然后 tsne, 再然后 pca;也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot(); cols,pt.size分别调整分组颜色和点的大小;
tsneplot<-TSNEPlot(pbmc,label = TRUE,pt.size = 1.5)+ NoLegend()
tsneplot



绘制Marker 基因的tsne图;

FeaturePlot(pbmc, features =c(\"MS4A1\", \"CD14\"))



6为分群重新指定细胞类型

new.cluster.ids <- c(\"Naive CD4T\", \"Memory CD4 T\", \"CD14+ Mono\", \"B\", \"CD8T\", \"FCGR3A+ Mono\", \"NK\", \"DC\",\"Platelet\")
将pbmc的水平属性赋值给new.cluster.ids的names属性;
names(new.cluster.ids)
NULL


levels(pbmc)
[1]\"0\" \"1\" \"2\" \"3\" \"4\"\"5\" \"6\" \"7\" \"8\" \\


names(new.cluster.ids) <- levels(pbmc)

names(new.cluster.ids)
[1]\"0\" \"1\" \"2\" \"3\" \"4\"\"5\" \"6\" \"7\" \"8\"
pbmc <- RenameIdents(pbmc,new.cluster.ids)

绘制tsne图(修改标签后的);
tsneplot2<-TSNEPlot(pbmc,label = TRUE,pt.size = 0.5)+ NoLegend()
tsneplot2



本文所用的练习数据已上传到omicshare论坛,欢迎大家登陆下载~
下载链接:http://www.omicshare.com/forum/thread-4982-1-1.html

参考资料
https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html

本文作者:基迪奥莫北


来源: 单细胞转录组学习笔记之Seurat 3.0(一)
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
100 小时

发表于 2021.7.5 11:30:29 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2019.7.25
在线时间
174 小时

发表于 2022.2.23 09:26:53 | 显示全部楼层
代码里面的///////之类的有点干扰
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2018.5.23
在线时间
12 小时

发表于 2022.3.27 15:29:43 | 显示全部楼层
厉害厉害
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表