單細胞交響樂31-實戰(zhàn)十四 10X HCA計劃的38萬骨髓細胞

劉小澤寫于2020.7.22
為何取名叫“交響樂”涡贱？因為單細胞分析就像一個大樂團促煮，需要各個流程的協(xié)同配合
單細胞交響樂1-常用的數(shù)據(jù)結構SingleCellExperiment
單細胞交響樂2-scRNAseq從實驗到下游簡介
 單細胞交響樂3-細胞質控
 單細胞交響樂4-歸一化
 單細胞交響樂5-挑選高變化基因
 單細胞交響樂6-降維
 單細胞交響樂7-聚類分群
 單細胞交響樂8-marker基因檢測
 單細胞交響樂9-細胞類型注釋
 單細胞交響樂9-細胞類型注釋
 單細胞交響樂10-數(shù)據(jù)集整合后的批次矯正
 單細胞交響樂11-多樣本間差異分析
 單細胞交響樂12-檢測Doublet
單細胞交響樂13-細胞周期推斷
 單細胞交響樂14-細胞軌跡推斷
 單細胞交響樂15-scRNA與蛋白豐度信息結合
 單細胞交響樂16-處理大型數(shù)據(jù)
單細胞交響樂17-不同單細胞R包的數(shù)據(jù)格式相互轉換
 單細胞交響樂18-實戰(zhàn)一 Smart-seq2
單細胞交響樂19-實戰(zhàn)二 STRT-Seq
單細胞交響樂20-實戰(zhàn)三 10X 未過濾的PBMC數(shù)據(jù)
單細胞交響樂21-實戰(zhàn)三批量處理并整合多個10X PBMC數(shù)據(jù)
單細胞交響樂22-實戰(zhàn)五 CEL-seq2
單細胞交響樂23-實戰(zhàn)六 CEL-seq
單細胞交響樂24-實戰(zhàn)七 SMARTer 胰腺細胞
 單細胞交響樂25-實戰(zhàn)八 Smart-seq2 胰腺細胞
 單細胞交響樂26-實戰(zhàn)九胰腺細胞數(shù)據(jù)整合
 單細胞交響樂27-實戰(zhàn)十 CEL-seq-小鼠造血干細胞
 單細胞交響樂28-實戰(zhàn)十一 Smart-seq2-小鼠造血干細胞
 單細胞交響樂29-實戰(zhàn)十二 10X 小鼠嵌合體胚胎
 單細胞交響樂30-實戰(zhàn)十三 10X 小鼠乳腺上皮細胞

1 前言

前面的種種都是作為知識儲備殃饿，但是不實戰(zhàn)還是記不住前面的知識
這是第十四個實戰(zhàn)練習，不過這次的練習對電腦要求比較高

單細胞領域不得不提的HCA計劃

參考：
官網(wǎng)：https://data.humancellatlas.org/
http://www.casisd.cn/zkcg/ydkb/kjqykb/2016/201612/201707/t20170703_4821952.html
https://baike.baidu.com/tashuo/browse/content?id=6937242830be0d08a0384772
http://med.china.com.cn/content/pid/134504/tid/1015

細胞類型多種多樣

細胞是生命最基本的單位毁靶，一般人至少有37兆2000億個細胞辜限，但現(xiàn)有的研究一般是組織層面，對細胞內部機制理解的還不是很深入穷吮。人體內的細胞有多少種類逻翁？這是一個非常基礎的問題捡鱼。美國國立衛(wèi)生研究院給出了一個答案：200種八回，包括像神經元細胞、心臟細胞驾诈、肌肉細胞等幾大類缠诅。然而，如果你去問一個免疫學家乍迄，他會告訴你光免疫細胞就至少有200種管引；如果你去問一個專門研究T細胞的免疫學家，他會告訴你光T細胞就至少有200種闯两。Regev說她自己光講人體細胞種類就要花15分鐘褥伴，并且一個細胞的種類可以不斷地往下劃分，亞型下面還有亞型漾狼，僅視網(wǎng)膜組織就至少包括100種不同種類的神經元重慢。更棘手的是，有些種類的細胞會在一定條件下轉化為其他種類的細胞逊躁，并且每個亞型的細胞會根據(jù)不同環(huán)境呈現(xiàn)出不同的狀態(tài)伤锚。

細胞圖譜的重要性

一個完整的人體細胞圖譜將賦予我們每種細胞類型唯一的“身份證”，可以幫助解釋不同類型的細胞如何協(xié)作并形成組織的志衣。另外與疾病相關的基因在哪些細胞更活躍屯援，不同細胞類型產生的機理又是怎樣的呢？

HCA計劃的開啟

2016年Chan Zuckerberg Initiative（CZI）基金會計劃在未來十年內投入30億美元支持科學研究念脯，其中最矚目的就是“人類細胞圖譜計劃”（the Human Cell Atlas）【可與“人類基因組計劃”相媲美】狞洋。2017年10月，HCA正式公布了首批擬資助的38個項目绿店，將為40萬億個細胞繪制圖譜吉懊，其中清華大學教授張學工負責的項目作為其中唯一一個由中國科學家承擔的項目庐橙，從全球范圍內征集的近500個項目中脫穎而出。

2018年3月8日借嗽，Sanger研究所在其官網(wǎng)宣布了人類發(fā)育細胞圖譜計劃（HDCA）的最新成果研究團隊從捐贈的人類發(fā)育組織（包括肝臟态鳖、皮膚、腎臟恶导、胎盤）中分離出超25萬個細胞浆竭，并利用強大的單細胞基因組測序工具獲取了人類早期發(fā)育，以及影響健康或者導致疾病的信息惨寿。HDCA是HCA的重要一部分邦泄，旨在構建參與人類發(fā)育的所有重要細胞的基因組參考圖譜。該項目聚焦的其他主要領域包括裂垦，改善對血細胞如何形成顺囊，以及免疫系統(tǒng)如何發(fā)揮功能的理解。

另外蕉拢，2018年特碳，美國費城兒童醫(yī)院和辛辛那提兒童醫(yī)院領導的國際研究團隊建議在HCA聯(lián)盟中建立一個縱向的兒童研究分支，并在HCA白皮書中概述了兒童細胞圖譜計劃（Pediatric Cell Atlas晕换，PCA）测萎，以在兒童健康和人類發(fā)展背景下，展開針對兒童獨特生物學特征的跨學科研究届巩。

HCA的主要內容

編目所有人體細胞類型（如免疫細胞硅瞧、腦細胞）及其子類型；
繪制不同細胞類型在組織和人體中的分布圖恕汇；
區(qū)分細胞狀態(tài)（如免疫細胞在被病原體激活前后的狀態(tài)）腕唧；
捕捉細胞轉換過程的關鍵特征（如干細胞細胞激活和分化）；
追蹤細胞譜系歷史（如骨髓前體干細胞到功能性的紅細胞）

數(shù)據(jù)準備

我們這里使用的數(shù)據(jù)是：10X技術得到的人類約380,000個骨髓細胞

數(shù)據(jù)鏈接：https://share.weiyun.com/JLvQvpHS

# 如果要自己下載
# library(HCAData)
# sce.bone <- HCAData('ica_bone_marrow')
load('sce.bone.new.RData')
sce.bone
# class: SingleCellExperiment 
# dim: 33694 378000 
# metadata(0):
#   assays(1): counts
# rownames(33694): ENSG00000243485
# ENSG00000237613 ... ENSG00000277475
# ENSG00000268674
# rowData names(2): ID Symbol
# colnames(378000):
#   MantonBM1_HiSeq_1-AAACCTGAGCAGGTCA-1
# MantonBM1_HiSeq_1-AAACCTGCACACTGCG-1 ...
# MantonBM8_HiSeq_8-TTTGTCATCTGCCAGG-1
# MantonBM8_HiSeq_8-TTTGTCATCTTGAGAC-1
# colData names(1): Barcode
# reducedDimNames(0):
#   altExpNames(0):

數(shù)據(jù)初探

# 這里的數(shù)據(jù)大小41M
object.size(counts(sce.bone))
# 41741200 bytes

# 真實文件大小700多M
file.info(path(counts(sce.bone)))$size
# [1] 769046295

# 其中包含供體信息
head(sce.bone$Barcode)
# [1] "MantonBM1_HiSeq_1-AAACCTGAGCAGGTCA-1"
# [2] "MantonBM1_HiSeq_1-AAACCTGCACACTGCG-1"
# [3] "MantonBM1_HiSeq_1-AAACCTGCACCGGAAA-1"
# [4] "MantonBM1_HiSeq_1-AAACCTGCATAGACTC-1"
# [5] "MantonBM1_HiSeq_1-AAACCTGCATCGATGT-1"
# [6] "MantonBM1_HiSeq_1-AAACCTGCATTCCTGC-1"

# 提取供體信息
sce.bone$Donor <- sub("_.*", "", sce.bone$Barcode)
table(sce.bone$Donor)
# 
# MantonBM1 MantonBM2 MantonBM3 MantonBM4 MantonBM5 
# 48000     48000     48000     48000     48000 
# MantonBM6 MantonBM7 MantonBM8 
# 42000     48000     48000

ID轉換

還是要獲得染色體信息

library(EnsDb.Hsapiens.v86)
rowData(sce.bone)$Chr <- mapIds(EnsDb.Hsapiens.v86, keys=rownames(sce.bone),
    column="SEQNAME", keytype="GENEID")

# 其中包括13個線粒體信息
table(grepl('MT',rowData(sce.bone)$Chr))
# 
# FALSE  TRUE 
# 33681    13

整合行名

library(scater)
rownames(sce.bone) <- uniquifyFeatureNames(rowData(sce.bone)$ID,
    names = rowData(sce.bone)$Symbol)

2 質控

使用線粒體信息進行過濾瘾英，并且因為數(shù)據(jù)量很大枣接，可以調用多線程

library(BiocParallel)
# 調用8線程
bpp <- MulticoreParam(8)

start=Sys.time()
sce.bone <- unfiltered <- addPerCellQC(sce.bone, BPPARAM=bpp,
                                     subsets=list(Mito=which(rowData(sce.bone)$Chr=="MT")))
end=Sys.time()
(end-start)
# Time difference of 2.167978 mins

過濾掉了6萬多個細胞

colSums(as.matrix(qc))
# low_lib_size            low_n_features 
# 33997                     42756 
# high_subsets_Mito_percent                   discard 
# 44105                     61275

作圖

unfiltered$discard <- qc$discard

gridExtra::grid.arrange(
    plotColData(unfiltered, x="Donor", y="sum", colour_by="discard") +
        scale_y_log10() + ggtitle("Total count"),
    plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
        scale_y_log10() + ggtitle("Detected features"),
    plotColData(unfiltered, x="Donor", y="subsets_Mito_percent",
        colour_by="discard") + ggtitle("Mito percent"),
    ncol=2
)

image-20200722115024437

再看看線粒體含量與文庫大小的關系

plotColData(unfiltered, x="sum", y="subsets_Mito_percent", 
    colour_by="discard") + scale_x_log10()

image-20200722115105131

3 歸一化

這里為了減少計算量，使用原來計算好的文庫size factor缺谴，不需要重復計算

sce.bone <- logNormCounts(sce.bone, size_factors = sce.bone$sum)
summary(sizeFactors(sce.bone))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.05    0.47    0.65    1.00    0.89   42.38

4 找高變異基因

將供體作為批次信息加入到構建模型中但惶，這一步依然需要多線程處理

library(scran)
start=Sys.time()
dec.bone <- modelGeneVar(sce.bone, block=sce.bone$Donor, BPPARAM=bpp)
end=Sys.time()
(end-start) # Time difference of 10.58185 mins

top.bone <- getTopHVGs(dec.bone, n=5000)

5 矯正批次效應

library(batchelor)
library(BiocNeighbors)

set.seed(1010001)
merged.bone <- fastMNN(sce.bone, batch = sce.bone$Donor, subset.row = top.bone,
     BSPARAM=BiocSingular::RandomParam(deferred = TRUE), 
     BNPARAM=AnnoyParam(),
     BPPARAM=bpp)

reducedDim(sce.bone, 'MNN') <- reducedDim(merged.bone, 'corrected')

看下結果：lost.var 值越大表示丟失的真實生物異質性越多

##      MantonBM1 MantonBM2 MantonBM3 MantonBM4 MantonBM5 MantonBM6 MantonBM7
## [1,]  0.010616  0.008241  0.000000  0.000000  0.000000  0.000000  0.000000
## [2,]  0.007894  0.007741  0.023662  0.000000  0.000000  0.000000  0.000000
## [3,]  0.005834  0.003823  0.005423  0.025272  0.000000  0.000000  0.000000
## [4,]  0.003482  0.002868  0.002581  0.003067  0.027117  0.000000  0.000000
## [5,]  0.005287  0.003544  0.003173  0.005702  0.006551  0.031890  0.000000
## [6,]  0.004601  0.004535  0.004533  0.004067  0.004934  0.005547  0.034452
## [7,]  0.002610  0.002238  0.003135  0.002799  0.001963  0.002697  0.002404
##      MantonBM8
## [1,]   0.00000
## [2,]   0.00000
## [3,]   0.00000
## [4,]   0.00000
## [5,]   0.00000
## [6,]   0.00000
## [7,]   0.04033

6 降維聚類

降維

這么大的細胞數(shù)據(jù)，一般會使用UMAP

set.seed(01010100)
sce.bone <- runUMAP(sce.bone, dimred="MNN",
    BNPARAM=AnnoyParam(),
    BPPARAM=bpp,
    n_threads=bpnworkers(bpp))

另外關于tSNE和UMAP的對比：可以看看https://towardsdatascience.com/how-exactly-umap-works-13e3040e1668

SNE does not preserve global data structure, meaning that only within cluster distances are meaningful while between cluster similarities are not guaranteed. But UMAP Can Preserve Global Structure

tSNE does not scale well for rapidly increasing sample sizes in scRNAseq.

UMAP is faster than tSNE when it concerns a) large number of data points, b) number of embedding dimensions greater than 2 or 3, c) large number of ambient dimensions in the data set

UMAP的改進：UMAP overall follows the philosophy of tSNE, but introduces a number of improvements such as another cost function and the absence of normalization of high- and low-dimensional probabilities.

總結：Despite tSNE served the Single Cell research area for years, it has too many disadvantages such as speed and the lack of global distance preservation

聚類

因為這里的細胞數(shù)量實在太大湿蛔，因此可以先借用kmeans方法對它們聚類膀曾，比如先聚成1000個小類，然后再對每個小類進行graph-based聚類

set.seed(1000)
clust.bone <- clusterSNNGraph(sce.bone, use.dimred="MNN", 
    use.kmeans=TRUE, kmeans.centers=1000)
colLabels(sce.bone) <- factor(clust.bone)
table(colLabels(sce.bone))
## 
##     1     2     3     4     5     6     7     8     9    10    11    12 
## 21938 42850 19861 38064 38024 71978 19237 24627  7583 16361  3179 13023

做個熱圖

tab <- table(Cluster=colLabels(sce.bone), Donor=sce.bone$Donor)
library(pheatmap)
pheatmap(log10(tab+10), color=viridis::viridis(100))

看到每個cluster中都會有幾個不同的樣本阳啥，這也符合預期添谊，因為所有的樣本都是來自不同供體的重復，它們實際都是骨髓細胞

image-20200722160655442

7 找marker基因

markers.bone <- findMarkers(sce.bone, block = sce.bone$Donor, 
    direction = 'up', lfc = 1, BPPARAM=bpp)

然后檢查第5群細胞

top.markers <- markers.bone[["5"]]
best <- top.markers[top.markers$Top <= 10,]
lfcs <- getMarkerEffects(best)

library(pheatmap)
pheatmap(lfcs, breaks=seq(-5, 5, length.out=101))

看到LYZ察迟、S100A8斩狱、VCAN基因高表達耳高，說明cluster5可能是單核細胞（Monocyte）

image-20200722161432088

8 細胞類型注釋

使用參考數(shù)據(jù)集幫助注釋，只不過這里并不是對每個細胞單獨操作所踊，因為數(shù)量太太泌枪。而是基于12個已經分好的cluster的總體表達量進行注釋，提高了注釋速度秕岛，但同樣犧牲了研究cluster內部細胞異質性的機會碌燕。適用于快速注釋細胞類群

se.aggregated <- sumCountsAcrossCells(sce.bone, id=colLabels(sce.bone))

library(SingleR)
hpc <- HumanPrimaryCellAtlasData()
anno.single <- SingleR(se.aggregated, ref = hpc, labels = hpc$label.main,
    assay.type.test="sum")

看一下結果，包含12行也就是12個分群結果瓣蛀，看到利用參考注釋的方法，cluster5也是被注釋到了Monocyte這個細胞類型

anno.single
## DataFrame with 12 rows and 5 columns
##                             scores     first.labels      tuning.scores
##                           <matrix>      <character>        <DataFrame>
## 1   0.388013:0.636789:0.754052:...              CMP 0.554195:0.4238259
## 2   0.326563:0.692156:0.661226:...          NK_cell 0.579208:0.4459558
## 3   0.323785:0.668599:0.743138:... Pre-B_cell_CD34- 0.488346:0.0719803
## 4   0.346204:0.637353:0.615021:...          T_cells 0.628161:0.3531725
## 5   0.297369:0.638962:0.745831:... Pre-B_cell_CD34- 0.548995:0.2633979
## ...                            ...              ...                ...
## 8   0.313537:0.773897:0.672663:...           B_cell 0.722126:-0.288616
## 9   0.346503:0.703506:0.670777:... Pro-B_cell_CD34+ 0.774848: 0.710191
## 10  0.313591:0.774162:0.668809:...           B_cell 0.719013:-0.281128
## 11  0.369873:0.720021:0.654892:...           B_cell 0.523979: 0.326187
## 12  0.403994:0.603811:0.697425:...              MEP 0.430493: 0.286268
##               labels    pruned.labels
##          <character>      <character>
## 1                CMP              CMP
## 2            T_cells          T_cells
## 3           Monocyte         Monocyte
## 4            T_cells          T_cells
## 5           Monocyte         Monocyte
## ...              ...              ...
## 8             B_cell           B_cell
## 9   Pro-B_cell_CD34+ Pro-B_cell_CD34+
## 10            B_cell           B_cell
## 11            B_cell               NA
## 12        BM & Prog.       BM & Prog.

歡迎關注我們的公眾號～_～　　
我們是兩個農轉生信的小碩雷厂，打造生信星球惋增，想讓它成為一個不拽術語、通俗易懂的生信知識平臺改鲫。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!