【文献解读】Nature Communications：高通量且纠错的Nanopore单细胞转录组测序（方法解读）

                简介                 

标题：High throughput error corrected Nanopore single cell transcriptome sequencing 高通量且纠错的Nanopore单细胞转录组测序   

主题：方法学解读  

杂志：Nature Communication

影响因子：12.12

发表日期：2020年8月12日

解读：章小鱼zxy

编辑：很跩的土豆

导读：上周我们分享了Nature Communication杂志上一篇关于“高通量且纠错的Nanopore单细胞转录组测序”的文章，作者通过将Oxford Nanopore高通量测序与准确细胞条形码（cell barcode，cellBC）和唯一分子标识符（UMI）相结合，即ScNaUmi-seq方法，能够做到高精度的全长序列测序，并以胎鼠脑组织为例证实了ScNaUmi-seq能在单细胞水平上定义mRNA剪接和编辑。本次我们将详细阐述其具体过程，以期为关注三代测序的同道们带来一些帮助和启发。

具体步骤大致如下：分离小鼠脑组织并制备单细胞cDNA文库——Nanopore reads图形化——将细胞barcode和UMI分配给Nanopore reads——定义每个UMI的cDNA共有序列——分配Gencode转录isoform——定量单细胞基因表达并确定主要细胞类型、鉴别新型转录亚型、基因或转录本isoform表达的相关性——Gria2数据分析。

                正文                 

1.    分离小鼠脑组织并制备单细胞cDNA文库

将来自E18 C57BL/6小鼠的海马、皮层和室管膜区的组织与2mg/ml Papain在无钙Hibernate E培养基中37°C孵育20分钟，Hibernate-E/B27/GlutaMAX（HEB）培养基冲洗后，充分研磨并用40um FlowMi过滤器过滤，计数细胞量并计算存活率。

使用10X Genomics Chromium Single Cell 3’Library,Gel Bead &Multiplex kit和Chip kit（v2）将E18小鼠脑单细胞悬液转化为带有barcode的scRNA-seq库。PCR扩增并用，扩增条件为：变性，95°C 3分钟；循环（一般8个循环，对于缺少末端polyA/polyT的cDNA 5个循环）, 98°C 30s—64°C 30s—72°C 5min；最终延长，72°C 10分钟，引物浓度为1uM。

2.    Nanopore 序列比对

用minimap2 v2.17将纳米孔测序数据与Mus musculus Genome（mm10）进行比对 （命令：“ minimap2 -ax splice -uf -MD -sam-hit-only -junc-bed”）。使用minimap2中paftools.js从Gencode vM18 GTF生成bed文件。对于与已知基因匹配的reads，使用ScNaUmi-seq配套Java工具包Sicelore AddGeneNameTag方法将基因名称添加到相应的SAM记录（SAM tag：GE）中。使用Sicelore AddBamReadSequenceTag方法对SAM进行Nanopore 序列（SAM tag：US）和序列质量（SAM tag：UQ）进行注释并进行cellBC和UMI分配。

3.    将cellBC和UMI分配给Nanopore reads

利用Java执行以下分析步骤

（1） 解析Illumina数据：对于每个基因，从Illumina数据中提取了cellBC，并确定每种基因/cellBC组合的UMI；将基因组区域（窗口大小为500 nt）与cellBCs和UMIs关联起来，以解释与注释基因外部对应的序列。将解析后的Illumina数据作为Java对象存储在嵌套的Hash表中。

（2）检索polyA尾：在距核苷酸5'或3'末端分别读取100个核苷酸内搜索polyT和polyA序列（默认为85％A或T，≥20nt）。无poly（A或T）的reads和在两端都有poly（A或T）的reads不会进一步分析。

（3）检索3’端adaptor序列：为定义cellBC的位置，使用Needleman Wunsch比对检索read的末端以及上一步中鉴定的poly（A / T）序列之间的衔接子序列。使用最佳匹配（最小不匹配）的adapter对应的位置。

（4）检索序列内adapter和polyA：为标记与文库制备过程中产生的嵌合cDNA对应的reads，我们在富含A或T的序列附近检索内部嵌合序列，并将这些reads标记为输出文件中的嵌合序列。

（5）检索cellBC：每个碱基用二进制编码（A: 00, G:01, T:10, C:11），这样能够将整个16个核苷酸的cellBC编码为一个整数。提取每个基因组比对后的 Nanopore reads中adapter下游的序列（16nt）以及前后位置的barcode序列。比较（i）如果Nanopore reads与Illumina数据的已知基因匹配，则通过Illumina sequencing确定cellBC；（ii）如果Nanopore reads与未注释的基因组区域匹配，或者在Illumina数据中未找到基因名称，则使用Illumina cellBC 序列。

（6）检索UMI：为了识别10个核苷酸的UMI，在cellBC序列末端之后提取14个核苷酸。使用与cellBC分配相似的策略，在Illumina数据中搜索相同的cellBC或者在上一步中确定的基因组区域或基因的匹配UMI，并可根据搜索的复杂程度动态调整最大允许编辑距离。

（7）评估cellBC和UMI分配的准确性：为评估cellBC和UMI分配的准确性，首先扫描每个Nanopore SAM文件，以找到匹配的Illumina cellBC或UMI。然后，重复扫描，将从Nanopore读取的每个cellBC或UMI序列替换为一个随机序列，并对该随机序列进行突变，以允许与先前扫描中用于相同SAM文件的突变数（编辑距离）相同，在Illumina数据集中搜索匹配的cellBC或UMI。使用10x Genomics数据实现最大的细胞barcode和UMI分配效率。

（8）   使用10 x Genomics数据可实现最大的cellBC和UMI分配率

4.    定义每个UMI的cDNA一致性序列

筛出潜在的嵌合reads和低质量基因组比对的reads，每个细胞和基因的SAM文件按UMI分组。使用Sicelore ComputeConsensus提取TSO末端（SAM tag，TE）和polyA起始（SAM tag，PE）之间的cDNA序列以计算一致性序列。根据UMI可获取reads的数量，将以下序列指定为UMI一致性序列：（i）仅一个read时，将读段的cDNA序列分配给UMI；（ii）两个reads时，将mapping效果最佳read的cDNA片段定义为一致性序列；（iii）超过2个reads，在poa多重比对后，使用UMI所有cDNA序列，用racon对UMI的所有cDNA序列产生一个一致性序列。

5.    分配isoforms

使用minimap2 v2.17拼接比对模式将所有UMI的cDNA一致性序列与Mus musculus Genome（mm10）进行比对。分析与已知基因匹配的SAM文件，以找到与Gencode vM18转录本isoforms匹配的外显子（相同的外显子组成）。要将UMI分配给Gencode转录本，需要UMI和Gencode转录本外显子结构完全匹配。我们授权在外显子边界添加或缺少序列的两碱基裕度，以允许外显子交界处的插入缺失和minimap2的不精确mapping。遵循此策略（Sicelore IsoformMatrix），将63.6％的UMI分配给已知的Gencode转录本isoform，并生成Nanopore / Illumina gene的基因水平（中位数UMI /细胞= 6047）和isoform水平（中位数UMI /细胞= 3795）计数矩阵。

6.    定量单细胞基因表达并确定主要细胞类型

使用R / Bioconductor（版本3.5.2）和Seurat R包（版本3.1.4）处理由Cell Ranger生成的原始基因表达矩阵。共有190个细胞和951个细胞被检测到，两次重复均使用默认的Cell Ranger cutoff。去除超过95％ dropouts的细胞。从剩下的186和935个细胞（以下称为1121个细胞数据集）中，将基因表达矩阵的细胞水平标准化为10.000，并进行对数标准化。根据方差稳定化转化方法选择了前2000个高度可变的基因，并将其用于主成分分析。由于两个重复序列的测序深度不同，因此使用Seurat CCA方法对数据进行了整合。使用前11个canonical correlations来对1121个细胞进行聚类和t-SNE可视化。使用标准标记基因将t-SNE图中的簇分配给已知细胞类型。使用Seurat多模式功能，我们整合了Illumina和Nanopore基因水平和isoform水平的数据集，可以直接比较单个细胞中基因和isoform的表达。

7.    鉴别新型转录亚型

使用Sicelore CollapseModel方法，将1121个细胞数据集的UMI用于鉴定新的转录本isoforms。Gencode不支持的但具有相同外显子结构的cDNA一致性序列首先按基因分组，并使用具有相同外显子结构的序列定义潜在的新型转录本isoforms。放弃少于五个UMI支持的新型isoforms。具有相同外显子布局但在SNV或5'或3'末端不同的同工型被视为相同的isoforms。按照Tardaguila等人的建议对新型isoforms进行分类：（i）已知剪接连接的组合，仅由在Gencode转录本中发现的外显子-外显子junctions组成；（ii）已知剪接位点，单个供体和受体位点的组合是已知的，但所得的剪接点是新的；（iii）至少一个供体或受体位点在Gencode转录本中找不到。

接下来，按照以下要求过滤了新颖的转录本亚型：（i）Gencode中描述的所有外显子-外显子junctions或在E18皮质/中脑Illumina短读数据集（GSE69711）中得到确认；（ii）位于由CAGE识别的已知转录起始位点50个核苷酸内的5'末端（FANTOM5 mm9参考UCSC移至mm10，https：//fantom.gsc.riken.jp/5/datafiles/latest/extra/CAGE_peaks/）; （iii）聚腺苷酸化位点50个核苷酸内的3'端（GENCODE vM24 PolyA功能注释，https：//www.gencodegenes.org/mouse/ ）。

8.    基因或转录isoform表达的相关性

为了分析簇之间的表达相关性，将Illumina基因水平以及Nanopore基因水平和isoform水平数据下采样（R包DropletUtils downsampleMatrix方法）到1121细胞Nanopore转录本isoform水平数据集的中位UMI /细胞（3795 UMI /细胞）。然后为每个群集的每个重复单元分组，并使用簇中每个基因或isoform的平均表达来生成皮尔逊相关矩阵（R cor函数）。

以上6、7、8是功能分析的不同方面。

9.    Gria2数据分析

从951和190细胞数据集中提取了对应于Gria2（mm10：chr3：80,682,936-80,804,791）的9593个序列（2105 UMI）。计算每个分子的一致性序列，并重新mapping到mm10基因组，以使用Sicelore SNPMatrix进行SNP calling。

若要获取更为详尽的信息，请参考原文及补充材料。

                参考                 

[1] Lebrigand K，Magnone V，Barbry P，et.al. High throughput error corrected Nanopore single cell transcriptome sequencing. Nat Commun.2020.11:4025. doi: 10.1038/s41467-020-17800-6

                索引                 

【往期文献】

GUT: 新冠病毒病毒活性与COVID-19患者肠道菌群的关系

Microbiome: 化学法去除宿主DNA以提高唾液样本宏基因组测序质量（唾液样本）

                后记                 

随着测序技术的不断发展，科学研究进入了数据井喷的时代。然而，测序样本的处理流程、测序数据的分析流程甚至是数据分析过程中的数据库搭建问题，都给测序技术的普及化设置了壁垒，严重阻碍了该项技术向广大科研工作者中推广。此外，基于长读长的三代测序技术的发展更是引入了一套完全有别于二代测序数据处理的分析流程，为了让更多学者认识三代测序、在科学研究中用好三代测序，本公众号应运而生。期待与您一起学习、成长。

^_^ 边学习，边分享，每天进步一点点 ^_^