隐藏在PC轴中的秘密
分享是一种态度
作者 | 周运来
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树。
生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。
前言
正常情况下，按照目前主流的单细胞数据分析教程，是可以分析我们的数据的。但是，如果在分析过程中发现了不正常的现象，比如，batch这个幽灵真的在脑海里盘旋不去，我们就要检查batch的来源了。
Long long ago，在一个夜深人静的夜晚，看着手里的单细胞数据，几度分析都没有找到很好的落脚点，还惨杂着若有若无的批次（batch）。到底该如何处理批次，网上给出的方法几乎都是参考bulk的方法来的，天知道这样做是不是合适。我们必须找到batch的藏身之处。本着if you can't see it , you can't stop it 的理念，所有看不到的效应都不应以莫须有的罪名给予去除。
论单细胞数据分析我们还是比较信任Seurat的，这个2014就被开发来分析单细胞数据的R包。于是，我们开始遍历Rahul Satija（Seurat的主要作者之一）的文章，看看他是如何查看和处理单细胞batch的。不仅遍历他的文章摘要，还要读文章的源码。经过一番努力，我们找到一篇2017年预印2019年见刊NCB的文章：
文章摘要：
在脊椎动物中，位于咽部中胚层心肌细胞和鳃状头部肌肉的多能祖细胞，心肺多能和头部肌肉的命运选择仍然不清楚。在这里，我们使用单细胞基因组学在简单的脊索动物模型Ciona重建形成第一和第二心脏谱系和咽部肌肉前体的发育轨迹，并表征了心肺命运选择的分子基础。我们表明，FGF-MAPK信号维持多潜能并促进咽部肌肉的命运，而信号终止允许部署一个泛心脏计划，由第一和第二心脏谱系共享用以定义心脏同一性。在第二种心脏谱系中，Tbx1/10-Dach通路积极地抑制第一种心脏谱系程序，调节以后跳动心脏中的细胞多样性。最后，Ciona和小鼠的跨物种比较揭示了脊索动物的心咽网络的深层进化起源。
文章设计与思路：
Early cardiopharyngeal development in Ciona, and sampling stages and established lineage tree. Cardiopharyngeal lineage cells are shown for only one side and known cell-type-specific marker genes are indicated. st., FABA stage55; hpf, hours post-fertilization; TVC, trunk ventral cell; STVC, second trunk ventral cell; FHP, first heart precursor; ASMF, atrial siphon muscle founder cells; SHP, second heart precursor; iASMP, inner atrial siphon muscle precursor; oASMP, outer atrial siphon muscle precursor; LoM, longitudinal muscles; QC, quality control. Dotted line: midline.
这篇文章用的是Seurat v1.2，代码作者为Xiang Niu May 11, 2016 ，我们先来看一段去除batch的源码。
4. Remove Contamination and Batch Effect
```{r,warning=FALSE,message=FALSE,tidy=TRUE}
# PCA on all the genes
hpfall=pca(hpfall,pc.genes = rownames(hpfall@data),do.print = F)
pcScree(hpfall,rownames(hpfall@data),10)
abline(h=0.5)
# Top 7 PCs are selected (PVE>0.5%)
# PC1
# Population of potential contaminated cell type
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,1) ,cells.use = pcTopCells(hpfall,1),col.use = col)
pca.plot(hpfall,1,2)
# PC2
# Batch effect due to library preparation, it separates hpf12/14/16 with hpf18/20
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,2),cells.use = pcTopCells(hpfall,2),col.use = col)
pca.plot(hpfall,1,2)
# PC3
# Mesenchymal contamination
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,3),cells.use = pcTopCells(hpfall,3),col.use = col)
pca.plot(hpfall,1,3)
feature.plot(hpfall,c("RTP4","TWIST1"),reduction.use = "pca",dim.2 = 3)
# PC4
# Heart vs Muscle split
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,4),cells.use = pcTopCells(hpfall,4),col.use = col)
feature.plot(hpfall,c("NKX2-3","EBF1/2/3/4","TBX1/10","GATA4/5/6"),reduction.use = "pca",dim.2 = 4)
# PC5
# Another batch effect
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,5),cells.use = pcTopCells(hpfall,5),col.use = col)
pca.plot(hpfall,1,5)
# PC6
# Saperation by heart vs muscle
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,6),cells.use = pcTopCells(hpfall,6),col.use = col)
pca.plot(hpfall,1,6)
feature.plot(hpfall,c("NKX2-3","EBF1/2/3/4","TBX1/10","GATA4/5/6"),reduction.use = "pca",dim.2 = 6)
PC7
# Batch effect marked hpf14
doHeatMap(hpfall,remove.key = T,slim.col.label = T,genes.use = pcTopGenes(hpfall,7),cells.use = pcTopCells(hpfall,7),col.use = col)
pca.plot(hpfall,1,7)
我们看到作者在这里用doHeatMap和pca.plot查看了PVE>0.5%的七个PC轴的基因，并判断出每个PC轴潜在的生物学意义，如PC5 作者写道：Another batch effect。然后，有batch的PCs用RegressOut回归掉（这个函数在V3中放到了 ScaleData的参数vars.to.regress 中，在R中?Seurat::ScaleData）。
# Remove batch effect
hpfall.remv1=RegressOut(hpfall,c("PC2","PC5","PC7"),do.scale = T)
这一波操作启示我们，在去批次之前需要判断批次的有无以及在哪里。重要的是：
批次的直接反映是特定基因的表达。
举个例子，有八个单细胞小鼠数据雌雄各四例，而我们想看的不是雌雄之间的区别而是他们的共性，这时候数据表现为按雌雄分开了。检查我们的PCs所表征的基因，发现某PC的基因都是和性别相关的，这时候就可以回归掉(或者去掉)这个PCs。
这需要我们认识这些基因，不认识就需要做富集来看通路。
下面我们用Seurat V3+ 来做一个发现PC轴秘密的演示，首先我们还是清出我们的R包和老朋友pbmc3k的数据集。
library(Seurat)
library(SeuratData)
pbmc3k.final
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features, 2000 variable features)
2 dimensional reductions calculated: pca, umap
在标准流程中，我们均一化之后就可以执行ScaleData，但是一开始我们并不知道哪些变量需要vars.to.regress，所以我们使用默认参数。
# Seurat 标准流程，这里不需要运行，因为pbmc3k.final 已经做过这些计算了。
pbmc.counts <- Read10X(data.dir = "D:\\Users\\Administrator\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19")
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunUMAP(object = pbmc)
我们主要关心PCs 的Loading值。
t(Loadings(object = pbmc3k.final[["pca"]])[1:5,1:5])
PPBP LYZ S100A9 IGLL5 GNLY
PC_1 0.010990202 0.11623171 0.11541436 -0.007987473 -0.01523876
PC_2 0.011484260 0.01472515 0.01895146 0.054542385 -0.13375626
PC_3 -0.151760919 -0.01280613 -0.02368853 0.049015330 0.04101340
PC_4 0.104037367 -0.04414540 -0.05787777 0.066947217 0.06912322
PC_5 0.003299077 0.04990688 0.08538231 0.004603231 0.10455861
# Set the feature loadings for a given DimReduc
new.loadings <- Loadings(object = pbmc3k.final[["pca"]])
new.loadings <- new.loadings + 0.01
Loadings(object = pbmc3k.final[["pca"]]) <- new.loadings
VizDimLoadings(pbmc3k.final,dims = 1:4)
早些时候，我们看到这个图没有什么感觉，就是一堆字母，现在我们知道如果每个PC的基因大多是细胞类型的makergene，那么该PC就是细胞类型特异的，肯定是要保留的了。比如这里的PC2，PC3 ，PC4基本断定是和B细胞相关的。
# Examine and visualize PCA results a few different ways
print(pbmc3k.final[["pca"]], dims = 1:5, nfeatures = 5)
PC_ 1
Positive: CST3, TYROBP, LST1, AIF1, FTL
Negative: MALAT1, LTB, IL32, IL7R, CD2
PC_ 2
Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
Negative: NKG7, PRF1, CST7, GZMB, GZMA
PC_ 3
Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
Negative: PPBP, PF4, SDPR, SPARC, GNG11
PC_ 4
Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
Negative: VIM, IL7R, S100A6, IL32, S100A8
PC_ 5
Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
Negative: LTB, IL7R, CKB, VIM, MS4A7
DimPlot(pbmc3k.final, reduction = "pca",group.by = 'seurat_annotations') # group.by可以换成样本分组
我们也看看PC热图的情况
#DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc3k.final, dims = 1:15, cells = 500, balanced = TRUE)
结合碎石图，判断纳入多少PC（这里并没有看基因）
ElbowPlot(pbmc3k.final)
按照示例文章的做法，这里我们选择了PCs1:10，结合上面的基因情况，判断这10个PCs哪些留下哪些去掉。即可。
这里我们科普下pve的概念，PVE（有一种翻译是：percent variance explained ）按照这个概念我们在Seurat里面可以这样计算：
(pbmc3k.final[["pca"]]@stdev^2/sum(pbmc3k.final[["pca"]]@stdev^2)) *100
[1] 20.468928 8.209683 6.092213 5.709129 4.086683 2.631764 1.737523 1.536987 1.386379 1.367403 1.346245 1.299317 1.285963 1.254616 1.246295 1.238943
[17] 1.218781 1.215172 1.205959 1.202547 1.198533 1.195067 1.188333 1.181939 1.181202 1.177457 1.174989 1.168232 1.167059 1.161573 1.158020 1.148486
[33] 1.145992 1.144709 1.139288 1.139030 1.133048 1.129764 1.128323 1.124985 1.119834 1.116586 1.115314 1.111942 1.111581 1.105309 1.102191 1.100031
[49] 1.096533 1.094121
sum((pbmc3k.final[["pca"]]@stdev^2/sum(pbmc3k.final[["pca"]]@stdev^2)) *100 )
[1] 100
计算累计贡献度
(cumsum((pbmc3k.final[["pca"]]@stdev^2/sum(pbmc3k.final[["pca"]]@stdev^2)) *100))
最后，我们可以在umap空间中可视化PC信息；
# Defining the information in the seurat object of interest
columns <- c(paste0("PC_", 1:16),
"seurat_annotations",
"UMAP_1", "UMAP_2")
# Extracting this data from the seurat object
pc_data <- FetchData(pbmc3k.final,
vars = columns)
# Adding cluster label to center of cluster on UMAP
umap_label <- FetchData(pbmc3k.final,
vars = c("seurat_annotations", "UMAP_1", "UMAP_2")) %>%
group_by(seurat_annotations) %>%
summarise(x=mean(UMAP_1), y=mean(UMAP_2))
# Plotting a UMAP plot for each of the PCs
map(paste0("PC_", 1:16), function(pc){
ggplot(pc_data,
aes(UMAP_1, UMAP_2)) +
geom_point(aes_string(color=pc),
alpha = 0.7) +
scale_color_gradient(guide = FALSE,
low = "grey90",
high = "blue") +
geom_text(data=umap_label,
aes(label=seurat_annotations, x, y)) +
ggtitle(pc)+ theme_bw()
}) %>%
plot_grid(plotlist = .)
如果有样本分组信息，自然是可以把这里的细胞类型换成分组信息了，再看哪些样本在哪个PC中高亮，进而推断其是不是特殊的batch，对应的PC以及基因是否需要去掉（或回归掉）。
在节目的最后我们再次强调单细胞数据分析的非线性过程，许多时候，如果看不到就不能判断是否需要做某个处理，所以需要一个反馈的过程。在单细胞数据科学中PCA分析是属于特征选择的过程，即，哪些特征哪来分析，这当然是值得谨慎处理的。单细胞数据分析的默认参数（default parameters）时代已经一去不复返了。
References
[1] A single-cell transcriptional roadmap for cardiopharyngeal fate diversification: https://www.nature.com/articles/s41556-019-0336-z[2] https://satijalab.org/people/[3] https://as.nyu.edu/content/nyu-as/as/faculty/rahul-satija.html[4] https://tem11010.github.io/Plotting-PCAs/[5] https://cmdlinetips.com/2019/04/introduction-to-pca-with-r-using-prcomp/[6] https://eranraviv.com/understanding-variance-explained-in-pca/[7] https://hbctraining.github.io/scRNA-seq/lessons/08_SC_clustering_quality_control.html[8] https://github.com/stevexniu/single-cell-ciona/blob/master/Preprocess.Rmd
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注