【单细胞分析】P2.5、聚类，筛选marker基因，可视化

「生信技能树」单细胞数据挖掘_哔哩哔哩_bilibilihttps://www.bilibili.com/video/BV1pa4y1s76J?p=8

 #5.1 聚类
  pc.num=1:20
  #基于PCA数据
  scRNA <- FindNeighbors(scRNA, dims = pc.num) 
  # dims参数，需要指定哪些pc轴用于分析；这里利用上面的分析，选择20
  scRNA <- FindClusters(scRNA, resolution = 0.5)
  table(scRNA@meta.data$seurat_clusters)

这里构建pc.nmu这个数列，相当于选取20个元素进行后续的计算。

Seurat识别细胞类群的原理（FindNeighbors和FindClusters） - 简书众所周知，seurat在降维之后主要依据两个函数来进行细胞分类，这里我们来深入了解一下seurat如何进行细胞分类的。首先我们来看有关分类的两个函数 我们来一一解决其中的问题...https://www.jianshu.com/p/ad6e616db6d6

findneighbors和findclusters是两种计算集合距离的方法，本质是用以判断两个集合的相似性。

dims在这里是选取的细胞数，resolution是分辨率，如果值越高，分的cluster越多，越低则产生的cluster越少。

  # 进行线性降维处理
  scRNA = RunTSNE(scRNA, dims = pc.num)
  DimPlot(scRNA, reduction = "tsne",label=T)
  ?RunTSNE
  p3_1 <- DimPlot(scRNA, reduction = "tsne",label=T) +
    labs(tag = "E")
  p3_1

t-SNE

笔记 | 什么是TSNE - 知乎https://zhuanlan.zhihu.com/p/49073961里面有个英文教程，很适合理解t-SNE的基本原理。

本质是通过t分布转化高维数据的相似性信息并将之映射到低维的数据点中。

 5.2 marker gene

  #5.2 marker gene
  #进行差异分析，一般使用标准化数据
  scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize")
  #结果储存在"data"slot里 scRNA@assays$RNA@来查看这个对象中的内容
  GetAssayData(scRNA,slot="data",assay="RNA")[1:8,1:4]
  
  #if test.use is "negbinom", "poisson", or "DESeq2", slot will be set to "counts
  

这里主要是将每个簇的marker基因筛选出来。

如果是获得的是DESeq2这类的数据，本身其使用的是counts数据，来进行差异分析和聚类。包含一些文库、细胞因子等信息等。

那么这里其实是使用了标准化之后的表达矩阵，因此直接使用了wilcox检验。

 #if test.use is "negbinom", "poisson", or "DESeq2", slot will be set to "counts
  diff.wilcox = FindAllMarkers(scRNA)##默认使用wilcox方法挑选差异基因，大概4-5min
  # 这里有准备好的文件，根据位置可以修改
  # load("../../tmp/diff.wilcox.Rdata")
  head(diff.wilcox)
  dim(diff.wilcox)

得到的这里获得的data.frame就是diff.wilcox应该是每个基因在不同簇之间的差异p值。

基于wilcox进行差异检验获取基因。这里的%>%是管道函数，将前者直接作为后者的第一个元素进行计算。

在all.markers函数中就已经根据筛选要求将这些表达量中的基因筛选了出来。

  library(tidyverse)
  all.markers = diff.wilcox %>% select(gene, everything()) %>%
    subset(p_val<0.05 & abs(diff.wilcox$avg_log2FC) > 0.5)
  #An adjusted P value < 0.05and | log 2 [fold change (FC)] | > 0.5 
  #were considered the 2 cutoff criteria for identifying marker genes.
  dim(all.markers)
  summary(all.markers)

Tip：管道函数

就是把左件的值发送给右件的表达式，并作为右件表达式函数的第一个参数，就是管道函数。

例如：

anscombe_tidy <- anscombe %>%mutate(observation = seq_len(n()))

以上代码等价于：

anscombe_tidy=mutate(anscombe,observation = seq_len(n()))

扩展资料：

1、管道函数的作用

％>％来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤，可以大量减少内存中的对象，节省内存。

符号%>%，这是管道操作，其意思是将%>%左边的对象传递给右边的函数，作为第一个选项的设置（或剩下唯一一个选项的设置）

  top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
  # 所有基因先分组，再根据avg_log2FC进行排序。
  top10
  top10 = CaseMatch(search = as.vector(top10$gene), match = rownames(scRNA)) 
  top10
  length(top10)
  # 检查是否在多个cluster中都是差异基因。
  length(unique(sort(top10))) 
  
  p3_2 <- DoHeatmap(scRNA, features = top10, group.by = "seurat_clusters")
  # 将marker基因差异做个热图看一下每组的分布。
  p3_2
  p3_1 | p3_2 #下图

top_n()是取前多少以及后多少的一个函数，可根据tw进行选择所参考的参数。这里要注意dplyr包对该函数有个更新，可能后续会用新的语法。具体用法可以参见?top_n()

casematch()是将case中作为要给vector再在所有的scRNA中进行match，返回的就是match到的值，本例中就是基因名。

相当于再确认一次，在删除RNA中，并筛选相应的基因名，作为vector进入后续的分析和研究中。

而在length(unique(sort()))中就是将之排序后，再将重复的基因去掉，这样就能看有多少个独特的基因被筛选出来。

因为取了10个top基因，因此应该是获得cluster数*top基因数（10）的总的基因数目，进而再去看，去重的话剩多少个基因数目。

最后将marker基因作为差异进行作图，将结果拼图。

P 2.6 拼图，|为并列，/为换行。其余无殊。