全基因組 - 人類基因組變異分析（PacBio) （5）-- pbsv

以Pacific Biosciences （PacBio）公司為代表的第三代單分子實時熒光測序平臺為以測序為基礎(chǔ)的科學(xué)研究帶來了革命性的突破。目前該技術(shù)廣泛應(yīng)用于基因組Denovo組裝、全長轉(zhuǎn)錄本檢測、宏基因組，基因組重測序等多個方向，并且在染色體結(jié)構(gòu)變異（Structure Variation, SV）的檢測中有著不可替代的優(yōu)勢。前面我們講了PacBio三代測序數(shù)據(jù)的類型、預(yù)處理、比對和SNPs/INDELs變異檢測等基本內(nèi)容。本期我們就繼續(xù)沿著分析流程圖一起看看基于比對結(jié)果檢測染色體結(jié)構(gòu)變異（SV）分析方法。

一、染色體結(jié)構(gòu)變異

染色體結(jié)構(gòu)變異（Structure Variation, SV），指基因組上發(fā)生的長度大于50bp的大片段插入(Insertion, INS)、缺失(Deletion, DEL)、倒位(Inversion, INV)、易位(Translocation)、重復(fù)（Duplication, DUP）等類型的變異，其中占比最大的就是大片段的插入和缺失（圖1）。插入缺失很好理解就是，多了一段或者少了一段DNA序列；重復(fù)就是有一段區(qū)域的序列重復(fù)出現(xiàn)；倒位就是序列翻轉(zhuǎn)了一下，如本來那個位置該是AATTG的，結(jié)果變成了GTTAA；易位的話就是序列位置的變化，又進一步分為染色體內(nèi)易位和染色體間易位。據(jù)統(tǒng)計，基因組結(jié)構(gòu)變異可能導(dǎo)致的遺傳性疾病已經(jīng)超過1,000種，對于每個人來講其基因組都有至少20,000個的結(jié)構(gòu)變異，這些變異帶來的影響或許比SNVs或InDels帶來的影響更大。

圖1. 結(jié)構(gòu)變異的類型

二、三代測序在結(jié)構(gòu)變異檢測上的優(yōu)勢

如果是僅僅研究像SNPs/INDELs一樣的很小的變異，二代測序就能夠勝任；但是如果要研究很大的結(jié)構(gòu)變異（>50bp），則二代測序的短讀長很難識別變異位點，尤其是復(fù)雜的高CG含量的區(qū)域，而且僅僅依靠算法很難克服技術(shù)上的缺陷。三代測序的長讀長能夠很有效的跨越覆蓋識別出結(jié)構(gòu)變異位點，得到結(jié)構(gòu)變異的全貌，輕松測通基因組上的復(fù)雜重復(fù)區(qū)域。通過三代測序技術(shù)，在人類基因組中發(fā)現(xiàn)了數(shù)萬個結(jié)構(gòu)變異，而這些變異通常無法通過二代測序技術(shù)進行識別（圖2）。

圖2. 三代測序之于二代測序的優(yōu)勢

從Medhat Mahmoud等人 (5) 的研究中可以看出（圖3B），在人類中，三代長度長測序比二代短讀長測序在結(jié)構(gòu)變異檢測數(shù)量上平均高2-4倍。

圖3. 三代測序技術(shù)在結(jié)構(gòu)變異檢測方面明顯優(yōu)于二代測序

三、SV檢測相關(guān)軟件

目前，對于PacBio三代測序平臺結(jié)構(gòu)變異檢測的軟件有PBSV, Sniffles, cuteSV、PAV、PBHoney、SMRT-SV、SVIM 等 (圖4，圖5，圖6）。

圖4. 適合于PacBio數(shù)據(jù)的結(jié)構(gòu)變異檢測軟件

  1. PBHoney，文章表于2014年 (7)，軟件最后的更新停留在了2017年，所以已經(jīng)不推薦了。

  2. SMRT-SV, 第一個版本由Chaisson et al. 發(fā)表于2015年 (8)，早已經(jīng)不再更新。SMRT-SV2由Audano et al.發(fā)表于2019年 (9)，替代SMRT-SV?，F(xiàn)在SMRT-SV2也已經(jīng)停止更新了，并在其github網(wǎng)站上推薦以下工具：

  1. DeBreak，文章于2023年發(fā)表在Nature Communication上，github上更新到 2022年12月1號(version 1.2）。從文章中的數(shù)據(jù)看來，DeBreak 在模擬數(shù)據(jù)（圖5）和實測數(shù)據(jù)HG002（圖6）中recall, precision和 F1 score上優(yōu)于同類軟件。
圖5. DeBreak在模擬數(shù)據(jù)上的表現(xiàn)

圖6. DeBreak在HG002數(shù)據(jù)上的表現(xiàn)
  1. PAV，Phased Assembly Variant Caller，文章于2021年發(fā)表于Science (11)，github更新到2023年10月13號(version 2.3.4)。
  2. PBSV，PacBio官方開發(fā)的結(jié)構(gòu)變異軟件，github上更新至2023年3月14日（version 2.9.0）。
  3. Sniffles2，版本一于2018年發(fā)表于Nature methods（12），Sniffles2于2022年發(fā)表于BioRxiv（13），github更新到2023年7月14號（version 2.2）
  4. cuteSV，文章于2020年發(fā)表于Genome Biology (14)， github上更新至2023年11月14號（version 2.1）。

綜上所述，現(xiàn)在常用以及保持更新的有以下這些，PacBio官方的PBSV、Sniffles2、cuteSV，DeBreak 和 PAV了 （此處可能總結(jié)的不全面，歡迎留言補充）。

前面咱們是用PacBio官方的比對軟件pbmm2，由于pbmm2生成的 .bam文件MD tag 格式有些特別，無法用于后續(xù)的sniffles進行SV檢測。所以這里我們使用pbsv進行后續(xù)SV檢測。后面我會針對minimap2+Sniffles2組合 以及minimap2+cuteSV再出一期教程，有時間也會嘗試一下PAV和DeBreak。

四、pbsv具體安裝和使用

主頁網(wǎng)址：https://github.com/PacificBiosciences/pbsv

  • pbsv is a suite of tools to call and analyze structural variants in diploid genomes from PacBio single molecule real-time sequencing (SMRT) reads.
  • pbsv calls insertions, deletions, inversions, duplications, and translocations. Both single-sample calling and joint (multi-sample) calling are provided.
  • 軟件安裝
#pvsv v.2.9.0, conda 安裝
$ conda install -c bioconda pbsv
  • 工作流程

整個工作流程從序列比對到獲得結(jié)構(gòu)變異一共分三步（圖7）：

  1. PacBio reads和參考基因組進行比對，從.subreads.bam / .ccs.fastq到比對排序后的.bam文件。
  2. 尋找結(jié)構(gòu)變異的特征，.bam到.svsig.gz。
  3. 獲得單個或者所有樣本的結(jié)構(gòu)變異和基因型，.svsig.gz到.vcf
圖7.pbsv軟件工作流程
  • 具體分析命令

數(shù)據(jù)我們還是使用德系猶太人家系：HG002（子）、HG003（父）、HG004（母），具體參考全基因組 - 人類基因組變異分析（PacBio) （3）-- pbmm2

1. Align PacBio reads to a reference genome 序列比對。

# 以CCS bam input為示例
$ pbmm2 align ref.fa movie1.ccs.bam ref.movie1.bam --sort --preset CCS --sample sample1

#實際運行
#HG002_1
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG002_1.log --sample HG002_1 \
Human_ref/GRCh38.mmi m84011_220902_175841_s1.hifi_reads.bam HG002_1.align.bam 

#HG003
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG003.log --sample HG003 \
Human_ref/GRCh38.mmi m84010_220919_235306_s2.hifi_reads.bam HG003.align.bam 

#HG004
$ pbmm2 align --preset CCS --sort -j 24 -J 4 \
 --log-level INFO --log-file pbmm2.HG004.log --sample HG004 \
Human_ref/GRCh38.mmi m84010_220919_232145_s1.hifi_reads.bam HG004.align.bam

2. Discover signatures of structural variation 尋找結(jié)構(gòu)變異的特征。

舉例：

$ pbsv discover ref.movie1.bam ref.sample1.svsig.gz
$ pbsv discover ref.movie2.bam ref.sample2.svsig.gz

# optionally index svsig.gz to allow random access via `pbsv call -r`
tabix -c '#' -s 3 -b 4 -e 4 ref.sample1.svsig.gz
tabix -c '#' -s 3 -b 4 -e 4 ref.sample2.svsig.gz

It is highly recommended to provide one tandem repeat annotation .bed file of your reference to pbsv discover via --tandem-repeats. This increases sensitivity and recall. Feel free to use the following for human SV calling: GRCh38 or hs37d5/hg19. 加上提供的重復(fù)串聯(lián)序列區(qū)域注釋文件可以提高敏感度和召回率。

實際運行：

#HG002_1
$ nohup pbsv discover --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed HG002_1.align.bam HG002_1.svsig.gz &

#HG003
$ nohup pbsv discover --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed HG003.align.bam HG003.svsig.gz &

#HG004
$ nohup pbsv discover --tandem-repeats Human_ref/human_GRCh38_no_alt_analysis_set.trf.bed HG004.align.bam HG004.svsig.gz &

3. Call structural variants and assign genotypes 獲得單個或者所有樣本的結(jié)構(gòu)變異和基因型。

#示例，如果輸入序列是CCS序列，加上--ccs選項
$ pbsv call ref.fa ref.sample1.svsig.gz ref.sample2.svsig.gz ref.var.vcf

#實際運行，Joint calling

nohup pbsv call -j 24 --ccs Human_ref/GRCh38.fa HG002_1.svsig.gz HG003.svsig.gz HG004.svsig.gz HG.SV.vcf &

至此，我們已經(jīng)獲得了三個樣本包含結(jié)構(gòu)變異信息的vcf格式文件了。

pbsv call 幫助文檔， 除上述示例中的參數(shù)，PBSV還有諸多篩選參數(shù)，包括輸出的SV類型、長度、基因型等等，需要大家根據(jù)自己的實際項目情況，進行探索，選擇合適參數(shù)。比如用 -t 可以指定特定結(jié)構(gòu)變異類型，-r 可以指定基因組區(qū)域。全部參數(shù)說明如下圖所示：

pbsv call - Call structural variants from SV signatures and assign genotypes (SVSIG to VCF).

Usage:
  pbsv call [options] <ref.fa|xml> <ref.in.svsig.gz|fofn> <ref.out.vcf>

  ref.fa|xml                                  STR   Reference genome assembly against which to call variants.
  ref.in.svsig.gz|fofn                        STR   SV signatures from one or more samples.
  ref.out.vcf                                 STR   Variant call format (VCF) output.

Basic Options:
  -r,--region                                 STR   Limit discovery to this reference region: CHR|CHR:START-END.
  --hifi,--ccs                                      Use options optimized for HiFi reads: -S 0 -P 10.

Variant Options:
  -t,--types                                  STR   Call these SV types: "DEL", "INS", "INV", "DUP", "BND".
                                                    [DEL,INS,INV,DUP,BND]
  -m,--min-sv-length                          STR   Ignore variants with length < N bp. [20]
  --max-ins-length                            STR   Ignore insertions with length > N bp. [15K]
  --max-dup-length                            STR   Ignore duplications with length > N bp. [1M]

SV Signature Cluster Options:
  --cluster-max-length-perc-diff              INT   Do not cluster signatures with difference in length > P%. [25]
  --cluster-max-ref-pos-diff                  STR   Do not cluster signatures > N bp apart in reference. [200]
  --cluster-min-basepair-perc-id              INT   Do not cluster signatures with basepair identity < P%. [10]

Consensus Options:
  -x,--max-consensus-coverage                 INT   Limit to N reads for variant consensus. [20]
  -s,--poa-scores                             STR   Score POA alignment with triplet match,mismatch,gap. [1,-2,-2]
  --min-realign-length                        STR   Consider segments with > N length for re-alignment. [100]

Call Options:
  -A,--call-min-reads-all-samples             INT   Ignore calls supported by < N reads total across samples. [3]
  -O,--call-min-reads-one-sample              INT   Ignore calls supported by < N reads in every sample. [3]
  -S,--call-min-reads-per-strand-all-samples  INT   Ignore calls supported by < N reads per strand total across samples
                                                    [1]
  -B,--call-min-bnd-reads-all-samples         INT   Ignore BND calls supported by < N reads total across samples [2]
  -P,--call-min-read-perc-one-sample          INT   Ignore calls supported by < P% of reads in every sample. [20]
  --preserve-non-acgt                               Preserve non-ACGT in REF allele instead of replacing with N.

Genotyping:
  --gt-min-reads                              INT   Minimum supporting reads to assign a sample a non-reference
                                                    genotype. [1]

Annotations:
  --annotations                               FILE  Annotate variants by comparing with sequences in fasta.Default
                                                    annotations are ALU, L1, SVA.
  --annotation-min-perc-sim                   INT   Annotate variant if sequence similarity > P%. [60]

Variant Filtering Options:
  --min-N-in-gap                              STR   Consider >= N consecutive "N" bp as a reference gap. [50]
  --filter-near-reference-gap                 STR   Flag variants < N bp from a gap as "NearReferenceGap". [1K]
  --filter-near-contig-end                    STR   Flag variants < N bp from a contig end as "NearContigEnd". [1K]

  -h,--help                                         Show this help and exit.
  --version                                         Show application version and exit.
  -j,--num-threads                            INT   Number of threads to use, 0 means autodetection. [0]
  --log-level                                 STR   Set log level. Valid choices: (TRACE, DEBUG, INFO, WARN, FATAL).
                                                    [WARN]
  --log-file                                  FILE  Log to a file, instead of stderr.

Copyright (C) 2004-2023     Pacific Biosciences of California, Inc.
This program comes with ABSOLUTELY NO WARRANTY; it is intended for
Research Use Only and not for use in diagnostic procedures.

五、結(jié)構(gòu)變異結(jié)果說明

以50bp為臨界，不論是小突變（SNV/INDEL）還是大的染色體結(jié)構(gòu)變異（SV），通常情況下生成的結(jié)果文件均為.vcf文件格式，也是信息分析中比較常見的一種格式，百度里也有詳盡的結(jié)果文件格式說明。圖8為HG002，HG003和HG004三個樣本joint calling結(jié)構(gòu)變異后的部分結(jié)果展示。

圖8. vcf文件展示

參考文獻

  1. 神燈寶典之PB三代重測序分析實錄（一）
  2. 神燈寶典之三代重測序分析實錄（二）
  3. 三代測序時代的臨床科研
  4. 三代重測序到底能干什么？
  5. Mahmoud, M., Gobet, N., Cruz-Dávalos, D. I., Mounier, N., Dessimoz, C., & Sedlazeck, F. J. (2019). Structural variant calling: the long and the short of it. Genome Biology.
  6. giab_tools_methods
  7. English, Adam C., William J. Salerno, and Jeffrey G. Reid. "PBHoney: identifying genomic variants via long-read discordance and interrupted mapping." BMC Bioinformatics 15 (2014).
  8. Chaisson, Mark JP, John Huddleston, Megan Y. Dennis, Peter H. Sudmant, Maika Malig, Fereydoun Hormozdiari, Francesca Antonacci et al. "Resolving the complexity of the human genome using single-molecule sequencing." Nature (2015). SMRT-sv
  9. Audano, P. A., Sulovari, A., Graves-Lindsay, T. A., Cantsilieris, S., Sorensen, M., Welch, A. E.,
    Dougherty, M. L., Nelson, B. J., Shah, A., Dutcher, S. K., et al. (2019). Characterizing the major structural variant alleles of the human genome. Cell.
  10. Chen, Yu, Amy Y. Wang, Courtney A. Barkley, Yixin Zhang, Xinyang Zhao, Min Gao, Mick D. Edmonds, and Zechen Chong. "Deciphering the exact breakpoints of structural variations using long sequencing reads with DeBreak." Nature Communications (2023).
  11. Ebert et al., “Haplotype-Resolved Diverse Human Genomes and Integrated Analysis of Structural Variation,” Science. 2021.
  12. Sedlazeck, Fritz J., Philipp Rescheneder, Moritz Smolka, Han Fang, Maria Nattestad, Arndt Von Haeseler, and Michael C. Schatz. "Accurate detection of complex structural variations using single-molecule sequencing." Nature methods (2018).
  13. Smolka, Moritz, Luis F. Paulin, Christopher M. Grochowski, Dominic W. Horner, Medhat Mahmoud, Sairam Behera, Ester Kalef-Ezra et al. "Comprehensive structural variant detection: from mosaic to population-level." BioRxiv (2022).
  14. Jiang, Tao, Yongzhuang Liu, Yue Jiang, Junyi Li, Yan Gao, Zhe Cui, Yadong Liu, Bo Liu, and Yadong Wang. "Long-read-based human genomic structural variation detection with cuteSV." Genome biology (2020).
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容