search

ixxmu / mp_duty

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

单细胞基地 | Cell Ranger 助你轻松搞定上游分析 | 10X Genomics 平台下机数据比对质控流程 #4810

Closed ixxmu closed 7 months ago

ixxmu commented 7 months ago

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

ixxmu commented 7 months ago

单细胞基地 | Cell Ranger 助你轻松搞定上游分析 | 10X Genomics 平台下机数据比对质控流程 by 生信小白要知道

噔噔！填坑来了填坑来了，今天来填坑来了！！！

上次我们在单细胞10x Genomics平台测序流程整理这个笔记里面介绍了单细胞10X Genomics平台的建库流程，相信大家对基础的文库构建过程有了一个比较系统的了解。在当时的笔记中，我们也提及到了10X平台的单细胞整体流程，流程共计分为5步，我们说会主要介绍第2步（文库构建）和第4步（数据比对质控），当时的笔记整理了第2步的内容，还剩第4步的坑还没有填，今天我们就来填它!

为什么要整理比对质控的流程？

每一篇单细胞笔记的开头，我们都会做这样一个疑问，给阅读这篇笔记的大家讲清楚我们写这篇笔记的目的，一般不会写得很深刻，但却是我们作为用户比较关心的一个角度，那么本篇笔记也不例外。

大家可以回忆一下自己选择单细胞测序服务的经历，一般是向公司提交订单，公司会根据你的需求帮你完成建库测序以及测序下机数据的比对质控，把质控报告和过滤矩阵这些输出结果交付给你。关于数据的比对质控这一步，主要是出于测序下机数据的数据量大，对计算资源的需求也大，所以绝大多数的人会选择交付给公司完成，这是很常见的现象。

但是呢，这个时候我们就面临到另一个问题，那就是当我们想要去复现或者借鉴一个感兴趣的文章的时候，刚好这个文章只提供了测序的下机数据，而我们刚好也没有接触过上游比对质控分析的流程 ... 那么在这种情况下，我们这篇笔记就可以发挥作用啦！

如何了解学习比对质控的流程？

虽然我们了解到了学习比对质控流程很有必要，但是具体该怎么学呢？那么我想当我们想要去了解一个生物软件或者方法的时候，不妨去对应的官网看一看，为了方便使用，官网一般都会有系统的使用指南。10X平台也不例外，在其官网(https://www.10xgenomics.com/)上面我们也能找到对应的信息(https://www.10xgenomics.com/support)。由于我们的笔记暂时关注的是最常规的单细胞基因表达谱，所以我们选择“Single Cell Gene Expression”这个选项进入，界面如下：

接着，我们进入这个“Software Analysis”的界面，这里面提供了3种工具，我们在图中分别简单地描述了其功能。其中，“Cell Ranger”是我们本篇笔记重点关注的内容，那么就让我们进入“Cell Ranger”的界面一探究竟吧。

由于Cell Ranger是面向于10X平台的多款单细胞测序产品开发的一款工具，所以它能提供的分析功能非常多样化，导致这个界面里面信息非常多。我们今天主要只关注3'端基因表达谱数据的分析，我们就采取用到哪里介绍哪里的策略，就不一一解析了，大家感兴趣的话可以点击对应菜单栏查看相应的信息。

配置Cell Ranger分析环境

接下来呢，我们就进入本文的正题，开始为我们的分析做准备了。

01 Cell Ranger的下载和配置

首先呢，我们需要从官网下载Cell Ranger。我们可以点击下图中的“Download Center”这一栏，进入下载页面，下载最新的Cell Ranger 8.0.0的压缩包。

下载命令如下：

# 使用curl命令下载
curl -o cellranger-8.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.0.tar.gz?Expires=1712352870&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=lAp0t6wo-oNLmAhcfNGu6Nqie6OOcWcWoZWx9Lj2JxE8Gj1Do8p2rVYH~PE4f0wCemkrQKilRUDa42FRKskUugvKT59MrZDj9fkWy4lxlPlMOBFOLJG~eZOQJxwKgTzkXdYGuBr~J~k0KGTF0rCYBH0p6s5BB5Bff5Rx033~sXxr-YU8tCgHtx7YZOg67ov-NEMqh3y-S~CcyRV2RcQkv1nrcFMrLoC4SxoV8U4LjlzGqTPxFzSMZ2-xQYOs5wWw0N2TnrlhV~ofY6xHvPmSfEzhnWt3UHBTc5zjBSoP3bKpsFaFWFInyrZXUMILV6iXzeiEvp4T7hq5g9S2TfmaUQ__"

# 使用wget命令下载
wget -O cellranger-8.0.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.0.tar.gz?Expires=1712352870&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=lAp0t6wo-oNLmAhcfNGu6Nqie6OOcWcWoZWx9Lj2JxE8Gj1Do8p2rVYH~PE4f0wCemkrQKilRUDa42FRKskUugvKT59MrZDj9fkWy4lxlPlMOBFOLJG~eZOQJxwKgTzkXdYGuBr~J~k0KGTF0rCYBH0p6s5BB5Bff5Rx033~sXxr-YU8tCgHtx7YZOg67ov-NEMqh3y-S~CcyRV2RcQkv1nrcFMrLoC4SxoV8U4LjlzGqTPxFzSMZ2-xQYOs5wWw0N2TnrlhV~ofY6xHvPmSfEzhnWt3UHBTc5zjBSoP3bKpsFaFWFInyrZXUMILV6iXzeiEvp4T7hq5g9S2TfmaUQ__"

这个下载速度视大家的网速而定，如果时间过长，大家可以先忙别的事情。

由于Cell Ranger解压即用，所以我们只需要执行解压命令，并将其配置到bashrc中即可使用。

# 解压下载的压缩包
tar -zxvf cellranger-8.0.0.tar.gz

# 进入解压后的文件夹中，然后输入下面的命令，将cellranger的配置写进bashrc中
echo "export PATH=${PWD}:\$PATH" >> ~/.bashrc

# 然后调用cellranger空运行，弹出了帮助文档，说明安装成功
$ cellranger
## cellranger cellranger-8.0.0
## 
## Process 10x Genomics Gene Expression, Feature Barcode, and Immune Profiling data
## 
## Usage: cellranger <COMMAND>
## 
## Commands:
##   count           Count gene expression and/or feature barcode reads from a single sample and GEM well
##   multi           Analyze multiplexed data or combined gene expression/immune profiling/feature barcode data
##   multi-template  Output a multi config CSV template
##   vdj             Assembles single-cell VDJ receptor sequences from 10x Immune Profiling libraries
##   aggr            Aggregate data from multiple Cell Ranger runs
##   reanalyze       Re-run secondary analysis (dimensionality reduction, clustering, etc)
##   mkvdjref        Prepare a reference for use with CellRanger VDJ
##   mkfastq         Run Illumina demultiplexer on sample sheets that contain 10x-specific sample index sets
##   testrun         Execute the 'count' pipeline on a small test dataset
##   mat2csv         Convert a feature-barcode matrix to CSV format
##   mkref           Prepare a reference for use with 10x analysis software. Requires a GTF and FASTA
##   mkgtf           Filter a GTF file by attribute prior to creating a 10x reference
##   upload          Upload analysis logs to 10x Genomics support
##   sitecheck       Collect Linux system configuration information
##   help            Print this message or the help of the given subcommand(s)
## 
## Options:
##   -h, --help     Print help
##   -V, --version  Print version

02 参考基因组文件的下载

接着呢，我们下载依据小鼠的参考基因组创建的比对索引文件，为后续比对分析做好准备（之所以下载小鼠的是因为小鼠的基因组小一点，下载起来会快一些，哈哈哈）。进入下载界面，运行下载命令，如下图所示：

下载命令如下：

# 使用curl命令下载
curl -O "https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCm39-2024-A.tar.gz"

# 使用wget命令下载
wget "https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCm39-2024-A.tar.gz"

这个下载速度也是视大家的网速而定，如果时间过长，大家也可以先忙别的事情。

下载好之后解压待用，命令如下：

tar -zxvf refdata-gex-GRCm39-2024-A.tar.gz

03 下载测试数据

我们的软件安装好了，参考基因组也下载完毕了，接下来就摩拳擦掌准备实践了，但是我们手头上还缺少测试的数据，那么这个测试数据从哪来呢？不用担心，10X的官网上给我们提供了大量可供尝试的数据，我们从这里下载，操作如下图所示：

我从中挑选了一个E18时期小鼠的心脏数据，细胞数大概在10K左右，fastq文件21.6G，数据路径在这https://www.10xgenomics.com/datasets/10-k-heart-cells-from-an-e-18-mouse-v-3-chemistry-3-standard-3-0-0

同样地，也是下载之后进行解压待用，命令如下：

wget -bc https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.0/heart_10k_v3/heart_10k_v3_fastqs.tar

tar -xvf heart_10k_v3_fastqs.tar

Cell Ranger参数解析

上面的部分我们已经准备好了测试Cell Ranger所需要的文件，但是要真正地使用它，我们还需要系统地了解一下它能提供的功能。我们空运行一下弹出帮助文档，然后对它的子命令的功能进行一个简单的解析。至于它的子命令的使用细节，我们就不做过多补充，大家可以根据需要去这个网址(https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials)查看，我们在介绍过程中同时也会提供每个子命令的教程网址。

$ cellranger
## cellranger cellranger-8.0.0

## Process 10x Genomics Gene Expression, Feature Barcode, and Immune Profiling data

## Usage: cellranger <COMMAND>

## Commands:
##   count           Count gene expression and/or feature barcode reads from a single sample and GEM well
##   multi           Analyze multiplexed data or combined gene expression/immune profiling/feature barcode data
##   multi-template  Output a multi config CSV template
##   vdj             Assembles single-cell VDJ receptor sequences from 10x Immune Profiling libraries 
##   aggr            Aggregate data from multiple Cell Ranger runs
##   reanalyze       Re-run secondary analysis (dimensionality reduction, clustering, etc)
##   mkvdjref        Prepare a reference for use with CellRanger VDJ
##   mkfastq         Run Illumina demultiplexer on sample sheets that contain 10x-specific sample index sets
##   testrun         Execute the 'count' pipeline on a small test dataset
##   mat2csv         Convert a feature-barcode matrix to CSV format
##   mkref           Prepare a reference for use with 10x analysis software. Requires a GTF and FASTA
##   mkgtf           Filter a GTF file by attribute prior to creating a 10x reference
##   upload          Upload analysis logs to 10x Genomics support
##   sitecheck       Collect Linux system configuration information
##   help            Print this message or the help of the given subcommand(s)

## Options:
##   -h, --help     Print help
##   -V, --version  Print version
  • cellranger count: 这个子命令是将测序产生的FASTQ文件中的测序reads与参考基因组进行比对质控产生可用于后续分析的一系列文件，我们主要关注质控报告web_summary.html和矩阵文件filtered_feature_bc_matrix，我们下面的教程也是主要围绕这个子命令进行的，详细信息可以参考下述路径(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/running-pipelines/cr-gex-count), (https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-ct)
  • cellranger multi: 这个子命令是用于处理多样本混样数据的。这个命令处理的数据在实验操作上也有不同，在实验过程中，多个样本使用Cell Multiplexing Oligos (CMOs)进行唯一标记，使多个样本可以集中在单个GEM孔中，就可以同时测多个样本的表达数据。这个子命令的详细信息大家可以参考这个路径(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/running-pipelines/cr-3p-multi), (https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-cp)
  • cellranger multi-template: 这个子命令可以理解为对cellranger multi命令的补充，可以通过cellranger multi-template --output=/path/to/FILE.csv这个命令生成cellranger multi需要的配置文件模板，以及通过cellranger multi-template --parameters生成cellranger multi的所有可配置参数的信息。好像并没有其他什么功能了。
  • cellranger vdj: 这个子命令是用于分析单细胞5'V(D)J测序产生的数据的。该子命令的描述信息是“组装来自10X免疫表达谱文库的单细胞VDJ受体序列”，是关于T细胞和B细胞受体的免疫表达谱测序，我们不做过多涉及，详细信息可以参考下述路径(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/running-pipelines/cr-5p-vdj), (https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-vdj)
  • cellranger aggr: 这个子命令是用于整合数据的。我们可以理解为，当多样本的数据在多个GEM中建库的时候，我们就需要分别运用cellranger count对其进行分析，这样每个GEM就会产生一个输出。cellranger aggr就可以将这些输出的结果进行整合生成一个单一的矩阵。按照我的理解，这个整合和我们在下游分析的时候进行的多样本整合是有异曲同工的作用的，我们一般会选择在下游对不同样本的数据进行整合，所以这个我们也不做过多讲解，详细信息可以参考下述路径(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/running-pipelines/cr-3p-aggr), (https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-ag)
  • cellranger reanalyze: 这个子命令是对数据进行降维聚类处理的，依赖于上文cellranger count产生的矩阵h5文件，作用其实类似于seurat的常规分析。详细信息可以参考下述路径(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/running-pipelines/cr-3p-reanalyze), (https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-ra)
  • cellranger mkvdjref: 这个子命令是依据参考基因组为cellranger vdj分析创建比对索引文件的。详细信息可以参考下述路径(https://www.10xgenomics.com/support/software/cell-ranger/latest/analysis/inputs/cr-5p-references)
  • cellranger mkfastq: 这个子命令是将Illumina测序的碱基信号文件(BCL, base call files)转换成序列FASTQ文件，实际上只是封装了Illumina的bcl2fastq2的程序。详细信息可以参考下述路径(https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-fq)
  • cellranger testrun: 这个子命令是对一些小型的测试数据集执行cellranger count命令，我看官网上的作用就是拿它来测试cellranger安装是否成功的，好像确实并没有什么实际作用。
  • cellranger mat2csv: 这个子命令是将feature-barcode矩阵转换成CSV表格格式，参数很简单，cellranger mat2csv --help随查随用。
  • cellranger mkref: 这个子命令是依据参考基因组为后续分析创建比对索引文件的，需要物种的参考基因组GTF文件和FASTA文件，我们稍后也会介绍这一部分的内容，详细信息可以参考下述路径(https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-mr)
  • cellranger mkgtf: 这个子命令是用来过滤GTF文件的，过滤一些不需要的biotype，详细信息可以参考下述路径(https://www.10xgenomics.com/support/software/cell-ranger/latest/analysis/inputs/cr-3p-references#filter-with-mkgtf-f4a9b6)
  • cellranger upload: 这个子命令是用于反馈问题的，使用方法如下: cellranger upload your@email.edu genome_id.mri.tgz，留下你的邮箱和报错信息(扩展名是这样.mri.tgz的一个文件.mri.tgz)，等10X官方给你解答报错。
  • cellranger sitecheck: 这个子命令是用于获取你的Linux系统的配置信息的。
  • cellranger help: 这个子命令用于弹出帮助文档信息。

Cell Ranger分析流程

来来来，接下来本文的核心就来了！！！铺垫了这么久，我们总算进入到了数据分析部分了，这应该也是大家最关心的一部分，我们接下来就看看它到底是怎么一回事。我将从以下三个部分进行阐述：依据参考基因组创建比对索引文件；数据比对质控；结果解读。

01 依据参考基因组创建比对索引文件

首先，创建比对索引文件这一步，主要是针对除了人和小鼠之外的10X官方没有提供索引文件的物种而言的，如果你的数据是人和小鼠的数据，那就可以直接跳转到 "02 数据比对质控" 这一步了。我们接下来就以官网提供的斑马鱼为例，为大家整理一下利用参考基因组创建比对索引文件的步骤，参考路径如下(https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-mr)

1) 下载对应物种参考基因组的FASTA文件和GTF文件

斑马鱼的参考基因组文件下载命令如下：

# 下载斑马鱼参考基因组的GTF文件并解压
wget http://ftp.ensembl.org/pub/release-105/gtf/danio_rerio/Danio_rerio.GRCz11.105.gtf.gz
gzip -d Danio_rerio.GRCz11.105.gtf.gz

# 下载斑马鱼参考基因组的FASTA文件并解压
wget -bc http://ftp.ensembl.org/pub/release-105/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
gzip -d Danio_rerio.GRCz11.dna.primary_assembly.fa.gz

2) 过滤GTF文件

这一步的目的是什么呢？这是因为在GTF文件中包含有一些非PolyA转录本的序列信息，这些序列可能与具有蛋白编码功能的序列有overlap，这些overlap可能会导致一些测序的reads被映射到多个基因上，这种多重映射的reads就不会被计数，造成数据损失。通过在cellranger mkgtf这个命令运行时添加--attribute=gene_biotype:protein_coding这个过滤参数，就可以只保留gtf文件中gene_biotype为protein_coding的信息。

cellranger mkgtf \
    Danio_rerio.GRCz11.105.gtf \
    Danio_rerio.GRCz11.105.filtered.gtf \
    --attribute=gene_biotype:protein_coding

这个过程很快，基本上分分钟的事情。

3) 执行cellranger mkref构建比对索引文件

上面对gtf文件过滤完成之后，我们就可以构建索引文件了，设置索引文件名字以及参考基因组的FASTA和GTF文件路径即可。

gtf_path=/path/to/Danio_rerio.GRCz11.105.filtered.gtf
fa_path=/path/to/Danio_rerio.GRCz11.dna.primary_assembly.fa

cellranger mkref \
    --genome=Drerio_genome \    # 输出索引文件名字
    --fasta=${fa_path} \        # 参考基因组的FASTA文件路径
    --genes=${gtf_path} \       # 参考基因组的GTF文件路径
    --localcores=35 \           # 预设计算核心数，根据自己服务器的资源酌情设置
    --localmem=150              # 预设内存占用，根据自己服务器的资源酌情设置

该过程基本上要运行数小时，大家投递任务后可以先忙别的事情。

我这边资源给的比较多，我们可以看到大概半小时就跑完了

runtime3

4) 补充说明根据上面的操作，我们对斑马鱼参考基因组的GTF文件进行了过滤，保留了斑马鱼的gtf中gene_biotype为protein_coding的信息构建了斑马鱼的比对索引文件。这个索引文件就可以比对到测序数据中编码蛋白的基因的信息，其他的信息由于被过滤掉了，我们就比对不到这些信息了。那么当你想要保留其他信息时需要怎么操作呢？在官网教程(https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-mr)的"Reference build instructions for Rhesus macaque"部分有提到，你可以给--attribute这个参数不断添加信息，恒河猴的GTF过滤就如下所示：

cellranger mkgtf \
Macaca_mulatta.Mmul_10.105.gtf Macaca_mulatta.Mmul_10.105.filtered.gtf \
--attribute=gene_biotype:protein_coding \
--attribute=gene_biotype:lncRNA \
--attribute=gene_biotype:antisense \
--attribute=gene_biotype:IG_LV_gene \
--attribute=gene_biotype:IG_V_gene \
--attribute=gene_biotype:IG_V_pseudogene \
--attribute=gene_biotype:IG_D_gene \
--attribute=gene_biotype:IG_J_gene \
--attribute=gene_biotype:IG_J_pseudogene \
--attribute=gene_biotype:IG_C_gene \
--attribute=gene_biotype:IG_C_pseudogene \
--attribute=gene_biotype:TR_V_gene \
--attribute=gene_biotype:TR_V_pseudogene \
--attribute=gene_biotype:TR_D_gene \
--attribute=gene_biotype:TR_J_gene \
--attribute=gene_biotype:TR_J_pseudogene \
--attribute=gene_biotype:TR_C_gene

可以看到，这个过滤结果保留了gene_biotype为编码蛋白基因(protein_coding)，长链非编码RNA(lncRNA)，反义基因(antisense)，免疫球蛋白基因(IG_*_gene)，T细胞受体基因(TR_*_gene)等信息。

当然，如果你的gtf文件中没有gene_biotype这个信息，你也可以选择其他的筛选指标，只要能保留你需要的序列即可。但是注意，过滤后最小的GTF文件中至少要包含蛋白编码基因的外显子区域信息。

02 数据比对质控

由于中间插入了参数解析这一部分内容，这里我们捋一下比对质控需要的数据啊。首先是依据小鼠的参考基因组创建的比对索引文件"refdata-gex-GRCm39-2024-A"，这个文件我们已经下载并进行了解压，它的数据结构如下所示：

$ tree refdata-gex-GRCm39-2024-A
refdata-gex-GRCm39-2024-A
├── fasta
│   ├── genome.fa
│   └── genome.fa.fai
├── genes
│   └── genes.gtf.gz
├── reference.json
└── star
    ├── chrLength.txt
    ├── chrNameLength.txt
    ├── chrName.txt
    ├── chrStart.txt
    ├── exonGeTrInfo.tab
    ├── exonInfo.tab
    ├── geneInfo.tab
    ├── Genome
    ├── genomeParameters.txt
    ├── SA
    ├── SAindex
    ├── sjdbInfo.txt
    ├── sjdbList.fromGTF.out.tab
    ├── sjdbList.out.tab
    └── transcriptInfo.tab

3 directories, 19 files

接着是我们在10X官网上下载的小鼠心脏的测试数据"heart_10k_v3_fastqs"，我们也已经进行了下载解压，数据结构如下：

$ tree heart_10k_v3_fastqs
heart_10k_v3_fastqs
├── heart_10k_v3_S1_L001_I1_001.fastq.gz
├── heart_10k_v3_S1_L001_R1_001.fastq.gz
├── heart_10k_v3_S1_L001_R2_001.fastq.gz
├── heart_10k_v3_S1_L002_I1_001.fastq.gz
├── heart_10k_v3_S1_L002_R1_001.fastq.gz
└── heart_10k_v3_S1_L002_R2_001.fastq.gz

0 directories, 6 files

这里我们需要对测试数据中各文件的命名规则进行一个介绍，详细信息大家可以参考这个网址(https://www.10xgenomics.com/cn/support/software/cell-ranger/latest/analysis/inputs/cr-specifying-fastqs)的 "FASTQ file naming convention" 部分，这里我们直接借鉴他的说法：

文中说，要用作 cellranger 的输入，FASTQ 文件应符合下面的命名规则：[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz

[Sample Name]: 样本名，我们这个数据的样本名就是 "heart_10k_v3";

[Lane Number]: 测序泳道编号，当数据量过大的时候，就要切换到不同的泳道上进行测序，基本对后续分析无影响，不过可以选择对相同泳道的数据进行分析；

[Read Type]: 序列类型，主要分为以下四种。不过为了方便大家的理解，我还准备了V3的测序文库结构图，结合图片能够帮助大家更好地理解这些概念：

  • I1: Sample index read (optional) (样本序列标签，在文库结构图上是i7 index那一部分，长度为8bp，用于区分样本，不需要区分的时候就没有这个信息，所以是optional)
  • I2: Sample index read (optional) (这个信息在现在的数据上看起来并没有，我猜应该只是多加一个样本序列标签而已，位置应该和i7 index的位置对应，应该是个i5 index)
  • R1: Read 1 (在文库结构图上是10X Barcode和UMI两部分，总长28bp，其中10X Barcode长为16bp，UMI长为12bp)
  • R2: Read 2 (在文库上是91bp长的一个序列，这个就是我们测序要测的那一段序列)
好，那么我们的数据都准备好了，接下来，我们就真的来跑这个Cell Ranger了，从上面可以看出来，我们是只有一个样本的数据，所以跑起来相对比较简单。关于多样本的数据，我们会在后续专门更新讲解 (又挖个坑，哈哈哈)。

运行代码如下：

reference=/path/to/refdata-gex-GRCm39-2024-A
fq_path=/path/to/heart_10k_v3_fastqs

cellranger count --id=heart_10k_v3 \    # 样本名，需要根据自己样本信息更改
--transcriptome=${reference} \          # 索引文件路径，需要更改为对应物种的索引文件
--fastqs=${fq_path} \                   # 测序文件路径，需要更改为自己的数据路径
--sample=heart_10k_v3 \                 # 和样本名同名
--create-bam=true \                     # 是否创建bam文件
--localcores=48 \                       # 预设计算核心数，根据自己服务器的资源酌情设置
--localmem=150                          # 预设内存占用，根据自己服务器的资源酌情设置

该过程基本上要运行数小时，大家投递任务后可以先忙别的事情。

由于我这边投递任务给到的资源比较多，我们可以看到不到一个小时就跑完了。

03 结果解读

1) 数据结构

接下来，就让我们看看这个结果的数据结构是什么样子的吧。由于整个输出结果层级嵌套，过于冗余，因此我们就先展示个输出结果第一层级的截图吧，输出结果第一层级的结构如下图所示：

我们主要关注的是outs中的结果信息，所以就重点展示一下它的数据结构。其中，我们重点关注的是涉及数据质量的质控报告web_summary.html以及用于后续分析的矩阵文件filtered_feature_bc_matrix。矩阵文件就不用提了，我们下游分析的输入文件，关于其解读，我们在这个笔记单细胞分析 | Seurat基础流程 | 保姆级教程的开头部分有介绍，大家可以跳转阅读。所以接下来，我们将主要对质控报告中的信息进行解读。

$ tree outs/ -L 2
outs/
├── analysis
│   ├── clustering
│   ├── diffexp
│   ├── pca
│   ├── tsne
│   └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
└── web_summary.html

8 directories, 14 files

2) 质控报告解读

接下来，就是我们本次笔记的最后一部分了，我们通过质控报告来看看我们的数据质量如何。我在官网对应数据页面(https://cf.10xgenomics.com/samples/cell-exp/3.0.0/heart_10k_v3/heart_10k_v3_web_summary.html)也找到了基于Cell Ranger 3.0.0版本的质控报告信息，我们刚好拿我们的Cell Ranger 8.0.0版本的结果和这个比较一下，看看两者有什么区别。另外，10X官方也提供了质控报告结果的解读文档，大家可以在这个网址获取(https://cdn.10xgenomics.com/image/upload/v1660261286/support-documents/CG000329_TechnicalNote_InterpretingCellRangerWebSummaryFiles_RevA.pdf)

① 样本信息数据

样本信息数据都比较直观，我们这里就不做讲解了，大家可以自行了解。

② 细胞层面数据

细胞层面的数据是评估测序质量的重要指标，我们重点讲解这一信息，先把图摆上来：Barcode Rank Plot: 这个图是我们首先要讲解的内容，该图显示的是每个Barcode对应的GEM内的UMI计数的降序分布情况。就可以理解为根据Barcode的UMI计数来确定Barcode标记的GEM是否与细胞相关，是否能判定为一个有效细胞。因为含有细胞的GEM比不含细胞的GEM含有更多的转录本捕获(即UMI counts)。

我们来看10X官方文档中对于“标准样本”的描述(结合下图)：在绿色箭头前面这一段“急剧的下降”表示与细胞相关的Barcode和空胞的Barcode的良好分离。在这个骤降的中间位置之前的Barcode对应的细胞基本上就是可以参与后续分析的高质量细胞。

同时，文中也提到，一个理想的Barcode Rank Plot应该包含一个Cliff和一个knee，Cliff就是下图中绿色箭头所指的这段骤降的区域，knee就是下面蓝色箭头所指的这个拐角区域。

结合上述信息，回看我们的数据，我们的Barcode Rank Plot截取的高质量细胞数为7889，大致在骤降区域的中间位置，合理；同时也具备Cliff和knee，是一个合格的数据。

  • Estimated Number of Cells: 捕获到的高质量细胞的数量，主要指标
  • Fraction Reads in Cells: 捕获到的胞内Reads的比例
  • Mean Reads per Cell: 每个细胞中的Reads数的平均数
  • Median UMI Counts per Cell: 每个细胞中的UMI计数的中位数
  • Median Genes per Cell: 每个细胞中的基因数的中位数，主要指标
  • Total Genes Detected: 检测到的基因的总数，主要指标

③ 测序层面数据

测序层面的数据主要用于评估测序质量，信息如下：

  • Number of Reads: 测序文库中捕获到的Reads总数。
  • Number of Short Reads Skipped: 略
  • Valid Barcodes: 有效Barcodes的比例，有效Barcodes低可能表明存在测序问题，一般这个指标>75%被认为是合理范围。
  • Valid UMIs: 含有效UMI的Reads的比例，有效UMI低可能表明测序或文库质量存在问题，一般这个指标>75%被认为是合理范围。
  • Sequencing Saturation: 测序饱和度，取决于文库复杂性、测序深度和实验分析目标。较低的测序饱和度表明文库的大部分复杂性尚未通过测序捕获。
  • Q30 Bases in Barcode: Q-score≥30的细胞Barcodes的比例。
  • Q30 Bases in RNA Read: Q-score≥30的RNA reads的比例
  • Q30 Bases in UMI: Q-score≥30的UMI的比例

Q30碱基百分比低可能表明存在测序问题，例如上样浓度不理想。

关于测序的这些指标大家感兴趣可以自行检索，我对这些指标的理解也很有限，就做简单的翻译了，重点是了解各个指标的一个标准，以符合我们的应用需求。

④ 比对层面数据

比对层面的数据一定程度上决定了数据的可用性，如果比对质量太差，则数据不一定能正常使用，主要与参考基因组的注释质量和测序质量有关，如果能排除测序质量的问题，则主要确定为参考基因组的注释质量问题。

  • Reads Mapped to Genome: 能够与基因组匹配的Reads的比例，主要指标
  • Reads Mapped Confidently to Genome: 能够与基因组唯一匹配的Reads的比例，主要指标
  • Reads Mapped Confidently to Intergenic Regions: 匹配到基因间区的Reads的比例
  • Reads Mapped Confidently to Intronic Regions: 匹配到基因内含子区的Reads的比例
  • Reads Mapped Confidently to Exonic Regions: 匹配到基因外显子区的Reads的比例，主要指标
  • Reads Mapped Confidently to Transcriptome: 匹配到转录组区域的Reads的比例
  • Reads Mapped Antisense to Gene: 匹配到基因反义链区域的Reads的比例

基本上我们关注比对结果主要关注和基因组匹配的Reads的比例以及和外显子区域匹配的Reads的比例，比例过低会丢失关键基因的表达，对后续分析造成影响。当这两个指标出问题的时候，就需要回头评估参考基因组的质量以及检查测序实验流程是否出现纰漏。

至于第二页的聚类和差异分析，我觉得大家可以简单看一看，毕竟这一部分的内容是我们在seurat中还要专门处理的，大家可以针对seurat的笔记单细胞分析 | Seurat基础流程 | 保姆级教程专门学习。

好了，到这里，我们本期的笔记就整理完了，10X平台的信息量太大了，单单是Cell Ranger的基础知识我们就写了几千字。关于多样本数据的处理，我们放到下一期的笔记中再做整理吧，到这里就算把Cell Ranger初步整理完了，收工了收工了。

最后，再给自己的帐号打个广告，哈哈哈，我的账号“杨小杨的科研笔记”，上面会更新一些分析过程中踩过的坑以及一些细枝末节可能有用的技巧，大家可以关注起来，说不定什么时候就用的上了呢，哈哈哈。

那今天的分享就到这里啦！我们下期再见哟！

推荐完小杨，最后顺便也给自己推荐一下嘿嘿嘿！

如果我们的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！

蟹蟹你们的喜欢和支持！！！