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注释图