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

image.png
已经被我下载到了，原来的网页目录在https://www.sciencedirect.com/science/article/pii/S2211124719300634?via%3Dihub
但是数据在GEO数据库中，因为前面的project的细胞注释中已经讲述过哪里下载数据，此处不在赘述，，，，因为这篇文章中后续使用了拟时间序列分析，所以在此处尝试一下，以后还有的单细胞富集分析，和`蛋白互作网络、受体 - 配体信号相互作用等等还需好好学习s（这些以后可以找点文章再去做做。
此外这篇文章简书进行了很好的描述
https://www.jianshu.com/p/8549dbc6d4e4
对于分化发育等细胞，可以利用monocle包进行时间序列分析。根据基因表达情况推测细胞当前所处的每个分化阶段位置。

1.新建文件夹

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

image.png
但是突然发现自己直接做时间序列分析就行，，教程先用的
https://www.jianshu.com/p/e79ab1cc0a67
后来发现要用的是10X，所以重新mv一下数
image.png

1、输入文件

image.png
教程是这样的，但是我却报错，Error in load_cellranger_matrix(D2) : could not find function "load_cellranger_matrix",然后去查了这个函数
后来找到一个网页教程也很好https://www.haomeiwen.com/subject/yiigbctx.html

生成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
新教程链接https://mp.weixin.qq.com/s?__biz=MzI1Njk4ODE0MQ==&mid=2247488352&idx=2&sn=c08863668cbde5d66fd5d3e592f3aac9&chksm=ea1f15e2dd689cf4c0fb80add7625bc2c60d074475491e8116131ae7368b52fd5f451597d553&mpshare=1&scene=1&srcid=09256G4F36KTNXpUGJVIZ8KE&sharer_sharetime=1600998206903&sharer_shareid=aac470b1ad84e1cbcb4ee2353b49492c&key=80efbdacc232da71ba8c5e67343ce054031a98182bc589275e939669132b557258112b1062758dce3a36eca03979183ab62c973da6e391ced5fb8defacf529da97fce5f9d73782b6d1f4d629bd93fbc02f577cc83c38be9dff3e711f8f27a84f796dbf150dcdbfdb008d19af66a378d3967c3e1486f7d9d342206e92ac04ce9c&ascene=1&uin=MjU0MDE0Njc5&devicetype=Windows+10+x64&version=62090529&lang=zh_CN&exportkey=Ay2euVhX0Yo0Ilzg28PK4hw%3D&pass_ticket=W%2B7dpO6yPmBPi6BSUQks6DDJw0deCcUhtCSocFC7Pqr%2FbTCA3Ab%2BZ89q7Sv%2FCOPb&wx_header=0

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