macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

様々な構造変化を検出する TIDDIT

2021 6/6 インストール追記

 

 ゲノム構造変異(SV)は大きなゲノムの再編成と定義され、逆位、転座、ならびに欠失および重複からなる(preprintより ref.1)。SVは、多くの異なるヒト遺伝的障害における直接的原因および寄与因子の両方であることが示されており、また、より一般的な表現型形質(ref.2-4)にも同様にも関係している。

遺伝子診断では、FISH(ref.5)やマイクロアレイ研究(ref.6)などの現在の技術では分解能が限られており、得られた情報は正しい解釈のための追加の方法で補完する必要がある場合が多い(ref.5,6)。

イルミナHiSeqXのようなシステムの出現により、研究者は、全ヒトゲノム(WHG)を高い(すなわち30倍)カバレッジで比較的低コスト(すなわち、WHGあたり$ 1000 )でシーケンスすることが可能になる(ref.8)。前例のないスピードで大量のデータが出ることで、構造変異および/または染色体切断点(すなわち、SVが起こる染色体中の正確な位置)を同定(コール)することができる計算ツールの繁栄を促した。これらのツールは、一般にバリアントコーラーと呼ばれる。バリアント者コーラーは、一般に、(BAM形式の)アラインメントファイルを入力として要求し、リファレンスゲノムとドナー/患者ゲノムとの間の差異を特定しようと試みる。構造的バリアントを検出するために、コーラーは一般にWGSデータ内の異なる信号に基づくヒューリスティックを使用する。これらの信号には、discordant pairs(ref.9,10)、 read-depth(ref.11)、およびsplit-reads(ref.12)が含まれる。一部のコーラーは、ローカル(ref.13)またはグローバルな(ref.14)デノボアセンブリ技術のいずれかを適用することによって患者シーケンスを再構成しようとする。配列データのサイズ、バリアントタイプ、および特性に応じて、バリアントを検出する最も適切な方法は異なってくる(ref.15)。

 比較的低コストで高品質のシーケンシングデータを生成する能力と、単一の実験から任意の変異を検出する可能性のおかげで、エクソームシーケンシング(ref.16,17)および全ゲノムシーケンシング(ref.7,18)は、臨床診断において非常に有用であり得る。特に変異を引き起こす希少疾患を研究することが期待される。しかし、高い検証コストを避けるためには、高感度だが、検出精度も高いコーラーがが必要となる。豊富なシーケンス決定プラットフォームおよびライブラリー調製方法が利用可能である。これらの異なる供給源から生成される配列決定データは、リード長およびカバレージなど異なる特性を有する(ref.20)。(一部略)

 ここでは、新しいバリアントコーラーTIDDITを示す。この名前は、さまざまな種類のSVを検出する機能を強調しており、転座、逆位、欠失、散在重複、挿入およびタンデム重複が挙げられるが、これらに限定されない。TIDDITは、構造変異のゲノム上での位置、変異の分類および品質評価のためのリードのデプス情報を検出するために、discordant pairsおよびsplit-readsを利用する。これらのWGS信号を統合することにより、TIDDITは平衡型と不平衡型の両方を検出することができる。最後に、TIDDITは、インサートサイズとペアの向きが異なる複数のペアエンドシーケンシングライブラリタイプをサポートしている。

 まれな疾患の原因となる変異の検索を簡略化するため、TIDDITはSVDB(Structural Variant DataBase)と呼ばれるデータベース機能とともに配布される。 SVDBは構造変異の頻度データベース作成に使用される。このデータベースは、変異を引き起こす珍しい病気、および特定の継承パターンに従った変異のために照会される。データベースの機能性を利用して珍しい変異の解析を優先させることで、変異を引き起こす希少疾患の診断を迅速化することができる。著者らの知る限りでは、構造変化を引き起こす珍しい病気を呼び出して評価するための枠組みを提供しているようなコーラーは存在しない。

 

 

インストール

cent OSに導入した。

依存

  • c++/c libraries,
  • python 2.7
  • Numpy
  • scipy
  • cmake

Github

https://github.com/SciLifeLab/TIDDIT

#conda (link)
mamba install -c bioconda tiddit

#from source
git clone https://github.com/SciLifeLab/TIDDIT.git
cd TIDDIT
#TIDDITのビルド (numpy cythonもpipで導入される)
./INSTALL.sh

> python TIDDIT.py --sv -h

$ python TIDDIT.py --sv -h

usage: TIDDIT --sv --bam inputfile [-o prefix] --ref ref.fasta

       [-h] [--sv] --bam BAM [-o O] [-i I] [-d D] [-p P] [-r R] [-q Q] [-Q Q]

       [-n N] [-m M] [-e E] [-l L] [-z Z] [--force_ploidy] [--n_mask N_MASK]

       --ref REF

 

optional arguments:

  -h, --help       show this help message and exit

  --sv             call structural variation

  --bam BAM        coordinate sorted bam file(required)

  -o O             output prefix(default=output)

  -i I             paired reads maximum allowed insert size. Pairs aligning on

                   the same chr at a distance higher than this are considered

                   candidates for SV (default=3std + mean_insert_size)

  -d D             expected reads orientations, possible values "innie" (->

                   <-) or "outtie" (<- ->). Default: major orientation within

                   the dataset

  -p P             Minimum number of supporting pairs in order to call a

                   variation event (default 5)

  -r R             Minimum number of supporting split reads to call a small

                   variant (default 5)

  -q Q             Minimum mapping quality to consider an alignment (default

                   10)

  -Q Q             Minimum regional mapping quality (default 20)

  -n N             the ploidy of the organism,(default = 2)

  -m M             minimum variant size,(default = 100)

  -e E             clustering distance parameter, discordant pairs closer than

                   this distance are considered to belong to the same

                   variant(default = sqrt(insert-size*2)*12)

  -l L             min-pts parameter (default=3),must be set > 1

  -z Z             minimum variant size (default=100)

  --force_ploidy   force the ploidy to be set to -n across the entire genome

                   (i.e skip normalisation based on the -s list of

                   chromosomes)

  --n_mask N_MASK  exclude regions from coverage calculation if they contain

                   more than this fraction of N (default = 0.5)

  --ref REF        reference fasta

python TIDDIT.py --cov -h

$ python TIDDIT.py --cov -h

usage: TIDDIT --cov --bam inputfile [-o prefix] [-h] [--cov] --bam BAM [-o O]

                                                [-z Z]

 

optional arguments:

  -h, --help  show this help message and exit

  --cov       generate a coverage bed file

  --bam BAM   coordinate sorted bam file(required)

  -o O        output prefix(default=output)

  -z Z        use bins of specified size(default = 500bp) to measure the

              coverage of the entire bam file, set output to stdout to print

              to stdout

 

 

ラン

 bamとreferenceを指定してランする。

python TIDDIT.py --sv --bam input.bam --ref ref.fa

 4つのファイルが出力される。

1、バリアント情報を記載したvcfファイル。intrachromosomal variantsと interchromosomal variantsに分けて出力される。

2、100-bpのウィンドウサイズでゲノム全体のカバレッジを示したファイル。

3、カバレッジ情報とそこから推測した各クロモソームのploidyを示したファイル。

出力はVCFフォーマットに準拠している。 INFOフィールドについてはGIthubの説明を参照。

4、データベースファイル(binary)

VCFファイル。

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Z

chr     1       SV_1_1  N       ]chr:3571103]N  .       PASS    SVTYPE=BND;CIPOS=0,119;CIEND=-148,0;COVM=31.9323711771;COVA=27.853062701;COVB=34.8007151875;LFA=37;LFB=37;LTE=37;E1=0;E2=0;OR=0,0,16,0;ORSR=0,21;QUALA=60;QUALB=60      GT:CN:DV:RV     1/1:.:16:21

chr     386408  SV_2_1  N       <DUP>   .       MinSize SVTYPE=DUP;CIPOS=-1,0;CIEND=0,1;END=386411;SVLEN=4;COVM=63.6749484193;COVA=63.6749484193;COVB=63.6749484193;LFA=32;LFB=32;LTE=16;E1=31.8374742097;E2=31.8374742097;OR=0,0,0,0;ORSR=0,16;QUALA=60;QUALB=60       GT:CN:DV:RV     0/1:4:0:16

chr     852461  SV_3_1  N       ]chr:1204006]N  .       PASS    SVTYPE=BND;CIPOS=0,147;CIEND=-289,0;COVM=32.0590581331;COVA=69.3346325494;COVB=36.146842485;LFA=44;LFB=68;LTE=44;E1=1;E2=0;OR=0,0,12,1;ORSR=0,31;QUALA=47;QUALB=54      GT:CN:DV:RV     1/1:.:13:31

chr     853371  SV_4_1  N       N[chr:1203881[  .       FewLinks        SVTYPE=BND;CIPOS=-1,0;CIEND=0,32;COVM=31.9775727391;COVA=33.1334091555;COVB=37.4980640557;LFA=21;LFB=66;LTE=21;E1=0;E2=0;OR=0,0,0,6;ORSR=0,15;QUALA=50;QUALB=51 GT:CN:DV:RV     1/1:.:6:15

 

 

カバレッジの分布だけ調べたいなら、TIDDIT.pyも使える。

 2、500-bpのウィンドウサイズでゲノム全体のカバレッジを計算。

python TIDDIT.py --cov -z 500 --bam input.bam -o output
  • -o   the prefix of the output files
  • -z    compute the coverage within bins of a specified size across the entire genome, default bin size is 500

output.tabが出力される。

> head output.tab

#CHR    start   end     coverage        quality

chr     0       500     24.748  60

chr     500     1000    24.296  60

chr     1000    1500    28.628  60

chr     1500    2000    27.13   60

chr     2000    2500    33.76   60

chr     2500    3000    31.386  60

chr     3000    3500    29.052  60

chr     3500    4000    25.58   60

chr     4000    4500    30.846  60

 

filterをPASSした変異には#FILTERのカラムにPASSがつきます。PASSの定義はGithubに書いてありますが、PASSした変異だけ抽出することで、検出の精度が大きく向上するとのことです。

 

 

引用

TIDDIT, an efficient and comprehensive structural variant caller for massive parallel sequencing data

Version 2. F1000Res. 2017; 6: 664.

Published online 2017 Jun 30. doi: 10.12688/f1000research.11168.2