單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析|| 隱藏在PC軸中的秘密

正常情況下婆芦，按照目前主流的單細(xì)胞數(shù)據(jù)分析教程蜂怎，是可以分析我們的數(shù)據(jù)的穆刻。但是，如果在分析過程中發(fā)現(xiàn)了不正常的現(xiàn)象杠步，比如氢伟，batch這個(gè)幽靈真的在腦海里盤旋不去，我們就要檢查batch的來源了幽歼。

Long long ago朵锣，在一個(gè)夜深人靜的夜晚，看著手里的單細(xì)胞數(shù)據(jù)甸私，幾度分析都沒有找到很好的落腳點(diǎn)诚些，還慘雜著若有若無的批次（batch）。到底該如何處理批次皇型，網(wǎng)上給出的方法幾乎都是參考bulk的方法來的诬烹，天知道這樣做是不是合適。我們必須找到batch的藏身之處弃鸦。本著if you can't see it , you can't stop it 的理念绞吁，所有看不到的效應(yīng)都不應(yīng)以莫須有的罪名給予去除。

論單細(xì)胞數(shù)據(jù)分析我們還是比較信任Seurat的唬格，這個(gè)2014就被開發(fā)來分析單細(xì)胞數(shù)據(jù)的R包家破。于是颜说，我們開始遍歷Rahul Satija（Seurat的主要作者之一）的文章，看看他是如何查看和處理單細(xì)胞batch的员舵。不僅遍歷他的文章摘要脑沿，還要讀文章的源碼。經(jīng)過一番努力马僻，我們找到一篇2017年預(yù)印2019年見刊NCB的文章：

文章摘要：

在脊椎動(dòng)物中，位于咽部中胚層心肌細(xì)胞和鰓狀頭部肌肉的多能祖細(xì)胞注服，心肺多能和頭部肌肉的命運(yùn)選擇仍然不清楚韭邓。在這里，我們使用單細(xì)胞基因組學(xué)在簡單的脊索動(dòng)物模型Ciona重建形成第一和第二心臟譜系和咽部肌肉前體的發(fā)育軌跡溶弟，并表征了心肺命運(yùn)選擇的分子基礎(chǔ)女淑。我們表明，F(xiàn)GF-MAPK信號維持多潛能并促進(jìn)咽部肌肉的命運(yùn)辜御，而信號終止允許部署一個(gè)泛心臟計(jì)劃鸭你，由第一和第二心臟譜系共享用以定義心臟同一性。在第二種心臟譜系中擒权，Tbx1/10-Dach通路積極地抑制第一種心臟譜系程序袱巨，調(diào)節(jié)以后跳動(dòng)心臟中的細(xì)胞多樣性。最后碳抄，Ciona和小鼠的跨物種比較揭示了脊索動(dòng)物的心咽網(wǎng)絡(luò)的深層進(jìn)化起源愉老。

文章設(shè)計(jì)與思路：

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%的七個(gè)PC軸的基因，并判斷出每個(gè)PC軸潛在的生物學(xué)意義璧尸，如PC5 作者寫道：Another batch effect咒林。然后，有batch的PCs用RegressOut回歸掉（這個(gè)函數(shù)在V3中放到了 ScaleData的參數(shù)vars.to.regress 中爷光，在R中?Seurat::ScaleData）垫竞。

# Remove batch effect
hpfall.remv1=RegressOut(hpfall,c("PC2","PC5","PC7"),do.scale = T)

這一波操作啟示我們，在去批次之前需要判斷批次的有無以及在哪里瞎颗。重要的是：

批次的直接反映是特定基因的表達(dá)件甥。

舉個(gè)例子，有八個(gè)單細(xì)胞小鼠數(shù)據(jù)雌雄各四例哼拔，而我們想看的不是雌雄之間的區(qū)別而是他們的共性引有，這時(shí)候數(shù)據(jù)表現(xiàn)為按雌雄分開了。檢查我們的PCs所表征的基因倦逐，發(fā)現(xiàn)某PC的基因都是和性別相關(guān)的譬正，這時(shí)候就可以回歸掉(或者去掉)這個(gè)PCs宫补。

這需要我們認(rèn)識這些基因，不認(rèn)識就需要做富集來看通路曾我。

下面我們用Seurat V3+ 來做一個(gè)發(fā)現(xiàn)PC軸秘密的演示粉怕，首先我們還是清出我們的R包和老朋友pbmc3k的數(shù)據(jù)集。

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

在標(biāo)準(zhǔn)流程中抒巢，我們均一化之后就可以執(zhí)行ScaleData贫贝，但是一開始我們并不知道哪些變量需要vars.to.regress，所以我們使用默認(rèn)參數(shù)蛉谜。

# Seurat 標(biāo)準(zhǔn)流程稚晚，這里不需要運(yùn)行，因?yàn)閜bmc3k.final 已經(jīng)做過這些計(jì)算了型诚。
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)

我們主要管關(guān)心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)

早些時(shí)候，我們看到這個(gè)圖沒有什么感覺狰贯，就是一堆字母也搓，現(xiàn)在我們知道如果每個(gè)PC的基因大多是細(xì)胞類型的makergene，那么該P(yáng)C就是細(xì)胞類型特異的涵紊，肯定是要保留的了傍妒。比如這里的PC2，PC3 栖袋，PC4基本斷定是和B細(xì)胞相關(guān)的拍顷。

# 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)

結(jié)合碎石圖，判斷納入多少PC（這里并沒有看基因）

 ElbowPlot(pbmc3k.final)

按照示例文章的做法塘幅，這里我們選擇了PCs1:10昔案，結(jié)合上面的基因情況，判斷這10個(gè)PCs哪些留下哪些去掉电媳。即可踏揣。

這里我們科普下pve的概念，PVE（有一種翻譯是：percent variance explained ）按照這個(gè)概念我們在Seurat里面可以這樣計(jì)算：

(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

計(jì)算累計(jì)貢獻(xiàn)度

(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 = .)

如果有樣本分組信息，自然是可以把這里的細(xì)胞類型換成分組信息了拼缝，開看哪些樣本在哪個(gè)PC中高亮娱局，進(jìn)而推斷其是不是特殊的batch，對應(yīng)的PC以及基因是否需要去掉（或回歸掉）咧七。

在節(jié)目的最后我們再次強(qiáng)調(diào)單細(xì)胞數(shù)據(jù)分析的非線性過程衰齐，許多時(shí)候，如果看不到就不能判斷是否需要做某個(gè)處理继阻，所以需要一個(gè)反饋的過程耻涛。在單細(xì)胞數(shù)據(jù)科學(xué)中PCA分析是屬于特征選擇的過程废酷，即，哪些特征哪來分析抹缕，這當(dāng)然是值得謹(jǐn)慎處理的澈蟆。單細(xì)胞數(shù)據(jù)分析的默認(rèn)參數(shù)（default parameters）時(shí)代已經(jīng)一去不復(fù)返了。

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