单细胞转录组数据分析|| 隐藏在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 stage⁵⁵; 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）时代已经一去不复返了。

A single-cell transcriptional roadmap for cardiopharyngeal fate diversification
https://satijalab.org/people/
https://as.nyu.edu/content/nyu-as/as/faculty/rahul-satija.html
https://tem11010.github.io/Plotting-PCAs/
https://cmdlinetips.com/2019/04/introduction-to-pca-with-r-using-prcomp/
https://eranraviv.com/understanding-variance-explained-in-pca/
https://hbctraining.github.io/scRNA-seq/lessons/08_SC_clustering_quality_control.html
https://github.com/stevexniu/single-cell-ciona/blob/master/Preprocess.Rmd