scCancer包：自动分析肿瘤单细胞转录组利器
分享是一种态度
scCancer包是2020年 清华大学Jin Gu团队开发的一个专门用于分析肿瘤单细胞转录组的工具包。它集成了单细胞数据分析的基本流程(Seurat)，增添若干适合肿瘤数据特征的进阶分析；并且操作简单，两行代码，即可跑完全部分析流程；同时分析结果html格式简洁明了、可读性高，并且提供最后保存的Seurat对象给使用者执行个性化的分析。
官方教程：https://github.com/wguo-research/scCancer 安装命令
checkPkg <- function(pkg){
return(requireNamespace(pkg, quietly = TRUE))
}
if(!checkPkg("BiocManager")) install.packages("BiocManager")
if(!checkPkg("devtools")) install.packages("devtools")
library(devtools)
if(!checkPkg("harmony")) install_github("immunogenomics/harmony")
if(!checkPkg("RcppArmadillo")) install.packages("RcppArmadillo")
if(!checkPkg("RcppProgress")) install.packages("RcppProgress")
if(!checkPkg("NNLM")) install_github("linxihui/NNLM")
if(!checkPkg("liger")) install_github("MacoskoLab/liger")
install_github("wguo-research/scCancer")
1、分析流程（只需要两步）
（1）准备数据[单样本]
示例数据：http://lifeome.net/software/sccancer/KC-example.tar.gz
经cellranger处理得到的10X单细胞表达数据。
sampleFolder即为代表单个样本。raw_feature_bc_matrix与filtered_feature_bc_matrix分别代表初过滤empty droplet前后的单细胞表达数据。
list.files("./data",recursive = T)
#[1] "sample1/filtered_feature_bc_matrix/barcodes.tsv.gz"
#[2] "sample1/filtered_feature_bc_matrix/features.tsv.gz"
#[3] "sample1/filtered_feature_bc_matrix/matrix.mtx.gz"
#[4] "sample1/raw_feature_bc_matrix/barcodes.tsv.gz"
#[5] "sample1/raw_feature_bc_matrix/features.tsv.gz"
#[6] "sample1/raw_feature_bc_matrix/matrix.mtx.gz"
list.files("./results",recursive = T,include.dirs = T)
#[1] "sample1"
（2）实战分析[只需两步]
library(scCancer)
第一步：计算QC过滤指标 runScStatistics
dataPath <- "./data/sample1" # The path of cell ranger processed data
savePath <- "./results/sample1" # A path to save the results
sampleName <- "sample1" # The sample name
authorName <- "Xiaobei" # The author name to mark the report
stat.results <- runScStatistics(
dataPath = dataPath,
savePath = savePath,
sampleName = sampleName,
authorName = authorName
)
# 大概2-3 min
查看分析报表，即计算、可视化的过滤指标
第二步：全部分析流程 runScAnnotation
dataPath <- "./data/sample1" # The path of cell ranger processed data
statPath <- "./results/sample1" # The path of the scStatistics results
savePath <- "./results/sample1" # A path to save the results
sampleName <- "sample1" # The sample name
authorName <- "Xiaobei" # The author name to mark the report
anno.results <- runScAnnotation(
dataPath = dataPath,
statPath = statPath,
savePath = savePath,
authorName = authorName,
sampleName = sampleName,
geneSet.method = "average" # or "GSVA"
)
# 大概25min
查看分析报表，即全部分析结果。根据报表里的解释也能很清楚做了哪些分析，得到哪些结论。下面会有具体解读~
（3）其它输入文件情况
我们处理自己或公共的数据集有时单细胞数据可能并不会提供完全的两套数据，而是
情况1：只有filtered_feature_bc_matrix(~1w cells)
由于raw_feature_bc_matrix只起辅助作用，所以没有太大影响；而且也没有解决办法，如果作者没提供的话。
情况2：只有raw_feature_bc_matrix(~50w cells)
此时我们可以提前加载DropletUtils包，runScStatistics()函数会自动调用该包，过滤raw_feature_bc_matrix里的empty droplet，再进行计算过滤指标操作。
情况3：只有单细胞表达矩阵(一般是过滤empty droplet后的count矩阵)信息
我们可以使用scCancer包提供的generate10Xdata()函数根据表达矩阵生成filtered_feature_bc_matrix文件夹，然后进行上述的分析。
情况4：多个样本的合并
首先对每个样本进行上述的两部分析；然后使用scCombination()函数进行合并校正，涉及多种校正算法。示例代码如下
raw：直接合并，不做校正SeuratMNN：调用Seurat包(V3)的MNN合并算法NormalMNN：考虑到肿瘤恶性细胞异质性，故结合runScAnnotation()注释结果，仅对非肿瘤恶性细胞进行MNN合并校正【默认校正方法】Regression：对假设的系统性的批次效应进行回归校正LIGER：调用liger包的合并算法Harmony：调用harmony包的合并算法
#因为暂时没有合适示例数据，先假设有sample1，sample2，sample3三个样本
# The paths of all sample's "runScAnnotation" results
single.savePaths <- c("./results/sample1", "./results/sample1", "./results/sample1")
sampleNames <- c("sample1", "sample2", "sample3") # The labels for all samples
savePath <- "./results/sample123-comb" # A path to save the results
combName <- "sample123-comb" # A label of the combined samples
authorName <- "Xiaobei" # The author name to mark the report
comb.method <- "NormalMNN" # Integration methods ("NormalMNN", "SeuratMNN", "Harmony", "Raw", "Regression", "LIGER")
comb.results <- runScCombination(
single.savePaths = single.savePaths,
sampleNames = sampleNames,
savePath = savePath,
combName = combName,
authorName = authorName,
comb.method = comb.method
)
2、报表结果解读
2.1 runScStatistics报表解读
这一步如上所说主要计算了低质量细胞的过滤指标，产生了cell.QC.thres.txt筛选阈值结果文件。过滤细胞主要以下两个层面，计算思路均基于outliner。
（1）UMI：细胞当前状态捕获的cDNA总数，即细胞文库大小
首先UMI过低的细胞，会被cellranger算法识别为empty droplet过滤掉(raw→empty，不需要我们计算)
其次UMI过高，则可能是混合了多个细胞的测序结果。类似的，表达基因种类数过低则表示捕获cDNA效率低，不能反映细胞的真实状态；过高则表示潜在的多细胞混杂。
最后是根据特征基因的含量 高线粒体/核糖体基因表达比例该细胞在分离、被捕获的过程中出现了溶解、老化等不正常状态，需要去除。
此外作者考虑到肿瘤细胞在分离过程中，容易产生 stress response。为避免干扰，所以需要去除 dissociation-associated gene异常高表达的细胞
综上基于outliner角度，计算了6个过滤细胞的指标；如若按照此阈值的过滤流程，最后会保留8980个细胞。Raw : 434012
->cellranger3 : 10227->nUMI<12496 : 9876->nGene>=200 : 9828->nGene<3268 : 9806->mito.percent<0.229 : 9507->ribo.percent<0.371 : 9454->diss.percent<0.055 : 8980
我们也可以自定义修改cell.QC.thres.txt阈值文件，设置成认为自己认为更合适的阈值，提供给runScAnnotation过滤。
此外报表中还提到了ambient RNA contamination的因素，具体就不多介绍了。在之前的OSCA笔记中也有相关的学习，详见[OSCA单细胞数据分析笔记-13、Multi-sample comparison]
2.2 runScAnnotation报表解读
(1) QC质控过滤
首先按照cell.QC.thres.txt阈值指标过滤细胞；然后过滤低表达的基因(expressed in too less (nCell < 3)
(2) 创建Seurat对象，按照Seurat标准流程分析
标准化/log转换/归一化、高变基因；
降维(线性/非线性)，默认选择的PCA主成分数为Top30；聚类，默认的resolution分辨率为0.8
按照cluster间差异分析[optional]，默认展示Top5 差异基因
以上为单细胞分析的最基础部分。下面的分析
runScAnnotation()默认均执行，但可以设置对应的bool逻辑值参数选择部分执行
(3) Doublet打分预测
调用了scdsR包的算法，计算每个细胞为Doublet(双细胞的磁珠)的可能性。
(4) 肿瘤微环境细胞类型预测
根据参考数据集，采用OCLR(one-class logistic regression)算法，预测每个细胞分别属于各个细胞类型的可能性，取最高的作为该细胞的类型标签。如果对于给定的细胞类型可能性均很低，则会贴上
Unknown的标签。scCancer包根据肿瘤微环境组成，提供的参考数据集包含：endothelial cells, fibroblasts, immune cells (CD4+ T cells, CD8+ T cells, B cells, nature killer cells, and myeloid cells)
我们也可以提供自己的细胞类型参考数据集，供该包自动注释。详见https://github.com/wguo-research/scCancer/wiki/5.-Other-personalized-settings#train-new-cell-type-templates
(5) 恶性肿瘤细胞预测
调用
infercnvR包，在肿瘤恶性细胞存在较多 copy number alterations (CNV)变异背景下，计算每个细胞的Malignant score。同时再计算内置的正常细胞的Malignant score，观察分布，确定Malignant score判定阈值。
(6) 恶性肿瘤细胞间异质性评价
如上得到恶性肿瘤细胞的分类之后，可以进一步从以下4个方面探究恶性肿瘤细胞之间的异质性。
A：细胞周期：主要使用Seurat包的AddModuleScore() 函数计算细胞周期相关基因(G2/M and S phase markers)在每个细胞的相对表达程度。
B：细胞干性：主要使用参考数据集(stem/progenitor cells)，类似上述细胞注释的方法，使用OCLR算法，计算每个细胞在表达水平上的相关性。
C：特定(通路)基因集表达活性：scCancer包默认会计算来自MSigDB数据库的 50 hallmark gene sets。算法有GSVA与Seurat的AddModuleScore()[默认]两种。此外我们也可以提供自己感兴趣的基因集。
D: 潜在特征基因表达谱分析：上一种方法是已知基因集具有确定的生物学意义。另一方面，我们可以通过NMF非负矩阵分解的方法：首先根据表达矩阵分为若干个主成分(默认为50，每个主成分代表一种未知的表达模式)，然后计算每个细胞对于这些主成分的表达情况。最后可根据感兴趣的cluster表达分布特征，赋予选定的主成分以某个生物学意义。
(7) 基于cluster的细胞间相互作用分析
主要根据参考ligand–receptor(配体-受体) gene pairs，判断cluster两两之间是否存在配对基因平均高表达的基因之一。
we estimated the global interaction strength between any two cell clusters by the number of the remaining gene pairs and the sum of their scores. 如下图所展示的：The size of point means the number of ligand-receptor paris with scores larger than 0.1. The color of point means the sum of the ligand-receptor pairs scores. To compare conveniently, the bottom subplot shows the predicted cell type fraction of each cell set.
以上就是 runScAnnotation()报表结果的基本组成。里面有很多参数可以调整，或者加载结果提供的Seurat对象进行感兴趣的个性分析。
最后的最后，关于这个scCancer包，我目前的感受是集成了很多单细胞数据包的分析手段；分析html结果对于我们理解单细胞数据或者说肿瘤单细胞数据集分析流程会有很大帮助。而且其中一些分析也适用于其它单细胞数据集。实际使用时，可以根据预期以及结果，灵活调整细节参数，以得到满意的结果。
Rethinking batch effect removing methods—MNN
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注