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 markers
pbmc.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 marker
VlnPlot(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注释图
网友评论