这篇文章描述了如何使用pbsv软件从Pacbio测序数据中检测基因组结构变异（SV）。
- 作者：Aaron Wenger， Yuan Li
2.1. 安装PacBio pbsv软件
- pbsv软件是PacBio SMRTAnalysis系统软件的一部分，下载地址： SMRT Analysis Software - PacBio
- 输入以下命令以检测pbsv是否可从bash终端（bash terminal）直接调用
which pbsv
2.2. 准备测试 测试之前准备工作包括
2.1 下载测试数据包
这个测试数据包文件大小为5MB，包含从12个SMRTCells中挑选的3个人类样本的数据，分别为
- HG00733 (女儿样本)：6条reads
- HG00731 (父亲样本)：3条reads
- HG00732 (母亲样本)：3条reads
这些reads都来自人类基因组中大约100KB的区域。
在bash终端输入以下命令以下载并解压测试数据包
wget https://dl.dnanex.us/F/D/8XK8Ff42Kf1zfqVqVQ6ggXf5V20zxF3V34014P2Q/subreads.tar.gz
tar -xzvf subreads.tar.gz
文件夹包含三个子文件夹HG00731,HG00732,HG00733，每个子文件夹对应一个人类样本。
2.2 下载人类参照基因组
如果服务器已经包含人类参照基因组hg38的FASTA文件，不需要重新下载；亦可使用其他版本的人类参照基因组。 人类参照基因组hg38的FASTA文件大小是3 GB。
在bash终端输入：
wget https://dl.dnanex.us/F/D/Bp35xf0kqQGpqk0zybgg340f1Kfv3xGf9Zz75ZP2/hg38.fa
使用pbsv软件对多样本测序数据进行联合检测基因组结构变异 执行时间：此测试通常在一小时内完成。
3.1 给每个人类样本分别生成一份subreadset.xml文件
- subreadset.xml — PacBio Dataset文件，可通过PacBio软件接口访问原始测序数据reads
在bash终端输入以下命令以生成三份subreadset.xml文件
for sample in HG00731 HG00732 HG00733; do
find ${sample} -name "*.subreads.bam" |
xargs dataset create --type SubreadSet --name "${sample} subreads" ${sample}.subreadset.xml
done
3.2 生成movie2sample.json文件
- movie2sample.json 是用户自定义JSON文件，提供从SMRT cells文件到人类标本的映射。
- 如要生成自定义的movie2sample.json文件， 有三种方法，任选一种即可。
3.2.1 手动编辑Json文件
movie2sample.json文件包含一个列表，每个元素是一对SMRTCell到标本的映射，元素之间用,逗号分隔（注意不是中文的逗号，）。
比如以下文件表示SMRTCell m54006_170929_182745, m54086_171004_235550, m54086_171005_100202测序样本为HG00731， m54006_170930_184310, m54006_171001_065147, m54086_171006_062139, 测序样本为HG00732，其他SMRTCell测序样本为HG00733。
[["m54006_170929_182745","HG00731"],["m54086_171004_235550","HG00731"],["m54086_171005_100202","HG00731"], ["m54006_170930_184310","HG00732"],["m54006_171001_065147","HG00732"],["m54086_171006_062139","HG00732"], ["m54006_171002_070909","HG00733"],["m54006_171025_222253","HG00733"],["m54007_171026_232716","HG00733"], ["m54007_171027_113626","HG00733"],["m54007_171027_234711","HG00733"],["m54043_171109_220647","HG00733"]]
3.2.2 bash终端awk命令
在bash终端输下命令，手动生成movie2sample.json文件
find HG00731 HG00732 HG00733 -name "*.subreads.bam" | cut -d'.' -f1 | awk -F'/' '{ print "[\"" $2 "\",\"" $1 "\"]"; }' | tr '\n' ',' | sed -e 's/,$/\n/' | awk '{ print "[" $0 "]"; }' > movie2sample.json
3.2.3 使用pbsvutil make-samples
如果你事先在SMRTLink Data Management模块手动编辑每个SMRTCell的样本名，并合并多个SMRTCells得到一个包含多个样本的SubreadSet，你也可以通过pbsvutil make-samples来生成movie2sample.json文件。
pbsvutil make-samples merged.subreadset.xml out_samples.json movie2sample.json
3.3 pbsv align比较和匹配测试数据与参照基因组
使用pbsv align，把每个样本的测序数据与参照基因进行比较和匹配，共生成三个Alignment BAM文件。pbsv align程序从 movie2sample.json文件中得到测序数据和样本的映射关系，从而确定结构变异在哪些样本中存在。
for sample in HG00731 HG00732 HG00733; do
pbsv align hg38.fa ${sample}.subreadset.xml hg38.${sample}.bam --movienames2samples_json movie2sample.json
done
3.4 生成FOFN文件
FOFN即 file of file names，FOFN文件中的每一行都是一个Alignment BAM文件路径。
find . -name "hg38.HG*.bam" > trio.fofn
3.5 pbsv call以检测结构变异
pbsv call hg38.fa trio.fofn hg38.trioSv.vcf
输入文件：hg38.fa 以及trio.fofn
输出文件：hg38.trioSv.vcf，其中包括两个Insertion结构变异。第一个位于chr1:37,967,108 ，在HG00731样本中为纯合子，在其余两个样本中为杂合子。第二个结构变异位于chr1:37,982,030， 在所有三个样本中都是杂合子。
4.1 pbsv align
测试步骤3.3中的三个pbsv align的进程相互独立，可以在不同服务器上同时并行处理。比如，可以通过qsub把三个pbsv align进程分配给不同服务器运行。
hg38.fa-enc.2.ngm, hg38.fa-ht-13-2.2.ngm。
touch empty.fa
pbsvutil ngmlr hg38.fa empty.fa empty.out.bam
4.2 pbsv call
在不同染色体的结构变异检测相互独立，可以在多个服务器上并行。比如以下每个命令都可以分别qsub到不同服务器执行。
for chromosome in $(seq 1 1 22) X Y; do
echo pbsv call hg38.fa trio.fofn out.chr${chromosome}.vcf --reference_regions chr${chromosome}
done
所有染色体的结构变异检测进程结束后，执行以下命令以合并所有染色体的结构变异。
cat out.chr1.vcf > out.merged.vcf
for chromosome in $(seq 2 1 22) X Y; do
cat out.chr${chromosome}.vcf >> out.merged.vcf
done