cytofWorkflow之基本质量控制（三）
上一讲我们构造好了SingleCellExperiment对象，后续全部的分析都会以这个SingleCellExperiment对象为准，大家务必熟悉SingleCellExperiment对象的各种结构，有了这个SingleCellExperiment对象，第一步是对它做一些质量控制操作。其实前面教程里面的看这些抗体在我们的不同病人的表达量分布情况：
require(cytofWorkflow)
p <- plotExprs(sce, color_by = "condition")
p$facet$params$ncol <- 6
p密度图显示，信号值被归一化。这已经是一种质量控制的啦。
再次回顾一下SingleCellExperiment对象构建的代码
全部的可以复制粘贴就执行的代码如下：
require(cytofWorkflow)
if(!file.exists(file = 'input_for_cytofWorkflow.Rdata')){
library(readxl)
url <- "http://imlspenticton.uzh.ch/robinson_lab/cytofWorkflow"
md <- "PBMC8_metadata.xlsx"
download.file(file.path(url, md), destfile = md, mode = "wb")
md <- read_excel(md)
head(data.frame(md))
table(md[,3:4])
# 样本的表型信息
## 真正的表达矩阵
library(HDCytoData)
fs <- Bodenmiller_BCR_XL_flowSet()
fs
# 如果网络不好，也可以自行下载
# 然后：loaded into R as a flowSet using read.flowSet() from the flowCore package
# 这个cytof的panel的抗体信息表格：
panel <- "PBMC8_panel_v3.xlsx"
download.file(file.path(url, panel), destfile = panel, mode = "wb")
panel <- read_excel(panel)
head(data.frame(panel))
# spot check that all panel columns are in the flowSet object
all(panel$fcs_colname %in% colnames(fs))
# 有了样本的表型信息，panel的抗体信息，以及表达量矩阵，就可以构建对象：
# specify levels for conditions & sample IDs to assure desired ordering
md$condition <- factor(md$condition, levels = c("Ref", "BCRXL"))
md$sample_id <- factor(md$sample_id,
levels = md$sample_id[order(md$condition)])
save(fs, panel, md,file = 'input_for_cytofWorkflow.Rdata')
}
load(file = 'input_for_cytofWorkflow.Rdata')
# expression values, 第一个样品是2838个细胞，39列（其中6列并不是抗体）
dim(exprs(fs[[1]]))
exprs(fs[[1]])[1:6, 1:5]
colnames(exprs(fs[[1]]))
# construct SingleCellExperiment
panel
md
library(CATALYST)
sce <- prepData(fs, panel, md, features = panel$fcs_colname)
pro='basic_qc'
p <- plotExprs(sce, color_by = "condition")
p$facet$params$ncol <- 4
p
ggsave2(filename = paste0(pro,'_all_markers_density.pdf'))
第一个图如下所示：
密度图看抗体信号强度
可以看到，某些抗体在两个不同组的表达量差异很明显， 他们将会是主要的细胞亚群差异来源。
然后看看2个分组的各自8个样品的细胞数量
代码也超级简单，这就是构建SingleCellExperiment对象的方便之处，如下：
n_cells(sce)
plotCounts(sce, group_by = "sample_id", color_by = "condition")
ggsave2(filename = paste0(pro,'_plotCounts.pdf'))
可以看到细胞数量差异是比较大的，这个是cytof数据的客观规律，不用过于在意。
每个样品的细胞数量
当然了，也可以自己写代码去看细胞数量，然后自己绘制这个条形图：
df=data.frame(n_cells=unlist(lapply(1:length(fs), function(i){
dim(exprs(fs[[i]]))[1]
})),
condition=md$condition,sample_id=md$sample_id)
library(ggplot2)
ggplot2::ggplot(df,aes(x=sample_id,y=n_cells,fill=condition))+ geom_bar(stat = "identity")
当然了，需要调整很多细节啦，其实看看作者自己的 plotCounts 函数就明白了。绝大部分情况下，初学者都无需自己写太多函数，各种工具调用的溜起来就行了。
然后看看样品之间的距离关系
代码也超级简单，这就是构建SingleCellExperiment对象的方便之处，如下：
# pbMDS: Pseudobulk-level MDS plot.
# A multi-dimensional scaling (MDS) plot on aggregated measurement values may be rendered with pbMDS.
pbMDS(sce, color_by = "condition", label_by = "sample_id")
ggsave2(filename = paste0(pro,'_pbMDS.pdf'))
可以看到cytofWorkflow给出来的示例数据超级正常，不同分组的病人距离很远，泾渭分明，如下：
样品之间的距离关系
这样的图表在我一直讲解的表达矩阵分析也是如此。很多人仅仅是学习了代码，不明白后面的统计学意义，PCA图的意义，我一直强调，做表达矩阵分析一定要有三张图，见：你确定你的差异基因找对了吗？ 。
这里面的MDS看样本距离关系，跟我们一直强调的PCA看样本距离关系，大同小异。
表达量热图看差异
代码也超级简单，这就是构建SingleCellExperiment对象的方便之处，如下：
# 10 lineage markers and 14 functional markers across all cells measured for each sample in the PBMC dataset
pdf( paste0(pro,'_plotExprHeatmap.pdf'))
plotExprHeatmap(sce, scale = "last",
hm_pal = rev(hcl.colors(10, "YlGnBu")))
dev.off()
出图如下：
抗体信号强度在所有病人的热图
可以看到，在Ref组里面非常多的CD4表达，在BCRXL组却低，所以最后对这些样本的细胞进行聚类分群后，是可以看到CD4阳性T细胞的比例差异的。
另外一种展现表达量差异的条形图
代码也超级简单，这就是构建SingleCellExperiment对象的方便之处，如下：
# PCA-based non-redundancy score (NRS)
plotNRS(sce, features = "type", color_by = "condition")
ggsave2(filename = paste0(pro,'_plotNRS.pdf'))
出图如下：
PCA-based non-redundancy score 比较难理解，反正绝大部分情况下我是忽略这张图的，虽然说cytofWorkflow教程里面展现了它。
其实这个时候的分析，跟纯粹的单细胞转录组就非常类似了。单细胞转录组数据分析的细节，以及背景我就不赘述了，看我在《单细胞天地》的单细胞基础10讲。