注: docker イメージのリンクも紹介してますが、テストするとエラーを吐きました。condaを使いlinuxマシンでに導入するのが無難なようです。
SNVsの分析はヒト遺伝学を研究するための標準的な技術となっているが[論文より ref.1]。、DNA配列(indels)の挿入と欠失は確実に検出することはできない[ref.2,3]。 Indelsはヒトゲノムで最も2番目に一般的な変異であり、構造変異中では最も多い[ref.4]。マイクロサテライト(単純配列反復、SSR、1〜6bpモチーフ)内で、indelsはリピートモチーフの長さを変え、40以上の神経学的疾患に関連している[ref.5]。 Indelsもまた、自閉症において重要な遺伝的要素を担っている。コードされたタンパク質を破壊する可能性のあるde novo indelsは、影響を受けていない兄弟よりも2倍近くも豊富である[ref.6]。
indel検出は、いくつかの理由から困難である。(1)indel配列とオーバーラップするリードはアライメントが難しく、gapではなく複数のミスマッチとして扱われることがある。 (2)エキソームシーケンシングのキャプチャ効率のばらつきおよび不均一なリード分布は、偽陽性の数を増加させる。 (3)エラー率増加は、マイクロサテライト内での検出を非常に困難にする。この研究で示されているように、(4)局在化、ほぼ同一の反復配列は、高い陽性率をもたらす可能性がある。これらの理由から、利用可能なソフトウェアツールで検出可能なindelサイズは比較的小さく、数十塩基を超えるものは少ない[ref.8]。
現在、indels検出には2つの主要なパラダイムが使用されている。最も一般的なアプローチは、リードマッパー(BWA、Bowtie、Novoalignなど)を使用してすべてのリードをリファレンスゲノムにマッピングすることだが、利用可能なアルゴリズムは数塩基以上のindel間のマッピングには有効ではない。先進的なアプローチではより長い変異を検出するためにペアエンド情報を使い local realignments を行うが(例えば、GATK UnifiedGenotyper[ref.1]およびDindel[ref.9])、実際には、より長い変異(≧20bp)ではその感度が大幅に低下する。 Split-read methods(例えば、Pindel[ref.10]およびSplitread[ref.11])は、理論的には任意のサイズの欠失を検出できるが、現在のシーケンス技術ではリード長が短いために(論文執筆時点)挿入を検出する能力は限られている。第2のパラダイムは、デノボ全ゲノムアセンブリを行い、組み立てられたコンティグとリファレンスゲノムとの間の変異を検出することからなる[ref.12,13]。より大きな突然変異を検出する可能性を有する一方で、実際には、このパラダイムは、ホモ接合型およびヘテロ接合型突然変異を正確に報告するために、細かくかつ局在化した分析が必要である。最近では、de novo aasemblyを使ったGATK HaplotypeCaller、SOAPindel[ref.14]、およびCortex[ref.15]の3種類のアプローチが開発されている。他の最近のアプローチであるTIGRA[ref.16]も、ローカルアセンブリを使用するが、ブレークポイントのみ検出するよう調整されており、indelsの配列は報告しない。
著者らは、exome-captureデータ内のindelsを検出するマイクロアセンブリパイプラインScalpelを提示する(論文より 図1)。マッピングとアセンブリの力を組み合わせることにより、Scalpelはde Bruijn graphを慎重に検索し、各エキソンにまたがるシーケンスパス(コンティグ)を探す。このアルゴリズムには、各エキソンのオンザフライリピート組成分析と、セルフチューニングのk-mer戦略が含まれる。
公式HP
http://scalpel.sourceforge.net/manual.html
マニュアル1
http://scalpel.sourceforge.net/manual.html
マニュアル2
https://sourceforge.net/p/scalpel/wiki/Manual/
Scalpelに関するツイート。
インストール
ubuntu18.04のAnaconda2.4.2でテストした。
#Anaconda環境ならcondaを使う(linux only)
conda install -c bioconda scalpel
#dockerイメージも提供されている。
docker pull hanfang/scalpel:0.5.3
docker imagesでIDを調べてから
> scalpel-discovery -h
$ scalpel-discovery -h
Local date and time: Sat Jul 7 10:13:11 2018
Program: scalpel-discovery (micro-assembly variant detection)
Version: 0.5.3 (beta), January 25 2016
Contact: Giuseppe Narzisi <gnarzisi@nygenome.org>
usage: scalpel-discovery <COMMAND> [OPTIONS]
COMMAND:
--help : this (help) message
--verbose : verbose mode
--single : single exome study
--denovo : family study (mom,dad,affected,sibling)
--somatic : normal/tumor study
> scalpel-discovery --single
$ scalpel-discovery --single
Local date and time: Sat Jul 7 10:14:20 2018
Program: scalpel-discovery (micro-assembly variant detection)
Version: 0.5.3 (beta), January 25 2016
Contact: Giuseppe Narzisi <gnarzisi@nygenome.org>
usage: scalpel-discovery --single --bam <BAM file> --bed <BED file> --ref <FASTA file> [OPTIONS]
Detect indels in one single dataset (e.g., one individual).
OPTIONS:
--help : this (help) message
--verbose : verbose mode
Required:
--bam <BAM file> : BAM file with the reference-aligned reads
--bed <BED file> : file with list of regions (BED format) in sorted order or single region in format chr:start-end (example: 1:31656613-31656883)
--ref <FASTA file> : reference genome in FASTA format (same one that was used to create the BAM file)
Optional:
--kmer <int> : k-mer size [default 25]
--covthr <int> : threshold used to select source and sink [default 5]
--lowcov <int> : threshold used to remove low-coverage nodes [default 2]
--covratio <float> : minimum coverage ratio for sequencing errors (default: 0.01)
--radius <int> : left and right extension (in base-pairs) [default 100]
--window <int> : window-size of the region to assemble (in base-pairs) [default 400]
--maxregcov <int> : maximum average coverage allowed per region [default 10000]
--step <int> : delta shift for the sliding window (in base-pairs) [default 100]
--mapscore <int> : minimum mapping quality for selecting reads to assemble [default 1]
--pathlimit <int> : limit number of sequence paths to [default 1000000]
--mismatches <int> : max number of mismatches in near-perfect repeat detection [default 3]
--dir <directory> : output directory [default ./outdir]
--numprocs <int> : number of parallel jobs (1 for no parallelization) [default 1]
--sample <string> : only process reads/fragments in sample [default ALL]
--coords <file> : file with list of selected locations to examine [default null]
Output:
--format : export mutations in selected format (annovar | vcf) [default vcf]
--intarget : export mutations only inside the target regions from the BED file
--logs : keep log files
Note 1: the list of detected indels is saved in file: OUTDIR/variants.indel.*
where OUTDIR is the output directory selected with option "--dir" [default ./outdir]
Note 2: use the export tool (scalpel-export) to export mutations using different filtering criteria
> scalpel-discovery --somatic
$ Local date and time: Sat Jul 7 10:14:53 2018
Program: scalpel-discovery (micro-assembly variant detection)
Version: 0.5.3 (beta), January 25 2016
Contact: Giuseppe Narzisi <gnarzisi@nygenome.org>
usage: scalpel-discovery --somatic --normal <BAM file> --tumor <BAM file> --bed <BED file> --ref <FASTA file> [OPTIONS]
Detect somatic indels in a tumor/normal pair
OPTIONS:
--help : this (help) message
--verbose : verbose mode
Required:
--normal <BAM file> : normal BAM file
--tumor <BAM file> : tumor BAM file
--bed <BED file> : file with list of regions (BED format) in sorted order or single region in format chr:start-end (example: 1:31656613-31656883)
--ref <FASTA file> : reference genome in FASTA format (same one that was used to create the BAM file)
Optional:
--kmer <int> : k-mer size [default 25]
--covthr <int> : threshold used to select source and sink [default 5]
--lowcov <int> : threshold used to remove low-coverage nodes [default 2]
--covratio <float> : minimum coverage ratio for sequencing errors (default: 0.01)
--radius <int> : left and right extension (in base-pairs) [default 100]
--window <int> : window-size of the region to assemble (in base-pairs) [default 400]
--maxregcov <int> : maximum average coverage allowed per region [default 10000]
--step <int> : delta shift for the sliding window (in base-pairs) [default 100]
--mapscore <int> : minimum mapping quality for selecting reads to assemble [default 1]
--pathlimit <int> : limit number of sequence paths to [default 1000000]
--mismatches <int> : max number of mismatches in near-perfect repeat detection [default 3]
--dir <directory> : output directory [default ./outdir]
--numprocs <int> : number of parallel jobs (1 for no parallelization) [default 1]
--coords <file> : file with list of selected coordinates to examine [default null]
--two-pass : perform second pass of analysis to confirm candidate calls
Output:
--format : export mutations in selected format (annovar | vcf) [default vcf]
--intarget : export mutations only inside the target regions from the BED file
--logs : keep log files
Note 1: the list of somatic indels is saved in file: OUTDIR/somatic.indel.*
where OUTDIR is the output directory selected with option "--dir" [default ./outdir].
Note 2: use the export tool (scalpel-export) to export mutations using different filtering criteria
> scalpel-discovery --denovo
$ Local date and time: Sat Jul 7 10:15:27 2018
Program: scalpel-discovery (micro-assembly variant detection)
Version: 0.5.3 (beta), January 25 2016
Contact: Giuseppe Narzisi <gnarzisi@nygenome.org>
usage: scalpel-discovery --denovo --dad <BAM file> --mom <BAM file> --aff <BAM file> --sib <BAM file> --bed <BED file> --ref <FASTA file> [OPTIONS]
Detect de novo indels in a family of four individuals (mom, dad, aff, sib).
OPTIONS:
--help : this (help) message
--verbose : verbose mode
Required:
--dad <BAM file> : father BAM file
--mom <BAM file> : mother BAM file
--aff <BAM file> : affected child BAM file
--sib <BAM file> : sibling BAM file
--bed <BED file> : file with list of regions (BED format) in sorted order or single region in format chr:start-end (example: 1:31656613-31656883)
--ref <FASTA file> : reference genome in FASTA format (same one that was used to create the BAM file)
Optional:
--kmer <int> : k-mer size [default 25]
--covthr <int> : threshold used to select source and sink [default 5]
--lowcov <int> : threshold used to remove low-coverage nodes [default 2]
--covratio <float> : minimum coverage ratio for sequencing errors (default: 0.01)
--radius <int> : left and right extension (in base-pairs) [default 100]
--window <int> : window-size of the region to assemble (in base-pairs) [default 400]
--maxregcov <int> : maximum average coverage allowed per region [default 10000]
--step <int> : delta shift for the sliding window (in base-pairs) [default 100]
--mapscore <int> : minimum mapping quality for selecting reads to assemble [default 1]
--pathlimit <int> : limit number of sequence paths to [default 1000000]
--mismatches <int> : max number of mismatches in near-perfect repeat detection [default 3]
--dir <directory> : output directory [default ./outdir]
--numprocs <int> : number of parallel jobs (1 for no parallelization) [default 1]
--coords <file> : file with list of selected coordinates to examine [default null]
--two-pass : perform second pass of analysis to confirm candidate calls
Output:
--format : export mutations in selected format (annovar | vcf) [default vcf]
--intarget : export mutations only inside the target regions from the BED file
--logs : keep log files
Note 1: the list of de novo indels is saved in file: OUTDIR/denovos.indel.*
where OUTDIR is the output directory selected with option "--dir" [default ./outdir]
Note 2: use the export tool (scalpel-export) to export mutations using different filtering criteria
ラン
1、シングルサンプル。
bedでエキソームキャプチャ領域を指定する。
scalpel-discovery --single --bam file.bam --bed regions.bed --ref genome.fa
デフォルトではoutdir/にvariants.indel.*というファイル名で出力される。scalpel-exportを使うと、そこから特定のタイプのバリアントのみ選抜できる。
scalpel-export --single --db --bed regions.bed --ref genome.fa > single.vcf outdir/variants.db.dir
2、 somatic変異。
tumor/normal比較から体細胞変異を検出する。正常組織由来のbamと腫瘍由来のbamを指定する。
scalpel-discovery --somatic --normal normal.bam --tumor tumor.bam --bed regions.bed --ref genome.fa --two-pass
#somaticをexport。1の例と同様、dbファイルを指定する。
scalpel-export --somatic --db DBfile --bed regions.bed --ref genome.fa> somatic.vcf
バリアント検出時、--two-passオプションをつけることで、ランタイムは増えるがfalse positiveがへりprecisionが改善する。
ビューアで確認する。
igv -g human_g1k_v37_decoy.fa normal.bam/tumor.bam, somatic.vcf
tumor(下半分)特異的な変異が出力されている。
3、de novo変異。
例えばdad、mom、sib、affectedのFamily quad解析を行いaffのde novo変異を検出する。
scalpel-discovery --denovo --dad dad.bam --mom mom.bam --aff aff.bam --sib sib.bam --bed regions.bed --ref genome.fa --two-pass
#de novoをexport。-1の例と同様、dbファイルを指定する。
scalpel-export --denovo --db DBfile --bed regions.bed --ref genome.fa
scalpel-discoveryが終わると、outdir/main/下に、dad、mam、aff、sibの4つのサブディレクトリができ、その各々にindels.vcfとdatabaseファイルが保存される。family trio解析に使う場合、--affと--sibに同じサンプルbamを指定する。他の pedigree typesについては今後対応予定らしい。
bed以外にchr:start-endで1つの領域を直接指定することで、関心のある領域だけ素早く解析することができるようになっている。
scalpel-discovery --single --bam file.bam --ref genome.fa --bed chr1:1000-4000
exportコマンドについては、以下のオプションを考慮するようアドバイスされている。
- mininum altertnative count:
- Change this number according to the type of study (germline, somatic, denovo) and the expected coverage.
- Smaller values will give more sensitivity but increase the number of false-positive calls.
- minimum variant allele frequency VAF (altCoverage/totalCoverage)
- Similarly to the minimum alternative count parameter, smaller values will increase sensitivity but reduce specificity.
- maximum chi-squared score:
- Chi square test statistic computed using the reference and alternative coverage for the mutation.
- Larger values will give more sensitivity but produce a large number of false-positives.
- For germline and denovo discovery we recommend using chi-square score ≤ 20 to select high confidence indels.
- minimum fisher exact test score:
- Fisher exact test statistic computed using the reference and alternative counts in tumor and normal samples.
- Goal is to test the independence between the allele balances in the tumor and the normal.
- We recommend using a fisher score > 10 to select high confidence somatic indels.
ヒトのwhole exomeをモデルに開発されているが、whole genomeにも使用できる。WGS に使用する時は、メモリ使用量の関係から染色体ごとに実施するよう推奨されている。またウィンドウサイズをdefaultの400から"--window 600"に変更することが提案されている(--window <int> : window-size of the region to assemble (in base-pairs) [default 400] )。 "--output-format annovar"をつけると、annnovarフォーマットで出力できる(デファルトはvcf)。
引用
Accurate detection of de novo and transmitted indels within exome-capture data using micro-assembly
Giuseppe Narzisi, Jason A. O'Rawe, Ivan Iossifov, Han Fang, Yoon-ha Lee, Zihua Wang, Yiyang Wu, Gholson J. Lyon, Michael Wigler, and Michael C. Schatz
Nat Methods. 2014 Oct; 11(10): 1033–1036. Published online 2014 Aug 17.
Reducing INDEL calling errors in whole-genome and exome sequencing data
Fang H, Wu Y, Narzisi G, O’Rawe JA, Jimenez Barron LT, Rosenbaum J, Ronemus M, Iossifov I, Schatz MC, Lyon GJ
Genome Medicine (2014) doi:10.1186/s13073-014-0089-z
Indel variant analysis of short-read sequencing data with Scalpel
Fang H, Grabowska EA, Arora K, Vacic V, Zody M, Iossifov I, O’Rawe JA, Y Wu, Jimenez-Barron LT, Rosenbaum J, Ronemus M, Lee Y, Wang Z, Dikoglu E, Jobanputra V, Lyon GJ, Wigler M, Schatz MC, Narzisi G*,
Preprint in bioRxiv (2016) doi: dx.doi.org/10.1101/028050