系列回顧：
ArchR官網(wǎng)教程學(xué)習(xí)筆記1：Getting Started with ArchR
ArchR官網(wǎng)教程學(xué)習(xí)筆記2:基于ArchR推測(cè)Doublet
ArchR官網(wǎng)教程學(xué)習(xí)筆記3:創(chuàng)建ArchRProject

這一章作者用了大量的篇幅來(lái)解釋ArchR是如何進(jìn)行降維的童漩。其中有很多專業(yè)術(shù)語(yǔ)边翁，我看著也困難，翻譯出來(lái)感覺(jué)也不是人話。谊娇。贮配。個(gè)人覺(jué)得看得懂更好掘托，看不懂也不用太糾結(jié)茅诱，主要是會(huì)運(yùn)行代碼，知道有這么個(gè)步驟就行了檩赢。

由于數(shù)據(jù)的稀疏性吕嘀，scATAC-seq的降維非常具有挑戰(zhàn)性。在scATAC-seq中贞瞒，一個(gè)特定的位點(diǎn)的可接近性可以在一個(gè)等位基因上偶房，兩個(gè)等位基因上，或者沒(méi)有等位基因军浆。即使在高質(zhì)量的scATAC-seq數(shù)據(jù)中棕洋，大多數(shù)可接近的區(qū)域沒(méi)有調(diào)換（transposed），導(dǎo)致許多l(xiāng)oci的可接近的等位基因是0瘾敢。此外拍冠，當(dāng)我們看到(例如)在單個(gè)細(xì)胞中三個(gè)Tn5插入一個(gè)peak里尿这，數(shù)據(jù)的稀疏性使我們懷疑在這個(gè)位點(diǎn)上是否真的有3倍高的可接近性（相比于其他細(xì)胞）簇抵。因此，許多分析策略都在二進(jìn)制的scATAC-seq數(shù)據(jù)矩陣上工作射众。這個(gè)二進(jìn)制的矩陣仍然大部分是0碟摆。然而，需要注意的是叨橱，scATAC-seq中的0可能意味著“不可接近”或“未被采樣”典蜕，而且從生物學(xué)角度來(lái)看，這兩種推斷是非常不同的罗洗。正因?yàn)槿绱耍?代表有信息而0沒(méi)有愉舔。這種低信息量使得scATAC-seq數(shù)據(jù)變得稀疏。

如果你要在這個(gè)稀疏的矩陣上執(zhí)行一個(gè)標(biāo)準(zhǔn)的降維操作伙菜，如主成分分析（PCA）轩缤，并繪制前兩個(gè)主成分，你不會(huì)得到想要的結(jié)果，因?yàn)橄∈栊詫?dǎo)致在所有的0位置上的細(xì)胞間高度相似火的。為了解決這個(gè)問(wèn)題壶愤，我們使用分層降維方法。首先馏鹤，我們使用Latent Semantic Indexing(LSI)征椒，一種來(lái)自處理自然語(yǔ)言的方法，最初設(shè)計(jì)用于基于字?jǐn)?shù)相似度來(lái)評(píng)估文檔的相似性湃累。這個(gè)解決方案是為自然語(yǔ)言處理創(chuàng)建的勃救，因?yàn)閿?shù)據(jù)是稀疏的和noisy的(許多不同的單詞和許多低頻率的單詞)。Cusanovich等人Cusanovich et al. (Science 2015)首次將LSI引入scATAC-seq中治力。在scATAC-seq中剪芥，不同的樣本是documents，不同的區(qū)域/峰是words琴许。首先税肪，我們通過(guò)每個(gè)單細(xì)胞的深度歸一化來(lái)計(jì)算頻率。然后榜田，這些值按照反向文檔頻率進(jìn)行規(guī)范化（文檔頻率根據(jù)特征出現(xiàn)的頻率對(duì)特征進(jìn)行衡量益兄，以識(shí)別更“特異的”而非commonly可訪問(wèn)的特征）。由此得到TF-IDF矩陣箭券，這個(gè)矩陣反映了一個(gè)word(即區(qū)域/峰)對(duì)document(即樣品)的重要性净捅。然后，通過(guò)奇異值分解(SVD)技術(shù)辩块，識(shí)別樣本間最有價(jià)值的信息蛔六，并在較低維空間中表示出來(lái)。LSI允許你將稀疏矩陣的維度從幾千減少到幾十或幾百废亭。在此基礎(chǔ)上国章，采用UMAP或t-SNE等傳統(tǒng)的降維技術(shù)來(lái)實(shí)現(xiàn)數(shù)據(jù)的可視化。在ArchR中豆村，這些可視化方法稱為嵌入（embeddings）液兽。

（一）ArchR的LSI的實(shí)現(xiàn)

ArchR實(shí)現(xiàn)了一些不同的LSI方法，并且我們已經(jīng)在多個(gè)不同的測(cè)試數(shù)據(jù)集上對(duì)其中的許多方法進(jìn)行了基準(zhǔn)測(cè)試掌动。ArchR的默認(rèn)的LSI實(shí)現(xiàn)與Timothy Stuart在Signac中引入的方法有關(guān)四啰，該方法使用一個(gè)term頻率，該頻率被深度歸一化為常數(shù)(10,000)粗恢，然后使用逆文檔頻率歸一化柑晒，再對(duì)結(jié)果矩陣進(jìn)行對(duì)數(shù)變換(又名log(TF-IDF))。

LSI降維的關(guān)鍵輸入之一是初始矩陣眷射。到目前為止匙赞，scATAC-seq的兩種主要策略是(1)使用峰區(qū)或(2)使用全基因組恋追。然而，在峰區(qū)域使用LSI本質(zhì)上是具有挑戰(zhàn)性的罚屋，因?yàn)樵诮稻S之前苦囱，我們沒(méi)有進(jìn)行聚類（cluster）或集群特定的峰（cluster-specific peaks）。此外脾猛，在clustering之前call peaks會(huì)模糊細(xì)胞特異類型的峰撕彤。而且，在實(shí)驗(yàn)中加入新樣本時(shí)猛拴，任何union峰set都會(huì)發(fā)生變化羹铅，使得該策略不太穩(wěn)定。第二種策略是使用全基因組分片愉昆，通過(guò)使用一致且無(wú)偏差的特征集(全基因組分片)來(lái)解決這些問(wèn)題职员。然而，對(duì)所有細(xì)胞的全基因組分片會(huì)使矩陣過(guò)大跛溉。由于這個(gè)原因焊切，大多數(shù)LSI實(shí)現(xiàn)使用size大于或等于5kb的“片”。這大大降低了該方法的分辨率芳室，因?yàn)榇蠖鄶?shù)可接近區(qū)域只有幾百個(gè)bp長(zhǎng)度专肪。

由于Arrow files的設(shè)計(jì)方式，ArchR能夠使用全基因組500bp分片非晨昂睿快速地執(zhí)行LSI嚎尤。這解決了分辨率的問(wèn)題，并允許在call peaks之前鑒定clusters伍宦。問(wèn)題是500bp的bin會(huì)生成大約600萬(wàn)個(gè)特征芽死。雖然ArchR能夠通過(guò)將矩陣分組將大量數(shù)據(jù)讀取到R中，但我們也實(shí)現(xiàn)了一種“estimated LSI”方法次洼，該方法對(duì)總細(xì)胞的一個(gè)子集執(zhí)行初始維數(shù)縮減关贵。這種estimated LSI方法有兩個(gè)主要的用途- (i)加速降維，(ii)隨著細(xì)胞數(shù)的減少滓玖，將減少數(shù)據(jù)的粒度（granularity）坪哄。這種粒度的減少可以用于減少數(shù)據(jù)中的批次效應(yīng)。然而势篡，它也可能掩蓋真實(shí)的生物學(xué)效應(yīng)，所以estimated LSI方法應(yīng)該在人工監(jiān)督下使用模暗。

（二）迭代LSI

在scRNA-seq中禁悠，識(shí)別variable基因是計(jì)算降維的常用方法(如PCA)。之所以這樣做兑宇，是因?yàn)檫@些高度可變的基因在生物學(xué)上有可能是很重要的碍侦，這就減少了實(shí)驗(yàn)噪音。在 scATAC-seq 中，數(shù)據(jù)是二進(jìn)制的瓷产，因此你不能識(shí)別variable峰來(lái)進(jìn)行降維站玄。我們嘗試使用最易接近的特性作為L(zhǎng)SI的輸入。然而濒旦，當(dāng)運(yùn)行多個(gè)樣本時(shí)株旷，結(jié)果顯示噪音高，重現(xiàn)性低尔邓。為了解決這個(gè)問(wèn)題晾剖，我們引入了“迭代LSI”方法（Satpathy, Granja et al. Nature Biotechnology 2019和Granja, Klemm and McGinnis* et al. Nature Biotechnology 2019）。這種方法在最易接近的分片上計(jì)算初始LSI轉(zhuǎn)換梯嗽，并識(shí)別不存在批處理混淆的低分辨率clusters齿尽。例如，當(dāng)對(duì)外周血單核細(xì)胞進(jìn)行檢測(cè)時(shí)灯节，將識(shí)別出與主要細(xì)胞類型(T細(xì)胞循头、B細(xì)胞和單核細(xì)胞)相對(duì)應(yīng)的簇。然后ArchR計(jì)算每個(gè)cluster所有特征的平均可接近性炎疆。之后贷岸，ArchR識(shí)別出這些cluster中變化最高的峰，并再次將這些特性用于LSI磷雇。在第二次迭代中偿警，變化最高的峰更類似于在scRNA-seq LSI實(shí)現(xiàn)中使用的可變基因。用戶可以設(shè)置應(yīng)該執(zhí)行多少次LSI迭代唯笙。我們發(fā)現(xiàn)這種方法可以使觀察到的批次效應(yīng)最小化螟蒸，并允許在一個(gè)合理大小的特征矩陣上進(jìn)行降維操作。

為了在ArchR中執(zhí)行迭代LSI崩掘，我們使用addIterativeLSI()函數(shù)七嫌。默認(rèn)參數(shù)應(yīng)該涵蓋大多數(shù)情況，但我們鼓勵(lì)你探索可用的參數(shù)苞慢，以及它們?nèi)绾斡绊懩愕臄?shù)據(jù)集诵原。使用?addIterativeLSI查看更詳細(xì)的說(shuō)明。最常見(jiàn)的調(diào)整參數(shù)是iterations挽放、varFeatures和resolution绍赛。重要的是要注意到LSI是不確定性的。這意味著即使你以完全相同的方式使用完全相同的參數(shù)運(yùn)行LSI辑畦，你也不會(huì)得到完全相同的結(jié)果吗蚌。當(dāng)然，它們高度相似纯出，但不完全相同蚯妇。因此敷燎，一旦確定了理想的降維，請(qǐng)確保保存你的ArchRProject或相關(guān)的LSI信息箩言。

為了練習(xí)使用硬贯，我們要?jiǎng)?chuàng)建一個(gè)reducedDims對(duì)象，稱為“IterativeLSI”：

> projHeme2 <- addIterativeLSI(
  ArchRProj = projHeme2,
  useMatrix = "TileMatrix", 
  name = "IterativeLSI", 
  iterations = 2, 
  clusterParams = list( #See Seurat::FindClusters
    resolution = c(0.2), 
    sampleCells = 10000, 
    n.start = 10
  ), 
  varFeatures = 25000, 
  dimsToUse = 1:30
)

Checking Inputs...
ArchR logging to : ArchRLogs\ArchR-addIterativeLSI-28e836691ea-Date-2020-11-20_Time-00-32-15.log
If there is an issue, please report to github with logFile!
2020-11-20 00:32:17 : Computing Total Across All Features, 0.004 mins elapsed.
......
2020-11-20 00:35:38 : Finished Running IterativeLSI, 3.357 mins elapsed.

如果您看到下游有細(xì)微的批次效應(yīng)陨收，可做的是添加更多的LSI迭代饭豹，并從較低的初始clustering分辨率開(kāi)始，如下所示畏吓。此外墨状，可變特征的數(shù)量可以降低，以增加對(duì)更多的可變特征的關(guān)注菲饼。

出于舉例的目的肾砂，我們將把這個(gè)reducedDims對(duì)象命名為“IterativeLSI2”，但是我們不會(huì)在下游分析里使用它：

> projHeme2 <- addIterativeLSI(
    ArchRProj = projHeme2,
    useMatrix = "TileMatrix", 
    name = "IterativeLSI2", 
    iterations = 4, 
    clusterParams = list( #See Seurat::FindClusters
        resolution = c(0.1, 0.2, 0.4), 
        sampleCells = 10000, 
        n.start = 10
    ), 
    varFeatures = 15000, 
    dimsToUse = 1:30
)

（三）Estimated LSI

對(duì)于非常大的scATAC-seq數(shù)據(jù)集宏悦，ArchR可以利用LSI投射來(lái)估計(jì)LSI降維镐确。這個(gè)過(guò)程類似于迭代LSI的工作流，但是LSI過(guò)程有所不同饼煞。首先源葫，使用隨機(jī)選擇的“l(fā)andmark”細(xì)胞的子集進(jìn)行LSI降維。其次砖瞧，使用從“標(biāo)志”細(xì)胞確定的反向文檔頻率對(duì)其余細(xì)胞進(jìn)行TF-IDF規(guī)范化息堂。第三，將這些歸一化的細(xì)胞投影到由“標(biāo)志”細(xì)胞定義的SVD子空間中块促。這導(dǎo)致了基于小細(xì)胞群的LSI轉(zhuǎn)換為剩余細(xì)胞的“標(biāo)志”荣堰。這種estimated LSI過(guò)程在ArchR中是有效的，因?yàn)楫?dāng)將新的細(xì)胞投射到LSI的“標(biāo)志”細(xì)胞時(shí)竭翠，ArchR迭代地從每個(gè)樣本中讀取細(xì)胞振坚，而LSI不將它們?nèi)看鎯?chǔ)在內(nèi)存中。這種優(yōu)化使得內(nèi)存消耗最小斋扰，并進(jìn)一步提高了超大數(shù)據(jù)集的可伸縮性渡八。重要的是，所需的“標(biāo)志”細(xì)胞集大小取決于數(shù)據(jù)集中不同細(xì)胞的比例传货。

ArchR通過(guò)addIterativeLSI()函數(shù)利用設(shè)置sampleCellsFinal和projectCellsPre參數(shù)來(lái)訪問(wèn)Estimated LSI屎鳍。samplesCellsFinal指定“標(biāo)志”細(xì)胞群的大小，而projectCellsPre告訴ArchR使用該“標(biāo)志”細(xì)胞群對(duì)其余細(xì)胞進(jìn)行投影损离。

（四）使用Harmony進(jìn)行批次效應(yīng)矯正

有時(shí)迭代LSI方法不足以校正強(qiáng)的批次效應(yīng)哥艇。因此，ArchR實(shí)現(xiàn)了一個(gè)常用的批次效應(yīng)修正工具Harmony僻澎，它最初是為scRNA-seq設(shè)計(jì)的貌踏。我們提供了一個(gè)封裝器，它將降維對(duì)象從ArchR直接傳遞給HarmonyMatrix()函數(shù)窟勃。其他參數(shù)可以通過(guò)附加參數(shù)(…)直接傳遞給函數(shù)中的HarmonyMatrix()祖乳。更多細(xì)節(jié)請(qǐng)參見(jiàn)?addHarmony()。

> projHeme2 <- addHarmony(
  ArchRProj = projHeme2,
  reducedDims = "IterativeLSI",
  name = "Harmony",
  groupBy = "Sample"
)

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 3 iterations

這個(gè)過(guò)程會(huì)在我們的projHeme2對(duì)象里創(chuàng)建一個(gè)新的reducedDims對(duì)象秉氧，叫“Harmony”眷昆。