https://www.jianshu.com/p/abc49787422e
前面一个文章介绍了,cellranger对于这个拟南芥例子的分析,发现鉴定的基因数目和文章中说的还是有出入,不确定是不是用的filter参数不一致的原因。
A single-cell RNA sequencing profiles the developmental landscape of Arabidopsis root分析结果anyway,我们基于我们的结果,接下来对其进行聚类和轨迹分析。
=====用seurat进行聚类分析=======
输入矩阵文件从输入矩阵结果可以看出,有33602个基因(其实鉴定到的基因是24060个),14539个cell,20990986个非0的表达值
pbmc.data <- Read10X(data.dir = "filtered_feature_bc_matrix/") //cellranger的输出
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200) //进行了初步质控,过滤掉了小于200个基因的cell和小于3个cell的基因
pbmc基本信息注:经初步过滤后,可以发现有22,229个基因,14539个cell。
==查看线粒体和叶绿体比率来进行过滤===
feature文件信息注:这里通过查看基因name实现,feature第二列
pbmc[["percent.mg"]] <- PercentageFeatureSet(pbmc, pattern ="^ATMG") //计算每个cell鉴定到了线粒体比例
pbmc[["percent.cg"]] <- PercentageFeatureSet(pbmc, pattern ="^ATCG") //计算每个cell鉴定到的叶绿体比例
VlnPlot(pbmc, features = c("nFeature_RNA","nCount_RNA","percent.mg","percent.cg"), ncol =4)
pbmc主要统计量小提琴图注:#nFeature_RNA:代表的是在该细胞中共检测到的表达量大于0的基因个数,nCount_RNA:代表的是该细胞中所有基因的表达量之和,percent.mg:代表的是线粒体基因表达量的百分比;percent.cg:代表的是叶绿体体基因表达量的百分比;
指标相关性分析从nCount和gene之间的关系图可以看到非常明显的一个相关性,当gene个数为5000时对应的umi在50000左右。以nFeature_RNA为例,可以看到数值在5000以上的细胞是非常少的,可以看做是离群值,所以在筛选时,如果一个细胞中检测到的基因个数大于5000,就可以进行过滤。
所以我们的过滤标准如下:
pbmc <- subset(pbmc, subset = nFeature_RNA >500& nFeature_RNA <5000& percent.mg <1 & percent.mg <3)
过滤之后结果pbmc <- NormalizeData(pbmc, normalization.method ="LogNormalize", scale.factor =10000) //归一化
pbmc <- FindVariableFeatures(pbmc, selection.method ="vst", nfeatures =5000) //提取细胞间变异系数较大的基因用于下游分析,我们这个例子中提取的前5000
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes) //这两步scale,保持同一个基因在不同cell之间的方差为1,均值为0,此步骤在下游分析中具有相同的权重,因此高表达的基因不会占主导地位
注:聚类分析用于识别细胞亚型,在Seurat中,不是直接对所有细胞进行聚类分析,而是首先进行PCA主成分分析,然后挑选贡献量最大的几个主成分,用挑选出的主成分的值来进行聚类分析
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) //PCA分析
VizDimLoadings(pbmc, dims =1:2, reduction ="pca")
各个基因在PC1和PC2中值DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1:9, cells = 500, balanced = TRUE) //主要用来查看数据集中的异质性的主要来源,并且可以确定哪些PC维度可以用于下一步的下游分析。
下面我们需要确定哪些主成分所代表的基因可以进入下游分析,这里可以使用JackStraw做重抽样分析。用JackStrawPlot可视化看看哪些主成分可以进行下游分析。
“ElbowPlot”:基于每个分量所解释的方差百分比对主要成分进行排名。 在此示例中,我们可以在PC20周围观察到“elbow ”,这表明大多数真实信号是在前20台PC中捕获的。
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)
head(Idents(pbmc), 5) //聚类细胞
==非线性维度约化(UMAP/TSNE)========
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "umap", label = TRUE)
UMAP结果 tSNE结果=====发现差异表达特征======
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
cluster1中top5 markerspbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
cluster0 top5 markerVlnPlot(pbmc, features = c("AT3G53980", "AT4G23670"),slot = "counts", log = TRUE)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
====给每个cluster根据marker基因命名====
在这个数据集的情况下,我们可以使用 canonical markers 轻松地将无偏聚类与已知的细胞类型相匹配。比如,我们可以利用已报道的基因,对于相应的cluster进行注释:cluster 17为lateral root
marker基因列表比如根据上述几个cluster,我们可以确定下面几个cluster的命名(PS:当然我还没分完)
cluster命名new.cluster.ids <- c("0","1","2","root cap","4","lateral root","6","7","root stele","9","phloem","endodermis","12","13","14","15","16","lateral root")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5)+ NoLegend() //同理可以把别的cluster也标注出来
UMAP注释图
网友评论