學(xué)習(xí)生信代碼的朋友可以直接跳轉(zhuǎn)到下面2.7 實(shí)戰(zhàn)案例煌寇，有完整和詳盡的代碼和分析流程轮锥。

2. 原始數(shù)據(jù)處理

在本篇中兵多，我們將介紹單細(xì)胞RNA測(cè)序(scRNA-seq)數(shù)據(jù)的“預(yù)處理preprocessing”步驟裁替。盡管這是常見的術(shù)語(yǔ)粘衬，但似乎有點(diǎn)用詞不當(dāng)荞估，因?yàn)榇诉^(guò)程涉及幾個(gè)步驟，這些步驟在開始下游分析之前至關(guān)重要稚新。 在這里勘伺，我們將主要將此處理階段稱為“原始數(shù)據(jù)處理raw data processing”，我們的重點(diǎn)將放在數(shù)據(jù)分析階段褂删，該階段從lane-demultiplexed的FASTQ文件開始飞醉，最后得到一個(gè)計(jì)數(shù)矩陣，表示每個(gè)量化細(xì)胞內(nèi)每個(gè)基因產(chǎn)生的不同分子的估計(jì)數(shù)量（圖 2.1）屯阀。

然后缅帘，該計(jì)數(shù)矩陣可作為多種方法的輸入噪裕，這些方法已開發(fā)用于使用 scRNA-seq 數(shù)據(jù)進(jìn)行的各種分析，從歸一化股毫、積分和過(guò)濾方法到推斷細(xì)胞類型的方法膳音， 分化軌跡和表達(dá)動(dòng)態(tài)。鑒于它是所有這些分析的起點(diǎn)铃诬，對(duì)該矩陣的穩(wěn)健而準(zhǔn)確的估計(jì)是支持和實(shí)現(xiàn)準(zhǔn)確而可靠的后續(xù)分析的基礎(chǔ)和關(guān)鍵步驟祭陷。原始數(shù)據(jù)處理中的根本錯(cuò)誤估計(jì)可能會(huì)導(dǎo)致更高級(jí)別分析中的無(wú)效推論，并可能妨礙基于原始數(shù)據(jù)中存在的信號(hào)的發(fā)現(xiàn)趣席，而原始數(shù)據(jù)處理中的信號(hào)已被遺漏兵志、減少或扭曲。盡管我們尋求輸入和輸出的直觀性質(zhì)宣肚，但在原始數(shù)據(jù)處理過(guò)程中出現(xiàn)了一些重要且困難的挑戰(zhàn)想罕。在這里，我們將介紹原始數(shù)據(jù)處理的基本步驟霉涨，包括讀取比對(duì)/映射read alignment/mapping按价、細(xì)胞條形碼(cell barcode)識(shí)別和校正，以及通過(guò)唯一分子標(biāo)識(shí)符(UMI)的解析來(lái)估計(jì)分子計(jì)數(shù)笙瑟。我們還將提到執(zhí)行這些處理步驟時(shí)出現(xiàn)的選擇和挑戰(zhàn)楼镐。

2.1 原始數(shù)據(jù)質(zhì)控

獲得原始FASTQ文件后，可以通過(guò)運(yùn)行質(zhì)量控制(QC)工具（例如FastQC）來(lái)快速診斷讀數(shù)本身的讀數(shù)質(zhì)量往枷、序列內(nèi)容等框产。如果運(yùn)行成功，FastQC會(huì)為每個(gè)輸入的FASTQ文件生成QC報(bào)告错洁”蓿總體而言囱晴，這份報(bào)告總結(jié)了質(zhì)量得分回梧、堿基內(nèi)容和某些相關(guān)統(tǒng)計(jì)數(shù)據(jù)精堕，以發(fā)現(xiàn)文庫(kù)制備或測(cè)序中產(chǎn)生的潛在問(wèn)題竹握。
盡管如今單細(xì)胞原始數(shù)據(jù)處理工具內(nèi)部評(píng)估了一些對(duì)單細(xì)胞數(shù)據(jù)很重要的質(zhì)量檢查诚隙，例如序列的 N 含量和映射讀數(shù)的比例枢赔，但運(yùn)行快速質(zhì)量檢查仍然是一個(gè)好習(xí)慣鉴吹，因?yàn)樗鼤?huì)評(píng)估其他有用的QC指標(biāo)驰坊。通常嗡载，這一部分由合作的測(cè)序公司完成窑多，但是熟悉基本流程和原理是大有裨益的。

2.2 比對(duì)和映射Alignment and mapping

映射或比對(duì)是單細(xì)胞原始數(shù)據(jù)處理的基本步驟洼滚。它是指確定每個(gè)測(cè)序片段的潛在起源基因座loci（例如埂息，與讀數(shù)相似的基因組或轉(zhuǎn)錄組基因座集）的過(guò)程。根據(jù)測(cè)序方案，生成的原始序列文件包含細(xì)胞水平信息（通常稱為細(xì)胞條形碼 (CB)）千康、唯一分子標(biāo)識(shí)符 (UMI) 以及從分子生成的原始cDNA序列（讀取序列）享幽。作為單細(xì)胞樣本原始數(shù)據(jù)處理的第一步（圖 2.1），準(zhǔn)確執(zhí)行映射或比對(duì)對(duì)于所有下游分析都至關(guān)重要拾弃，因?yàn)椴徽_的讀取到轉(zhuǎn)錄本/基因映射可能會(huì)導(dǎo)致錯(cuò)誤或不準(zhǔn)確的矩陣 值桩。
雖然將讀取序列映射到參考序列的時(shí)間遠(yuǎn)遠(yuǎn)早于scRNA-seq出現(xiàn)的時(shí)間，但當(dāng)前scRNA-seq樣本規(guī)模的快速增長(zhǎng)（通常為數(shù)億至數(shù)十億個(gè)讀群来弧）使得該階段的計(jì)算量特別大奔坟。此外，使用RNA-seq對(duì)齊器需要使用單獨(dú)的工具來(lái)執(zhí)行其他步驟搭盾，例如解復(fù)用demultiplexing和UMI分解咳秉。這種增加的步驟在計(jì)算上可能會(huì)很麻煩。此外鸯隅，它通常需要大量的磁盤空間來(lái)存儲(chǔ)中間文件澜建，并且處理此類中間文件所需的額外輸入和輸出操作進(jìn)一步增加了運(yùn)行時(shí)間要求。
為此蝌以，專門構(gòu)建了幾個(gè)用于對(duì)齊或映射scRNA-seq數(shù)據(jù)的專用工具炕舵，這些工具可以自動(dòng)和/或內(nèi)部處理這些額外的處理要求。Cell Ranger（10x Genomics 的商業(yè)軟件）饼灿、zUMIs幕侠、alevin帝美、RainDrop碍彭、 kallisto|bustools、STARsolo和alevin-fry提供了專門的處理方法悼潭。 雖然所有工具都為用戶提供了相對(duì)簡(jiǎn)化的界面庇忌，但這些工具在內(nèi)部使用各種不同的方法，其中一些工具生成傳統(tǒng)的中間文件（例如BAM文件）并隨后對(duì)其進(jìn)行處理舰褪，而其他工具則完全在內(nèi)存中工作或使用簡(jiǎn)化的中間步驟以減少輸入/輸出負(fù)擔(dān)皆疹。

2.2.1 比對(duì)的種類

在這里，我們介紹scRNA-seq數(shù)據(jù)映射中最常應(yīng)用的三種主要映射算法：剪接比對(duì)占拍、連續(xù)比對(duì)和輕量級(jí)映射變體略就。
首先，讓我們區(qū)分基于比對(duì)的方法和基于輕量級(jí)映射的方法（圖 2.2）晃酒。比對(duì)方法應(yīng)用一系列不同的啟發(fā)法來(lái)確定可能產(chǎn)生讀數(shù)的潛在基因座表牢，然后隨后在每個(gè)假定基因座處對(duì)讀數(shù)和參考之間的最佳核苷酸水平比對(duì)進(jìn)行評(píng)分，通常使用基于動(dòng)態(tài)編程的方法贝次。使用動(dòng)態(tài)規(guī)劃算法來(lái)解決對(duì)齊問(wèn)題有著悠久而豐富的歷史崔兴。所使用的動(dòng)態(tài)規(guī)劃算法的確切類型取決于所尋求的對(duì)齊類型。全局比對(duì)適用于要比對(duì)整個(gè)查詢序列和參考序列的情況。另一方面敲茄，當(dāng)可能僅期望查詢的子序列與引用的子序列匹配時(shí)位谋，可以應(yīng)用局部對(duì)齊。通常堰燎，最適用于短讀對(duì)齊的模型既不是完全全局的也不是完全局部的掏父，而是屬于半全局對(duì)齊的類別，其中大多數(shù)查詢預(yù)計(jì)與引用的某些子字符串對(duì)齊（這通常稱為 “擬合”對(duì)齊）秆剪。此外损同，有時(shí)通常允許對(duì)比對(duì)進(jìn)行軟剪裁，以便忽略或降低對(duì)出現(xiàn)在讀取開始或結(jié)束處的不匹配鸟款、插入或刪除的懲罰膏燃。這通常是使用“擴(kuò)展”對(duì)齊來(lái)完成的。
除了這些通用的比對(duì)技術(shù)之外何什，還設(shè)計(jì)了許多更復(fù)雜的修改和啟發(fā)式方法组哩，以使比對(duì)過(guò)程在比對(duì)基因組測(cè)序讀數(shù)的背景下更加實(shí)用和高效。例如处渣，banded alignment是一種流行的啟發(fā)式方法伶贰，包含在許多流行工具的對(duì)齊實(shí)現(xiàn)中，如果用戶對(duì)低于某個(gè)閾值的對(duì)齊分?jǐn)?shù)不感興趣罐栈，它可以避免計(jì)算動(dòng)態(tài)規(guī)劃表的大部分內(nèi)容黍衙。同樣，其他啟發(fā)式方法荠诬，如X-drop和Z-drop也很流行用于盡早修剪沒(méi)有希望的對(duì)齊琅翻。
基于比對(duì)的方法知道read的每個(gè)潛在映射的質(zhì)量，然后可以使用這些讀數(shù)來(lái)過(guò)濾與映射中的參考對(duì)齊良好的讀數(shù)柑贞》阶担基于比對(duì)的方法包括傳統(tǒng)的“完全比對(duì)”方法，如STAR和STARsolo等工具中實(shí)施的方法钧嘶，以及在salmon和alevin中實(shí)現(xiàn)的選擇性對(duì)齊等方法棠众，這些方法對(duì)映射進(jìn)行評(píng)分，但省略了最佳對(duì)齊回溯的計(jì)算有决。

基于比對(duì)的方法可以進(jìn)一步分為剪接比對(duì)和連續(xù)比對(duì)方法闸拿。剪接比對(duì)方法能夠在參考的幾個(gè)不同片段上比對(duì)序列讀數(shù)，從而允許與讀數(shù)的相應(yīng)部分良好對(duì)準(zhǔn)的參考區(qū)域之間存在潛在的大間隙书幕。這些比對(duì)方法通常在將RNA-seq讀數(shù)與基因組比對(duì)時(shí)應(yīng)用新荤。剪接比對(duì)是一個(gè)具有挑戰(zhàn)性的問(wèn)題，特別是在只有一小部分讀數(shù)跨越剪接點(diǎn)的情況下按咒。另一方面迟隅，連續(xù)比對(duì)尋找與讀取對(duì)齊良好的參考的連續(xù)子串但骨。在此類比對(duì)問(wèn)題中，盡管可以允許小的插入和缺失智袭，但通常不能容忍大的間隙奔缠，例如在剪接比對(duì)中觀察到的間隙。

2.3 細(xì)胞條形碼校正

基于液滴的單細(xì)胞分離系統(tǒng)吼野，例如10x Genomics提供的系統(tǒng)校哎，已成為研究細(xì)胞異質(zhì)性原因和后果的重要工具。在該分離系統(tǒng)中瞳步，每個(gè)捕獲細(xì)胞的RNA材料與帶有條形碼的珠子一起在水基液滴封裝內(nèi)被提取闷哆。這些珠子用獨(dú)特的寡核苷酸（稱為細(xì)胞條形碼 (CB)）標(biāo)記單個(gè)細(xì)胞的 RNA 內(nèi)容，隨后將其與從 RNA 內(nèi)容反轉(zhuǎn)錄的cDNA片段一起進(jìn)行測(cè)序单起。這些珠子含有高多樣性 DNA 條形碼抱怔，能夠?qū)?xì)胞分子內(nèi)容進(jìn)行并行條形碼編碼，并在計(jì)算機(jī)中將測(cè)序讀數(shù)解復(fù)用到各個(gè)細(xì)胞箱中嘀倒。

2.3.1 條形碼的錯(cuò)誤種類

用于單細(xì)胞分析的標(biāo)簽屈留、序列和解復(fù)用方法通常效果很好。然而测蘑，在基于液滴的文庫(kù)中觀察到的細(xì)胞條形碼 (CB) 數(shù)量可能與最初封裝的細(xì)胞數(shù)量顯著差異數(shù)倍灌危。一些主要的錯(cuò)誤來(lái)源可能導(dǎo)致這種觀察：

• 雙峰/多重峰：?jiǎn)蝹€(gè)條形碼可能與兩個(gè)或多個(gè)細(xì)胞相關(guān)聯(lián)，并導(dǎo)致細(xì)胞計(jì)數(shù)不足碳胳。
• 空液滴：可能捕獲沒(méi)有封裝細(xì)胞的液滴勇蝙，這會(huì)導(dǎo)致細(xì)胞計(jì)數(shù)過(guò)多。
• 序列錯(cuò)誤：PCR擴(kuò)增或測(cè)序過(guò)程中可能會(huì)出現(xiàn)錯(cuò)誤挨约，并可能導(dǎo)致細(xì)胞計(jì)數(shù)不足或過(guò)多味混。

用于將 RNA-seq 讀數(shù)解復(fù)用到細(xì)胞特異性 bin 中的計(jì)算工具使用各種診斷指標(biāo)來(lái)過(guò)濾掉人工數(shù)據(jù)或低質(zhì)量數(shù)據(jù)。 例如烫罩，存在多種去除環(huán)境 RNA 污染的方法 [Lun 等人惜傲，2019，Muskovic 和 Powell贝攒，2021，Young 和 Behjati时甚，2020]隘弊、雙峰檢測(cè) [Bais 和 Kostka，2019荒适，DePasquale 等人梨熙，2019， McGinnis 等人刀诬，2019咽扇，Wolock 等人，2019] 和基于核苷酸序列相似性的細(xì)胞條形碼校正序列錯(cuò)誤。

2.4 UMI分解

細(xì)胞條形碼(CB)校正后质欲，讀數(shù)要么被丟棄树埠，要么分配給校正后的CB。隨后嘶伟，我們希望量化每個(gè)校正的CB中每個(gè)基因的豐度怎憋。
由于擴(kuò)增偏差的存在，必須根據(jù)UMI對(duì)讀數(shù)進(jìn)行重復(fù)數(shù)據(jù)刪除九昧，以評(píng)估采樣分子的真實(shí)計(jì)數(shù)绊袋。此外，在嘗試執(zhí)行此估計(jì)時(shí)铸鹰，其他幾個(gè)復(fù)雜因素也需要考慮在內(nèi)癌别。
UMI重復(fù)步驟旨在識(shí)別實(shí)驗(yàn)中捕獲和測(cè)序的每個(gè)細(xì)胞中的每個(gè)原始PCR前分子的讀數(shù)和UMI集。該過(guò)程的結(jié)果是將分子計(jì)數(shù)分配給每個(gè)細(xì)胞中的每個(gè)基因蹋笼，隨后在下游分析中用作該基因的原始表達(dá)规个。我們將查看觀察到的UMI及其相關(guān)映射讀數(shù)的集合，并嘗試推斷每個(gè)基因產(chǎn)生的觀察到的分子的原始數(shù)量的過(guò)程姓建，稱為UMI解析诞仓。
為了簡(jiǎn)化說(shuō)明，在下文中速兔，比對(duì)到參考基因組（例如基因的基因組位點(diǎn)）的reads被稱為該參考基因組的reads墅拭，它們的UMI標(biāo)簽被稱為該參考基因組的UMI；由UMI標(biāo)記的讀數(shù)集稱為該UMI的讀數(shù)涣狗。一次讀取只能由一個(gè)UMI標(biāo)記谍婉，但如果它映射到多個(gè)引用，則可以屬于多個(gè)引用镀钓。此外穗熬，由于scRNA-seq中每個(gè)細(xì)胞的分子條形碼過(guò)程通常是孤立和獨(dú)立的，UMI分解將針對(duì)特定細(xì)胞進(jìn)行解釋丁溅。相同的過(guò)程通常獨(dú)立地應(yīng)用于所有細(xì)胞唤蔗。

2.4.1 UMI分解的必要性

在理想的情況下，在正確的（未改變的）UMI標(biāo)簽讀取的情況下窟赏，每個(gè)UMI的讀取唯一地映射到公共參考基因妓柜，并且UMI和預(yù)PCR分子之間存在雙射。因此涯穷，UMI重復(fù)數(shù)據(jù)刪除過(guò)程在概念上很簡(jiǎn)單：UMI的讀數(shù)是來(lái)自單個(gè)PCR前分子的PCR重復(fù)棍掐。每個(gè)基因的捕獲和測(cè)序分子的數(shù)量是針對(duì)該基因觀察到的不同UMI的數(shù)量。
然而拷况，實(shí)踐中遇到的問(wèn)題使得上述簡(jiǎn)單規(guī)則不足以識(shí)別一般UMIs的基因起源作煌，因此需要開發(fā)更復(fù)雜的模型掘殴。在這里，我們主要關(guān)注兩個(gè)挑戰(zhàn)粟誓。

• 我們討論的第一個(gè)問(wèn)題是UMI中的錯(cuò)誤奏寨。當(dāng)讀段的測(cè)序UMI標(biāo)簽包含PCR或測(cè)序過(guò)程中引入的錯(cuò)誤時(shí)，就會(huì)發(fā)生這種情況努酸。常見的UMI錯(cuò)誤包括PCR期間的核苷酸替換和測(cè)序期間的讀取錯(cuò)誤服爷。未能解決此類UMI錯(cuò)誤可能會(huì)導(dǎo)致估計(jì)的分子數(shù)量增加。
• 第二個(gè)問(wèn)題是多重映射获诈，它出現(xiàn)在讀取或UMI屬于多個(gè)引用的情況下仍源，例如多基因reads/UMI。如果UMI的不同讀取映射到不同的基因舔涎、一個(gè)讀取映射到多個(gè)基因或兩者兼而有之笼踩，就會(huì)發(fā)生這種情況。這個(gè)問(wèn)題的后果是多基因reads/UMIs的基因起源不明確亡嫌，這導(dǎo)致這些基因的采樣前PCR分子計(jì)數(shù)不確定嚎于。簡(jiǎn)單地丟棄多基因reads/UMIs可能會(huì)導(dǎo)致數(shù)據(jù)丟失或傾向于產(chǎn)生多重映射讀數(shù)的基因之間的偏差估計(jì)，例如序列相似的基因家族挟冠。

2.5 計(jì)數(shù)矩陣質(zhì)量控制

生成計(jì)數(shù)矩陣后于购，執(zhí)行質(zhì)量控制 (QC) 評(píng)估非常重要。有幾種不同的評(píng)估通常屬于質(zhì)量控制的范疇知染。通常會(huì)記錄和報(bào)告基本的全局指標(biāo)肋僧，以幫助評(píng)估測(cè)序測(cè)量本身的整體質(zhì)量。這些指標(biāo)由諸如映射reads的總分?jǐn)?shù)控淡、每個(gè)細(xì)胞觀察到的不同UMI的分布嫌吠、UMI 重復(fù)數(shù)據(jù)刪除率的分布、每個(gè)細(xì)胞檢測(cè)到的基因的分布等數(shù)量組成掺炭。這些類似的指標(biāo)通常由量化工具本身記錄辫诅，因?yàn)樗鼈冏匀划a(chǎn)生，并且可以在讀取映射涧狮、細(xì)胞條形碼校正和UMI解析過(guò)程中計(jì)算炕矮。同樣，有多種工具可以幫助組織和可視化這些基本指標(biāo)勋篓，例如 Loupe browser吧享、alevinQC，或者kb_python report譬嚣，報(bào)告，具體取決于所使用的量化方法钞它。 除了這些基本的全局指標(biāo)之外拜银，在這個(gè)分析階段殊鞭，QC指標(biāo)主要旨在幫助確定哪些細(xì)胞 (CB) 已“成功”測(cè)序，以及哪些細(xì)胞表現(xiàn)出需要過(guò)濾或糾正尼桶。

2.5.1 空滴檢測(cè)

第一個(gè)QC步驟之一是確定哪些細(xì)胞條形碼對(duì)應(yīng)于“高置信度”測(cè)序細(xì)胞操灿。在基于液滴的實(shí)驗(yàn)方案中，某些條形碼與環(huán)境RNA而不是捕獲細(xì)胞的RNA相關(guān)泵督。當(dāng)液滴無(wú)法捕獲細(xì)胞時(shí)就會(huì)發(fā)生這種情況趾盐。這些空滴仍然傾向于產(chǎn)生測(cè)序讀數(shù)，盡管這些讀數(shù)的特征看起來(lái)與與正確捕獲的細(xì)胞相對(duì)應(yīng)的條形碼相關(guān)的特征明顯不同小腊。存在許多方法來(lái)評(píng)估條形碼是否可能對(duì)應(yīng)于空液滴救鲤。一種簡(jiǎn)單的方法是檢查細(xì)胞條形碼的累積頻率圖，其中條形碼按照與其關(guān)聯(lián)的不同UMI數(shù)量的降序排列秩冈。
這導(dǎo)致了幾種專門設(shè)計(jì)用于檢測(cè)空的或損壞的液滴或通常被認(rèn)為“低質(zhì)量”的細(xì)胞的工具的開發(fā)本缠。這些工具結(jié)合了各種不同的細(xì)胞質(zhì)量測(cè)量方法，包括不同UMI的頻率入问、檢測(cè)到的基因數(shù)量以及線粒體RNA的分?jǐn)?shù)丹锹，并且通常通過(guò)對(duì)這些特征應(yīng)用統(tǒng)計(jì)模型來(lái)對(duì)高質(zhì)量細(xì)胞進(jìn)行分類 假定的空滴或損壞的細(xì)胞。這意味著通撤沂В可以對(duì)單元進(jìn)行評(píng)分楣黍，并且可以基于單元不為空或受損的估計(jì)后驗(yàn)概率來(lái)選擇最終過(guò)濾。雖然這些模型通常適用于單細(xì)胞RNA-seq數(shù)據(jù)棱烂，但可能必須應(yīng)用幾個(gè)額外的過(guò)濾器或啟發(fā)式方法才能在單核RNA-seq數(shù)據(jù)中獲得穩(wěn)健的過(guò)濾租漂，就像DropletUtils的emptyDropsCellRanger函數(shù)中那樣。

2.5.2 雙聯(lián)體檢測(cè)

除了確定哪些細(xì)胞條形碼對(duì)應(yīng)于空滴或受損細(xì)胞之外垢啼，人們還可能希望識(shí)別那些對(duì)應(yīng)于雙聯(lián)體或多聯(lián)體的細(xì)胞條形碼窜锯。當(dāng)給定的液滴捕獲兩個(gè)（雙聯(lián)體）或更多（多聯(lián)體）細(xì)胞時(shí)，這可能會(huì)導(dǎo)致這些細(xì)胞條形碼在數(shù)量上出現(xiàn)偏態(tài)分布芭析，例如它們代表的讀數(shù)和UMI數(shù)量锚扎，以及它們顯示的基因表達(dá)譜。還開發(fā)了許多工具來(lái)預(yù)測(cè)細(xì)胞條形碼的雙聯(lián)態(tài)狀態(tài)馁启。一旦檢測(cè)到驾孔，被確定為可能是雙聯(lián)體和多重聯(lián)體的細(xì)胞可以被去除或以其他方式在后續(xù)分析中進(jìn)行調(diào)整。

2.6 計(jì)數(shù)數(shù)據(jù)表示

當(dāng)完成最初的原始數(shù)據(jù)處理和質(zhì)量控制并繼續(xù)進(jìn)行后續(xù)分析時(shí)惯疙，重要的是要承認(rèn)并記住基因計(jì)數(shù)矩陣充其量只是原始樣本中測(cè)序分子的近似值翠勉。在原始數(shù)據(jù)處理流程的幾個(gè)階段，進(jìn)行簡(jiǎn)化以生成該計(jì)數(shù)矩陣霉颠。例如对碌，reads映射并不完美，細(xì)胞條形碼的修正也存在問(wèn)題蒿偎。準(zhǔn)確解析UMI尤其具有挑戰(zhàn)性朽们，與UMI相關(guān)的多重映射讀取問(wèn)題經(jīng)常被忽視怀读。

2.7 實(shí)戰(zhàn)案例

鑒于我們已經(jīng)涵蓋了各種原始數(shù)據(jù)處理方法背后的概念，我們現(xiàn)在將注意力轉(zhuǎn)向演示如何使用特定工具（在本例中為alevin-fry）來(lái)處理小型示例數(shù)據(jù)集骑脱。首先菜枷，我們需要來(lái)自單細(xì)胞實(shí)驗(yàn)的FASTQ格式的測(cè)序讀數(shù)以及讀數(shù)將被比對(duì)的參考（例如轉(zhuǎn)錄組）。通常叁丧，參考基因組分別包括FASTA和GTF格式的已測(cè)序物種的基因組序列和相應(yīng)的基因注釋啤誊。
在此示例中，我們將使用人類基因組的5號(hào)染色體及其相關(guān)基因注釋作為參考拥娄，該參考是來(lái)自10x Genomics構(gòu)建的人類參考基因組GRCh38 (GENCODE v32/Ensembl 98) reference的子集蚊锹。相應(yīng)地，我們從10x Genomics的人腦腫瘤數(shù)據(jù)集中提取映射到生成參考的讀數(shù)子集条舔。
Alevin-fry[He et al., 2022]是一種快速枫耳、準(zhǔn)確且節(jié)省內(nèi)存的單細(xì)胞和單核數(shù)據(jù)處理工具。Simpleaf是一個(gè)用rust編寫的程序孟抗，它公開了一個(gè)統(tǒng)一且簡(jiǎn)化的接口迁杨，用于使用alevin-fry處理流程的一些最常見協(xié)議和數(shù)據(jù)類型。這里我們首先展示如何使用兩個(gè)simpleaf命令處理單細(xì)胞原始數(shù)據(jù)凄硼。然后铅协，我們描述了這些 simpleaf命令對(duì)應(yīng)的完整的salmon alevin和alevin-fry命令集，以概述本節(jié)中描述的步驟發(fā)生的位置摊沉，并傳達(dá)可能的不同處理選項(xiàng)狐史。這些命令將從命令行運(yùn)行，conda將用于安裝運(yùn)行此示例所需的所有軟件说墨。

2.7.1 準(zhǔn)備

在開始之前骏全，我們?cè)诮K端中創(chuàng)建一個(gè)conda環(huán)境并安裝所需的包。Simpleaf依賴于 alevin-fry尼斧、salmon和pyroe姜贡。 它們都在bioconda上可用，并且在安裝simpleaf時(shí)會(huì)自動(dòng)安裝棺棵。

conda create -n af -y -c bioconda simpleaf
conda activate af

接下來(lái)楼咳，我們創(chuàng)建一個(gè)工作目錄af_xmpl_run，并從遠(yuǎn)程主機(jī)下載并解壓縮示例數(shù)據(jù)集烛恤。

# Create a working dir and go to the working directory
## The && operator helps execute two commands using a single line of code.
mkdir af_xmpl_run && cd af_xmpl_run

# Fetch the example dataset and CB permit list and decompress them
## The pipe operator (|) passes the output of the wget command to the tar command.
## The dash operator (-) after `tar xzf` captures the output of the first command.
## - example dataset
wget -qO- https://umd.box.com/shared/static/lx2xownlrhz3us8496tyu9c4dgade814.gz | tar xzf - --strip-components=1 -C .
## The fetched folder containing the fastq files are called toy_read_fastq.
fastq_dir="toy_read_fastq"
## The fetched folder containing the human ref files is called toy_human_ref.
ref_dir="toy_human_ref"

# Fetch CB permit list
## the right chevron (>) redirects the STDOUT to a file.
wget -qO- https://raw.githubusercontent.com/10XGenomics/cellranger/master/lib/python/cellranger/barcodes/3M-february-2018.txt.gz | gunzip - > 3M-february-2018.txt

參考文件（基因組FASTA文件和基因注釋GTF文件）和讀取記錄（FASTQ文件）準(zhǔn)備就緒后母怜，我們現(xiàn)在可以應(yīng)用上面討論的原始數(shù)據(jù)處理流程來(lái)生成基因計(jì)數(shù)矩陣。

2.7.2 簡(jiǎn)化的原始數(shù)據(jù)處理流程

Simpleaf旨在簡(jiǎn)化單細(xì)胞原始數(shù)據(jù)處理的alevin-fry接口缚柏。 它將整個(gè)處理流程封裝為兩個(gè)步驟：

• 1.simpleaf index索引所提供的參考或創(chuàng)建剪接參考（剪接轉(zhuǎn)錄本+內(nèi)含子）并對(duì)其進(jìn)行索引苹熏。

• 2.simpleaf quant將測(cè)序讀數(shù)與索引參考進(jìn)行映射，并對(duì)映射記錄進(jìn)行量化以生成基因計(jì)數(shù)矩陣。
  可以在此處找到使用simpleaf進(jìn)行映射的更多高級(jí)用法和選項(xiàng)柜裸。
  運(yùn)行simpleaf index時(shí)缕陕，如果提供了基因組FASTA文件（-f）和基因注釋GTF文件（-g）粱锐，則會(huì)生成splici引用并對(duì)其進(jìn)行索引疙挺；如果僅提供轉(zhuǎn)錄組FASTA文件（--refseq），它將直接對(duì)其進(jìn)行索引怜浅。目前铐然，我們推薦使用 splici 索引。

# simpleaf needs the environment variable ALEVIN_FRY_HOME to store configuration and data.
# For example, the paths to the underlying programs it uses and the CB permit list
mkdir alevin_fry_home & export ALEVIN_FRY_HOME='alevin_fry_home'

# the simpleaf set-paths command finds the path to the required tools and write a configuration JSON file in the ALEVIN_FRY_HOME folder.
simpleaf set-paths

# simpleaf index
# Usage: simpleaf index -o out_dir [-f genome_fasta -g gene_annotation_GTF|--refseq transcriptome_fasta] -r read_length -t number_of_threads
## The -r read_lengh is the number of sequencing cycles performed by the sequencer to generate biological reads (read2 in Illumina).
## Publicly available datasets usually have the read length in the description. Sometimes they are called the number of cycles.
simpleaf index \
-o simpleaf_index \
-f toy_human_ref/fasta/genome.fa \
-g toy_human_ref/genes/genes.gtf \
-r 90 \
-t 8

在輸出目錄simpleaf_index中恶座，ref文件夾包含splici引用搀暑；索引文件夾包含基于拼接參考構(gòu)建的salmon index。
下一步跨琳，simpleaf quant自点，使用索引目錄和映射記錄FASTQ文件來(lái)生成基因計(jì)數(shù)矩陣。該命令封裝了本節(jié)中討論的所有主要步驟脉让，包括映射桂敛、單元格條形碼校正和 UMI 解析。

# Collecting sequencing read files
## The reads1 and reads2 variables are defined by finding the filenames with the pattern "_R1_" and "_R2_" from the toy_read_fastq directory.
reads1_pat="_R1_"
reads2_pat="_R2_"

## The read files must be sorted and separated by a comma.
### The find command finds the files in the fastq_dir with the name pattern
### The sort command sorts the file names
### The awk command and the paste command together convert the file names into a comma-separated string.
reads1="$(find -L ${fastq_dir} -name "*$reads1_pat*" -type f | sort | awk -v OFS=, '{$1=$1;print}' | paste -sd,)"
reads2="$(find -L ${fastq_dir} -name "*$reads2_pat*" -type f | sort | awk -v OFS=, '{$1=$1;print}' | paste -sd,)"

# simpleaf quant
## Usage: simpleaf quant -c chemistry -t threads -1 reads1 -2 reads2 -i index -u [unspliced permit list] -r resolution -m t2g_3col -o output_dir
simpleaf quant \
-c 10xv3 -t 8 \
-1 $reads1 -2 $reads2 \
-i simpleaf_index/index \
-u -r cr-like \
-m simpleaf_index/index/t2g_3col.tsv \
-o simpleaf_quant

運(yùn)行這些命令后溅潜，可以在simpleaf_quant/af_quant/alevin文件夾中找到生成的量化信息术唬。在該目錄中，存在三個(gè)文件：quants_mat.mtx滚澜、quants_mat_cols.txt和quants_mat_rows.txt粗仓，它們分別對(duì)應(yīng)于計(jì)數(shù)矩陣、該矩陣每列的基因名稱以及校正设捐、過(guò)濾后的該矩陣的每一行細(xì)胞條形碼借浊。這些文件的尾行如下所示。這里值得注意的是alevin-fry在USA模式（未剪接萝招、剪接和模糊模式）下運(yùn)行蚂斤，因此對(duì)每個(gè)基因的剪接和未剪接狀態(tài)進(jìn)行了量化——生成的quants_mat_cols.txt文件的行數(shù)將等于已注釋基因數(shù)的三倍，對(duì)應(yīng)于每個(gè)基因的剪接（S）即寒、未剪接（U）和剪接不明確變體（A）使用的名稱橡淆。

# Each line in `quants_mat.mtx` represents
# a non-zero entry in the format row column entry
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat.mtx
138 58 1
139 9 1
139 37 1

# Each line in `quants_mat_cols.txt` is a splice status
# of a gene in the format (gene name)-(splice status)
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat_cols.txt
ENSG00000120705-A
ENSG00000198961-A
ENSG00000245526-A

# Each line in `quants_mat_rows.txt` is a corrected
# (and, potentially, filtered) cell barcode
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat_rows.txt
TTCGATTTCTGAATCG
TGCTCGTGTTCGAAGG
ACTGTGAAGAAATTGC

我們可以使用pyroe中的load_fry函數(shù)將計(jì)數(shù)矩陣作為AnnData對(duì)象加載到Python中。fishpondR包中函數(shù)loadFry已經(jīng)實(shí)現(xiàn)了類似的功能母赵。

import pyroe

quant_dir = 'simpleaf_quant/af_quant'
adata_sa = pyroe.load_fry(quant_dir)

默認(rèn)將Anndata對(duì)象的x層加載為每個(gè)基因的剪接計(jì)數(shù)和模糊計(jì)數(shù)的總和逸爵。然而，最近的工作和更新的實(shí)踐表明凹嘲，即使在單細(xì)胞RNA-seq數(shù)據(jù)中包含內(nèi)含子計(jì)數(shù)师倔，也可能會(huì)提高靈敏度并有利于下游分析。雖然利用這些信息的最佳方法是正在進(jìn)行的研究的主題周蹭，但由于alevin-fry自動(dòng)量化每個(gè)樣本中的剪接趋艘、未剪接和模糊讀段疲恢，因此可以簡(jiǎn)單地獲得包含每個(gè)基因總計(jì)數(shù)的計(jì)數(shù)矩陣，如下所示 :

import pyroe

quant_dir = 'simpleaf_quant/af_quant'
adata_usa = pyroe.load_fry(quant_dir, output_format={'X' : ['U','S','A']})

2.7.3. 完整的alevin-fry分析流程

Simpleaf使得通過(guò)幾個(gè)命令以“標(biāo)準(zhǔn)”方式處理單細(xì)胞原始數(shù)據(jù)成為可能瓷胧。接下來(lái)显拳，我們將展示如何通過(guò)調(diào)用pyroe、salmon和alevin-fry命令來(lái)生成相同的量化結(jié)果搓萧。除了教學(xué)價(jià)值之外杂数，如果僅需要重新運(yùn)行管道的一部分或者需要指定simpleaf當(dāng)前未公開的某些參數(shù)，那么了解每個(gè)步驟的確切命令將很有幫助瘸洛。
請(qǐng)注意揍移，應(yīng)提前執(zhí)行準(zhǔn)備部分中的命令。以下命令中調(diào)用的所有工具反肋，pyroe那伐、salmon和alevin-fry，在安裝simpleaf時(shí)已經(jīng)安裝完畢石蔗。

2.7.3.1 建立索引

首先罕邀，我們處理基因組FASTA文件和基因注釋GTF文件以獲得剪接索引。以下代碼塊中的命令類似于上面討論的simpleaf index命令抓督。這包括兩個(gè)步驟：

• 1.使用基因組和基因注釋文件燃少，通過(guò)調(diào)用pyroe make-splici構(gòu)建剪接參考（剪接轉(zhuǎn)錄本 + 內(nèi)含子）。
• 2.通過(guò)調(diào)用salmon index對(duì)splici引用建立索引铃在。

# make splici reference
## Usage: pyroe make-splici genome_file gtf_file read_length out_dir
## The read_lengh is the number of sequencing cycles performed by the sequencer. Ask your technician if you are not sure about it.
## Publicly available datasets usually have the read length in the description.
pyroe make-splici \
${ref_dir}/fasta/genome.fa \
${ref_dir}/genes/genes.gtf \
90 \
splici_rl90_ref

# Index the reference
## Usage: salmon index -t extend_txome.fa -i idx_out_dir -p num_threads
## The $() expression runs the command inside and puts the output in place.
## Please ensure that only one file ends with ".fa" in the `splici_ref` folder.
salmon index \
-t $(ls splici_rl90_ref/*\.fa) \
-i salmon_index \
-p 8

splici索引可以在salmon_index目錄中找到阵具。

2.7.3.2. 比對(duì)和量化

接下來(lái)，我們將通過(guò)調(diào)用salmon alevin來(lái)映射根據(jù)剪接索引記錄的測(cè)序讀數(shù)定铜。這將生成一個(gè)名為almon_alevin的輸出文件夾阳液，其中包含我們使用alevin-fry處理映射讀取所需的所有信息。

# Collect FASTQ files
## The filenames are sorted and separated by space.
reads1="$(find -L $fastq_dir -name "*$reads1_pat*" -type f | sort | awk '{$1=$1;print}' | paste -sd' ')"
reads2="$(find -L $fastq_dir -name "*$reads2_pat*" -type f | sort | awk '{$1=$1;print}' | paste -sd' ')"

# Mapping
## Usage: salmon alevin -i index_dir -l library_type -1 reads1_files -2 reads2_files -p num_threads -o output_dir
## The variable reads1 and reads2 defined above are passed in using ${}.
salmon alevin \
-i salmon_index \
-l ISR \
-1 ${reads1} \
-2 ${reads2} \
-p 8 \
-o salmon_alevin \
--chromiumV3 \
--sketch

然后揣炕，我們使用alevin-fry執(zhí)行細(xì)胞條形碼校正和UMI解析步驟帘皿。此過(guò)程涉及三個(gè) alevin-fry命令：

• 1.generate-permit-list命令用于細(xì)胞條形碼校正。
• 2.collate命令過(guò)濾掉無(wú)效的映射記錄畸陡，更正細(xì)胞條形碼并整理源自相同已更正細(xì)胞條形碼的映射記錄鹰溜。
• 3.quant命令執(zhí)行UMI解析和量化。

# Cell barcode correction
## Usage: alevin-fry generate-permit-list -u CB_permit_list -d expected_orientation -o gpl_out_dir
## Here, the reads that map to the reverse complement strand of transcripts are filtered out by specifying `-d fw`.
alevin-fry generate-permit-list \
-u 3M-february-2018.txt \
-d fw \
-i salmon_alevin \
-o alevin_fry_gpl

# Filter mapping information
## Usage: alevin-fry collate -i gpl_out_dir -r alevin_map_dir -t num_threads
alevin-fry collate \
-i alevin_fry_gpl \
-r salmon_alevin \
-t 8

# UMI resolution + quantification
## Usage: alevin-fry quant -r resolution -m txp_to_gene_mapping -i gpl_out_dir -o quant_out_dir -t num_threads
## The file ends with `3col.tsv` in the splici_ref folder will be passed to the -m argument.
## Please ensure that there is only one such file in the `splici_ref` folder.
alevin-fry quant -r cr-like \
-m $(ls splici_rl90_ref/*3col.tsv) \
-i alevin_fry_gpl \
-o alevin_fry_quant \
-t 8

運(yùn)行這些命令后丁恭，生成的量化信息可以在alevin_fry_quant/alevin中找到曹动。有關(guān)映射、CB校正和UMI分解步驟的其他相關(guān)信息可分別在salmon_alevin牲览、alevin_fry_gpl和alevin_fry_quant文件夾中找到墓陈。
在此處給出的示例中，我們演示了使用simpleaf和alevin-fry處理10x Chromium 3’ v3 數(shù)據(jù)集。Alevin-fry和simpleaf提供了許多其他選項(xiàng)來(lái)處理不同的單細(xì)胞方案贡必，包括但不限于Dropseq兔港、sci-RNA-seq3和其他10x Chromium平臺(tái) 。有關(guān)不同處理階段的可用選項(xiàng)的更全面的列表和描述可以在alevin-fry和simpleaf文檔中找到仔拟。
當(dāng)然衫樊，還有許多其他原始數(shù)據(jù)處理工具也存在類似的資源，包括 zUMIs理逊、alevin橡伞、kallisto|bustools、 STARsolo和CellRanger晋被。