干货 | 使用PBmm2+PBSV基于pacbio三代数据比对结果检测染色体结构变异的分析流程简介

随着基因组研究的不断推进，以PacBio测序为代表的第三代基因测序技术在生命科学相关学科的科研领域中发挥了不可替代的作用。目前，三代测序技术相比于传统的二代测序技术而言，由于测得的序列长度长，广泛应用于基因组Denovo、全长转录本检测、重测序等多个方向，并且在染色体结构变异（SV）的检测中有着不可替代的优势。在科学研究中，我们经常会对物种测序序列与参考基因组之间作比较，研究测序物种与参考基因组之间的变异类型，从而更好的揭示物种的多样性特征。下面，就跟随本文一起学习如何对三代测序数据做染色体结构变异分析吧。

结构变异简介

通常，基因组上的变异种类，按照目前业界的看法可以分为如下三个大类：一是单核苷酸多态性(SNP)；二是基因组上很短的Insertion和Deletion，也常被我们合并起来称为Indel。三是基因组结构性变异（Structure Variantions，简称SVs），也是本文进行检测的重点变异类型，一般是指基因组上大长度的序列变化和位置关系变化。有很多种类型（图1），通常定义是长度大于50bp的插入（Insertion）、缺失（Deletion）、串联重复（Tandem repeate）、染色体倒位（Inversion）、染色体内部或染色体之间的序列易位（Translocation）、拷贝数变异（CNV）以及形式更为复杂的嵌合性变异。其中，占比最大的就是Indel。

图1 | 结构性变异的不同种类

SV检测软件的介绍

PBmm2：pbmm2是PB官方基于minimap2进行优化的版本，优化的目的是支持原生PacBio数据的输入和输出。它提供了一组推荐参数，针对Pacbio的三代测序数据适配性好，测试表明，pbmm2在序列识别、比对上参考基因组的碱基数量，尤其是运行时间等方面优于BLASR，是Pacbio公司官方推荐的BLASR的替代比对软件。

PBSV：PBSV是一套call结构变异的软件，它通过对Pacbio公司的单分子实时测序数据进行分析，从而得到二倍体基因组中的大片段结构变异，工作流程如图所示（图2）。

PBSV可以有效的鉴定如下的结构变异类型。

· 20 bp到10 kb的插入

· 20 bp到100 kb的缺失

· 200 bp到10 kb的倒位

· 20 bp到10 kb的重复

· 不同染色体之间的易位或单个染色体上相隔100kb以上的易位

图2 | PBSV的工作流程

我们以人类hg38版本基因组序列为例，介绍SV分析流程。PBmm2与PBSV均可以通过bioconda获取。

conda install -c bioconda pbmm2

conda install -c bioconda pbsv

1.将PacBio测序数据比对到参考基因组上

由于是PB官方推出的比对软件，整体和下机数据格式的契合度比较高，几乎可以实现一行命令从下机bam文件到比对后的bam文件，且无需再进行比对上的reads提取和排序。整个软件的速度也比较快。

PBmm2提供了多种数据格式的输入，根据输入文件的格式运行对应的命令。

Subreads BAM input:

pbmm2 align hg38.fa movie1.subreads.bam hg38.movie1.bam --sort --median-filter --sample sample1

CCS BAM input:

pbmm2 align hg38.fa movie1.ccs.bam hg38.movie1.bam --sort --preset CCS --sample sample1

CCS FASTQ input:

pbmm2 align hg38.fa movie1.Q20.fastq hg38.movie1.bam --sort --preset CCS --sample sample1 --rg ''@RG\\tID:movie1''

常用参数整理

--sample指定样品名称

-j指定比对线程

-J指定排序步骤线程

--sort输出排序后的bam

--preset选择参数集

subreads数据可以使用--preset SUBREADS

如果担心比对结果不佳，需要输出未比对的reads进行排查，可以加上--unmapped参数。

2.发现结构变异的特征

对已经比对完的BAM文件，识别结构变异的特征，结构变异特征信息将会被存储在.svsig.gz文件中：

pbsv discover hg38.movie1.bam hg38.sample1.svsig.gz

3.Call结构变异并分配基因型

从样本的结构变异特征文件中call结构变异，此步骤可以同时输入多个.svsig.gz文件，如果输入为CCS reads，请将--CCS添加到下面的参数中：

pbsv call hg38.fa hg38.sample1.svsig.gz hg38.sample2.svsig.gz hg38.var.vcf

所有样本的结构变异信息，将会被call出，并写入到vcf格式的文件中。

小技巧：由于call sv这一步骤比较耗时，尤其是对大基因组物种以及高测序深度的数据。因此，推荐此步骤拆分进行。可以将参考基因组分解为多个染色体，对每一个单独的染色体并行运行，以提高运行速度。

分割.svsig.gz文件为每一条单独的染色体：

for i in $(samtools view -H hg38.movie1.bam | grep ''^@SQ'' | cut -f2 | cut -d'':'' -f2); do pbsv discover --region $i hg38.movie1.bam hg38.sample1.$i.svsig.gz done

Call SVs

pbsv call -j 8 hg38.fa hg38.sample1.*.svsig.gz hg38.sample1.vcf

-j为运行线程数

除上述示例中的参数，PBSV还有诸多筛选参数，包括输出的SV类型、长度、基因型等等，需要大家根据自己的实际项目情况，进行探索，选择合适参数。全部参数说明如下图所示：

图3 | pbsv参数说明

至此，我们已经获得了包含变异信息的vcf格式文件啦，你学会了嘛~

本文介绍了SV变异的基本知识，并且从生物信息分析的角度介绍了三代序列SV检测的分析过程和注意事项。结合本文，相信读者可以对SV的分析过程有一定的了解，如果各位小伙伴手上恰好有相关数据，不妨通过本文的流程实操一遍，相信你一定可以独立的完成从测序原始数据到结构变异的vcf文件的获取。一起动手实践起来吧~~

1、干货 | 秘籍宝典：免疫浸润分析-TIMER2.0介绍

2、干货 | UCSC数据库如何预测转录因子

3、干货 | 转录组建库起始量及不同建库试剂盒的系统比较

4、干货 | 单细胞测序送样指南

END

Shuai  撰文

本文系欧易生物原创

欢迎转发到朋友圈

转载请注明文本转自欧易生物