跟着Theis Lab学单细胞: alevin-fry实战

· 6 min read · Tutorials

跟着Theis Lab学单细胞: alevin-fry实战

跟着Theis Lab学单细胞: alevin-fry实战

单细胞最佳实践

原始数据处理实战

本章导读 SECTION PREVIEW

1. 了解alevin-fry这个轻量级单细胞分析工具

2. 手把手教你如何配置环境和准备数据

3. 掌握两种运行方式：简单版（simpleaf）和完整版（逐步运行）

4. 学习如何读取和解释结果数据

本节对应原书 3.8 A real-world example

01

说明

在之前的文章中，我们详细介绍了单细胞数据分析的上游处理流程。在本章中，我们将通过一个实际案例，展示如何使用alevin-fry处理一个小规模数据集。

由于教程需要，本例中仅使用人类染色体 5 及其相关注释作为参考基因组（拆分自10x Genomics人类参考基因组）。示例数据集是10x Genomics的 人类脑瘤数据集 中来源于chr 5的部分。

（本教程中用到的所有数据，代码中都包含了下载链接，本期内容包中也含有所有的代码文件和数据文件！获取方式见文末）

本例中使用到的工具如下：

• alevin‐fry
一个高效、准确且内存占用较低的单细胞与单核数据处理工具，可自动对剪接型（spliced）、未剪接型（unspliced，即内含子区域）以及模糊（ambiguous）读段进行定量。
• simpleaf
一个用 Rust 语言开发的工具，用于简化alevin-fry的处理流程。simpleaf将整个过程封装为两个步骤：

1. 构建参考索引：将基因组和基因注释文件组合生成一个“splici”参考基因组（即将剪接和内含子信息合并的参考基因组）。

2. 定量分析：将测序读段比对到这个参考上，并最终生成基因表达计数矩阵。

02

准备工作

环境配置

# 创建并激活conda环境
conda create -n af -y -c bioconda simpleaf pyroe
conda activate af

⚠️ 注意：如果你使用的是M芯片的Mac，则需要按照下面这个方式配置成普通64位环境（仅在当前环境生效，不影响其他conda环境）：

CONDA_SUBDIR=osx-64 conda create -n af -y -c bioconda simpleaf pyroe
conda activate af
conda env config vars set CONDA_SUBDIR=osx-64

数据准备

# 创建并进入工作目录
mkdir af_xmpl_run && cd af_xmpl_run

# 下载示例数据并解压
wget -qO- https://umd.box.com/shared/static/lx2xownlrhz3us8496tyu9c4dgade814.gz | tar xzf - --strip-components=1 -C .
# 定义fastq目录和参考基因组目录
fastq_dir="toy_read_fastq"
ref_dir="toy_human_ref"

# 下载barcode白名单并解压（原文链接有误，此处已更正）
wget -qO- https://raw.githubusercontent.com/f0t1h/3M-february-2018/master/3M-february-2018.txt.gz | gunzip - > 3M-february-2018.txt

03

开始实战

方式一：使用simpleaf（简单版）

Step 1. 构建参考基因组

# 配置环境变量
export ALEVIN_FRY_HOME=$(pwd)/alevin_fry_home
mkdir -p $ALEVIN_FRY_HOME
simpleaf set-paths

simpleaf set-paths会自动检测所需软件是否都安装好了。你可能会看到Could not find macs3 executable 这条信息，但是因为本教学中不涉及ChIP-seq，因此不影响后续使用。有需要的读者可以自行安装conda install macs3。

# 建立索引
simpleaf index \
-o simpleaf_index \
-f $ref_dir/fasta/genome.fa \
-g $ref_dir/genes/genes.gtf \
-r 90 \
-t 8 \
--no-piscem

如果出错，记得检查ref_dir是否已定义（应在准备工作部分就已设置好）。此外，目前默认的piscem有兼容性问题，所以使用时要加上--no-piscem这个选项，用salmon构建索引。实际使用的时候，可以增加-t 的线程数，以加快速度。额外需要注意的是，-r 需要是实际的r2读段长度（测序仪的循环数），这个长度根据测序的协议版本可能有所不同。 一般的10x数据，v3是91bp，v4是90bp，具体可以查阅此处，也可以自己检视fastq文件确定。

Step 2. 定量分析

# 获取reads文件路径
reads1=$(ls -1 "${fastq_dir}"/*_R1_* | sort | tr '\n' ',' | sed 's/,$//')
reads2=$(ls -1 "${fastq_dir}"/*_R2_* | sort | tr '\n' ',' | sed 's/,$//')

# 运行定量
simpleaf quant \
-c 10xv3 -t 8 \
-1 $reads1 -2 $reads2 \
-i simpleaf_index/index \
-u -r cr-like \
-m simpleaf_index/index/t2g_3col.tsv \
-o simpleaf_quant

成功运行以上命令后，定量文件可以在simpleaf_quant/af_quant/alevin中找到。这个文件夹有三个文件：

quants_mat.mtx: 计数矩阵，类似cellranger的matrix.mtx
quants_mat_cols.txt: 基因名，类似cellranger的features.tsv
quants_mat_rows.txt: 细胞条码，类似cellranger的barcodes.tsv

方式二：逐步运行alevin（完整版）

Step 1. 构建参考基因组

⚠️开始前记得检查ref_dir和fastq_dir这几个变量有没有设置好哦～

#构建splici参考组
pyroe make-splici \
${ref_dir}/fasta/genome.fa \
${ref_dir}/genes/genes.gtf \
90 \
splici_rl90_ref

# 构建salmon索引
salmon index \
-t $(ls splici_rl90_ref/*\.fa) \
-i salmon_index \
-p 8

Step 2. 序列比对

# 获取reads文件路径
reads1="$(ls -1 "${fastq_dir}"/*_R1_* | sort | tr '\n' ' ' | sed 's/ $//')"
reads2="$(ls -1 "${fastq_dir}"/*_R2_* | sort | tr '\n' ' ' | sed 's/ $//')"

# 进行比对
salmon alevin \
-i salmon_index \
-l ISR \
-1 ${reads1} \
-2 ${reads2} \
-p 8 \
-o salmon_alevin \
--chromiumV3 \
--sketch

Step 3. 细胞条码解析、UMI定量

# 生成barcode白名单
alevin-fry generate-permit-list \
-u 3M-february-2018.txt \
-d fw \
-i salmon_alevin \
-o alevin_fry_gpl

# 整理比对结果
alevin-fry collate \
-i alevin_fry_gpl \
-r salmon_alevin \
-t 8

# UMI处理和定量
alevin-fry quant \
-r cr-like \
-m $(ls splici_rl90_ref/*3col.tsv) \
-i alevin_fry_gpl \
-o alevin_fry_quant \
-t 8

完成定量后，结果储存在alevin_fry_quant/alevin中，包含跟方法一同样的quants_mat.mtx、quants_mat_cols.txt和quants_mat_rows.txt三个文件。

04

输出解读

输出文件

虽然输出的三个文件乍一看跟cellranger非常类似，但实际内容有不小的区别。

# 计数矩阵
# 每一行三个数字代表 行-列-非零计数
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat.mtx
138 58 1
139 9 1
139 37 1

# 基因名
# 格式为基因名-剪切状态（A-ambiguous，S-spliced，U-unspliced）
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat_cols.txt
ENSG00000120705-A
ENSG00000198961-A
ENSG00000245526-A

# 细胞条码
# 没啥特殊的
$ tail -3 simpleaf_quant/af_quant/alevin/quants_mat_rows.txt
TTCGATTTCTGAATCG
TGCTCGTGTTCGAAGG
ACTGTGAAGAAATTGC

读取方式

读取时也不能用Seurat和scanpy默认的函数。对于alevin-fry处理的数据，正确的打开方式是使用Python的pyroe包或R的fishpond包（其实很好记！大家看出来名字里的梗了吗？salmon-roe-fishpond，三文鱼-鱼子酱-鱼塘）。

# python
import pyroe
# 只读取spliced和ambiguous counts
adata_sa = pyroe.load_fry('simpleaf_quant/af_quant')
# 读取所有counts（包括unspliced）
adata_usa = pyroe.load_fry('simpleaf_quant/af_quant',
                          output_format={'X' : ['U','S','A']})

# R
library(fishpond)
# 读取数据
sce <- loadFry('alevin_fry_quant')

05

后记

第三章的实战，大家感受如何？本书关于上游的部分就结束了，之后我们将会正式进入下游流程。首先介绍的是基于python的anndata-mudata-scanpy生态系统，包括具体的环境配置、anndata和mudata的结构解剖，以及和其他单细胞生态相互转换的方式等等。敬请期待！

本章英文原文、扩展阅读，以及用到的示例数据、完整代码，可以在后台回复scbp03获取。

----------

觉得本系列对你有用的话，欢迎加🌟，不错过最新推送！

Single Cell Tutorial