美文网首页单细胞学习
单细胞分析-拟南芥例子体验2

单细胞分析-拟南芥例子体验2

作者: jjjscuedu | 来源:发表于2021-05-17 20:33 被阅读0次

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

    相关文章

      网友评论

        本文标题:单细胞分析-拟南芥例子体验2

        本文链接:https://www.haomeiwen.com/subject/imxjjltx.html