单细胞数据挖掘

生信
5.3k 词

案例里用了《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基因过多细胞

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

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

筛选双细胞

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

周期判断

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

1
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实现一下。 おわり