全长转录组 | Iso-Seq 三代测序数据分析流程 (PacBio) （1）

我也认为长读长测序是 RNA 测序的未来！随着价格的降低和碱基质量的提升，传统的二代RNA-seq会被逐渐取代。

很多物种的转录本非常多样和复杂，绝大多数真核生物基因不符合“一基因一转录本”的模式，这些基因往往存在多种可变剪切（Alternative splicing，AS）形式。目前，基于第二代测序技术的RNA测序（RNA-seq）技术已被广泛用于各种转录组研究。但其测序的序列读长较短（50-300bp），大多只能覆盖转录本的一小部分，导致难以精确重构同一转录本的同源异构体（isoform），因此使得二代RNA测序对于全长转录本的重构是不准确的，片面的。

全长转录组（Full-length transcriptome）测序和分析是基于PacBio和Oxford Nanopore三代测序平台，利用其长读长的特性，建库测序时无需对RNA进行打断，如直接获得包含5’UTR、3’UTR、polyA尾的mRNA全长序列及完整结构信息，从而准确分析有参考基因组物种的可变剪接及融合基因等结构信息，克服无参考基因组物种转录本拼接组装较短、信息不完整的难题 （图1）。

图1. 长度长能精确重构mRNA同源异构体

最早PacBio关于全长转录组的产品命名为Iso-Seq, 全称叫做 Isoform-Sequencing, 是对自家开发的转录本测序技术的规范化命名。现在使用最新的试剂盒SMRTbell prep kit 3.0进行测序文库的构建。2023年10月初推出的Kinnex全长RNA建库试剂盒，能将5个转录本串联成一条测序read，充分利用其读长长的优势，提高测序通量，配合Revio系统一起使用，使得对全长转录本进行定量变得更为实际。

Iso-Seq 方法可对整个 cDNA 分子（长达 10 kb 或更长）进行测序，无需进行生物信息学转录本组装，因此可以对批量（bulk）和单细胞转录本组中的新基因和异构体进行表征，并进一步：

• 鉴定可变剪接 (AS) 事件，包括可变起始位点、终止位点、内含子保留和外显子跳跃事件。
• 通过开放阅读框 (ORF) 预测新型同源异构体的功能影响。
• 检测差异表达的同源异构体和同源异构体的转换事件。
• 发现肿瘤样本中的基因融合事件。
• 识别等位基因同源异构体。

一、PacBio Iso-Seq 的实验流程

通过 PacBio SMRTbell prep kit 3.0 构建Iso-Seq测序文库，适用于 PacBio Sequel II 和 Revio 仪器型号 （图2）。

图2. Iso-Seq的建库实验流程

（1）通过带有Oligo-dT的引物富集含有polyA尾的mRNA。

（2）使用 Iso-Seq RT enzyme对mRNA进行反转录。

（3）加入模板转换oligo (Template Switch Oligo, TSO)

（4）通过PCR扩增富集合成的cDNA，此步骤可加入barcode，用于混样。

（5）对全长cDNA进行损伤修复、末端修复、末端加A尾。

（6）连接SMRT哑铃型测序接头，最后结合测序引物，绑定DNA聚合酶，形成完整的SMRT-bell测序文库。

二、Iso-Seq 全长转录组分析流程

1. Iso-seq基础概念 （1）

• ROI：reads of insert

ROI , 全称 reads of insert，可以理解为插入片段，首先看下三代测序文库构建阶段的reads示意图，如图3：

图3. 全长转录组文库的结构

对于上述的文库片段，测序产生的reads 示意图，如图4：

图4. 全长转录组reads的结构

由于是一个环状分子， 随着测序反应的进行，会循环测序；如果把插入片段的正负链都测了一次，就做1个full pass。对于CCS 而言，要求至少有2个full pass , 才能去生成CCS reads。三代测序的特点就是读长很长，可以达到十几kb, 对于短的插入片段而言，CCS这样定义当然没有问题，但是对于全长转录本而言，转录本长度很长，比如转录本长度1kb, 读长3kb, 此时在一个零模波导孔（ZMW）中测序的reads 就不可能达到2个full pass , 也就产生不了CCS reads, 为了解决这个问题，提高reads的利用率，提出了ROI 的概念，ROI 指的就是插入片段，图4测序reads 产生的ROI 如图5：

图5. ROI的结构

ROI 不要求满足2个full pass， 相对CCS 而言，更加适合全长转录本的分析。

• Artifacts, 文库构建过程中可能产生的非正常转录本

可以理解为，共有两种来源：

Artificial Concatemer

图6. Artificial Concatemer的结构

这种序列是由于文库制备阶段，adapter 序列错误的将两条转录本的序列链接构成了一个环状分子，这个和adapter 浓度有关，通常这种reads 产生的比例很少，小于0.5%， 在后续的分析中，这部分reads 需要去除。

PCR Chimera

图7. PCR Chimera

在PCR 反应中，由于不完全延伸的产物作为了下次扩增反应的引物，导致出现嵌合体序列，直观上看，就是PCR产物来源于两条或者多条reads。PCR 产生的嵌合体序列，在PCR 反应体系中，这种序列是不可避免的，大约有3%的比例，在后续的分析过程中，可以借助软件去除这部分reads。

• FL Reads

FL , Full-length reads, 全长转录本。从raw data 到 ROI , 在从ROI 去除 artifacts reads 之后，我们就得到了用于后续分析的clean reads。clean reads 就已经是转录本的序列了，我们首先看一下clean reads 当中，哪些是全长转录本；哪些不是全长转录本。

全长转录本的示意图，如图8：

图8. 全长转录本的定义示意图

对于全长转录本而言，其ROI reads 中包含5‘ primer 和 3‘ primer; 而且会出现polyA 为结构；（polyA 针对mRNA和部分lncRNA）。

2. Iso-Seq软件

Iso-Seq是PacBio官方开发的一款用PacBio subreads 或 HiFi 数据进行全长转录组分析的一款软件，最终输出高质量转录本的全长序列。截止到2023年6月7号，最新版本为4.0.0。

Github主页：https://github.com/PacificBiosciences/IsoSeq

• 软件安装：isoseq和lima

#使用conda安装isoseq，v4.0.0. 
$ conda install -c bioconda isoseq

#lima,用于拆解barcode.
$ conda install -c bioconda lima

• Iso-Seq分析流程

整个Iso-Seq的分析流程如图9：

（1）原始下机数据 subreads.bam 通过CCS软件获得HiFi Reads, hifi.reads.bam。

（2）cDNA两端引物的去除，拆解barcode，调整转录本5' - 3' (polyA)方向。

（3）3' polyA尾和嵌合体（concatemer）序列的去除。

（4）转录本聚类。

（5）Consensus的转录本序列以.fasta格式输出。

图9. Iso-seq分析流程

• 单个样本官方示例数据演示

（1）示例数据下载

# Download toy S-read dataset
# This is a toy dataset consisting of ~80k segmented reads (S-reads) from a Kinnex full-length RNA library

$ wget https://downloads.pacbcloud.com/public/dataset/IsoSeq_sandbox/human_80k_Sreads.segmented.bam

# Download the Iso-Seq v2 cDNA primers (from Iso-Seq express 2.0 kit)

$ wget https://downloads.pacbcloud.com/public/dataset/Kinnex-full-length-RNA/REF-primers/IsoSeq_v2_primers_12.fasta

图10. Hifi reads作为输入序列

Hifi reads两端含有 5' 和 3' 引物序列。

（2）demultiplex，使用lima拆解barcode

$ lima --version
lima 2.9.0

$ lima human_80k_Sreads.segmented.bam IsoSeq_v2_primers_12.fasta human_80k.bam --isoseq --peek-guess

对于iso-seq数据，使用--isoseq加--peek-guess参数来降低假阳性率。可以使用 --biosample-csv input.csv添加样本名称, bio sample name。

图10. 运行lima后产生的文件

图11.Demultplex和引物去除

Demultplex和 5' - 3' 引物去除后，得到含有polyA尾序列的 Full-Length reads (FL reads)。

（3）refine，使用isoseq refine去除poly(A)和嵌合体（concatemer）序列

• 输入文件为：<primer--pair>.fl.bam和primers.fasta。
• 输出文件主要为：<movie>.flnc.bam。

$ isoseq refine human_80K.IsoSeqX_bc10_5p--IsoSeqX_3p.bam IsoSeq_v2_primers_12.fasta human_80K.flnc.bam --require-polya

$ ls human_80K.flnc.*
human_80K.flnc.bam
human_80K.flnc.bam.pbi
human_80K.flnc.consensusreadset.xml
human_80K.flnc.filter_summary.report.json
human_80K.flnc.report.csv

--require-polya ：只有有polyA尾序列才会被认作full length，并且去除polyA序列。

-j，--num-threads： CPU线程数。

--min-polya-length： 最小polyA尾长度，默认值为20。

3' polyA尾和嵌合体（concatemer）序列的去除后得到 Full-Length Non-Concatemer (FLNC) reads。

（4）cluster，转录本聚类

• 输入文件：<movie>.flnc.bam or flnc.fofn。
• 输出文件：<prefix>.bam

$ isoseq cluster2 human_80K.flnc.bam human_80K.transcripts.bam

$ ls human_80K.transcripts.*
human_80K.transcripts.bam
human_80K.transcripts.bam.pbi
human_80K.transcripts.cluster_report.csv

输出转录本的同源异构体（isoforms）至少有两个或两个以上的FLNC（full-length non-concatemer）序列支持。如果想包含singletons，可以加入--singletons。

关于满足什么条件isoseq cluster算法会将两条序列聚类为同一个转录本序列（图13）：

图13. 转录本聚类原则

（A）5'端差异小于100bp以内。

（B）3'端差异小于30bp以内。

（C）小于10bp以内的gaps，gaps的数量没有上限。

图14. Cluster Isoform

• 混样数据官方示例演示

（1）示例数据下载

# This is a 12-plex regular Iso-Seq (non-Kinex) run on Sequel II system consisting of ~3 million HiFi reads.
# Download HiFi reads from a non-Kinnex (regular Iso-Seq) BAM file
$ wget https://downloads.pacbcloud.com/public/dataset/Kinnex-full-length-RNA/DATA-SQ2-UHRR-Monomer/1-CCS/m64307e_230628_025302.hifi_reads.bam
$ wget https://downloads.pacbcloud.com/public/dataset/Kinnex-full-length-RNA/DATA-SQ2-UHRR-Monomer/1-CCS/m64307e_230628_025302.hifi_reads.bam.pbi

# Download the Iso-Seq v2 cDNA primers (from Iso-Seq express 2.0 kit)
$ wget https://downloads.pacbcloud.com/public/dataset/Kinnex-full-length-RNA/REF-primers/IsoSeq_v2_primers_12.fasta

（2）demultiplex，使用lima拆解barcode

$ lima --version
lima 2.9.0

# Demux and primer removal
$ lima --isoseq --peek-guess m64307e_230628_025302.hifi_reads.bam IsoSeq_v2_primers_12.fasta UHRR.bam

每个barcode对输出一个.bam文件，一共12个.bam文件对应12个样本。

图15. lima运行后的输出文件

（3）合并拆分样本的文件名

# Combine inputs
$ ls UHRR.IsoSeqX*bam > all.fofn

$ cat all.fofn

UHRR.IsoSeqX_bc01_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc02_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc03_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc04_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc05_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc06_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc07_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc08_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc09_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc10_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc11_5p--IsoSeqX_3p.bam
UHRR.IsoSeqX_bc12_5p--IsoSeqX_3p.bam

fofn: "files-of-file-names"的缩写。

（4）refine，使用isoseq refine去除poly(A)和嵌合体（concatemer）序列

# Remove poly(A) tails and concatemer
$ isoseq refine all.fofn IsoSeq_v2_primers_12.fasta UHRR.flnc.bam --require-polya


$ ls UHRR.flnc.*
UHRR.flnc.bam
UHRR.flnc.bam.pbi
UHRR.flnc.consensusreadset.xml
UHRR.flnc.filter_summary.report.json
UHRR_80K.flnc.report.csv

（5）cluster，转录本聚类

$ isoseq cluster2 UHRR.flnc.bam UHRR.transcripts.bam

（6）Polish （可选择）

$ isoseq cluster flnc.fofn clustered.bam --verbose --use-qvs

这里使用isoseq cluster，而不是isoseq cluster2, cluster相比于cluster2比较耗时。Polish会略微提升数据质量，不是必须步骤。 运行完成以后获得以下文件：

• <prefix>.bam
• <prefix>.hq.fasta.gz with predicted accuracy ≥ 0.99
• <prefix>.lq.fasta.gz with predicted accuracy < 0.99
• <prefix>.bam.pbi
• <prefix>.transcriptset.xml

综上所展示，iso-seq分析流程及每一步输出文件总览，如图16。

图16. iso-seq总体流程及每一步产生文件总览

参考文献

1. Iso-seq 必备基础-blog.csdn
2. pacbio 三代全长转录组数据分析流程
3. PacBio Iso-Seq Workshop Online