各位童鞋耐齐，大家好，又到周一了蒋情，今天我們來(lái)分享一下一個(gè)批次矯正的方法-----liger埠况，在這個(gè)之前呢，我分享過(guò)Seurat多樣本整合去批次的原理棵癣，文章在Seurat包其中的FindIntegrationAnchors函數(shù)解析辕翰，分享了一些去批次的軟件，文章在批次效應(yīng),單細(xì)胞數(shù)據(jù)用Harmony算法進(jìn)行批次矯正, 今天我們來(lái)分享另外一個(gè)批次去除的方法狈谊，liger喜命，文章在Single-Cell Multi-omic Integration Compares and Contrasts Features of Brain Cell Identity, 2019年6月發(fā)表與cell（頂級(jí)期刊了），而同時(shí)在另一篇文獻(xiàn)中A benchmark of batch-effect correction methods for single-cell RNA sequencing data比較了很多種批次去除批次方法的畴，對(duì)Seurat渊抄，harmony，liger三種方法的評(píng)價(jià)最高丧裁，但是其中Seurat矯正存在嚴(yán)重的過(guò)矯正問(wèn)題，而harmony目前已經(jīng)普遍運(yùn)用含衔，我們今天分享的liger煎娇，它的獨(dú)到之處又是什么呢？為什么可以發(fā)到cell贪染，今天我們來(lái)參透它缓呛，還是原來(lái)的方法，先分享文章杭隙，后示例代碼哟绊。
SUMMARY
Defining cell types requires integrating diverse single-cell measurements from multiple experiments and biological contexts（這個(gè)不用多介紹了，一個(gè)樣本發(fā)文章的時(shí)代早就過(guò)去了）. To flexibly model singlecell datasets, we developed LIGER, an algorithm that delineates shared and dataset-specific features of cell identity. We applied it to four diverse and challenging analyses of human and mouse brain cells.
（1） First, we defined region-specific and sexually dimorphic gene expression in the mouse bed nucleus of the stria terminalis.（這個(gè)地方用到了形態(tài)學(xué)方法方面的輔助痰憎，以檢驗(yàn)整合結(jié)果的優(yōu)劣）
（2）Second, we analyzed expression in the human substantia nigra, comparing cell states in specific donors and relating cell types to those in the mouse.（跨物種之間的整合結(jié)果檢驗(yàn)）
（3）Third, we integrated in situ and singlecell expression data to spatially locate fine subtypes of cells present in the mouse frontal cortex.（原位和單細(xì)胞共同的分析檢驗(yàn)）票髓。
Finally, we jointly defined mouse cortical cell types using single-cell RNA-seq and DNA methylation profiles（DNA甲基化攀涵，這個(gè)不是我們今天的重點(diǎn)）, revealing putative mechanisms of cell-type-specific epigenomic regulation（表觀調(diào)控）. Integrative analyses using LIGER promise to accelerate investigations of celltype definition, gene regulation, and disease states（讓我們拭目以待）。
INTRODUCTION
The function of the mammalian brain is dependent upon the coordinated activity of highly specialized cell types.（第一句話就很重要洽沟，強(qiáng)調(diào)了細(xì)胞空間位置的重要性以故，這也是為什么現(xiàn)在推出10X空間轉(zhuǎn)錄組的原因）。單細(xì)胞技術(shù)have provided an unprecedented opportunity to systematically identify these cellular specializations裆操，across multiple regions怒详，in the context of perturbations，and in related species（每次讀到這里踪区，都會(huì)想空間轉(zhuǎn)錄組如果也是單細(xì)胞精度就非常完美了）昆烁，F(xiàn)urthermore, new technologies can now measure DNA methylation（甲基化的結(jié)果也是非常的重要，大家可以深入的學(xué)習(xí)缎岗，這個(gè)方面你的大牛是湯富籌（不知道名字打?qū)α藳](méi)））静尼，chromatin accessibility（這個(gè)就是ATAC），and in situ expression（原位雜交）密强，in thousands to millions of cells.（龐大的單細(xì)胞數(shù)據(jù)目前也是一個(gè)大問(wèn)題茅郎，其中張澤民團(tuán)隊(duì)研究的新冠文章細(xì)胞數(shù)量達(dá)到恐怖的百萬(wàn)級(jí)）Each of these experimental contexts and measurement modalities provides a different glimpse into cellular identity.
Integrative computational tools that can flexibly combine individual single-cell datasets into a unified, shared analysis offer many exciting biological opportunities.（整合分析的必要性），The major challenge of
integrative analysis lies in reconciling the immense heterogeneity observed across individual datasets.（現(xiàn)在不止免疫的個(gè)體異質(zhì)性了或渤，很多都設(shè)及到批次）系冗。However, in many kinds of analysis, both dataset similarities and differences are biologically important, such as when we seek to compare and contrast scRNA-seq data from healthy and disease-affected individuals。
To address these challenges, we developed a new computational method called LIGER (linked inference of genomic experimental relationships). We show here that LIGER enables the identification of shared cell types across individuals, species, and multiple modalities (gene expression, epigenetic, or spatial data), as well as dataset-specific features, offering a unified analysis of heterogeneous single-cell datasets.（在這里我們只關(guān)注樣本的差異去除薪鹦，至于物種可以了解一下）掌敬。
Result1 Comparing and Contrasting Single-Cell Datasets with Shared and Dataset-Specific Factors
LIGER takes as input multiple single-cell datasets, which may be scRNA-seq experiments from different individuals, time points, or species—or measurements from different molecular modalities, such as single-cell epigenome data or spatial gene expression data（個(gè)體，物種池磁，技術(shù)）
LIGER then employs integrative non-negative matrix factorization (iNMF)（不知道大家對(duì)非負(fù)矩陣分解有多少了解） to learn a low-dimensional space in which each cell is defined by one set of dataset-specific factors, or metagenes, and another set of shared metagenes奔害。
這個(gè)圖片應(yīng)該不陌生吧，降維的時(shí)候用到的地熄。
Each factor often corresponds to a biologically interpretable signal—like the genes that define a particular cell type.（不知道大家對(duì)降維后每個(gè)主成分的研究有多少）华临。A tuning parameter; λ; allows adjusting the
size of dataset-specific effects to reflect the divergence of the datasets being analyzed. We found that iNMF performs comparably to both NMF and principal-component analysis (PCA) in reconstructing the original data（肯定是差不多的）
，After performing iNMF, we use a novel strategy that increases robustness of joint clustering.We first assign each cell a label based on themaximum factor loading and then build a shared factor neighborhood graph 端考，in which we connect cells that have similar factor loading patterns（這個(gè)很智慧雅潭，需要格外注意）。
We derived a novel algorithm for iNMF optimization, which scales well with the size of large single-cell datasets却特，To aid selection of the key parameters—the number of factors k and the tuning parameter λ—we developed heuristics based on factor entropy and dataset alignment（看看如何啟發(fā)的）
Overall, these heuristics performed well across different analyses 扶供，though we have observed that manual tuning can sometimes improve the results.（下圖）：
Additionally, we derived novel algorithms for rapidly updating the factorization to incorporate new data or change parameters，We anticipate that this capability will be useful for leveraging a rapidly growing corpus of single-cell data.（其實(shí)一但涉及人工調(diào)整參數(shù)裂明，就需要我們有很強(qiáng)的生物學(xué)背景）椿浓。
Result2 Liger Shows Robust Performance on Highly Divergent Datasets
We assessed the performance of LIGER through the use of two metrics: alignment and agreement（這里應(yīng)該理解為指標(biāo)）。Alignment measures the uniformity of mixing for two or more samples in the aligned latent space.（衡量對(duì)齊的潛在空間中兩個(gè)或多個(gè)樣本的混合均勻性。）扳碍，This metric should be high when datasets share underlying cell types, and low when datasets do not share cognate populations.（我們暫且記住這個(gè)注釋）提岔。The second metric,agreement, quantifies the similarity of each cell’s neighborhood when a dataset is analyzed separately versus jointly with other datasets（量化相似性）。High agreement indicates that cell-type relationships are preserved with minimal distortion in the joint analysis.（高度aggrement表明左腔，在聯(lián)合分析中唧垦，細(xì)胞類型關(guān)系得以保留，并且失真最小）液样。
We calculated alignment and agreement metrics using published datasets振亮，comparing the LIGER analyses to those generated by the Seurat package（和Seurat相比較），We first ran our analyses on a pair of scRNA-seq datasets from human blood cells that show primarily technical differences（技術(shù)上帶來(lái)的批次）鞭莽，and should thus yield a high degree of alignment. Indeed, LIGER and Seurat show similarly high alignment statistics坊秸，and LIGER’s joint clusters match the published cluster assignments for the individual datasets
這個(gè)地方其實(shí)很容易實(shí)現(xiàn)，Seurat過(guò)矯正澎怒，所以對(duì)于擁有相同細(xì)胞數(shù)據(jù)類型的整合當(dāng)然效果很好褒搔，而liger采用主成分連接式的方法，當(dāng)然結(jié)果也很不錯(cuò)喷面。
LIGER and Seurat also performed similarly when integrating human and mouse pancreatic data, with LIGER showing slightly higher alignment（跨物種整合星瘾，這個(gè)不太建議大家做）。
In both analyses, LIGER produced considerably higher agreement than Seurat惧辈，suggesting better preservation of the underlying cell-type architectures in the integrated space.We expected this advantage should be especially beneficial when analyzing very divergent datasets that share few or no common cell populations.（這個(gè)地方才是我們需要關(guān)注的重點(diǎn)琳状，看看結(jié)果），LIGER generated minimal false alignment between these classes and demonstrated a good preservation of complex internal substructure（這個(gè)時(shí)候就要體現(xiàn)數(shù)據(jù)之間的差異性了）盒齿。
even across considerable changes in dataset proportion念逞，
In each of the three analyses described above, the LIGER joint clustering result closely matched the published cluster assignments for the individual datasets（匹配性很好）。
Together, these analyses indicate that LIGER sensitively detects common populations without spurious alignment and preserves complex substructure, even when applied across divergent datasets（這個(gè)地方我們最希望看到的就是不同樣本之間的生物學(xué)差異）边翁。
Result 3 Interpretable Factors Unravel Complex and Dimorphic Expression Patterns in the Bed Nucleus
An important application of integrative single-cell analysis in neuroscience is to quantify cell-type variation across different brain regions and different members of the same species. To examine LIGER’s performance in these tasks, we analyzed the bed nucleus of the stria terminalis (BNST), a subcortical region composed of multiple subnuclei翎承，implicated in social, stress-related, and reward behaviors，To date, scRNA-seq has not been performed on BNST, providing an opportunity to clarify how cell types are shared between this structure and datasets generated from related tissues.
We isolated, sequenced, and analyzed 204,737 nuclei enriched for the BNST region符匾。Initial clustering identified 106,728 neurons, of which 70.2% were localized to BNST by examination of marker expression in the Allen Brain Atlas (ABA)叨咖，Clustering analysis revealed 41 transcriptionally distinct populations of BNST-localized neurons（單純的聚類分析）
In agreement with previous estimates(跟之前估計(jì)的一致)，85.9% of the cells were inhibitory
(expressing Gad1 and Gad2), while the remaining 14% were excitatory (expressing Slc17a6 [9.4%] or Slc17a8 [4.7%])
Examination of cluster markers in the ABA showed that many cell types localized to specific BNST substructures, including the principal, oval, and anterior commissure nuclei啊胶。For example, we identified two molecularly distinct subpopulations in the oval nucleus of the anterior BNST (ovBNST)芒澜，a structure known to regulate anxiety and to manifest（顯現(xiàn)） a robust circadian rhythm of expression of Per2, similar to the superchiasmatic nucleus (SCN) of the hypothalamus.
Two clusters, BNST_Vip and BNSTp_Cplx3, expressed markers of caudal ganglionic eminence (CGE)-derived interneurons found in cortex and hippocampus. Part of the BNST has embryonic origins in the CGE，suggesting that this structure may harbor such cell types. To examine this possibility, we integrated the 352 nuclei from the BNST_Vip and BNST_Cplx3 clusters with 330 CGE interneuron cell profiles
sampled from our recent adult mouse frontal cortex dataset（我們來(lái)看看整合結(jié)果）创淡， Four clusters in the LIGER analysis showed meaningful alignment between BNST nuclei and cortical CGE cells。
One population (cluster 1), which was Vip-negative and likely localized to the posterior BNST , expressed Id2, Lamp5, Cplx3, and Npy, all markers known to be present in cortical neurogliaform (NG) cells南吮。
A second population (cluster 2) expressed Vip, Htr3a, Cck, and Cnr1, likely corresponding to VIP+ basket cells琳彩。
這個(gè)地方設(shè)及到因子分析，不知道大家沒(méi)有過(guò)多的分析過(guò)，我們下一篇文章分享這個(gè)露乏，但是這里的差異分析大家要關(guān)注一下碧浊，不知道大家知不知道這個(gè)差異基因排序的原理以及什么軟件做的，知道的話瘟仿，恭喜你箱锐，趕緊嘗試一下吧。
Result 4 Integration of Substantia Nigra Profiles across Different Human Postmortem Donors and Species
這里對(duì)不同技術(shù)進(jìn)行整合劳较，我們簡(jiǎn)單看一下
我們提煉一下其中的信息：
（1）LIGER accurately integrated each of the cell-type substituents of the SN across datasets
（2）這部分也包括多人和小鼠跨物種的數(shù)據(jù)整合驹止，總之效果不錯(cuò)。
Result5 Integrating scRNA-Seq and In Situ Transcriptomic Data Locates Frontal Cortex Cell Types in Space（這部分結(jié)果涉及到空間信息观蜗，我們看一下）
Result 6 LIGER Defines Cell Types Using Both Single-Cell Transcriptome and Single-Cell DNA Methylation Profiles（這部分涉及到甲基化的分析結(jié)果臊恋，我們簡(jiǎn)單帶過(guò)）
Method（我們來(lái)看一下其中關(guān)鍵的方法）
（1）Performing iNMF Using Block Coordinate Descent
哇，高深的數(shù)學(xué)讓大腦停止運(yùn)轉(zhuǎn)了墓捻。
（2）Efficient Updating for New K, New Data, and New λ
（3）Calculating alignment and agreement metrics
最后抖仅，我們看一下示例代碼
To perform online iNMF, we need to install the latest Liger package from GitHub. Please see the instruction below.
library(devtools)
install_github("welch-lab/liger")
We first create a Liger object by passing the filenames of HDF5 files containing the raw count data. The data can be downloaded here. Liger assumes by default that the HDF5 files are formatted by the 10X CellRanger pipeline. Large datasets are often generated over multiple 10X runs (for example, multiple biological replicates). In such cases it may be necessary to merge the HDF5 files from each run into a single HDF5 file. We provide the mergeH5 function for this purpose (see below for details).
library(rliger)
pbmcs = createLiger(list(stim = "pbmcs_stim.h5", ctrl = "pbmcs_ctrl.h5"))
We then perform the normalization, gene selection, and gene scaling in an online fashion, reading the data from disk in small batches.
pbmcs = normalize(pbmcs)
pbmcs = selectGenes(pbmcs, var.thresh = 0.2, do.plot = F)
pbmcs = scaleNotCenter(pbmcs)
Online Integrative Nonnegative Matrix Factorization
Now we can use online iNMF to factorize the data, again using only minibatches that we read from the HDF5 files on demand (default mini-batch size = 5000). Sufficient number of iterations is crucial for obtaining ideal factorization result. If the size of the mini-batch is set to be close to the size of the whole dataset (i.e. an epoch only contains one iteration), max.epochs needs to be increased accordingly for more iterations.
pbmcs = online_iNMF(pbmcs, k = 20, miniBatch_size = 5000, max.epochs = 5)
Quantile Normalization and Downstream Analysis
After performing the factorization, we can perform quantile normalization to align the datasets.
pbmcs = quantile_norm(pbmcs)
We can also visualize the cell factor loadings in two dimensions using t-SNE or UMAP.
pbmcs = runUMAP(pbmcs)
plotByDatasetAndCluster(pbmcs, axis.labels = c("UMAP1","UMAP2"))
Let’s first evaluate the quality of data alignment. The alignment score ranges from 0 (no alignment) to 1 (perfect alignment).
calcAlignment(pbmcs)
## [1] 0.9149238
With HDF5 files as input, we need to sample the raw, normalized, or scaled data from the full dataset on disk and load them in memory. Some plotting functions and downstream analyses are designed to operate on a subset of cells sampled from the full dataset. This enables rapid analysis using limited memory. The readSubset function allows either uniform random sampling or sampling balanced by cluster. Here we extract the normalized count data of 5000 sampled cells.
pbmcs = readSubset(pbmcs, slot.use = "norm.data", max.cells = 5000)
Using the sampled data stored in memory, we can now compare clusters or datasets (within each cluster) to identify differentially expressed genes. The runWilcoxon function performs differential expression analysis by performing an in-memory Wilcoxon rank-sum test on this subset. Thus, users can still analyze large datasets with a fixed amount of memory.
de_genes = runWilcoxon(pbmcs, compare.method = "datasets")
Here we show the top 10 genes in cluster 1 whose expression level significantly differ between two dataset.
de_genes = de_genes[order(de_genes$padj), ]
head(de_genes[de_genes$group == "1-stim",], n = 10)
## feature group avgExpr logFC statistic auc pval
## 14061 ISG15 1-stim -5.430433 13.820208 121239.0 0.9898839 6.419043e-114
## 14324 IFI6 1-stim -6.843238 14.763583 118624.5 0.9685372 2.942152e-108
## 23989 ISG20 1-stim -5.779989 9.286306 117816.5 0.9619401 4.421875e-99
## 21242 IFIT3 1-stim -8.320794 14.539736 114625.0 0.9358824 1.133535e-98
## 21244 IFIT1 1-stim -9.094809 13.669177 112456.0 0.9181731 2.367054e-92
## 27532 MX1 1-stim -8.935794 12.875135 111536.5 0.9106656 5.884017e-87
## 20364 LY6E 1-stim -8.721669 12.810795 111435.0 0.9098369 9.506436e-86
## 23754 B2M 1-stim -3.188800 0.461932 108547.5 0.8862612 3.907049e-69
## 21241 IFIT2 1-stim -11.928663 10.782880 102299.0 0.8352439 3.361002e-66
## 14634 IFI44L 1-stim -12.345872 9.720498 98779.5 0.8065081 3.120570e-55
## padj pct_in pct_out
## 14061 9.020681e-110 100 100
## 14324 2.067303e-104 100 100
## 23989 2.071354e-95 100 100
## 21242 3.982393e-95 100 100
## 21244 6.652841e-89 100 100
## 27532 1.378135e-83 100 100
## 20364 1.908485e-82 100 100
## 23754 6.863219e-66 100 100
## 21241 5.248019e-63 100 100
## 14634 4.385337e-52 100 100
We can show UMAP coordinates of sampled cells by their loadings on each factor (Factor 1 as an example). Underneath it displays the most highly loading shared and dataset-specific genes, with the size of the marker indicating the magnitude of the loading.
p_wordClouds = plotWordClouds(pbmcs, num.genes = 5, return.plots = T)
p_wordClouds[[1]]
We can generate plots of dimensional reduction coordinates colored by expression of specified gene. The first two UMAP dimensions and gene ISG15 (identified by Wilcoxon test in the previous step) is used as an example here.
plotGene(pbmcs, "ISG15", return.plots = F)
Furthermore, we can make violin plots of expression of specified gene for each dataset (ISG15 as an example).
plotGeneViolin(pbmcs, "ISG15", return.plots = F)
The online algorithm can be implemented on datasets loaded in memory as well. The same analysis is performed on the PBMCs, shown below.
stim = readRDS("pbmcs_stim.RDS")
ctrl = readRDS("pbmcs_ctrl.RDS")
pbmcs_mem = createLiger(list(stim = stim, ctrl = ctrl), remove.missing = F)
pbmcs_mem = normalize(pbmcs_mem)
pbmcs_mem = selectGenes(pbmcs_mem, var.thresh = 0.2, do.plot = F)
pbmcs_mem = scaleNotCenter(pbmcs_mem)
pbmcs_mem = online_iNMF(pbmcs_mem, k = 20, miniBatch_size = 5000, max.epochs = 5)
pbmcs_mem = quantile_norm(pbmcs_mem)
pbmcs_mem = runUMAP(pbmcs_mem)
plotByDatasetAndCluster(pbmcs_mem, axis.labels = c("UMAP1","UMAP2"))
如果有條件的話，不妨試一下砖第，如何靈活運(yùn)用這個(gè)軟件撤卢，就看大家的需求了
生活很好，有你更好
