Single-cell best practices

1. 现有书籍

1.1 Bioconductor OSCA 和 OSTA

OSCA和OSTA使用Bioconductor分别分析了单细胞转录组和空间转录组数据，这两本书籍是为Bioconductor 生态系统设计的，它们仅使用 Bioconductor 上提供的工具，所以不一定是最佳的分析。

1.2 当前单细胞 RNA-seq 分析的最佳实践：教程

这本书介绍了单细胞RNA测序分析的最佳实践。与 Bioconductor OSCA 相比，本书不受分析工具的影响，内容更加完整详细并且提供了该领域的最新观点。

2. 单细胞RNA测序

所有测序分析方法都有自己的优点和局限性，通过了解这些特性从而更好地了解数据中可能存在的偏差。

2.1 生命的基本要素

细胞是生命活动的基本单位，是生命的基本组成部分，生物体进行自我复制DNA中的遗传信息以保持生命周期的持续。测序是破译DNA序列的过程，得到其中基因和调控原件的位置和功能，以揭示遗传特征应用于突变与疾病。

eg：镰刀状红细胞病，血红蛋白异常，因β-肽链第6位氨基酸谷氨酸被缬氨酸所代替。

cell

2.2 测序发展

seq

2.2.1 第一代测序

弗里德里希·桑格 (Friedrich Sanger) 开发了一种使用放射性标记、快速测定DNA 测序的方法，称为“链终止法”，通常称为“桑格测序”。
Sanger：加热，扩增，加引物，DNA聚合酶、dNTP、ddNTP，随机结合ddNTP，长短不一双链DNA，电泳
放射性同位素标记的ddNTP不含羟基 3’OH，不能形成磷酸二脂键

优点：操作简单，准确率高

局限性：通量低，成本高
分析仪：ABI 370 荧光染料代替放射性分子
—

2.2.2 第二代测序

也称为焦磷酸测序、下一代测序NGS、高通量测序，合成测序。二代实现了边合成边测序，通过捕获新添加碱基所携带的特殊标记来确定DNA序列。
扩增，测序仪：Illumina、Roche 454、SOLiD

Illumina测序流程：样本准备文库制备（DNA片段化，加接头，PCR扩增），桥式PCR生成簇，测序，数据分析。

优点：通量高、试剂便宜

局限：测序仪贵、准确性不如第一代

2.2.3 第三代测序

三代测序技术最大的特点是可以实现长读长测序，而且有实时测序能力，便携可以现场工作。

测序仪：PacBio、GridION

优点：长读长、便携、速度快
局限性：错误率高、试剂贵

测序发展

2.3 NGS流程

存在多种NGS流程，测序一般步骤相同。
1、样品准备和文库制备
2、扩增和测序
3、数据输出和分析

2.4 RNA测序

中心法则，编码mRNA非编码ncRNA，可变剪切，mRNA前体通过不同位点外显子的组合，一个基因可以编码多种蛋白。
RNA测序可获得基因的表达水平即得到基因表达谱，能进一步能够检测基因亚型、基因融合、单核苷酸变异和许多其他有趣的特性。

2.5 单细胞RNA测序

RNA测序可主要分为bulk RNA-seq（细胞平均表达水平）和scRNA-seq（异质性）。
eg：肿瘤学中，可能存在导致复发的罕见耐药肿瘤细胞，即使在培养细胞上，也很难通过简单的批量 RNA 测序来识别。

2.5.2 转录本定量

转录本定量是对已测序转录本与基因序列比对的命中数进行计数的过程，最直接的方式就是统计mapping到这个转录本上的reads的个数，将reads数作为表达量，我们称这种表达量为raw count。

1个基因可能有多个转录本，同时不同转录本对应不同的外显子区域，是一对多的关系。

主要有两种方法：全长和基于标签
全长：力图捕获并均匀测序整条转录本，对每个转录本都试图获得一致的read覆盖度。
基于标签：通过捕获3’或者5‘末端，使用UMI标识符

2.5.3 单细胞测序方案

三种类型：基于微流体装置、基于孔板、基于 Fluidigm C1 微流控芯片

2.5.3.1 基于微流体装置的方案

将细胞捕获在水凝胶液滴内，液滴被设计为可以同时捕获beads和cell。

三种基于微流体装置的方法（beads不同）： Drop-seq beads 开源更喜欢具有较低 GC 含量的基因， InDrop beads 完全开源，10X Genomics beads 有利于捕获和扩增较短的基因和具有较高 GC 含量的基因。

10X Genomics：样品制备分离细胞（移液器，显微操纵术，流式细胞仪等），beads准备，形成油包水结构，细胞mRNA释放，beads上核苷酸捕获，反转录为cDNA形成cDNA文库,液滴破碎释放，测序。

优势：通量高、测序周期短、单细胞捕获率高以及适用范围广局限：存在一定doublet or empty droplet 的情况，导致部分数据不可用。
转录本检出率低，仅捕获3’而不是完整转录本。

2.5.3.2 基于板孔

将细胞物理分离到微孔板中，需要通过例如荧光激活细胞分选（FACS）进行细胞分选，其中根据特定的细胞表面marker对细胞进行分选。

基于板的测序方案：SMART-seq2、MARS-seq、QUARTZ-seq 、 SRCB-seq

优点：可以在文库制备之前收集信息，例如通过 FACS分选将细胞大小和任何使用的标记的强度等信息与孔坐标相关联。
局限：劳动密集，需要许多移液步骤，从而导致潜在的技术噪音和批次效应。

2.5.3.3 Fluidigm C1

Fluidigm C1 是一种微流控芯片，它以自动方式将细胞装载并分离到小型反应室中。

优点：允许完整长度的转录本覆盖。
局限：价格昂贵，使用的阵列仅捕获特定的细胞大小，这可能会使捕获的转录本产生偏差。

2.5.3.4 纳米孔单细胞转录组测序

ScNaUmi-seq:将纳米孔测序与barcode 和 UMI 分配相结合，将reads鉴定的barcode和UMI分配给Nanopore reads。

scCOLOR-seq ：单细胞校正长阅读测序技术，它可以对barcode和UMI进行纠错，并允许对单个细胞进行独立的cDNA纳米孔测序。

2.5.4 单细胞与单核

单细胞测序在组织解离过程中，某些细胞类型脆弱、难以捕获。单细胞测序高度依赖新鲜组织，因此很难利用组织生物库。
细胞核对机械力的抵抗力更强，并且可以安全地从冷冻组织中分离出来，而无需使用组织解离酶。

单核中RNA数量还是明显低于整个细胞测序得到的分子数。

3. 原始数据处理

processing

3.1 原始数据质量控制

原始fastq文件，通过质量控制QC工具来快速诊断数据质量，生成QC报告。

FASTA和FASTQ：

把FASTA作为存储有顺序的序列数据的文件后缀,如参考基因组序列。
FASTA文件用一个空格把头信息分为两个部分：第一部分是序列名字，它和大于号（>）紧接在一起；第二部分是注释信息，这个可以没有，就看具体需要，紧接的下一行是具体的序列内容。

>gene_00284728 length=231;type=dna
GAGAACTGATTCTGTTACCGCAGGGCATTCGGATGTGCTAAGGTAGTAATCCATTATAAGTAACATG
CGCGGAATATCCGGGAGGTCATAGTCGTAATGCATAATTATTCCCTCCCTCAGAAGGACTCCCTTGC
GAGACGCCAATACCAAAGACTTTCGTAAGCTGGAACGATTGGACGGCCCAACCGGGGGGAGTCGGCT
ATACGTCTGATTGCTACGCCTGGACTTCTCTT

FASTQ存的则是产生自测序仪的原始测序数据，文件后缀通常都是.fastq，.fq或者.fq.gz（gz压缩）。

@DJB775P1:248:D0MDGACXX:7:1202:12362:49613
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
+
JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA
@DJB775P1:248:D0MDGACXX:7:1202:12782:49716
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
+
IIIIIIIIIIIIIIIHHHHHHFFFFFFEECCCCBCECCCCCCCCCCCCCCCC

每四行成为一个独立的单元，我们称之为read。
第一行：以‘@’开头，是这一条read的名字，这个字符串是根据测序时的状态信息转换过来的，中间不会有空格，它是每一条read的唯一标识符。
第二行：测序read的序列，N代表的是测序时那些无法被识别出来的碱基。
第三行：以‘+’开头。
第四行：测序read的质量值，这个和第二行的碱基信息一样重要，它描述的是每个测序碱基的可靠程度，用ASCII码表示。

3.2 对齐和映射

映射或对齐原始数据处理的基本步骤，它是指确定每个测序片段的潜在起源位点的过程(例如，与reads相似的基因组或转录组位点集)。

3.2.1 映射类型

最常应用的三种主要映射算法：剪切比对、连续比对、轻量级映射的变体

3.2.2 针对不同参考序列的映射

主要三种类别：完整的参考基因组、带注释的转录组、增强的转录组

3.3 细胞条形码矫正

3.3.1 条形码错误类型

错误类型：
1、双峰/多重峰：单个条形码可以与两个或多个细胞相关联，并导致细胞计数不足。
2、空液滴：可以捕获没有封装细胞的液滴，其中环境 RNA 带有条形码标记并且可以进行测序；这会导致细胞计数过多。
3、序列错误：PCR 扩增或测序过程中可能会出现错误，并可能导致细胞计数不足或过多。

校正错误：
1、针对已知的潜在条形码列表进行校正，这种方法利用已知的许可列表来有效纠正许多可能已损坏的条形码。
2、基于膝盖或肘部的方法进行校正，构建条形码的累积频率图，其中条形码按照与其关联的不同reads或 UMI 的数量降序排列。频率图将包含拐点，作为正确捕获的细胞和空液滴之间可能的区分点。
3、基于强制数量的有效单元格进行过滤。频率大于或等于所选索引处的频率的所有条形码均被视为有效并被视为构成许可列表。

3.4 UMI分辨率

UMI 重复数据删除步骤旨在识别源自实验中捕获和测序的每个细胞中的每个原始 PCR 前分子的读数和 UMI 集。该过程的结果是将分子计数分配给每个细胞中的每个基因，随后在下游分析中用作该基因的原始表达估计。

3.4.1 UMI解析

将查看观察到的 UMI 及其相关映射读数的集合并尝试推断每个基因产生的观察到分子的原始数量的过程称为UMI 解析过程。

1、UMI错误：PCR期间核苷酸替换，测序期间读取错误。
2、多重映射：一个read映射到多个基因，导致了这些基因PCR前分子reads计数不确定，多基因reads/UMIs 的基因起源不明确。

3.5 计数矩阵质量控制

生成计数矩阵后，执行质量控制 (QC) 评估非常重要。

3.5.1 空滴检测

确定哪些细胞条形码对应于“高置信度”测序细胞。检查条形码的累积频率图，其中条形码按照与其关联的不同 UMI 数量的降序排列。该图通常包含一个“膝盖”，可以将其识别为正确捕获的细胞和空液滴之间的可能区分点。

缺点：并非所有累积直方图都显示明显的拐点，与条形码相关的总 UMI 计数可能并不是单独确定条形码是否与空的或损坏的细胞相关的最佳信号。

3.5.2 双峰检测

识别那些对应于双联体或多联体的细胞条形码。

3.6 计数数据表示

3.8 例子

演示如何alevin-fry使用特定工具来处理小型示例数据集。

3.8.1 准备

环境准备：在终端创建一个conda环境，并且安装Simpleaf包。

conda create -n af  
conda env list  
source activate af  
conda install simpleaf=0.15.0

数据准备：创建一个工作目录af_xmpl_run，下载并解压缩示例数据集。

mkdir af_xmpl_run 
cd af_xmpl_run

wget -qO- https://umd.box.com/shared/static/lx2xownlrhz3us8496tyu9c4dgade814.gz | tar xzf - --strip-components=1 -C .
wget -qO- https://raw.githubusercontent.com/10XGenomics/cellranger/master/lib/python/cellranger/barcodes/3M-february-2018.txt.gz | gunzip - > 3M-february-2018.txt


"toy_read_fastq", "toy_human_ref", 3M-february-2018.txt

3.8.2 简化的原始数据处理管道

1、simpleaf index ：采用 FASTA 格式的参考基因组和 GTF 格式的基因注释，提取剪接转录本和内含子区域的序列后，制作spliced+intronic ( splici ) 参考索引。
2、simpleaf quant ：可直接从提供的索引和测序读数生成基因×条形码计数矩阵。在检测和纠正细胞条形码后，它会映射读取并解析相关的 UMI。

提供了基因组FASTA文件（）和基因注释GTF文件（），则会生成splici引用并对其进行索引。

fastq_dir="toy_read_fastq"
ref_dir="toy_human_ref"

mkdir alevin_fry_home && export ALEVIN_FRY_HOME='alevin_fry_home'

##simpleaf需要环境变量ALEVIN_FRY_HOME来存储配置和数据
simpleaf set-paths    ##simpleaf set-paths命令找到所需工具的路径，并在ALEVIN_FRY_HOME文件夹中编写一个配置JSON文件。

##Usage: simpleaf index -o out_dir [-f genome_fasta -g gene_annotation_GTF|--refseq transcriptome_fasta] -r read_length -t number_of_threads  
 ##-r read_length是测序仪为生成生物reads而进行的测序周期数。  
 
 
simpleaf index \
-o simpleaf_index \
-f toy_human_ref/fasta/genome.fa \
-g toy_human_ref/genes/genes.gtf \
-r 90 \
-t 8

使用索引目录和映射记录 FASTQ 文件来生成基因计数矩阵。该命令封装了本节中讨论的所有主要步骤，包括映射、单元格条形码校正和 UMI 解析。

 ##reads1和reads2变量是通过从toy_read_fastq目录中查找模式为“_R1_”和“_R2_”的文件名来定义的。

reads1_pat="_R1_"
reads2_pat="_R2_"   


##指定的目录（${fastq_dir}）中查找文件名符合指定模式（*$reads1_pat*）的文件，并将它们按照文件名排序后以逗号分隔的形式输出。  

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,)"

## 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/ref/t2g_3col.tsv \
-o simpleaf_quant

在此目录中，存在三个文件：quants_mat.mtx、quants_mat_cols.txt和quants_mat_rows.txt，分别对应于计数矩阵、该矩阵每列的基因名称以及该矩阵每行的校正、过滤的细胞条形码。

tail -3 simpleaf_quant/af_quant/alevin/quants_mat.mtx  

tail -3 simpleaf_quant/af_quant/alevin/quants_mat_cols.txt  

tail -3 simpleaf_quant/af_quant/alevin/quants_mat_rows.txt

使用load_fry中的函数将计数矩阵作为对象加载到 Python 中pyroe。


import pyroe

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

获得包含每个基因总计数的计数矩阵。


import pyroe

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