单细胞 RNA-seq 原始信分析 - 第二部分：原始数据处理

学习生信代码的朋友可以直接跳转到下面2.7 实战案例，有完整和详尽的代码和分析流程。

2. 原始数据处理

在本篇中，我们将介绍单细胞RNA测序(scRNA-seq)数据的“预处理preprocessing”步骤。尽管这是常见的术语，但似乎有点用词不当，因为此过程涉及几个步骤，这些步骤在开始下游分析之前至关重要。 在这里，我们将主要将此处理阶段称为“原始数据处理raw data processing”，我们的重点将放在数据分析阶段，该阶段从lane-demultiplexed的FASTQ文件开始，最后得到一个计数矩阵，表示每个量化细胞内每个基因产生的不同分子的估计数量（图 2.1）。

然后，该计数矩阵可作为多种方法的输入，这些方法已开发用于使用 scRNA-seq 数据进行的各种分析，从归一化、积分和过滤方法到推断细胞类型的方法， 分化轨迹和表达动态。鉴于它是所有这些分析的起点，对该矩阵的稳健而准确的估计是支持和实现准确而可靠的后续分析的基础和关键步骤。原始数据处理中的根本错误估计可能会导致更高级别分析中的无效推论，并可能妨碍基于原始数据中存在的信号的发现，而原始数据处理中的信号已被遗漏、减少或扭曲。尽管我们寻求输入和输出的直观性质，但在原始数据处理过程中出现了一些重要且困难的挑战。在这里，我们将介绍原始数据处理的基本步骤，包括读取比对/映射read alignment/mapping、细胞条形码(cell barcode)识别和校正，以及通过唯一分子标识符(UMI)的解析来估计分子计数。我们还将提到执行这些处理步骤时出现的选择和挑战。

2.1 原始数据质控

获得原始FASTQ文件后，可以通过运行质量控制(QC)工具（例如FastQC）来快速诊断读数本身的读数质量、序列内容等。如果运行成功，FastQC会为每个输入的FASTQ文件生成QC报告。总体而言，这份报告总结了质量得分、碱基内容和某些相关统计数据，以发现文库制备或测序中产生的潜在问题。
尽管如今单细胞原始数据处理工具内部评估了一些对单细胞数据很重要的质量检查，例如序列的 N 含量和映射读数的比例，但运行快速质量检查仍然是一个好习惯，因为它会评估其他有用的QC指标。通常，这一部分由合作的测序公司完成，但是熟悉基本流程和原理是大有裨益的。

2.2 比对和映射Alignment and mapping

映射或比对是单细胞原始数据处理的基本步骤。它是指确定每个测序片段的潜在起源基因座loci（例如，与读数相似的基因组或转录组基因座集）的过程。根据测序方案，生成的原始序列文件包含细胞水平信息（通常称为细胞条形码 (CB)）、唯一分子标识符 (UMI) 以及从分子生成的原始cDNA序列（读取序列）。作为单细胞样本原始数据处理的第一步（图 2.1），准确执行映射或比对对于所有下游分析都至关重要，因为不正确的读取到转录本/基因映射可能会导致错误或不准确的矩阵 。
虽然将读取序列映射到参考序列的时间远远早于scRNA-seq出现的时间，但当前scRNA-seq样本规模的快速增长（通常为数亿至数十亿个读取）使得该阶段的计算量特别大。此外，使用RNA-seq对齐器需要使用单独的工具来执行其他步骤，例如解复用demultiplexing和UMI分解。这种增加的步骤在计算上可能会很麻烦。此外，它通常需要大量的磁盘空间来存储中间文件，并且处理此类中间文件所需的额外输入和输出操作进一步增加了运行时间要求。
为此，专门构建了几个用于对齐或映射scRNA-seq数据的专用工具，这些工具可以自动和/或内部处理这些额外的处理要求。Cell Ranger（10x Genomics 的商业软件）、zUMIs、alevin、RainDrop、 kallisto|bustools、STARsolo和alevin-fry提供了专门的处理方法。 虽然所有工具都为用户提供了相对简化的界面，但这些工具在内部使用各种不同的方法，其中一些工具生成传统的中间文件（例如BAM文件）并随后对其进行处理，而其他工具则完全在内存中工作或使用简化的中间步骤以减少输入/输出负担。

2.2.1 比对的种类

在这里，我们介绍scRNA-seq数据映射中最常应用的三种主要映射算法：剪接比对、连续比对和轻量级映射变体。
首先，让我们区分基于比对的方法和基于轻量级映射的方法（图 2.2）。比对方法应用一系列不同的启发法来确定可能产生读数的潜在基因座，然后随后在每个假定基因座处对读数和参考之间的最佳核苷酸水平比对进行评分，通常使用基于动态编程的方法。使用动态规划算法来解决对齐问题有着悠久而丰富的历史。所使用的动态规划算法的确切类型取决于所寻求的对齐类型。全局比对适用于要比对整个查询序列和参考序列的情况。另一方面，当可能仅期望查询的子序列与引用的子序列匹配时，可以应用局部对齐。通常，最适用于短读对齐的模型既不是完全全局的也不是完全局部的，而是属于半全局对齐的类别，其中大多数查询预计与引用的某些子字符串对齐（这通常称为 “拟合”对齐）。此外，有时通常允许对比对进行软剪裁，以便忽略或降低对出现在读取开始或结束处的不匹配、插入或删除的惩罚。这通常是使用“扩展”对齐来完成的。
除了这些通用的比对技术之外，还设计了许多更复杂的修改和启发式方法，以使比对过程在比对基因组测序读数的背景下更加实用和高效。例如，banded alignment是一种流行的启发式方法，包含在许多流行工具的对齐实现中，如果用户对低于某个阈值的对齐分数不感兴趣，它可以避免计算动态规划表的大部分内容。同样，其他启发式方法，如X-drop和Z-drop也很流行用于尽早修剪没有希望的对齐。
基于比对的方法知道read的每个潜在映射的质量，然后可以使用这些读数来过滤与映射中的参考对齐良好的读数。基于比对的方法包括传统的“完全比对”方法，如STAR和STARsolo等工具中实施的方法，以及在salmon和alevin中实现的选择性对齐等方法，这些方法对映射进行评分，但省略了最佳对齐回溯的计算。

基于比对的方法可以进一步分为剪接比对和连续比对方法。剪接比对方法能够在参考的几个不同片段上比对序列读数，从而允许与读数的相应部分良好对准的参考区域之间存在潜在的大间隙。这些比对方法通常在将RNA-seq读数与基因组比对时应用。剪接比对是一个具有挑战性的问题，特别是在只有一小部分读数跨越剪接点的情况下。另一方面，连续比对寻找与读取对齐良好的参考的连续子串。在此类比对问题中，尽管可以允许小的插入和缺失，但通常不能容忍大的间隙，例如在剪接比对中观察到的间隙。

2.3 细胞条形码校正

基于液滴的单细胞分离系统，例如10x Genomics提供的系统，已成为研究细胞异质性原因和后果的重要工具。在该分离系统中，每个捕获细胞的RNA材料与带有条形码的珠子一起在水基液滴封装内被提取。这些珠子用独特的寡核苷酸（称为细胞条形码 (CB)）标记单个细胞的 RNA 内容，随后将其与从 RNA 内容反转录的cDNA片段一起进行测序。这些珠子含有高多样性 DNA 条形码，能够对细胞分子内容进行并行条形码编码，并在计算机中将测序读数解复用到各个细胞箱中。

2.3.1 条形码的错误种类

用于单细胞分析的标签、序列和解复用方法通常效果很好。然而，在基于液滴的文库中观察到的细胞条形码 (CB) 数量可能与最初封装的细胞数量显着差异数倍。一些主要的错误来源可能导致这种观察：

• 双峰/多重峰：单个条形码可能与两个或多个细胞相关联，并导致细胞计数不足。
• 空液滴：可能捕获没有封装细胞的液滴，这会导致细胞计数过多。
• 序列错误：PCR扩增或测序过程中可能会出现错误，并可能导致细胞计数不足或过多。

用于将 RNA-seq 读数解复用到细胞特异性 bin 中的计算工具使用各种诊断指标来过滤掉人工数据或低质量数据。 例如，存在多种去除环境 RNA 污染的方法 [Lun 等人，2019，Muskovic 和 Powell，2021，Young 和 Behjati，2020]、双峰检测 [Bais 和 Kostka，2019，DePasquale 等人，2019， McGinnis 等人，2019，Wolock 等人，2019] 和基于核苷酸序列相似性的细胞条形码校正序列错误。

2.4 UMI分解

细胞条形码(CB)校正后，读数要么被丢弃，要么分配给校正后的CB。随后，我们希望量化每个校正的CB中每个基因的丰度。
由于扩增偏差的存在，必须根据UMI对读数进行重复数据删除，以评估采样分子的真实计数。此外，在尝试执行此估计时，其他几个复杂因素也需要考虑在内。
UMI重复步骤旨在识别实验中捕获和测序的每个细胞中的每个原始PCR前分子的读数和UMI集。该过程的结果是将分子计数分配给每个细胞中的每个基因，随后在下游分析中用作该基因的原始表达。我们将查看观察到的UMI及其相关映射读数的集合，并尝试推断每个基因产生的观察到的分子的原始数量的过程，称为UMI解析。
为了简化说明，在下文中，比对到参考基因组（例如基因的基因组位点）的reads被称为该参考基因组的reads，它们的UMI标签被称为该参考基因组的UMI；由UMI标记的读数集称为该UMI的读数。一次读取只能由一个UMI标记，但如果它映射到多个引用，则可以属于多个引用。此外，由于scRNA-seq中每个细胞的分子条形码过程通常是孤立和独立的，UMI分解将针对特定细胞进行解释。相同的过程通常独立地应用于所有细胞。

2.4.1 UMI分解的必要性

在理想的情况下，在正确的（未改变的）UMI标签读取的情况下，每个UMI的读取唯一地映射到公共参考基因，并且UMI和预PCR分子之间存在双射。因此，UMI重复数据删除过程在概念上很简单：UMI的读数是来自单个PCR前分子的PCR重复。每个基因的捕获和测序分子的数量是针对该基因观察到的不同UMI的数量。
然而，实践中遇到的问题使得上述简单规则不足以识别一般UMIs的基因起源，因此需要开发更复杂的模型。在这里，我们主要关注两个挑战。

• 我们讨论的第一个问题是UMI中的错误。当读段的测序UMI标签包含PCR或测序过程中引入的错误时，就会发生这种情况。常见的UMI错误包括PCR期间的核苷酸替换和测序期间的读取错误。未能解决此类UMI错误可能会导致估计的分子数量增加。
• 第二个问题是多重映射，它出现在读取或UMI属于多个引用的情况下，例如多基因reads/UMI。如果UMI的不同读取映射到不同的基因、一个读取映射到多个基因或两者兼而有之，就会发生这种情况。这个问题的后果是多基因reads/UMIs的基因起源不明确，这导致这些基因的采样前PCR分子计数不确定。简单地丢弃多基因reads/UMIs可能会导致数据丢失或倾向于产生多重映射读数的基因之间的偏差估计，例如序列相似的基因家族。

2.5 计数矩阵质量控制

生成计数矩阵后，执行质量控制 (QC) 评估非常重要。有几种不同的评估通常属于质量控制的范畴。通常会记录和报告基本的全局指标，以帮助评估测序测量本身的整体质量。这些指标由诸如映射reads的总分数、每个细胞观察到的不同UMI的分布、UMI 重复数据删除率的分布、每个细胞检测到的基因的分布等数量组成。这些类似的指标通常由量化工具本身记录，因为它们自然产生，并且可以在读取映射、细胞条形码校正和UMI解析过程中计算。同样，有多种工具可以帮助组织和可视化这些基本指标，例如 Loupe browser、alevinQC，或者kb_python report，报告，具体取决于所使用的量化方法。 除了这些基本的全局指标之外，在这个分析阶段，QC指标主要旨在帮助确定哪些细胞 (CB) 已“成功”测序，以及哪些细胞表现出需要过滤或纠正。

2.5.1 空滴检测

第一个QC步骤之一是确定哪些细胞条形码对应于“高置信度”测序细胞。在基于液滴的实验方案中，某些条形码与环境RNA而不是捕获细胞的RNA相关。当液滴无法捕获细胞时就会发生这种情况。这些空滴仍然倾向于产生测序读数，尽管这些读数的特征看起来与与正确捕获的细胞相对应的条形码相关的特征明显不同。存在许多方法来评估条形码是否可能对应于空液滴。一种简单的方法是检查细胞条形码的累积频率图，其中条形码按照与其关联的不同UMI数量的降序排列。
这导致了几种专门设计用于检测空的或损坏的液滴或通常被认为“低质量”的细胞的工具的开发。这些工具结合了各种不同的细胞质量测量方法，包括不同UMI的频率、检测到的基因数量以及线粒体RNA的分数，并且通常通过对这些特征应用统计模型来对高质量细胞进行分类 假定的空滴或损坏的细胞。这意味着通常可以对单元进行评分，并且可以基于单元不为空或受损的估计后验概率来选择最终过滤。虽然这些模型通常适用于单细胞RNA-seq数据，但可能必须应用几个额外的过滤器或启发式方法才能在单核RNA-seq数据中获得稳健的过滤，就像DropletUtils的emptyDropsCellRanger函数中那样。

2.5.2 双联体检测

除了确定哪些细胞条形码对应于空滴或受损细胞之外，人们还可能希望识别那些对应于双联体或多联体的细胞条形码。当给定的液滴捕获两个（双联体）或更多（多联体）细胞时，这可能会导致这些细胞条形码在数量上出现偏态分布，例如它们代表的读数和UMI数量，以及它们显示的基因表达谱。还开发了许多工具来预测细胞条形码的双联态状态。一旦检测到，被确定为可能是双联体和多重联体的细胞可以被去除或以其他方式在后续分析中进行调整。

2.6 计数数据表示

当完成最初的原始数据处理和质量控制并继续进行后续分析时，重要的是要承认并记住基因计数矩阵充其量只是原始样本中测序分子的近似值。在原始数据处理流程的几个阶段，进行简化以生成该计数矩阵。例如，reads映射并不完美，细胞条形码的修正也存在问题。准确解析UMI尤其具有挑战性，与UMI相关的多重映射读取问题经常被忽视。

2.7 实战案例

鉴于我们已经涵盖了各种原始数据处理方法背后的概念，我们现在将注意力转向演示如何使用特定工具（在本例中为alevin-fry）来处理小型示例数据集。首先，我们需要来自单细胞实验的FASTQ格式的测序读数以及读数将被比对的参考（例如转录组）。通常，参考基因组分别包括FASTA和GTF格式的已测序物种的基因组序列和相应的基因注释。
在此示例中，我们将使用人类基因组的5号染色体及其相关基因注释作为参考，该参考是来自10x Genomics构建的人类参考基因组GRCh38 (GENCODE v32/Ensembl 98) reference的子集。相应地，我们从10x Genomics的人脑肿瘤数据集中提取映射到生成参考的读数子集。
Alevin-fry[He et al., 2022]是一种快速、准确且节省内存的单细胞和单核数据处理工具。Simpleaf是一个用rust编写的程序，它公开了一个统一且简化的接口，用于使用alevin-fry处理流程的一些最常见协议和数据类型。这里我们首先展示如何使用两个simpleaf命令处理单细胞原始数据。然后，我们描述了这些 simpleaf命令对应的完整的salmon alevin和alevin-fry命令集，以概述本节中描述的步骤发生的位置，并传达可能的不同处理选项。这些命令将从命令行运行，conda将用于安装运行此示例所需的所有软件。

2.7.1 准备

在开始之前，我们在终端中创建一个conda环境并安装所需的包。Simpleaf依赖于 alevin-fry、salmon和pyroe。 它们都在bioconda上可用，并且在安装simpleaf时会自动安装。

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

接下来，我们创建一个工作目录af_xmpl_run，并从远程主机下载并解压缩示例数据集。

# 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文件）准备就绪后，我们现在可以应用上面讨论的原始数据处理流程来生成基因计数矩阵。

2.7.2 简化的原始数据处理流程

Simpleaf旨在简化单细胞原始数据处理的alevin-fry接口。 它将整个处理流程封装为两个步骤：

• 1.simpleaf index索引所提供的参考或创建剪接参考（剪接转录本+内含子）并对其进行索引。

• 2.simpleaf quant将测序读数与索引参考进行映射，并对映射记录进行量化以生成基因计数矩阵。
  可以在此处找到使用simpleaf进行映射的更多高级用法和选项。
  运行simpleaf index时，如果提供了基因组FASTA文件（-f）和基因注释GTF文件（-g），则会生成splici引用并对其进行索引；如果仅提供转录组FASTA文件（--refseq），它将直接对其进行索引。目前，我们推荐使用 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引用；索引文件夹包含基于拼接参考构建的salmon index。
下一步，simpleaf quant，使用索引目录和映射记录FASTQ文件来生成基因计数矩阵。该命令封装了本节中讨论的所有主要步骤，包括映射、单元格条形码校正和 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

运行这些命令后，可以在simpleaf_quant/af_quant/alevin文件夹中找到生成的量化信息。在该目录中，存在三个文件：quants_mat.mtx、quants_mat_cols.txt和quants_mat_rows.txt，它们分别对应于计数矩阵、该矩阵每列的基因名称以及校正、过滤后的该矩阵的每一行细胞条形码。这些文件的尾行如下所示。这里值得注意的是alevin-fry在USA模式（未剪接、剪接和模糊模式）下运行，因此对每个基因的剪接和未剪接状态进行了量化——生成的quants_mat_cols.txt文件的行数将等于已注释基因数的三倍，对应于每个基因的剪接（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函数将计数矩阵作为AnnData对象加载到Python中。fishpondR包中函数loadFry已经实现了类似的功能。

import pyroe

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

默认将Anndata对象的x层加载为每个基因的剪接计数和模糊计数的总和。然而，最近的工作和更新的实践表明，即使在单细胞RNA-seq数据中包含内含子计数，也可能会提高灵敏度并有利于下游分析。虽然利用这些信息的最佳方法是正在进行的研究的主题，但由于alevin-fry自动量化每个样本中的剪接、未剪接和模糊读段，因此可以简单地获得包含每个基因总计数的计数矩阵，如下所示 :

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使得通过几个命令以“标准”方式处理单细胞原始数据成为可能。接下来，我们将展示如何通过调用pyroe、salmon和alevin-fry命令来生成相同的量化结果。除了教学价值之外，如果仅需要重新运行管道的一部分或者需要指定simpleaf当前未公开的某些参数，那么了解每个步骤的确切命令将很有帮助。
请注意，应提前执行准备部分中的命令。以下命令中调用的所有工具，pyroe、salmon和alevin-fry，在安装simpleaf时已经安装完毕。

2.7.3.1 建立索引

首先，我们处理基因组FASTA文件和基因注释GTF文件以获得剪接索引。以下代码块中的命令类似于上面讨论的simpleaf index命令。这包括两个步骤：

• 1.使用基因组和基因注释文件，通过调用pyroe make-splici构建剪接参考（剪接转录本 + 内含子）。
• 2.通过调用salmon index对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. 比对和量化

接下来，我们将通过调用salmon alevin来映射根据剪接索引记录的测序读数。这将生成一个名为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执行细胞条形码校正和UMI解析步骤。此过程涉及三个 alevin-fry命令：

• 1.generate-permit-list命令用于细胞条形码校正。
• 2.collate命令过滤掉无效的映射记录，更正细胞条形码并整理源自相同已更正细胞条形码的映射记录。
• 3.quant命令执行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

运行这些命令后，生成的量化信息可以在alevin_fry_quant/alevin中找到。有关映射、CB校正和UMI分解步骤的其他相关信息可分别在salmon_alevin、alevin_fry_gpl和alevin_fry_quant文件夹中找到。
在此处给出的示例中，我们演示了使用simpleaf和alevin-fry处理10x Chromium 3’ v3 数据集。Alevin-fry和simpleaf提供了许多其他选项来处理不同的单细胞方案，包括但不限于Dropseq、sci-RNA-seq3和其他10x Chromium平台 。有关不同处理阶段的可用选项的更全面的列表和描述可以在alevin-fry和simpleaf文档中找到。
当然，还有许多其他原始数据处理工具也存在类似的资源，包括 zUMIs、alevin、kallisto|bustools、 STARsolo和CellRanger。