search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
103 stars 30 forks source link

单细胞转录组 —— Cell Ranger 原始数据处理 #5160

Closed ixxmu closed 2 months ago

ixxmu commented 2 months ago

https://mp.weixin.qq.com/s/AyqdJp7abVtoxT5qdV8hag

ixxmu commented 2 months ago

单细胞转录组 —— Cell Ranger 原始数据处理 by 生信杂货铺

前言

我们前面介绍了单细胞转录组原始数据的处理步骤，以及每一步中涉及的各种问题及解决方案。下面我们将介绍几个常用的原始数据的处理流程。

参考基因组

首先是参考基因组的选择，大多数 scRNA-seq 实验都是利用人体或小鼠组织、器官组织或细胞培养物进行的。

这些基因组组装和注释的更新比较频繁。有两种常用的组装文件来源：

  • UCSC 的组装文件名为 hg19、hg38、mm10 等
  • GRC：如 GRCh37、GRCh38、GRCm38 等

UCSC 和 GRC 的主要版本在主要染色体上是一致的（如 hg38 中的 chr1 等于 GRCh38 中的 chr1），但在附加等位基因和所谓的 ALT 位点上有所不同，这些位点在次要版本（如 GRCh38.p13）之间会发生变化。

更多信息请查阅 Genome Reference Consortium[1] 和 Heng Li's blog[2]。

基因组注释文件定义了基因组（基因）的转录区，以及注释具有外显子内含子边界的确切转录本，并为新定义的特征分配类型--如蛋白质编码、非编码等。

人类和小鼠基因组注释的常用来源是 RefSeq、ENSEMBL 和 GENCODE。

其中 RefSeq 是这三个来源中最保守的一个，每个基因注释的转录本往往最少。RefSeq 转录本 ID 以 NM_ 或 NR 开头，如 NM_12345。

ENSEMBL 和 GENCODE 非常相似，可以互换使用。其中基因名称以 ENSG（人类）和 ENSMUSG（小鼠）开头；转录本分别以 ENST 和 ENSMUST 开头。

除了基因 ID，大多数基因还有一个常用的名称（gene symbol）；例如，人肌动蛋白 B 的 ID 为 ENSG00000075624，symbol 为 ACTB。

人类基因名称由 HGNC[3] 定期更新和定义。老鼠的基因名称是由一个类似的组织 MGI[4] 决定的。

目前的 ENSEMBL/GENCODE 数据库中人类基因组注释包含约 6 万个基因，其中 2 万个为蛋白质编码基因，以及 23.7 万个转录本。

大多数基因可按类型粗略地分为蛋白质编码基因、长非编码 RNA、短非编码 RNA 和假基因。

如图所示，不同注释版本之间基因也不一样

数据处理流程

scRNA-seq 原始数据处理流程包括四个主要步骤：

  • 将 cDNA 片段映射到参考文献
  • 将 reads 分配给基因
  • 将 reads 分配到细胞（根据 barcode）
  • 计算 RNA 分子的数量（无重复 UMI 数量）

处理的结果会产生一个基因与细胞的计数矩阵，用来估算每个细胞中每个基因的 RNA 分子数量。

Cell Ranger

Cell Ranger 是处理 10x Genomics Chromium scRNAseq 数据的默认工具。它使用 STAR 软件将 reads 与基因组进行剪接比对。

它根据 reads 在基因组中的比对位置，将其分为外显子、内含子和基因间三类。即如果 reads 至少有 50% 的位置与外显子相交，则分为外显子类；如果为非外显子且与内含子相交，则分为内含子类；否则为基因间（如下图所示）。

在对 reads 进行类型分配后，会调整比对质量。即对于与一个外显子位点对齐但也与一个或多个非外显子位点对齐的 reads，外显子位点会被优先考虑，reads 会被认为有把握比对到外显子位点，并被赋予最高比对质量。

通过检查外显子和内含子与转录组的兼容性，Cell Ranger 进一步将外显子和内含子定位到注释的转录本。

reads 的分类基于它们是正的还是反义的，基于它们是外显子的还是内含子的，或者它们的剪接模式是否与该基因相关的转录本注释兼容。

当输入的是细胞核数据时，会有很大一部分 reads 来自未剪接的转录本，比对到内含子上。为了计算这些内含子读数，在运行 cellranger count 和 cellranger multi 流程时，可以使用 include-introns 参数。

如果使用了该参数，则任何比对到基因上的 reads（包括上图中标注为转录组（蓝色）、外显子（浅蓝色）和内含子（红色））都会纳入 UMI 数量的统计。

如果一个 reads 只比对到一个基因，那么它就被认为是唯一比对的。只有唯一比对的 reads 才会被纳入 UMI 计数；多重比对的 reads 会被 Cell Ranger 丢弃。

软件安装

目前 Cell Ranger 不支持 conda 安装，可以通过容器（docker 或 singularity）安装。但是很多时候需要 root 权限，比较麻烦，好在下载就能用。

进入软件 下载页面[5]，可以看到下载方式。可以使用 curl 下载软件包

curl -o cellranger-8.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.1.tar.gz?Expires=1719433552&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=EFfBnQitglqkGEJroanvt-Y9kHHTyUk2NnbE~O6udqXJCYRE9SsNqe1akDBFcY2lhj0ywbgCdasdt87d4387Ra~z-S7uejBuAp4w6wq7dKkyKarczUQbMJQeL1vTizLBkUIIdcKxzNrLoY4E1hJgC1HDcFwNMRnDaDKIGJunkjurBW35grgskwAAMf6vkGq8adrO56WSXoLZA9PmPp2zx50Ih0GbKZ79CwvwRTC-f67DNAXCgHRtSYqA99jlVfqb-JAHPiqKrKcb7GlOtpmRYqLah63jlUnzKn0d18x1iC1-rI~lWVN1jQcISUN-GJvRZT4HNnxhap8TiCebyIP1gA__"

或者使用 wget

wget -O cellranger-8.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.1.tar.gz?Expires=1719433552&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=EFfBnQitglqkGEJroanvt-Y9kHHTyUk2NnbE~O6udqXJCYRE9SsNqe1akDBFcY2lhj0ywbgCdasdt87d4387Ra~z-S7uejBuAp4w6wq7dKkyKarczUQbMJQeL1vTizLBkUIIdcKxzNrLoY4E1hJgC1HDcFwNMRnDaDKIGJunkjurBW35grgskwAAMf6vkGq8adrO56WSXoLZA9PmPp2zx50Ih0GbKZ79CwvwRTC-f67DNAXCgHRtSYqA99jlVfqb-JAHPiqKrKcb7GlOtpmRYqLah63jlUnzKn0d18x1iC1-rI~lWVN1jQcISUN-GJvRZT4HNnxhap8TiCebyIP1gA__"

解压

tar -xf cellranger-8.0.1.tar.gz 

可以将软件环境添加到环境变量中，例如

export PATH=$PATH:/path/to/cellranger

如果下面的测试成功，说明安装没有问题了

cellranger testrun --id=tiny

构建索引

定量之前首先要构建索引，可以下载 Cell Ranger 构建好的人类和小鼠的索引，也可以使用比对好的参考基因组及注释文件来手动构建索引。

Cell Ranger 所构建的索引，都是使用 primary 基因组装配版本（即不包括 ALT 位点）。注释的 GTF 文件经过过滤，具体的过滤脚本可以在 RELEASE NOTES[6] 中找到。

主要保留了蛋白质编码、长链非编码 RNA、反义 RNA，以及属于 BCR/TCR（即 V/D/J）基因（旧版可能不包括后者）。所有假基因和小非编码 RNA 都被移除。

索引的版本大概有如下几个，2020-A 是最新的参考版本。

此外，Cell Ranger 还包含人类+小鼠组合参考索引，这对于涉及人类和小鼠细胞的实验非常有用。

下载索引

在下载页面中，可以找到对应构建好的索引的链接。例如，下载并解压人类 GRCh38 版的索引

wget "https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2024-A.tar.gz"
tar -xf refdata-cellranger-arc-GRCh38-2020-A-2.0.0.tar.gz

构建索引

例如，从 UCSC 中下载小鼠 mm10 版的参考基因组及注释文件

wget https://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/mm10.fa.gz
wget https://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/genes/mm10.refGene.gtf.gz

解压

gunzip mm10.fa.gz
gunzip mm10.refGene.gtf.gz

这里我们没有对基因进行过滤，直接构建参考基因组

cellranger mkref \
  --genome=mm10 \
  --fasta=mm10.fa \
  --genes=mm10.refGene.gtf \
  --output-dir mm10_index

定量

barcode

细胞条形码（cell barcode，CB）序列是连接到珠子上的合成序列，用于识别单个细胞。

独特序列库称为白名单，取决于 Chromium 库制备试剂盒的版本。白名单文件可从 Cell Ranger repository[7] 中获取。

Chromium 有三个白名单：

  • 737K-april-2014_rc.txt
  • 737K-august-2016.txt
  • 3M-february-2018.txt。

第一个列表中的 CB 长度为 14bp，另外两个为 16bp。下表提供了常用 10x 单细胞测序试剂盒的 CB 和 UMI 长度，以及相应的白名单文件：

Cell Ranger 使用以下算法根据白名单纠正假定的条形码序列

  1. 计算数据集中每个白名单上的条形码的频率
  2. 对于数据集中未出现在白名单上的条形码，找到与之汉明距离为 1 的白名单序列。对于每个这样的序列
  • 根据碱基的质量分数，计算观察到的条形码源自白名单条形码的后验概率
  • 后验概率超过 0.975 的白名单条形码替换观察到的条形码

在输出的 BAM 文件中，原始未校正的条形码在 CR 标记中（SAM 比对记录中的 TAG），校正后的条形码序列在 CB 标记中。无法分配校正条形码的 reads 将不会有 CB 标记。

UMI 计数

在进行 UMI 计数之前，需要根据 UMI 序列进行 PCR 去重。Cell Ranger 会尝试纠正 UMI 序列中的测序错误。

比对到转录组的 reads 会被放入具有相同 CB、UMI 和基因的分组中。如果两个分组具有相同的条形码和基因，但它们的 UMI 相差一个碱基（即汉明距离为 1），那么其中一个 UMI 很可能是由于测序中的置换错误造成的。在这种情况下，reads 数较低组的 UMI 会被修正为 reads 较高的 UMI。

Cell Ranger 会再次根据 CB、UMI（可能已修正）和基因对 reads 进行分组。如果两组或更多组读数的 CB 和 UMI 相同，但基因不同，则保留 reads 数最多组比对上的基因来计算 UMI，其他组则被丢弃。

在这两个过滤步骤之后，每个观察到的 CB、UMI 和基因组合都会作为 UMI 计数记录在未经过滤的计数矩阵中。

细胞过滤

未经过滤的原始矩阵包含许多实际上是空液滴的列，这些液滴中的基因表达计数并非为零，原因是存在技术噪音，例如存在来自破碎细胞的环境 RNA。不过，通常可以通过 RNA 的数量将它们与真正的细胞区分开来。

在 Cell Ranger 中，有两种过滤细胞的算法，分别为 Cell Ranger 2.2 和 Cell Ranger 3.0。

最初的算法（2.2）能识别条形码数量与每个条形码的 UMI 图中的第一个 拐点。而 3.0 引入了改进的细胞识别算法，能更好地识别低 RNA 含量细胞群，尤其是当低 RNA 含量细胞混入高 RNA 含量细胞群时。

例如，肿瘤样本中往往混有大的肿瘤细胞和较小的肿瘤浸润淋巴细胞（TIL），研究人员可能会对 TIL 群体特别感兴趣。新算法基于 EmptyDrops[8] 方法

计数矩阵

上面提到的功能都已经封装到软件中，只需要调用 count 子命令即可计算出技术矩阵。

首先，获取 10x 的测序原始数据，可以从网上下载。比如，从 SRA 数据库中下载 PRJNA666791 项目的数据

prefetch --option-file SRR_Acc_List.txt --output-directory sra

解压

ls sra/*/*.sra | xargs fastq-dump --split-files --gzip -O raw

解压后每个 sra 会生成三个文件，后缀分别带有 _1、_2 和 _3，分别修改为 R1、R2 和 I1 文件。可以先确认一下文件中序列的长度再进行修改。

其中样本索引文件（I1 或 I2）一般是不需要输入的，后面可以忽略。

这套数据中包含两个样本，每个样本 8 个文件，修改文件名时可以看作是 8 条 lane

GSM4812353_S1_L001_R1_001.fastq.gz
GSM4812353_S1_L001_R2_001.fastq.gz
...
GSM4812353_S1_L008_R1_001.fastq.gz
GSM4812353_S1_L008_R2_001.fastq.gz

基本命令为

cellranger count \
  --id run_count_GSM4812353 \
  --sample GSM4812353 \
  --fastqs input_dir \
  --transcriptome mm10_index \
  --localcores 20 \
  --output-dir GSM4812353 \
  --create-bam true

分析结果主要存储在 GSM4812353/outs 中。例如

结果解读

输出结果中包含三个文件夹

  • analysis：分析结果，包括降维、聚类和差异表达
  • filtered_feature_bc_matrix：过滤后的计数结果
  • raw_feature_bc_matrix：原始的计数结果

我们一般会选择过滤后的结果（filtered_feature_bc_matrix）用于下游分析。该文件夹下包含的三个文件

  • barcodes.tsv.gz：CB 信息
  • features.tsv.gz：基因信息，来自于小鼠基因组的注释文件
  • matrix.mtx.gz：计数矩阵，只存储了有值的行列索引及其对应的 count

其中，filtered_feature_bc_matrix 和 raw_feature_bc_matrix 分别都有一个对应的 h5 格式文件，可以直接作为下游软件的输入数据格式。

molecule_info.h5 文件是分子水平的信息，可以用于 cellranger agrr 合并多个样本

web_summary.html 网页文件有详细的样本分析报告，方便查看

multi

multi 主要用于分析 3' Cell Multiplexing 数据，其他情况下，一般推荐使用 count。

什么是 Cell Multiplexing 呢？简单来说，Cell Multiplexing 是一种技术，通过对细胞或细胞核样本进行分子标签标记，然后将这些样本混合在一起进行处理。在测序后，利用分子标签信息将混合样本中的细胞进行区分。

该技术的主要优点有：

  • 提高单次实验的样本通量
  • 增加检测细胞的数量
  • 提高实验的重复性
  • 在分析前检测并去除多胞（multiplets）

更多信息可以参考官网的介绍 What is Cell Multiplexing?[9]

运行 cellranger multi 需要一个 CSV 配置文件。例如

cellranger multi --id=sample345 --csv=/home/jdoe/sample345.csv

不知道怎么填写这个配置文件，可以下载示例文件，然后修改对应的值

cellranger multi-template --output=/path/to/FILE.csv

下面有几个官方的示例，具体参考 Cell Multiplexing with cellranger multi[10]，其中红色字体为可修改部分。

配置文件主要包括 4 个部分

  • [gene-expression]：两列，用于指定与基因表达数据分析相关的参数，如参考基因组和期望的细胞数，以及其他通用参数。
  • [feature]：两列，用于指定与分析 Feature Barcode 文库相关的参数。
  • [libraries]：三个必填列，用于指定输入 FASTQ 文件的位置。
  • [samples]：两个必填列，用于指定 Cell Multiplexing 的样本信息。

count 与 multi 的区别

count 流程适用于

  • 单个样本测序：每个样本都是独立测序，没有合并多个样本的情况。例如，一个简单的单细胞 RNA-seq 实验，单次运行只处理一个样本的文库。

  • 单一数据类型：测序数据只包含一种类型的数据，不涉及多种数据类型的整合分析。

multi 流程适用于

  • 多个样本合并测序：多个样本被混合在一个文库中进行测序。需要将多个样本的数据从一个测序数据中拆分并分别分析。

  • 多种数据类型整合分析：同时包含多种类型的数据，如基因表达数据（RNA-seq）和细胞表面蛋白（TCR/BCR）。

  • 多文库、多样本组合：复杂的实验设计，需要在一次测序中处理来自多个文库和样本的数据。例如，一个实验中有多个样本，每个样本都有多个文库（如不同的测序平台或不同的标记方法）。

参考资料
[1]

GRC: https://www.ncbi.nlm.nih.gov/grc/human

[2]

Which human reference genome to use?: https://lh3.github.io/2017/11/13/which-human-reference-genome-to-use

[3]

HGNC: https://www.genenames.org/

[4]

MGI: https://www.informatics.jax.org/mgihome/nomen/

[5]

downloads: https://www.10xgenomics.com/support/software/cell-ranger/downloads

[6]

RELEASE NOTES: https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build#header

[7]

Cell Ranger repository: https://github.com/10XGenomics/cellranger/tree/master/lib/python/cellranger/barcodes

[8]

EmptyDrops: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1662-y

[9]

[What is Cell Multiplexing?: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cellplex

[10]

Cell Multiplexing with cellranger multi: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/multi#cmoreference