美文网首页
monocle2 使用笔记

monocle2 使用笔记

作者: RedStones | 来源:发表于2022-03-18 17:56 被阅读0次

    1. 文献怎么用 monocle2 ?

    (1) ex1

    01.png 02.png

    作者处理2种细胞使用的策略还不一样!

      1. 中性粒细胞 neutrophil 是先subset出来,去除污染细胞,然后Seurat聚类,找类之间 top20 DEG 做 cell order(就是monocle2的排序);
        根据前体细胞的marker基因的表达,选择那个state为起点。
      1. 而 monocyte 是subset出cluster1,去除污染细胞,然后像上文一样cluster,然后DEG是 Monocle 的 differentialGeneTest() 找的,过滤条件为 qval<1e-10 and cell>20 这样的基因 order cell;
        起点是包含最多WT库的分支。

    (2) ex2

    https://www.sciencedirect.com/science/article/pii/S0092867417305962

    p1.png

    Figure 6. The CD8+ and CD4+ T Cell State Transition Analysis Based on Integrated Expression and TCR Clonality

    (A) The ordering of CD8+ T cells along pseudotime in a two-dimensional state-space defined by Monocle2. Cell orders are inferred from the expression of most dispersed genes across CD8+ T cell populations sans MAIT. Each point corresponds to a single cell, and each color represents a T cell cluster.
    // 使用方差最大的基因 dispersed genes,没说多少个。

    (B) The same pseudotime plot as in (A) for four clusters of CD4+ T helper cells.

    (3) ex3

    https://www.nature.com/articles/s41586-018-0694-x

    To characterize the potential process of T cell functional changes and determine the potential lineage differentiation between diverse T cells, we applied the Monocle (version 2) algorithm50 with the top 700 signature genes of CD8+ T cells excluding MAIT cells, based on the rank of F statistic generated by ANOVA (Supplementary Table 5).
    // 使用ANOVA的F统计量,top 700 个基因。

    (4) ex4

    https://www.nature.com/articles/s41467-020-20019-0

    First, the top 1000 HVGs were selected using the function ‘FindVariableFeatures’ in Seurat v3.

    (10) tips

    • 建议 过滤、分群、差异基因都在 Seurat 中完成,然后只用 monocle2 做 拟时序分析。
    • 不同分组间的细胞尽量不要放在一起做轨迹分析,同一组的生物学重复可以一起分析。明显没有生物学相关性的细胞也不要放在一起做轨迹分析。

    2. 从 Seurat 对象新建 Monocle 对象

    #' get Monocle2 obj from Seurat obj
    #' 
    #' version 0.2
    #'
    #' @param sce Seurat obj
    #'
    #' @return Monocle2 obj
    #' @export
    #'
    #' @examples
    #' cds = Seurat2Monocle2(sce)
    Seurat2Monocle2=function(sce, do.init=T){
      # get data
      #expr_matrix=as.matrix(sce@assays$RNA@counts)
      expr_matrix=GetAssayData(sce, assay = 'RNA', slot = 'counts')
      cell_data=sce@meta.data
      gene_data=data.frame(
        gene_short_name=row.names(sce), #must have this column
        row.names = row.names(sce)
      )
      # new obj
      pd <- new("AnnotatedDataFrame", data = cell_data)
      fd <- new("AnnotatedDataFrame", data = gene_data)
      cds <- newCellDataSet(expr_matrix, 
                            phenoData = pd, 
                            featureData = fd,
                            lowerDetectionLimit = 0.5,
                            expressionFamily=negbinomial.size())
      # init processing
      if(do.init){
        cds <- estimateSizeFactors(cds) #add column: Size_Factor
        cds <- estimateDispersions(cds) #Removing 323 outliers 
        cds=detectGenes(cds, min_expr = 0.1) #add column: num_genes_expressed
      }
      print( head(pData(cds)) )
      #
      return(cds)
    }
    

    示例:

    # subset of Seurat obj 
    sce=subset(scObj_nue, idents = c(9,10), invert=T,
               downsample=50)
    sce=subset(sce, subset= origin=="APC")
    sce #298
    DimPlot(sce, label=T, split.by = "origin")
    
    
    # 1.load data to monocle2
    cds = Seurat2Monocle2(sce)
    class(cds)
    

    3. 拟时序分析

    The ordering workflow
    Step 1: choosing genes that define progress
    Step 2: reducing the dimensionality of the data
    Step 3: ordering the cells in pseudotime

    除了第一步比较灵活外,剩下2步是固定的。最后是各种可视化。

    (1) step1 获取差异基因列表 DEG list

    这个gene list 对拟时序结构具有决定作用。如果结果不理想,这里是优化的重点。

    ## Trajectory Step 1: choosing genes that define progress----
    #第1步) 获取 输入用于排序的基因 symbol list
    
    Monocle官网教程提供了4个选择方法:
        * 选择发育差异表达基因
        * 选择clusters差异表达基因
        * 选择离散程度高的基因
        * 自定义发育marker基因
        
        前三种都是无监督分析方法,细胞发育轨迹生成完全不受人工干预;
        最后一种是半监督分析方法,可以使用先验知识辅助分析。
    
    ##1)使用clusters差异表达基因
    deg.cluster <- FindAllMarkers(scRNA.Osteoclastic)
    diff.genes <- subset(deg.cluster,p_val_adj<0.05)$gene
    HSMM <- setOrderingFilter(HSMM, diff.genes)
    plot_ordering_genes(HSMM)
    
    ##2)使用seurat选择的高变基因
    var.genes <- VariableFeatures(scRNA.Osteoclastic)
    HSMM <- setOrderingFilter(HSMM, var.genes)
    plot_ordering_genes(HSMM)
    
    ##3)使用monocle选择的高变基因
    disp_table <- dispersionTable(HSMM)
    disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
    HSMM <- setOrderingFilter(HSMM, disp.genes)
    plot_ordering_genes(HSMM)
    
    
    ##4) Monocle2 的 differentialGeneTest() 计算的差异基因,比如保留 qval<1e-10 且 cell number >10 的基因
    diff_test_res <- differentialGeneTest(cds[expressed_genes,],
                                          fullModelFormulaStr = "~percent.mt")
    # 20:35 --> 20:37 特别慢,仅仅1k个细胞就需要10min。
    
    ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
    cds <- setOrderingFilter(cds, ordering_genes)
    plot_ordering_genes(cds)
    

    (2). step2 降维

    ## Trajectory step 2: reduce data dimensionality----
    # cds <- setOrderingFilter(cds, ordering_genes)
    cds <- reduceDimension(cds, max_components = 2, method = 'DDRTree')
    

    (3). 按照轨迹排序细胞

    ## Trajectory step 3: order cells along the trajectory ----
    cds <- orderCells(cds)
    

    (4). 可视化

    plot_cell_trajectory(cds, color_by = "Pseudotime")
    plot_cell_trajectory(cds, color_by = "seurat_clusters")
    
    plot_cell_trajectory(cds, color_by = "State")
    # "State" is just Monocle's term for the segment of the tree.
    
    ## 如何重叠太厉害,还可以分面
    plot_cell_trajectory(HSMM, color_by="XXX") +
        facet_wrap(~XXX, nrow=1)
    

    如果自动指定的不对,还可以手动指定哪个 state 是起点。
    如果没有时间序列,就需要根据生物学知识找到root。
    比如: 高度增殖的祖细胞群产生两种有丝分裂后的细胞。所以根应该有表达高水平增殖标记的细胞。

    cds<- orderCells(cds, root_state = 1)
    plot_cell_trajectory(cds, color_by = "Pseudotime")
    

    更多可视化函数:

    plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75) ####以state进行着色
    plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)+facet_wrap(~State, nrow = Cluster_num) ##绘制state的分面图
    
    plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 0.75) ##根据拟时间值着色
    
    plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75) ##根据细胞聚类的进行着色
    plot_cell_trajectory(HSMM_myo, color_by = "Cluster",cell_size = 0.75)+facet_wrap(~Cluster, nrow =Cluster_num) ###绘制clusster的分面图
    
    plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75) +scale_color_manual(values=col) ##如果有Seurat生的rds文件的话,按照seurat中分的群进行着色,如果不想用ggplot的默认色,可以提供颜色列表col list。
    plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters", cell_size = 0.75) + facet_wrap (~seurat_clusters, nrow =Cluster_num) ###按照seurat中分的群绘制分面图。
    
    plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75) ###按照样本进行着色
    plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)+facet_wrap(~stim, ncol = 2)
    ##绘制样本着色的分面图。
    

    (5) 查看核心基因轨迹

    # 根据 cds子集对象 画伪时间上的基因表达
    cds_subset <- HSMM_filtered[my_genes,]
    plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")
    plot_genes_in_pseudotime(cds_subset, color_by =  "State")
    

    (6) 查看核心基因表达热图

    diff_test_res = differentialGeneTest(cds[marker_genes,], fullModelFormulaStr="~sm.ns(Pseudotime)")
    sig_gene_names = row.names(subset(diff_test_rtes, qval<0.01)) #原文例子是 0.1,是不是太大了?
    
    # 主要是这个基因列表
    plot_pseudotime_heatmap(cds[sig_gene_names,], num_clusters=6, cores=1, show_rownames=T)
    

    热图可以采用top50的基因进行展示,当然,也可以自己筛选。
    一行为每个基因对应的表达量。一列为一个拟时间点,一个拟时间值可能代表多个细胞。

    4. BEAM 分析: 单细胞轨迹的“分支”分析

    # 查看分支点:就是黑色圆圈中的白色数字
    plot_cell_trajectory(cds, color_by = "State")
    
    # BEAM进行统计分析:我们分析branch_point = 1这个分支处的细胞命名分叉是如何进行的。
    BEAM_res = BEAM(cds[expressed_genes, ], branch_point=1, cores=4)
    
    # 对影响分析的基因根据qvalue进行排序
    BEAM_res=BEAM_res[order(BEAM_res$qval),]
    BEAM_res=BEAM_res[, c("gene_short_name", "pval", "qval")]
    
    
    # 绘制分支热图
    my_branched_heatmap = plot_genes_branched_heatmap(
        cds[row.names(subset(BEAM_res, qval<1e-4)),], 
        branch_point=1,
        num_clusters=4, cores=2, 
        use_gene_short_name=T,
        show_rownames=T)
    

    软件作者是这么说的:Cell fate 1 corresponds to the state with small id while cell fate 2 corresponds to state with bigger id


    能自定义颜色吗?如果基因过多,就别不显示了,看不清。

    tmp1=plot_genes_branched_heatmap(test[row.names(subset(BEAM_res,qval<1e-4)),],
        branch_point = 1,
        num_clusters = 4, #这些基因被分成几个group
        cores = 1,
        branch_labels = c("Cell fate 1", "Cell fate 2"), #规定标签
        #hmcols = NULL, #默认值
        hmcols = colorRampPalette(rev(brewer.pal(9, "PRGn")))(62), #自定义热图颜色
        branch_colors = c("#979797", "#F05662", "#7990C8"), #pre-branch, Cell fate 1, Cell fate 2分别用什么颜色
        use_gene_short_name = T,
        show_rownames = F,
        return_heatmap = T #是否返回一些重要信息
    )
    
    pdf("branched_heatmap.pdf",width = 5,height = 6)
    tmp1$ph_res
    dev.off()
    

    返回的列表tmp1包含热图对象,热图的数值矩阵,热图主题颜色,每一行的注释(基因属于哪一个group)等信息。
    根据行注释提取出基因group,可以去做富集分析,后期加到热图的旁边。像这样:
    https://www.cnblogs.com/TOP-Bio/p/14856979.html

    gene_group=tmp1$annotation_row
    gene_group$gene=rownames(gene_group)
    
    library(clusterProfiler)
    library(org.Hs.eg.db)
    allcluster_go=data.frame()
    for (i in unique(gene_group$Cluster)) {
      small_gene_group=filter(gene_group,gene_group$Cluster==i)
      df_name=bitr(small_gene_group$gene, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")
      go <- enrichGO(gene         = unique(df_name$ENTREZID),
                     OrgDb         = org.Hs.eg.db,
                     keyType       = 'ENTREZID',
                     ont           = "BP",
                     pAdjustMethod = "BH",
                     pvalueCutoff  = 0.05,
                     qvalueCutoff  = 0.2,
                     readable      = TRUE)
      go_res=go@result
      if (dim(go_res)[1] != 0) {
        go_res$cluster=i
        allcluster_go=rbind(allcluster_go,go_res)
      }
    }
    head(allcluster_go[,c("ID","Description","qvalue","cluster")])
    

    (2) 单个基因的伪时间曲线,类似 plot_genes_in_pseudotime(),不同的是该函数显示2个发育路径。

    展示具体的基因,这些基因可以是:

    • 随拟时序、state而改变的基因(第5步)
    • 细胞亚群的marker基因(seurat得到的差异基因)
    • 分支点分析得到的分支特异的基因(第6步BEAM函数得到的基因)

    We can plot a couple of these genes, such as Pdpn and Sftpb (both known markers of fate in this system), using the plot_genes_branched_pseudotime function, which works a lot like the plot_genes_in_pseudotime function, except it shows two kinetic trends, one for each lineage, instead of one.

    另外,可以绘制基因在所有细胞中的表达量在不同分支上的表达量情况。我们可以根据分叉的state位置来判定这个基因对分支的影响。在拟时间分支树构建的过程中,其实是有很多细小的分支,最终被拟合成一条直线,便于可视化的友好性。下图中的Branch表示挑选了两个不同的未被拟合的实际分支。

    plot_genes_branched_pseudotime(
        HSMM[rownames(subset(BEAM_res, qval<1e-8)),], 
        branch_point=1,
        color_by="orig.ident",
        ncol=3)
    
    > str(HSMM)
    
    
    
    # 示例2
    test_genes=c("TFF3","GUCA2B")
    pdf("genes_branched_pseudotime.pdf",width = 9,height = 4)
    plot_genes_branched_pseudotime(test[test_genes,],
                                   branch_point = 1,
                                   color_by = "celltype",
                                   cell_size=2,
                                   ncol = 2)
    dev.off()
    

    参考资料

    相关文章

      网友评论

          本文标题:monocle2 使用笔记

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