Mapping Distinct Bone Marrow Niche Populations and Their Differentiation Paths
骨髓龕中不同細胞群體的關聯(lián)及其分化途徑
文章信息介紹的文章利用單細胞轉(zhuǎn)錄組分析骨髓龕中不同細胞類群間的關聯(lián)及其分化軌跡中不同的轉(zhuǎn)錄調(diào)控因子的功能。文章于2019年3月7日發(fā)表在Cell 雜志。于2019年7月7日發(fā)表在Cell reports 雜志空免。文章標題是：Mapping Distinct Bone Marrow Niche Populations and Their Differentiation Paths牢酵；DOI：https://doi.org/10.1016/j.celrep.2019.06.031。文章連接
Finally, we validated our findings using lineage-specific reporter strains and targeted knockdowns. Our analysis reveals differentiation hierarchies for maturing stromal cells, determines key transcription factors along these trajectories, and provides an understanding of the complexity of the bone marrow microenvironment.
摘要
骨髓微環(huán)境是由復雜表型和細胞成熟軌跡未知的異質(zhì)性的非造血細胞群體組成拿诸。在這些非造血細胞群體中，間充質(zhì)細胞維持產(chǎn)生間質(zhì)細胞，骨細胞旅东、脂肪細胞和軟骨細胞。闡明骨髓內(nèi)這些獨特的細胞亞群仍然具有挑戰(zhàn)性十艾。本研究中抵代，我們使用非造血骨髓細胞的單細胞RNA測序來定義特定的細胞亞群。此外忘嫉，通過結(jié)合細胞狀態(tài)層次的計算預測和已知的關鍵轉(zhuǎn)錄因子的表達荤牍，我們繪制了針對骨細胞，軟骨細胞和脂肪細胞譜系的分化途徑庆冕。 最后康吵，我們使用譜系特異性報告菌株和靶向敲除來驗證我們發(fā)現(xiàn)的基因。 我們的分析揭示了成熟基質(zhì)細胞的分化層次访递，確定了伴隨這些分化軌跡的關鍵轉(zhuǎn)錄因子涎才，并增加了對骨髓微環(huán)境復雜性的理解。
測序數(shù)據(jù)介紹
For scRNA-seq, we used inDrops (Klein et al., 2015) following a previously described protocol (Zilionis et al., 2017) with the following modifications: the sequence of the primer on the hydrogel beads was 5′-CGATGACGTAATACGACTCACTATAGGGTGTCGGGTGCAG[bc1,8nt]GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG[bc2,8nt]NNNNNNTTTTTTTTTTTTTTTTTTTV-3′; the sequence of the PE2-N6 primer (step 151 in (Zilionis et al., 2017) was 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGNNNNNN-3′; and the sequences of the PCR primers (steps 157 and 160 in (Zilionis et al., 2017)) were 5′-AATGATACGGCGACCACCGAGATCTACACXXXXXXXXTCGTCGGCAGCGTC-3′ and 5′-CAAGCAGAAGACGGCATACGAGATGGGTGTCGGGTGCAG-3′. Following droplet barcoding reverse transcription, emulsions were split into aliquots of approximately 1,000 (in vivo samples) or 3,000 (cultured samples) single-cell transcriptomes and frozen at ?80°C. For the in vivo samples, two libraries (n = 1,533 cells total) were prepared for mouse 1 and three libraries (n = 3,574 cells total) were prepared for mouse 2. For the three cultured samples, one library per sample was prepared (n = 2,837, n = 2,164, and n = 2,520 total cells, respectively). These cell numbers correspond to the final number of transcriptomes detected upon removal of background barcodes and stressed or dying cells (see section below).
數(shù)據(jù)分析情況
數(shù)據(jù)過濾：統(tǒng)計每個文庫數(shù)據(jù)所有獨特的12 bp的cell barcodes(CBs),排除出現(xiàn)次數(shù)小于100次的CBs力九，并過濾掉CBs包含八個相同核苷酸的片段耍铜。然后，排除CBs相關的隨機UMI(隨機相關定義為任何一個堿基的UMI不超過90%)
比對: hg38 genome使用STAR(2.5.3)軟件
基因定量：對轉(zhuǎn)錄本使用“-‘‘–quantMode TranscriptomeSAM”進行定量跌前。生成包含每個細胞和每個基因的UMI數(shù)量的表達矩陣棕兼。
細胞過濾：下游分析的要求細胞至少具有1,000個UMIs比對到至少500個獨特基因。我們還排除了含有超過20％線粒體或者核糖體基因的細胞抵乓。
聚類分析：質(zhì)控后獲得7698個正常的骨髓細胞（normal BM）伴挚。隨機過濾掉BM5 CD38+CD38-的1590個細胞中的783個細胞,來降低CD38+CD38-細胞對細胞群體的代表性靶衍。剩余的6915個正常的骨髓細胞被用來聚類形成細胞類型使用BackSPIN無監(jiān)督進行細胞類型聚類。BackSPIN clustering (優(yōu)點：可以克服細胞和基因同時聚類的困難).
可視化：使用KNN和t-SNE兩種方法對細胞進行可視化分析茎芋。
表達矩陣可以下載：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132151
Sequencing and read mapping
All libraries were sequenced on two runs of a NextSeq 500 (one for the in vivo samples and one for cultured samples). Raw sequencing data (FASTQ files) were processed using the previously described (
) inDrops.py bioinformatics pipeline (available at https://github.com/indrops/indrops), with a few modifications: Bowtie v.1.1.1 was used with parameter -e 80; all ambiguously mapped reads were excluded from analysis; and reads were aligned to the Ensembl release 85 Mus musculus GRCm38 reference.
Quantification and Statistical Analysis of scRNA-seq
Cell filtering and normalization
Each library was initially filtered to include only abundant barcodes (> 500 total counts for all in vivo libraries; > 800 counts for cultured sample 1; > 1,000 counts for cultured samples 2 and 3), based on visual inspection of the histograms of total transcript counts per cell barcode. Next, we excluded putatively stressed or dying cells with > 30% (in vivo samples) or > 20% (cultured samples) of their transcripts coming from mitochondrial genes.
The gene expression counts of each cell were then normalized using a variant of total-count normalization that avoids distortion from very highly expressed genes. Specifically, we calculated
, the normalized transcript counts for gene j in cell i, from the raw counts
and [圖片上傳失敗...(image-cd9a1-1564118858268)]
is the average of [圖片上傳失敗...(image-e1fc15-1564118858268)]
over all cells. To prevent very highly expressed genes from correspondingly decreasing the relative expression of other genes, we excluded genes comprising > 5% of the total counts of any cell when calculating
In order to focus on heterogeneity within the stromal cell population, we first clustered the data and excluded hematopoietic and endothelial cell clusters based on their expression of previously known marker genes, as well as putative cell doublets. Specifically, we identified genes that were highly variable (top 25% by v-score (
), a measure of above-Poisson noise) and expressed at reasonably high levels (at least 3 counts in at least 5 cells). The counts for these genes were z-score normalized and used to perform principal components analysis (PCA), keeping the top 35 dimensions. After PCA, a k-nearest-neighbor (kNN) graph (k = 4) was constructed by connecting each cell to its four nearest neighbors, using Euclidean distance in the principal component space. Finally, we applied spectral clustering (scikit-learn SpectralClustering function with assign_labels = ’discretize’) to the kNN graph and visualized the clustering by projecting the graph into two dimensions using a force-directed graph layout (SPRING(
)).
We then identified enriched genes in each cluster and assigned cell type labels based on well-characterized cell type-specific marker genes (Figure S1D). Using this approach, we excluded putative endothelial cells, granulocytes, lymphoid progenitors, megakaryocytes, and erythroid progenitors.
For the in vivo samples, we also used Scrublet(
) to identify two clusters of cell doublets that co-expressed marker genes of distinct cell types. 142 putative doublets were excluded.
Clustering and visualization of stromal cells
We repeated cell clustering and visualization using only the non-hematopoietic, non-endothelial clusters. Gene filtering, PCA, and kNN graph construction were performed as above, except only the top 25 principal components were used, and only seven spectral clusters were generated.
Permutation test for gene enrichment
To find significantly enriched genes in each cell cluster, we used a parameter-free permutation-based test to calculate p values, with the difference in means as the test statistic (
). We accounted for multiple hypotheses testing with a false discovery rate of 5% using the Benjamini-Hochberg procedure (
). To be considered for differential gene expression analysis, genes had to be expressed by at least 5 of the cells in the cluster of interest.
Gene set enrichment analysis (GSEA)
We used the online GSEA tool (http://software.broadinstitute.org/gsea/login.jsp)(
,
) to find terms enriched in cluster-specific genes. As input, we used significantly enriched genes with > 2-fold higher average expression (adding a pseudocount of 0.1 transcript counts) in the cluster of interest compared to the remaining cells. The following gene set collections were tested: H (hallmark), C2 (curated), and C5 (Gene Ontology).
Population balance analysis (PBA)
The PBA algorithm calculates a scalar “potential” for each cell that is analogous to a distance, or pseudotime, from an undifferentiated source, and a vector of fate probabilities that indicate the distance to fate branch points. These fate probabilities and temporal ordering were computed using the Python implementation of PBA (available online https://github.com/AllonKleinLab/PBA), as described(
).
The inputs to the PBA scripts are a set of comma-separated value (.csv) files encoding: the edge list of a kNN graph (k = 50) of the cell transcriptomes (A.csv); a vector assigning a net source/sink rate to each graph node (R.csv); and a lineage-specific binary matrix identifying the subset of graph nodes that reside at the tips of branches (S.csv). These files are provided online at http://kleintools.hms.harvard.edu/paper_websites/bone_marrow_stroma/. PBA is then run according to the following steps:
-
(1)
Apply the script ‘compute_Linv.py -e A.csv’, here inputting edges (flag ‘-e’) from the SPRING kNN graph (see above). This step outputs the random-walk graph Laplacian, Linv.npy.
-
(2)
Apply the script ‘compute_potential.py -L Linv.npy -R R.csv’, here inputting the inverse graph Laplacian (flag ‘-L’) computed in step (1) and the net source/sink rate to each graph node (flag ‘-R’). This step outputs a potential vector (V.npy) that is used for temporal ordering (cells ordered from high to low potential).
-
(3)
Apply the script ‘compute_fate_probabilities.py -S S.csv -V V.npy -e A.csv -D 1’, here inputting the lineage-specific exit rate matrix (flag ‘-S’), the potential (flag ‘-V’) computed in step (2), the same edges (flag ‘-e’) used in step (1) and a diffusion constant (flat ‘-D’) of 1. This step outputs fate probabilities for each cell.
Estimation of net source/sink rate vector R
A complete definition of the vector R in terms of biophysical quantities has been published previously(
). We assigned negative values to R for the five cells with the highest expression of marker genes for each of the three terminal lineages. Specifically, for each lineage, we identified genes enriched in the most mature cell cluster (cluster 1 for adipocytes, cluster 6 for osteoblasts, and cluster 7 for chondrocytes), keeping genes expressed in > 25% of cells with an average expression level of > 0.5 transcript counts and > 2-fold higher average expression within the cluster than in the rest of the cells. We then identified the five cells with the highest average z-score normalized expression of these marker genes. We used the same procedure to identify ten starting cells (cells with highest score of cluster 3 [MSC] genes). We assigned different exit rates to each of the three lineages using a fitting procedure that ensured that cells identified as the putative starting MSCs would have a uniform probability to become each fate. We assigned a single positive value to all remaining cells, with the value chosen to enforce the steady-state condition
. In the fitting procedure, all exit were initially set to one and iteratively incremented or decremented until the average fate probabilities of the putative starting MSCs were within 1% of uniform. The separate lineage exit rates were then used to form the lineage-specific exit rate matrix S.
Extracting and ordering cells for each lineage
To isolate the differentiation trajectory for each lineage (adipocyte, osteoblast, and chondrocyte), we ordered cells on the basis of their graph distance from the earliest predicted MSC progenitors, keeping only cells for which the probability of the given fate increased or remained constant with graph distance. Graph distance was measured by PBA potential, and starting with the cell closest to the MSC origin, we added the cell with next highest potential to the trajectory if the PBA-predicted lineage probability for cell i was at least 99.5% of the average lineage probability of the cell(s) already in the trajectory.
More formally, the procedure is as follows: order all N cells in the experiment from highest to lowest PBA potential V, with decreasing potential corresponding to increasing distance from MSCs. Let Ei be an indicator variable for the membership of ordered cell i in the erythroid trajectory (Ei = 1 if cell i is in the trajectory; otherwise, Ei = 0). If Pi is the PBA-predicted lineage probability for ordered cell i, then Ei = 1 if
Cells on a given lineage’s trajectory were then ordered by decreasing potential. Defining tj as the index of the jth cell on a given trajectory,
Throughout this paper, we report this cell order (akin to the “pseudotime” in other publications) as a percentage of ordered cells, with the first, least differentiated cell at 0% and the most mature cell at 100%.
Significant dynamic genes
To find genes with significant changes in expression across each lineage’s cell ordering, we used a modified version of permutation test described above (see “Permutation test for gene enrichment”). Specifically, we applied a sliding window (n = 50 cells) to the cell ordering and used the difference in means between the windows with the highest and lowest expression as the test statistic, comparing the observed difference to the differences obtained after permuting the cell ordering. To be considered for analysis, genes were required to have a mean expression of at least 0.01 transcript counts in the input cells.
Dynamic gene clustering
To find groups of genes with similar expression patterns along each lineage’s differentiation ordering, we clustered the smoothed expression traces for all significantly variable genes (see previous section) with at least two-fold change between the windows with minimum and maximum expression. In detail, we smoothed the gene expression traces using a Gaussian kernel (σ = 10% of cell ordering), z-score normalized the smoothed traces, and clustered the traces using k-means clustering.
Mapping cultured cell transcriptomes to freshly isolated cell data
For Figures 5A and S5B, cells from the cultured samples were projected into the same principal component space as the in vivo data, then mapped to their most similar in vivo neighbors. In detail, counts were first converted to TPM for all samples. Then, using only the in vivo cells, the top 25% most variable genes (measured by v-score) with at least three transcript counts in at least five cells were z-score normalized and used to find the top 35 principal components. Next, the cultured cells were z-score normalized using the gene expression means and s.d. from the in vivo data and transformed into the in vivo principal component space. Lastly, each cultured cell was mapped to its closest in vivo neighbor in principal component space (Euclidean distance). In the visualization in Figure 5A, the number of cultured cells mapping to each in vivo cell was smoothed over the kNN graph (see section “Smoothing over the kNN graph”). For Figure S5B, we compared cells in the in vivo MSC cluster to the cultured cells mapping to them.
Smoothing over the kNN graph
We smoothed data over the kNN graph for visualization of the density of cultured cells mapping to in vivo cells (Figure 5A). Smoothing was performed by diffusing the number of mapped cells over the graph, as described. In brief, if G is the kNN graph, then the smoothing operator S is
, where L is the Laplacian matrix of G,
, and expm is the matrix exponential. Then the smoothed vector
of a vector of raw values X (number of mapped cells) is
RNA velocity
In order to generate the input for Velocyto (v0.17.13) (
), which requires annotation of exons and introns for read alignments, the raw reads were reprocessed using dropEst (v0.8.5) (
). We first ran droptag with the default parameters, then aligned reads to the mouse genome (mm10) using STAR (v2.7.0a) (
), allowing unique alignments only (‘–outSAMmultNmax 1’). Then dropEst was run with default settings, aside from the following: ‘-m -V -b -F -L eiEIBA’. Cell barcodes were error-corrected using the Velocyto ‘dropest-bc-correct’ command, followed by generation of Velocyto loom files using ‘run-dropest’.
Velocyto.py was run following an example notebook (https://github.com/velocyto-team/velocyto-notebooks/blob/master/python/DentateGyrus.ipynb). Briefly, the loom files generated by dropEst were merged and then filtered to include cell barcodes used in this paper’s other analyses (see ‘Cell filtering and normalization’ section). Following gene filtering (3000 most variable genes with a minimum of 3 counts and detection in 3 cells), spliced and unspliced counts were normalized separately based on total counts per cell, with a target size of the mean total counts across cells. PCA was run using 33 components, followed by KNN imputation with 66 neighbors. Gamma fitting, RNA velocity calculations, and Markov process simulations were conducted as in the Dentate Gyrus example (code for the full analysis is available at kleintools.hms.harvard.edu/paper_websites/bone_marrow_stroma).
Monocle
Monocle (v2.10.1) (
) was run on our normalized counts matrix. Using the same gene filter as in the other analyses (see ‘Cell filtering and normalization’ section), we generated an embedding using the reduceDimension() function (‘max_components = 2, method = ”DDRTree” norm_method = ”none,” pseudo_expr = 0, relative_expr = FALSE, scaling = TRUE’). After generating an initial ordering with the orderCells() function, we identified the state corresponding to MSCs and re-ran orderCells() using this state as the root state. To compare the Monocle osteoblast cell ordering to that of PBA, we selected cells in the MSC and osteoblast states of the Monocle embedding and then ordered cells by Monocle pseudotime.
主要分群情況
識別正常骨髓樣品種的細胞類群
AML腫瘤生態(tài)系統(tǒng)的單細胞分析
臨床意義
在本報告中颅眶，我們提供了轉(zhuǎn)錄事件的關鍵見解，這些轉(zhuǎn)錄事件指導成骨細胞田弥，軟骨細胞和脂肪細胞從基質(zhì)細胞分化涛酗。這里產(chǎn)生的scRNA-seq基因表達譜允許實時描述與骨髓微環(huán)境內(nèi)的命運選擇相關的動態(tài)過程。該研究的主要結(jié)果是三種不同分化途徑的詳細表征偷厦。我們確定了每條路徑中的中間體商叹，這些中間體可以通過前瞻性選擇來測試其譜系潛力。最后只泼，我們發(fā)現(xiàn)這些種群與命運標記的譜系一致剖笙，并且它們預測的分化潛能在培養(yǎng)中被概括。
近年來在穩(wěn)態(tài)和疾病期間表征基質(zhì)群體方面取得了重大進展请唱，我們的研究為更好地理解調(diào)節(jié)骨髓微環(huán)境細胞分化的轉(zhuǎn)錄網(wǎng)絡提供了一個景觀（Hoggatt等弥咪，2016，Méndez- Ferrer等十绑，2010聚至，Mercier等，2011孽惰，Morrison和Scadden，2014鸥印，Tikhonova等勋功，2019，Baryawno等库说，2019）狂鞋。基質(zhì)細胞的失調(diào)與幾種病理生理過程有關潜的，例如肥胖骚揍，骨質(zhì)減少，骨質(zhì)疏松癥啰挪，癌癥信不，牙齒脫落和衰老（Engblom等，2017亡呵，Medyouf等抽活，2014，Mendelson和Frenette锰什，2014下硕，Raaijmakers等）丁逝。 al。梭姓，2010霜幼，Zambetti et al。誉尖，2016）罪既。因此，理解調(diào)節(jié)基質(zhì)細胞分化的機制可以提高對這些疾病的發(fā)病機理的理解释牺，并最終導致新的治療方法萝衩。
我們的偽時間結(jié)果以及轉(zhuǎn)錄本驗證支持亞群之間的關系，并允許我們探索基質(zhì)細胞表型的轉(zhuǎn)錄層次没咙。即使scRNA-seq工具及其許多應用的擴展猩谊，轉(zhuǎn)錄組快照也不能提供細胞狀態(tài)的完整圖像（Cie?lik和Chinnaiyan，2018祭刚，Kumar等牌捷，2017）。因此涡驮，采用多模態(tài)分析將有助于增強理解暗甥，因為這些技術具有測量多種分子表型的能力。此外捉捅，理解MSCs命運決定的進展將需要使用炎癥和疾病模型以及基于我們初步發(fā)現(xiàn)的患者樣本進行更深入的研究撤防。我們在這里顯示了使用命運映射和報告菌株驗證scRNA-seq數(shù)據(jù)的重要性“艨冢可以進一步利用相同的模型來發(fā)現(xiàn)特定轉(zhuǎn)錄因子在擾動期間譜系承諾決策中的重要性寄月。
總之，該研究提供了對骨髓微環(huán)境的細胞組成和沿成骨細胞无牵，軟骨細胞和脂肪細胞命運的分化途徑的轉(zhuǎn)錄中間體的重要見解漾肮。此外，盡管體內(nèi)和培養(yǎng)的基質(zhì)細胞之間存在差異茎毁，但體外分化實驗證明可用于測定不同基質(zhì)命運的轉(zhuǎn)錄因子相關性克懊。我們的數(shù)據(jù)集和分析（kleintools.hms.harvard.edu/paper_websites/bone_marrow_stroma）將作為研究基質(zhì)細胞分化的未來研究的資源。
(注：本文小編僅僅關注單細胞轉(zhuǎn)錄組七蜘，原文的單細胞基因分型也很有意義)
在本報告中谭溉，我們提供了關鍵的洞察轉(zhuǎn)錄事件，直接成骨細胞橡卤，軟骨細胞夜只，脂肪細胞分化為基質(zhì)細胞。這里生成的scrna-seq基因表達譜可以實時描述骨髓微環(huán)境中與命運選擇相關的動態(tài)過程蒜魄。這項研究的主要結(jié)果是詳細描述了三種不同的微分路徑扔亥。我們確定了沿著每一條路徑的中間產(chǎn)物场躯，這些中間產(chǎn)物可以被前瞻性地選擇來測試它們的譜系潛能。最后旅挤，我們發(fā)現(xiàn)這些群體與命運標記的譜系是一致的踢关，并且他們預測的分化潛能在文化中得到了概括。
近年來粘茄，在穩(wěn)定狀態(tài)和疾病期間基質(zhì)種群的特征化方面取得了重大進展签舞，我們的研究為更好地理解調(diào)節(jié)骨髓微環(huán)境細胞分化的轉(zhuǎn)錄網(wǎng)絡提供了一個前景（Hoggatt等人2016年，M_ndez-Ferrer等人柒瓣，2010年儒搭，Mercier等人，2011年芙贫，Morrison and Scadden搂鲫，2014年，Tikhonova等人磺平，2019年魂仍，Baryawno等人，2019年）拣挪〔磷茫基質(zhì)細胞的失調(diào)與一些病理生理過程有關，如肥胖菠劝、骨質(zhì)減少赊舶、骨質(zhì)疏松、癌癥赶诊、牙齒脫落和衰老（Engblom等人笼平，2017年，Medyouf等人甫何，2014年出吹，Mendelson和Frenette巍耗，2014年，Raaijmacers等人斋枢，2010年筒扒，Zambetti等人怯邪，2016年）。因此奔垦，了解調(diào)節(jié)間質(zhì)細胞分化的機制可能導致對這些疾病的發(fā)病機理的進一步了解摩疑，并最終獲得新的治療方法贼陶。
我們的假時間結(jié)果以及轉(zhuǎn)錄驗證支持亞群之間的關系刃泡，并允許我們探索基質(zhì)細胞表型的轉(zhuǎn)錄層次。即使scrna-seq工具及其許多應用程序的擴展碉怔，轉(zhuǎn)錄快照也無法提供細胞狀態(tài)的完整圖片（Cie_lik和Chinnaiyan烘贴，2018年，Kumar等人撮胧，2017年）桨踪。因此，采用多模分析將有助于獲得更深入的理解芹啥，因為這類技術具有測量多分子表型的能力锻离。此外，了解MSC命運決定的進展需要更深入的研究叁征，利用炎癥和疾病模型以及基于我們最初發(fā)現(xiàn)的患者樣本纳账。我們在這里展示了用命運圖和報告菌株驗證scrna-seq數(shù)據(jù)的重要性逛薇。同樣的模型也可以在微擾過程中進一步研究特定轉(zhuǎn)錄因子在譜系承諾決策中的意義捺疼。
總之，本研究對骨髓微環(huán)境的細胞組成和向成骨細胞永罚、軟骨細胞和脂肪細胞脂肪的分化途徑的轉(zhuǎn)錄中間產(chǎn)物提供了重要的見解啤呼。此外卧秘，盡管體內(nèi)和培養(yǎng)的基質(zhì)細胞之間存在差異，但體外分化實驗證明有助于分析轉(zhuǎn)錄因子與不同基質(zhì)命運的相關性官扣。我們的數(shù)據(jù)集和分析（kleintools.hms.harvard.edu/paper-websites/bone-mallow-stroma）將作為未來研究基質(zhì)細胞分化的資源翅敌。
