R语言中实现的包：Rtsne
来学习下核心函数 Rtsne() 的主要参数：
下面以常规模型构建过程中产生的risk和riskScore数据为例：
下面结合ggplot2进行绘图：
另一种降维方法：PCA，已在之前的笔记中写过。 【R>>PCA】主成分分析
参考链接：
t-SNE一种高效的降维算法
下游分析
cellranger count 计算的结果只能作为错略观测的结果，如果需要进一步分析聚类细胞，还需要进行下游分析，这里使用官方推荐 R 包（Seurat 3.0）
流程参考官方外周血分析标准流程（ https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html ）
Rstudio操作过程：
## 安装seurat
install.packages('Seurat')
## 载入seurat包
library(dplyr)
library(Seurat)
## 读入pbmc数据（文件夹路径不能包含中文，注意“/“的方向不能错误，这里读取的是10x处理的文件，也可以处理其它矩阵文件，具体怎样操作现在还不知道，文件夹中的3个文件分别是：barcodes.tsv，genes.tsv，matrix.mtx，文件的名字不能错，否则读取不到）
pbmc.data <- Read10X(data.dir = "D:/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
## 查看稀疏矩阵的维度,即基因数和细胞数
dim(pbmc.data)
pbmc.data[1:10,1:6]
## 创建Seurat对象与数据过滤，除去一些质量差的细胞（这里读取的是单细胞 count 结果中的矩阵目录；在对象生成的过程中，做了初步的过滤；留下所有在>=3 个细胞中表达的基因 min.cells = 3；留下所有检测到>=200 个基因的细胞 min.genes = 200。)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
##计算每个细胞的线粒体基因转录本数的百分比（%）,使用[[ ]] 操作符存放到metadata中，mit-开头的为线粒体基因
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
##展示基因及线粒体百分比（这里将其进行标记并统计其分布频率，"nFeature_RNA"为基因数，"nCount_RNA"为细胞数，"percent.mt"为线粒体占比）
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
## 过滤细胞：根据上面小提琴图中基因数"nFeature_RNA"和线粒体数"percent.mt"，分别设置过滤参数，这里基因数 200-2500，线粒体百分比为小于 5%，保留gene数大于200小于2500的细胞；目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据；而保留线粒体基因的转录本数低于5%的细胞,为了过滤掉死细胞等低质量的细胞数据。
pbmc <- subset(pbmc, subset = nFeature_RNA >200 &nFeature_RNA <2500 &percent.mt <5)
## 表达量数据标准化,LogNormalize的算法：A = log( 1 + ( UMIA ÷ UMITotal ) × 10000
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#pbmc <- NormalizeData(pbmc) 或者用默认的
## 鉴定表达高变基因(2000个）,用于下游分析,如PCA；
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
## 提取表达量变化最高的10个基因；
top10 <- head(VariableFeatures(pbmc), 10)
top10
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))
plot1<-VariableFeaturePlot(object=pbmc)
plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE)
CombinePlots(plots=list(plot1,plot2))
## PCA分析：
# PCA分析数据准备,使用ScaleData()进行数据归一化；默认只是标准化高变基因（2000个）,速度更快,不影响PCA和分群,但影响热图的绘制。
#pbmc <- ScaleData(pbmc,vars.to.regress ="percent.mt")
## 而对所有基因进行标准化的方法如下：
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
## 线性降维（PCA）,默认用高变基因集,但也可通过features参数自己指定；
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
## 展示 pca 结果（最简单的方法）
DimPlot(object=pbmc,reduction="pca")
## 检查PCA分群结果, 这里只展示前5个PC,每个PC只显示5个基因；
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
##PC_ 1
##Positive: RPS27, MALAT1, RPS6, RPS12, RPL13
##Negative: CSTA, FCN1, CST3, LYZ, LGALS2
##PC_ 2
##Positive: NKG7, GZMA, CST7, KLRD1, CCL5
##Negative: RPL34, RPL32, RPL13, RPL39, LTB
##PC_ 3
##Positive: MS4A1, CD79A, BANK1, IGHD, CD79B
##Negative: IL7R, RPL34, S100A12, VCAN, AIF1
##PC_ 4
##Positive: RPS18, RPL39, RPS27, MALAT1, RPS8
##Negative: PPBP, PF4, GNG11, SDPR, TUBB1
##PC_ 5
##Positive: PLD4, FCER1A, LILRA4, SERPINF1, LRRC26
##Negative: MS4A1, CD79A, LINC00926, IGHD, FCER2
## 展示主成分基因分值
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
## 绘制pca散点图
DimPlot(pbmc, reduction = "pca")
## 画第1个或15个主成分的热图；
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
## 确定数据集的分群个数
# 鉴定数据集的可用维度，方法1：Jackstraw置换检验算法；重复取样（原数据的1%）,重跑PCA,鉴定p-value较小的PC；计算‘null distribution’(即零假设成立时)时的基因scores。虚线以上的为可用维度，也可以调整 dims 参数，画出所有 pca 查看。
#pbmc <- JackStraw(pbmc, num.replicate = 100)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
#JackStrawPlot(pbmc, dims = 1:15)
# 方法2：肘部图（碎石图）,基于每个主成分对方差解释率的排名。
ElbowPlot(pbmc)
## 细胞聚类：分群个数这里选择10,建议尝试选择多个主成分个数做下游分析,对整体影响不大；在选择此参数时,建议选择偏高的数字（为了获取更多的稀有分群,“宁滥勿缺”）；有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。
## 非线性降维（UMAP/tSNE)基于PCA空间中的欧氏距离计算nearest neighbor graph,优化任意两个细胞间的距离权重（输入上一步得到的PC维数） 。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
## 接着优化模型,resolution参数决定下游聚类分析得到的分群数,对于3K左右的细胞,设为0.4-1.2 能得到较好的结果(官方说明)；如果数据量增大,该参数也应该适当增大。
pbmc <- FindClusters(pbmc, resolution = 0.5)
## 使用Idents（）函数可查看不同细胞的分群；
head(Idents(pbmc), 5)
## 结果：AAACCTGAGGTGCTAG AAACCTGCAGGTCCAC AAACCTGCATGGAATA AAACCTGCATGGTAGG AAACCTGCATTGGCGC
1 3 0 10 2
Levels: 0 1 2 3 4 5 6 7 8 9 10 11
## Seurat提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起）,比如UMAP和t-SNE,运行UMAP需要先安装'umap-learn'包，这里不做介绍，两种方法都可以使用，但不要混用，如果混用，后面的结算结果会将先前的聚类覆盖掉，只能保留一个。
## 这里采用基于TSNE的聚类方法。
pbmc <- RunTSNE(pbmc, dims = 1:10)
## 用DimPlot()函数绘制散点图,reduction = "tsne",指定绘制类型；如果不指定,默认先从搜索 umap,然后 tsne, 再然后 pca；也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot()； cols,pt.size分别调整分组颜色和点的大小；
DimPlot(pbmc,reduction = "tsne",label = TRUE,pt.size = 1.5)
## 这里采用基于图论的聚类方法
pbmc<-RunUMAP(object=pbmc,dims=1:10)
DimPlot(object=pbmc,reduction="umap")
## 细胞周期归类
pbmc<- CellCycleScoring(object = pbmc, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes)
head(x = [email protected])
DimPlot(pbmc,reduction = "tsne",label = TRUE,group.by="Phase",pt.size = 1.5)
## 存储结果
saveRDS(pbmc, file = "D:/pbmc_tutorial.rds")
save(pbmc,file="D:/res0.5.Robj")
## 寻找cluster 1的marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
## 结果: p_val avg_logFC pct.1 pct.2 p_val_adj
MT-CO1 0.000000e+00 -0.6977083 0.985 0.996 0.000000e+00
RPS27 2.182766e-282 0.3076454 1.000 0.999 3.480202e-278
MT-CO3 2.146399e-274 -0.4866429 0.995 0.997 3.422218e-270
DUSP1 2.080878e-247 -1.7621662 0.376 0.745 3.317752e-243
RPL34 8.647733e-244 0.3367755 1.000 0.997 1.378795e-239
##寻找每一cluster的marker
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_logFC)
# A tibble: 24 x 7
# Groups: cluster [12]
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl><dbl><dbl> <dbl><fct> <chr>
1 2.29e-123 0.636 0.344 0.097 3.65e-119 0 CD8B
2 7.62e-113 0.487 0.632 0.305 1.22e-108 0 LEF1
3 2.04e- 74 0.483 0.562 0.328 3.25e- 70 1 LEF1
4 1.39e- 61 0.462 0.598 0.39 2.22e- 57 1 ITM2A
5 0. 2.69 0.972 0.483 0. 2 GNLY
6 0. 2.40 0.964 0.164 0. 2 GZMB
7 1.31e-121 0.768 0.913 0.671 2.09e-117 3 JUNB
8 2.06e- 94 0.946 0.426 0.155 3.28e- 90 3 RGS1
9 2.05e-255 1.57 0.586 0.09 3.27e-251 4 GZMK
10 2.94e-140 1.57 0.69 0.253 4.68e-136 4 KLRB1
# ... with 14 more rows
## 存储marker
write.table(pbmc.markers,file="D:/allmarker.txt")
## 各种绘图
## 绘制Marker 基因的tsne图
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"),cols = c("gray", "red"))
## 绘制Marker 基因的小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
## 绘制分cluster的热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
剩下的便是寻找基因 marker 并对细胞类型进行注释(见下回分解)
近年来，由于细胞的异质性及发育分化等相关的问题越来越被研究者们所关注，单细胞转录组分析为研究异质细胞群的复杂生物学过程提供了方法和工具。每一个细胞进行转录组测序时就是细胞发育过程中的快照，单细胞拟时间分析软件Monocle2是基于R语言的安装包，其功能基于单细胞转录组的表达矩阵，通过无监督学习（Reversed Graph Embedding算法）的方式将细胞置于发育轨迹的不同分支上，从而模拟细胞群体生物学过程。也就是我们经常说的拟时序（pseudotime）分析，又称细胞轨迹（cell trajectory）分析。通过拟时分析可以推断出发育过程细胞的分化轨迹或细胞亚型的演化过程，在发育相关研究中使用频率较高。模拟细胞的分化轨迹的软件，最常用的软件为Monocle2。该算法不仅能模拟细胞的发育轨迹，同时也能对细胞进行聚类（t-SNE）。通过聚类获得不同状态下的差异基因，分析影响分支形成的关键基因及其功能，对研究相关生物学问题有指导性的作用。
Monocle2主要基于关键基因的表达模式，通过学习每个细胞必须经历的基因表达变化的序列，根据拟时间值中对单个细胞进行排序，模拟出时间发育过程的动态变化。而这个排序技术表现是一种在低维空间排布高维数据的降维技术。（具体描叙请参考： https://www.meiwen.com.cn/subject/zlgsvqtx.html ），目前monocle已经升级至monocle3，但在结果分析上monocle3的可读性不如monocle2，因此本研究中主推monocle2这个版本。
在做拟时序分析之前，先要明白几个研究目的：
先保证服务器上装有R ，我这里安装的版本是version 3.5.1。
对于输入文件，有3种类型的格式：
一种是单细胞运行获得的3个结果文件。格式如下表所示：
准备需要进行拟时间分析数据的三个文件：
1)表达量文件：exprs：基因在所有细胞中的count值矩阵。
格式示例：
注意：该文件的表头必须以“%%**\n%”的形式出现，否则就会出现报错。
2)表型文件 phenoData:
格式示例，该文件即为每个细胞的barcode信息。：
由于seurat软件包的升级，大量的函数命名发生了变化，对2和3这两种类型的版本的操作也是不一样。
Seurat2.4的版本：
Seurat3.0及更高版本：
如果你没有3个标准格式的文件，也没有生成rds文件，仅仅只有表达量矩阵，具体原理和有rds文件一致，可按照如下方法进行处理。
如果你采用的是第二种输入数据进行分析，即使你原先的Seurat对象已经归一化过了，官方推荐在monocle中重新进行归一化处理。
在进行降维聚类之前，我们先获得高表达的基因集作为每个聚类的用来order的Feature基因。当然我们可以使用所有的基因，但一些表达量特别低的基因提供的聚类信号往往会制造分析噪音，Feature基因的选择性很多，一种是我们可以根据基因的平均表达水平来进行筛选，另外我们也可以选择细胞间异常变异的基因。这些基因往往能较好地反映不同细胞的状态。以平均表达量高于0.1来进行筛选。
绘制用于order基因和非order的平均表达量的分布点图。这里不做展示。
对细胞进行聚类，在Seurat中采用的是分辨率来确定cluster的数目。而monocle中可以直接指定聚类数目。主要指出的这里所聚类获得的cluster数目比我们赋值的要少一个。即当num_clusters=3时，你只获得了2个cluster。具体解释有些难懂，在这里不做过多的解释。
拟时间分析不仅是要对分析的细胞群进行排序，还要获得觉得细胞排序的关键基因集。这种基因集有两种情况，在有先验知识的情况下，我们可以从系统生物学的角度获得一系列与细胞发育相关的基因集，从而对细胞进行排序，这种方式是最为保险的，但局限性是对未知的发育情况不能进行分析。另外一种情况就是使用无监督聚类方式计算关键基因集。接下来我们采用differentialGeneTest方式获得clustering_DEG_genes（与聚类相关的差异基因集）
上一步过程是对所有的细胞进行无监督训练，运行时间与细胞数和基因数相关，一般会花很长的时间。可以根据cores的数目进行并行。
differentialGeneTest这个函数测试每个基因的差异表达，取决于伪时间或根据指定的其他协变量。 “ differentialGeneTest”是Monocle的主要差异分析常规， 它接受一个CellDataSet和两个模型公式作为输入，指定由实现的广义谱系模型“ VGAM”包。 也就是说我们可以根据指定’~cluster’或者拟时间值来获得差异基因。差异基因的结果如下表所示：
在这个表格中，我们先看一下表头对应的关键列。第一列没有列名，为基因名称。pval,qval 为差异基因的显著检验打分。num_cells_expressed为表达这个基因的细胞数。use_for_ordering该基因是否是作为排序的差异基因。
根据排序好细胞进行结果可视化。
命令行如下所示：
上述部分结果如下图所示（不包括分面图）
上图表示在主成分中的细胞聚类分布的tsne图。不同颜色代表细胞群的不同细胞命运的分支。
上图表示依据seurat的cluster ID映射到拟时间的排序上。
接下来可视化展示细胞state相关的差异基因的表达量分布情况，可以根据这些基因来确定细胞的发育方向，下图仅展示qvalue值排在前6个基因。横纵左边意义见坐标轴的描述。
在monocle中的分析中发现，细胞群能从一个轨迹分叉成不同的分支点，在科研中，我们会比较关注发生分支的原因是什么，即分析分支点之间的差异。Monocle采用分支表达式分析建模，主要是BEAM函数，可以将分叉过程重构为一个分支轨迹，从而分析不同细胞命运下的差异。
命令行如下：
举个例子：我们分析branch_point = 1这个分支处的细胞命名分叉是如何进行的。
即下图中所示。下图只有1个分支点，即分析state1,2,3 这三个state的差异。
对影响分析的基因根据qvalue进行排序。
绘制与分支相关的基因热图。
关键参数为subset (BEAM_res, qval <1e-4))，挑选基因进行热图绘制，也可以设置成其他的阈值。
branch_point=1，分支点选为1。num_clusters=4，将基因根据表达相似性分成4个模块。结果如下图所示：
