文章来源于 The Neonatal and Adult Human Testis Defined at the Single-Cell Level

image.png

1.新建文件夹

首先先把之前/data/jiarf/10X/results/outs/filtered_feature_bc_matrix/Integrate/Rdata/artical_data做好的表达矩阵放到新建的那个文件夹中

image.png

image.png

1、输入文件

image.png

生成cellDataSet对象

Monocle把单细胞表达数据存放在CellDataSet对象里。

# 4. cell trajectory----------------------------------------------------------------------
rm(list=ls())
setwd('/data/jiarf/literiture')
library(monocle3)
library(Matrix)
HSMM_expre_matrix <- Matrix::readMM('./D2/GSM3526584_D2Total_matrix.mtx.gz')
HSMM_simple_sheet <- read.delim('./D2/GSM3526584_D2Total_barcodes.tsv.gz',header = F)
HSMM_gene_annotation <- read.delim('./D2/GSM3526584_D2Total_genes.tsv.gz',header = F)
#creat object of cellDataSet
pd <- new('AnnotatedDataFrame',data=HSMM_simple_sheet)
fd <- new('AnnotatedDataFrame',data=HSMM_gene_annotation)
HSMM <- newCellDataSet(as.matrix(HSMM_expre_matrix),
                       phenoData=pd,
                       featureData=fd,
                       expressionFamily = negbinomial.size())

给你的数据选择一个合适的distribution（分布）

image.png

image.png

2、估计size factor和分散度

HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)

3、过滤低质量的细胞

要看多少细胞里表达某个基因，或者看一个细胞里表达多少基因是非常容易的。Monocle提供了一个非常简单的函数来计算：

HSMM <- detectGenes(HSMM,min_expr = 0.1)
print(head(fData(HSMM)))
> HSMM <- detectGenes(HSMM,min_expr = 0.1)
> View(HSMM)
> fData(HSMM)
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘fData’ for signature ‘"CellDataSet"’

但是这个地方报错了，看了一下简书的教程，也只做了前面那一个，所以也就不纠结了
在CellDataSet对象里，存放每个细胞评分的data在phenoData里。评分属性作为一列。 你可以过滤那些不符合你要求的细胞。你也可以用FastQC来过滤细胞。这类软件可以鉴别RNA-seq文库里严重降解的RNA，也可以鉴别包含有核糖体，线粒体或者其他类型RNA污染的文库。

4、细胞分类和细胞计数

image.png

image.png

image.png

无监督：我们可以根据平均表达水平过滤基因，我们还可以另外选取那些在细胞之间不经常变动的基因。这些基因能够对细胞状态提供很有效的信息。

disp_table <- dispersionTable(HSMM)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) 
HSMM_myo <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
plot_ordering_genes(HSMM)

setOrderingFilter函数会标记基因（在后续clusterCells细胞聚类所用的基因）。 尽管我们也可以提供一个我们自己写的基因列表。plot_ordering_genes功能根据细胞之间的平局表达量来显示一个基因的变异性 。红线显示的是Monocle基于这种关系的预期值。我们把用于细胞分群的基因标记为黑点，其他基因标记为灰点。
可是我的图为啥没有黑点，都是灰的，

image.png

image.png

5.降维(时间还是蛮久的）

HSMM_myo <- reduceDimension(HSMM_myo, 
                            max_components= 2, 
                            reduction_method = 'DDRTree')

排序

HSMM_myo <- orderCells(HSMM_myo)

State轨迹分布图

plot1 <- plot_cell_trajectory(HSMM_myo, color_by = "State")
# ggsave("./State.pdf", plot = plot1, width = 6, height = 5)
ggsave("./State.png", plot = plot1, width = 6, height = 5)

Cluster轨迹分布图

因为我没有做seurat的cluster，所以没有这个图

Pseudotime轨迹图

plot3 <- plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime")
# ggsave("./Pseudotime.pdf", plot = plot3, width = 6, height = 5)
ggsave("./Pseudotime.png", plot = plot3, width = 6, height = 5)

合并作图

library(patchwork)
plotc <- plot1|plot3
# ggsave("./Combination.pdf", plot = plotc, width = 10, height = 3.5)
ggsave("./Combination.png", plot = plotc, width = 10, height = 3.5)

保存结果

write.csv(HSMM_myo@phenoData@data, "./pseudotime.csv")

Combination.png

pseudotime.csv.png

6 轨迹图分面显示

p1 <- plot_cell_trajectory(HSMM_myo, color_by = "State") + facet_wrap(~State, nrow = 1)
# p2 <- plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters") + facet_wrap(~seurat_clusters, nrow = 1)
# plotc <- p1/p2
ggsave("./trajectory_facet.png", plot = p1, width = 6, height = 5)

trajectory_facet.png

7 Monocle基因可视化

image.png

教程中的代码报错了，很奇怪

Error in orig[[nm]][i, , ..., drop = drop] : subscript out of bounds

以为是基因的问题，换了四个还是不行，说明是别的，，然后去百度谷歌解决
哦对因为我的pdata（）和fData（）函数不能用（我也不知道为啥，就是每次都报错，不认识这俩函数，查的时候发现可以代替HSMM_myo@featureData@data == fData（），，，HSMM_myo@phenoData@data == pdata（））
刚刚是p1报错了，解决如下

 blast_genes <- row.names(subset(HSMM_myo@featureData@data,
                                 V2 %in% c("LURAP1", "FOXE3", "RAD54L")))
p1 <- plot_genes_jitter(HSMM_myo[blast_genes,], grouping = "State", color_by = "State")
ggsave("./genes_visual.png", plot = p1, width = 8, height = 4.5)

genes_visual.png

但是p2又报错了

> p2 <- plot_genes_violin(HSMM_myo[blast_genes,], grouping = "State", color_by = "State")
Error in plot_genes_violin(HSMM_myo[blast_genes, ], grouping = "State",  : 
  unused arguments (grouping = "State", color_by = "State")

找错的过程中又发现了一个教程：https://www.plob.org/article/20919.html
上面那个错误原因是plot_genes_violin函数中早已没有grouping这个参数，所以看了一下

image.png

image.png

> HSMM_SUBSET <- HSMM_myo[row.names(subset(HSMM_myo@featureData@data,
+                  V2 %in% c("LURAP1", "FOXE3", "RAD54L"))),]
> p2 <- plot_genes_violin(HSMM_SUBSET, group_cells_by = "State")
Error: methods::is(object = cds_subset, class2 = "cell_data_set") is not TRUE

还找到了官网教程https://cole-trapnell-lab.github.io/monocle3/docs/differential/
算了找不到解决的办法，老说这个报错
略过，做下一个

拟时相关基因聚类热图

Monocle中differentialGeneTest()函数可以按条件进行差异分析，将相关参数设为fullModelFormulaStr = "~sm.ns(Pseudotime)"时，可以找到与拟时相关的差异基因。我们可以按一定的条件筛选基因后进行差异分析，全部基因都输入会耗费比较长的时间。建议使用cluster差异基因或高变基因输入函数计算。分析结果主要依据qval区分差异的显著性，筛选之后可以用plot_pseudotime_heatmap函数绘制成热图。
我们选用monocle的高可变基因

#高变基因
disp_table <- dispersionTable(HSMM_myo)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
diff_test <- differentialGeneTest(HSMM_myo[disp.genes,], cores = 4, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 1e-04))
p2 = plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,], num_clusters=5,
                             show_rownames=T, return_heatmap=T)
ggsave("./pseudotime_heatmap2.png", plot = p2, width = 5, height = 10)

pseudotime_heatmap2.png

BEAM分析

单细胞轨迹中通常包括分支，它们的出现是因为细胞的表达模式不同。当细胞做出命运选择时，或者遗传、化学或环境扰动时，就会表现出不同的基因表达模式。BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。

mycds_sub <- HSMM_myo[disp.genes,]
plot_cell_trajectory(mycds_sub, color_by = "State")
beam_res <- BEAM(mycds_sub, branch_point = 1, cores = 8)#long time
beam_res <- beam_res[order(beam_res$qval),]
beam_res <- beam_res[,c("V2", "pval", "qval")]#beam_res <- beam_res[,c("gene_short_name", "pval", "qval")]
mycds_sub_beam <- mycds_sub[row.names(subset(beam_res, qval < 1e-4)),]
plot_genes_branched_heatmap(mycds_sub_beam,  branch_point = 1, num_clusters = 3, show_rownames = T)

genes_branched_heatmap.png

end