单细胞数据挖掘

目录

1. 实验原理
   1. 数据集
   2. 质控
   3. 挑选高变异基因
   4. 降维与聚类
   5. 细胞类型注释
   6. 轨迹分析
2. R语言实现
   1. 提取数据
   2. 质控
      1. 创建Seurat对象，简单筛选细胞
      2. 筛选线粒体基因和ERCC基因过多细胞
      3. 筛选双细胞
      4. 周期判断
      5. 质控绘图
   3. 挑选高变异基因
   4. 降维与聚类，主要是找marker基因
      1. PCA
      2. 聚类还有一堆东西
   5. QC之二
   6. 细胞类型注释
   7. 轨迹分析
3. 总结

案例里用了《Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients》，主要是使用integrated scRNA-seq（综合（？）单细胞RNA测序）和bulk RNA-seq（整体（？）RNA测序）来找出GBM分化相关基因然后预测患者总生存期之类的。

换一个类似的做实验，也不懂各种数据有什么差别，就选《Cell differentiation trajectory predicts patient potential immunotherapy response and prognosis in gastric cancer》这篇文章，主要是使用integrated scRNA-seq和bulk RNA-seq来找出GC分化相关基因然后预测患者总生存期之类的。（干的事都差不多，真是原原又神神啊）
这个用的数据集是GSE112302和GSE84437，和案例类似的部分用的是GSE112302。

实验原理

主要就是干了质控、基因可视化、降维与聚类、周期判断、double droplet、细胞类型注释、轨迹分析。

数据集

这堆数据集就是一轴是细胞或一坨样本，一轴是基因种类，每个点代表某个细胞/样本的某个基因表达量。

质控

质控主要就是过滤：

• 低质量的细胞(比如基因少、线粒体基因多)
• 双细胞（两个细胞被捕获到同一个液滴中）
• 高ERCC基因（External RNA Controls Consortium），“是一组合成的RNA分子，它们的浓度已知，并且在测序实验中添加到样本中，以便可以用来校正技术变异和进行数据标准化”，在所有单元格中的浓度应该是一致的，所以如果ERCC基因百分比高，那这个单元格的RNA可能较少。

Seurat包可以傻瓜实现。

案例里还提到了double droplet和周期判断。
double droplet就是上面说的双细胞。
而周期判断就是说，一些细胞处于分裂期，它们会高表达某些基因，可能会影响聚类，所以要去掉。
scran包可以傻瓜实现。

挑选高变异基因

highly Variable gene是在一组样本或细胞中，表达水平变化大的基因，它们可能是驱动细胞异质性的关键。 Seurat包同样可以傻瓜实现。

降维与聚类

降维机器学习里也有，“将数据从原始的高维空间转换到低维空间，同时尽可能保留数据的变异性”，总的来说目的就是去掉冗杂信息。一个方法就是PCA（主成分分析），能使得数据的方差尽可能大（保留变异性），而均方误差尽可能小（和原数据的差距尽可能小），欸网上看证明去吧。
PCA对数据的尺度（即每个特征的平均值和标准差）非常敏感，如果不同的feature有不同的尺度（0-1vs0-1000），那么PCA的结果将主要受到尺度大的特征的影响，所以PCA之前要标准化。

聚类就是把数据分成几类，使得同一个类内数据的相似性尽可能大，同时不同类数据的差异性也尽可能地大。欸还是无监督学习，感兴趣看看吧。

案例里PCA之后聚类然后又用t-SNE降维，有人说是因为先用PCA去除噪声，利于准确聚类，然后用t-SNE更好地可视化（“在低维空间中保留高维空间中的邻域结构，可以有效地揭示数据中的群体结构和模式” ）。
另外PCA是线性的，而t-SNE是非线性的。

Seurat包还是可以傻瓜实现。that's right the square hole

细胞类型注释

字面意思，就比如把一些细胞标注成B细胞之类的甚至更精细一些，幼稚B细胞之类的。这靠的是一些特异性高的marker基因。 可以用SingleR包。

轨迹分析

就是“构建细胞的变化轨迹，重塑细胞随时间的变化过程”，目的是绘制祖细胞到分化细胞的进展情况。
可以用Monocle包。

---

本来这点东西一天下班时间就能搞完的，现在下来两三天了，不光是因为只有业余时间，还真是怠惰啊。

---

R语言实现

代码详见single-cell-data-mining。

提取数据

怎么两个数据集形式差那么多，你们能不能统一一下（骂骂咧咧）。

先看GSE84465学习一下。
案例里的GSE84465集合到了一个.csv里，列名是某样本的plate_id（微孔板）和Well（孔）（我还没做过湿实验啊），行名就是基因。
然后还有一个SraRunTable.txt，刚刚还奇怪，上面的文件都没有分组信息，这个文件里就有，就是每个样本的一些元信息，比如组织类型、患者编号、对应的样本编号等。这个就是SRA里的元数据，进SRA Run Selector里就能找到。
所以这个SraRunTable.txt就是用来给上面那些sc RNA-seq数据分组。

我们要做的就是从SRA生成GSE对应的列名，并使用SRA的信息筛选出肿瘤组织，借助对应的列名筛选出GSE中的肿瘤组织。
还有SRA记录着分组信息，也要一起筛选出来。

所以说拿到GSE112302也是一样的，就是GSE112302没有给全部合在一起的数据，需要我们自己合起来（附加列，也就是样本），然后用SRA筛选出肿瘤组织。

质控

创建Seurat对象，简单筛选细胞

要先创建Seurat对象。CreateSeuratObject()函数的文档如下：

• counts：这是一个矩阵类对象，包含未经标准化的数据，其中单元格作为列，特征作为行，或者是一个Assay派生对象。
• assay：初始测定的名称。
• meta.data：要添加到Seurat对象的额外单元格级元数据。应该是一个数据框，其中行是单元格名称，列是额外的元数据字段。元数据中的行名称需要与计数矩阵的列名称匹配。
• project：Seurat对象的项目名称。
• ...：传递给其他方法的参数。
• min.cells：至少在这么多细胞中检测到的特征。也会对计数矩阵进行子集。要重新引入排除的特征，使用较低的截止值创建一个新对象。
• min.features：至少在这么多特征中检测到的细胞。

（对了R语言的函数是这样调用的：function_name(arg1 = value1, arg2 = value2, ...)）

其实只要用到counts、meta.data、min.cells、min.features就行了。 这里counts就是数据，meta.data就是分组信息，用data.frame()上面说了是这样的：...是列（额外的元数据字段，这里就是病人编号，要用来分组嘛），row.names是行（样本编号）,要和counts的列名对应。
然后后两个前者用来筛选表达的细胞的数量低的基因，后者用来筛选基因少的细胞，顺便还能过滤空细胞。

筛选线粒体基因和ERCC基因过多细胞

然后还要筛选掉线粒体基因过多的细胞：

scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-") #先全部算了比较好，PercentageFeatureSet()计算满足给定模式的特征在每个单元格中的百分比，后面那个参数就是正则（真是先进啊（大嘘））
scRNA <- subset(scRNA, subset = percent.mt < 5) #案例取5%，subset()就是Array.filter()（确信）

还有ERCC基因一样的，选40%。

筛选双细胞

才发现这个要等到先聚类...别急。

周期判断

这个要先标准化...别急...看seurat#7942。 我说怎么把这两个放在后面。

db.test <- scDblFinder::findDoubletClusters(GetAssayData(scRNA,slot="counts",assay="RNA"), clusters=scRNA@meta.data$seurat_clusters)

质控绘图

最后就是绘图。
看他们都是 - 图A：观察不同组cell的counts、feature分布 - 图B：nCount_RNA与对应的nFeature_RNA关系

图A用小提琴图，复习一下，VlnPlot()：

• object：Seurat对象。
• features：要绘制的特征。
• group.by：用于分组的元数据。
• pt.size：点的大小。
• cols：颜色。可以用rainbow(n)生成。

然后可以用+连接修改，theme(legend.position = "none")隐藏图例，labs(tag = "A")添加全局标签A。

图B用FeatureScatter()，绘制两个特征之间的关系。

• object：Seurat对象。
• feature1：x轴特征。
• feature2：y轴特征。
• group.by：用于分组的元数据。
• 剩下通用的都差不多。

案例还对ERCC基因比例和nCount用了FeatureScatter()，这个图一看也验证了上面说的，ERCC基因百分比高，nCount就低。

挑选高变异基因

用FindVariableFeatures()：

• object：Seurat对象。
• selection.method：选择方法。可以是"vst"（variance stabilizing transformation）或"mean.var.plot"（均值方差图）。
• nfeatures：选择为顶部变异特征的特征数量；只在selection.method设置为'dispersion'或'vst'时使用。
• 还有一些其它的，还有使用vst法时的参数。

然后用VariableFeaturePlot()绘制变异特征图。
LabelPoints()可以标记点，head()就可以取前几个。

降维与聚类，主要是找marker基因

PCA

要先标准化，这里用ScaleData()，对每个基因的表达量进行缩放和中心化，可以消除不同基因之间由于表达量级别不同而产生的差异。而不止用NormalizeData()，因为它对每个细胞之间的某个基因表达量进行归一化，可以消除不同细胞之间由于测序深度不同而产生的差异，所以归一化后，不同的基因可能仍然有不同的表达量级别。

然后就是简单的scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA))，绘图可以用DimPlot()。

最后挑选主成分，jackStraw()评估每个基因的PCA得分的显著性（对于每个主成分，它会对数据进行PCA，然后对PCA的结果进行随机重排，然后再次进行PCA，然后计算原始PCA（实际）和重排PCA（随机下理论）的差异，这个差异越大，说明这个主成分越重要。），ScoreJackStraw用于分析jackStraw函数的结果，JackStrawPlot()绘制jackStraw图。

这个图，empirical是实际的PCA得分，theoretical是随机重排的PCA得分（这两个对应彩色点），虚线是"均匀分布的 p 值的累积分布函数"，如果主成分轴显着，则二者应显示出差异。不是很清楚，反正选不一样的就完了。
这里结果都挺好的，盗用一下某人的例子，像11、12这种就不应该选了： 《JackStraw Plot (QQ-plot) を理解する》

聚类还有一堆东西

烦了，调包就完了。

这里结果可能不太一样（甚至教程的复现和原文不一样），数据集应该是没变，重复运行一样应该也不是随机过程的问题，暂且按下不表就当是版本差异。

FindAllMarkers()找marker基因，真慢啊，有必要换RRO或者python。
完了再筛选出P<0.05的，还有在不同类型细胞中表达量差异大的基因。

QC之二

刚才说的要先聚类或标准化的筛选双细胞和周期判断。

emptyDrops()检测空液滴，而检测双细胞可以用doubletFinder()，案例里的doubletCluster()其实寻找的是表达介于两个Cluster之间的Cluster，但他解释的是“有无double droplet聚在一起的类”（而且好像已经废弃了）。
CellCycleScoring()评估细胞周期。

细胞类型注释

按下不表。

轨迹分析

选择特定基因代表细胞不同阶段的特征，然后用Monocle包，newCellDataSet()创建数据集，estimateSizeFactors()估计因子，estimateDispersions()估计离散度，reduceDimension()降维，clusterCells()聚类，plot_cell_trajectory()绘制轨迹图。

总结

前面这一半的最终目标就是找出分化相关基因。
抛开具体的实现来说，要做的事就是：

• 质控：提升数据质量
• 降维：减少冗余信息
• 聚类：找出数据的内在结构
• 找marker基因

还顺带做了轨迹分析利于解析分化的调控机制、找分化路径之类的。

有时间再用Python实现一下。 おわり