scRNA-seq表达矩阵的构建

希望大家能有所收获！

目录

⊙引言—关于课程

⊙scRNA-seq简介

⊙scRNA-seq原始数据的质控

⊙scRNA-seq数据处理—文件格式小结

⊙scRNA-seq数据处理—demultiplexing

⊙scRNA-seq数据的处理—STAR

⊙scRNA-seq数据处理—Kallisto

正文

表达矩阵的构建

scRNA-seq数据的许多分析以表达矩阵为起点。按照惯例，表达矩阵的每一行代表一个基因，每列代表一个细胞（尽管一些作者使用转置矩阵）。每个条目代表给定细胞中特定基因的表达水平。基因表达的测量单位取决于protocol和使用的一般方式。

4.1 reads QC

scRNA-seq实验的结果是大量的cDNA reads的集合。第一步是确保reads具有高质量。可以使用一些常规工具（如FastQC或Kraken）执行质量控制。

假设我们的reads位于experiment.bam中，我们将运行FastQC

$<path_to_fastQC>/fastQC experiment.bam

下面是FastQC输出125 bp reads数据集的示例。该图显示了一个技术错误，导致几个碱基无法在read中心无法正确读取。但是，由于读取的其余部分质量很高，因此该错误很可能对比对效率的影响可以忽略不计。

此外，使用Integrative Genomics Browser（IGV）或SeqMonk对数据进行可视化通常很有帮助。

4.2 reads比对

在从读数中trim掉低质量碱基后，可以将剩余序列比对到参考基因组。同样，不需要特殊的方法，因此我们可以使用STAR或TopHat比对工具。对于有良好注释基因组的生物（例如小鼠，人）的全转录数据集，伪比对方法（例如Kallisto，Salmon）可能优于常规比对。对于具有数十或数十万个reads的基于drop-seq的数据集，伪比对可能更有效，因为它们的运行时间可能比传统的比对工具低几个数量级。

使用STAR来比对reads.bam的例子：

$<path_to_STAR>/STAR --runThreadN 1 --runMode alignReads
--readFilesIn reads1.fq.gz reads2.fq.gz --readFilesCommand zcat --genomeDir <path>
--parametersFiles FileOfMoreParameters.txt --outFileNamePrefix <outpath>/output

注意，如果使用spike-ins，参考序列应该与所述的DNA序列spike-ins在比对之前的分子。

请注意，使用UMI时，应从读取序列中删除其barcode。通常的做法是将barcode添加到read名称上。

一旦将每个细胞的reads比对到参考基因组，我们需要确保每个细胞的足够数量的reads可以比对到参考基因组。根据我们的经验，小鼠或人类细胞的可比对的reads比例为60-70％。但是，此结果可能会因protocol，reads长度和reads比对的设置而异。一般来说，我们希望所有细胞都具有相似的比对的reads部分，因此应检查并可能删除任何异常值。低比例的可比对reads通常意味着污染。

如何使用Salmon量化表达的一个例子是

$<path_to_Salmon>/salmon quant -i salmon_transcript_index -1 reads1.fq.gz -2 reads2.fq.gz -p #threads -l A -g genome.gtf --seqBias --gcBias --posBias

注意 Salmon根据我们的经验产生估计reads数和估计的每百万转录数（tpm），后者用于校正scRNASeq的长基因的表达，因此我们建议使用reads数。

4.3 比对例子

下面的直方图显示了scRNA-seq实验中映射到每个细胞的读数总数。每个条形代表一个细胞，它们按照每个细胞的总读数按升序排序。三个红色箭头表示在覆盖范围方面为异常值的细胞，应将其从进一步分析中删除。两个黄色箭头指向具有令人惊讶的大量未映射读数的细胞。在该实施例中，我们在比对QC步骤期间保持细胞，但是由于核糖体RNA读取的高比例，它们随后在细胞QC期间被去除

4.4 对比QC

在将原始测序映射到基因组后，我们需要评估映射的质量。有许多方法可以测量绘图质量，包括：映射到rRNA / tRNA的读数量，独特映射读数的比例，跨剪接点的读数映射，沿着转录本的读取深度。为大量RNA-seq开发的方法，如RSeQC，适用于单细胞数据：

python <RSeQCpath>/geneBody_coverage.py -i input.bam -r genome.bed -o output.txt
python <RSeQCpath>/bam_stat.py -i input.bam -r genome.bed -o output.txt
python <RSeQCpath>/split_bam.py -i input.bam -r rRNAmask.bed -o output.txt

然而，预期结果将取决于实验方案，例如许多scRNA-seq方法使用poly-A选择以避免对rRNA进行测序，这导致跨基因的读取覆盖中的3'偏差（也称为基因体覆盖）。下图显示了这个3'偏差以及三个单元格，它们是异常值并从数据集中删除：

4.5 reads量化

下一步是量化每个细胞的每个基因的表达水平。对于mRNA数据，我们可以使用为大量RNA-seq数据开发的工具之一，例如HT-seq或FeatureCounts

# include multimapping
<featureCounts_path>/featureCounts -O -M -Q 30 -p -a genome.gtf -o outputfile input.bam
# exclude multimapping
<featureCounts_path>/featureCounts -Q 30 -p -a genome.gtf -o outputfile input.bam

独特的分子标识符（UMI）使得计算分子的绝对数量成为可能，并且它们已被证明在scRNA-seq中很受欢迎。我们将在下一章讨论如何处理UMI。

4.6 独特的分子标识符（UMI）

感谢EMBL Monterotondo的 Andreas Buness 在本节的合作。

4.6.1 简介

独特的分子标记是在反转录过程中添加到转录本中的短（4-10bp）随机条形码。它们使测序读数能够分配到单个转录物分子，从而从scRNASeq数据中去除扩增噪声和偏差。

当对包含UMI的数据进行排序时，使用技术仅对包含UMI（通常为3'末端）的转录本的末端进行特异性测序。

4.6.2 比对barcode

由于独特的条形码数量（4N,其中N是UMI的长度）远小于每个细胞的分子总数（~106），每个条形码通常被分配给多个转录本。因此，为了识别独特的分子，必须使用条形码和映射位置（转录物）。第一步是映射UMI读数，我们建议使用STAR，因为它很快并输出高质量的BAM比对。此外，映射位置可用于例如。识别记录不良的3'UTR成绩单。

UMI测序通常由配对末端读数组成，其中每对读取一个读取细胞和UMI条形码，而另一个读取包含来自转录物的外显子序列（图4.5）。注意，建议修剪和/或过滤以去除含有poly-A序列的读段，以避免由于这些读取映射到具有内部poly-A / poly-T序列的基因/转录物而导致的错误。

处理UMI实验中的读取后，通常使用以下约定：

1. UMI被添加到另一个配对读取的读取名称中。
2. 读取按单元条形码分类到单独的文件中

• 对于极大的浅数据集，可以将单元条形码添加到读取名称中以减少文件数量。

4.6.3 比对barcode

理论上，每个独特的UMI转录物对应代表源自单个RNA分子的所有读数。但是，实际情况往往并非如此，最常见的原因是：

1. 不同的UMI不一定意味着不同的分子

• 由于PCR或测序错误，碱基对取代事件可导致新的UMI序列。较长的UMI会产生更多出错的机会，并且根据单元条码的估计，我们预计7-10％的10bp UMI至少会包含一个错误。如果没有纠正，这种类型的错误将导致高估转录本的数量。

1. 不同的转录物不一定意味着不同的分子

• 映射错误和/或多映射读取可能导致某些UMI被分配给错误的基因/转录本。这种类型的错误也会导致高估转录本的数量。

1. 相同的UMI不一定意味着相同的分子

• UMI频率和短UMI的偏差可导致相同的UMI附着于来自相同基因的不同mRNA分子。因此，可能低估了转录本的数量。

4.6.4 纠正误差

如何最好地解释UMI中的错误仍然是一个活跃的研究领域。我们知道解决上述问题的最佳方法是：

1. UMI工具的定向邻接方法实现了一个过程，该过程考虑了不匹配的数量和类似UMI的相对频率，以识别可能的PCR /排序错误。
2. 目前是一个未决问题。通过删除具有少量读取的UMI来支持它们与特定转录本的关联，或者通过移除所有多映射读取，可以减轻该问题。
3. 由Grun，Kester和van Oudenaarden（2014）提出的简单饱和度（又名“碰撞概率”）校正来估计分子的真实数量M：
4. 其中N =唯一UMI条形码的总数，n =观察到的条形码数。

这种方法的一个重要警告是它假设所有UMI都是同等频繁的。在大多数情况下，这是不正确的，因为通常存在与GC内容相关的偏差。

确定如何最好地处理和使用UMI目前是生物信息学界的一个活跃的研究领域。我们知道最近开发的几种方法，包括：

• UMI工具
• PoissonUMIs
• zUMIs
• dropEst

4.6.5 下游分析

当前的UMI平台（DropSeq，InDrop，ICell8）具有低且高度可变的捕获效率，如下图所示。

这种可变性会引入强烈的偏差，需要在下游分析中加以考虑。最近的分析通常基于细胞类型或生物途径将细胞/基因聚集在一起以增加功率。对这些数据的稳健统计分析仍然是一个开放的研究问题，还有待确定如何最好地调整偏差。

练习1我们为您提供了由三个不同个体产生的诱导多能干细胞的UMI计数和读数（Tung et al.2017）（有关此数据集的详细信息，请参阅：第7.1章）。

umi_counts <- read.table("tung/molecules.txt", sep = "\t")
read_counts <- read.table("tung/reads.txt", sep = "\t")

使用此数据：

1. 绘制捕获效率的可变性
2. 确定扩增率：每个UMI的平均读数。