macでインフォマティクス

macでインフォマティクス

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

体細胞変異を検出する NeuSomatic

 

 体細胞変異はガンの発生、進行、および治療における重要なサインである。体細胞変異の正確な検出は、腫瘍とノーマルの交差汚染、腫瘍の異質性、シークエンシングアーティファクト、およびカバレッジのために困難である。一般に、前述の問題で発生する誤検出のコールを効果的にフィルタリングし、低い対立遺伝子頻度(AF)で発生する、または複雑さが低い領域で発生する可能性のコールを正確に維持することが正確な体細胞変異検出アルゴリズムには不可欠である。

 これまでに、MuTect2、MuSE、VarDict、VarScan2、Strelka2、SomaticSniperなど、体細胞変異検出問題に対処するさまざまなツールが開発されてきた。これらのツールは異なる統計的およびアルゴリズム的アプローチを採用しており、それはそれらが設計された特定のガンまたはサンプルタイプにおいてうまく機能する。しかしながら、それらは一般化においてより広い範囲のサンプルタイプおよびシーケンシングに限定されており、したがってそのようなシナリオにおいて最適以下の精度を示すことがある(ref.7,8,9)。本著者らの以前の仕事であるSomaticSeq(link)では、アルゴリズム的に直交する方法を統合することによって感度を最大にするためにアンサンブルアプローチを使用した。また、機械学習を使用して、精度を高く維持するために約100の機能を統合し、すべての個々の方法よりも精度が向上した。それにもかかわらず、SomaticSeqで使用されている機械学習のバックボーンは、変異の位置について抽出された一連の機能に依存している。その結果、体細胞変異のゲノムコンテキストで生の情報を完全に捉えて、真の体細胞変異をバックグラウンドエラーと区別することができず、複雑性の低い領域や低腫瘍純度などの困難な状況でのパフォーマンスが制限されていた。

 ここで著者らは、一般可能性の限界および腫瘍シーケンシングデータの統計的モデリングの複雑さを、ディープConvolutional Neural Networks (CNNs) を活用することによって扱う。 CNNは最近、生殖細胞系変異コーリング(ref.11-13)および皮膚ガン分類(ref.14) を含むさまざまな分野からの分類問題において著しい性能を示している。それでも、体細胞変異検出の困難な問題にCNNを適用することは検討されていない。以前の深層学習に基づく唯一の試み(ref.15)は、6層の完全に接続されたニューラルネットワークを手動で抽出された特徴セットに適用することであった。このアプローチは、ローカル領域で見られるパターンを使用して生のデータから直接特徴表現を学ぶことであるCNNアーキテクチャによって提供される力を欠いている。その上、完全に接続されたネットワークの複雑さのために、CNNで見られるようにそれはより一般化性とスケーラビリティを持たない。

 シーケンスアラインメントや体細胞変異を正確に同定するための他の方法から得られるシグナルを効果的に活用することができる体細胞変異検出のための最初のCNNベースのアプローチ、NeuSomaticを紹介する。生殖細胞系列のバリアントに焦点を当てている他のディープラーニングベースの方法とは異なり、NeuSomaticは腫瘍サンプルの複雑さのために精度の点でより大きな満たされていないニーズに取り組んでいる。生データから直接重要な変異シグナルを効果的に捉え、異なるシーケンシング技術、サンプル純度、全ゲノムvsターゲットエンリッチメントなどのシーケンシング戦略に対して一貫して高い精度を達成することができる。
 NeuSomaticネットワークへの入力は、一致した(matched)ノーマルサンプルと腫瘍サンプルのシーケンスアライメントをスキャンすることによって同定された候補体細胞変異である(論文図1;補足図1および2)。他の方法によって報告された体細胞変異もまたこの候補のリストに含めることができる。各候補遺伝子座について、その遺伝子座を中心とする領域からの信号を捕捉するために、それぞれサイズ5×32のk個のチャネルからなる(サイズk×5×32の)3次元特徴行列Mを構築する。マトリックスMの各チャネルは、4つのヌクレオチド塩基およびギャップのある塩基(「 - 」)を表す5行と、候補位置の周囲のアライメント列を表す32列とを有する。

 最初の3つのチャンネルは、それぞれ、候補遺伝子座周辺の基準塩基、およびその領域の異なる塩基の頻度を要約した基準、腫瘍頻度、および通常頻度の各チャンネルです。挿入を捕獲するために、リードアラインメントにおける挿入に対応するギャップを有する候補遺伝子座の周りの参照配列を増強する(図1a;補足図1および2)。したがって、腫瘍および正常頻度マトリックスの各列は、それぞれ腫瘍および正常サンプルの対応するマルティプルアラインメント(MSA)列におけるA / C / G / T /ギャップ塩基の頻度を表す。残りのチャンネルは、カバレッジ、ベースクオリティ、マッピングクオリティ、ストランドバイアス、そして異なるベースをサポートするリードのためのクリッピング情報のような他の機能を要約する。 NeuSomaticがアンサンブルモードで使用されている場合は、個々の体細胞変異検出方法によって報告された機能に追加のチャンネルも使用する。この簡潔で包括的な構造表現により、NeuSomaticは腫瘍、ノーマル、リファレンスの必要な情報を使用して、低AFで難しい生殖細胞系変異とシーケンシングエラーを区別して体細胞変異を捕獲できる。この設計はまた、CNNにおける畳み込みフィルタの使用を可能にして、行列のサブブロックにおける文脈パターンを捕捉する。(以下省略)

 

 

インストール

依存

Python 3.7 and the following Python packages must be installed:

  • pytorch 1.0.1
  • torchvision 0.2.1
  • pybedtools 0.8.0
  • pysam 0.15.2
  • zlib 1.2.11
  • numpy 1.15.4
  • scipy 1.2.0
  • imageio 2.5.0
  • biopython 1.73

It also depends on the following packages:

  • cudatoolkit 8.0 (if you want to use GPU)
  • tabix 0.2.6
  • bedtools 2.27.1
  • samtools 1.9

本体 Github

Docker

ここではdockerイメージ(CPU-only)を使いテストする。トレーニングする場合はGPUが必要になる。導入手順はGithubを参照。

docker pull msahraeian/neusomatic

> call.py -h

# call.py -h

usage: call.py [-h] --candidates_tsv [CANDIDATES_TSV [CANDIDATES_TSV ...]]

               --reference REFERENCE --out OUT --checkpoint CHECKPOINT

               [--ensemble] [--num_threads NUM_THREADS]

               [--batch_size BATCH_SIZE]

               [--max_load_candidates MAX_LOAD_CANDIDATES]

               [--pass_threshold PASS_THRESHOLD]

               [--lowqual_threshold LOWQUAL_THRESHOLD]

 

simple call variants from bam

 

optional arguments:

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

  --candidates_tsv [CANDIDATES_TSV [CANDIDATES_TSV ...]]

                        test candidate tsv files

  --reference REFERENCE

                        reference fasta filename

  --out OUT             output directory

  --checkpoint CHECKPOINT

                        network model checkpoint path

  --ensemble            Enable calling for ensemble mode

  --num_threads NUM_THREADS

                        number of threads

  --batch_size BATCH_SIZE

                        batch size

  --max_load_candidates MAX_LOAD_CANDIDATES

                        maximum candidates to load in memory

  --pass_threshold PASS_THRESHOLD

                        SCORE for PASS (PASS for score => pass_threshold)

  --lowqual_threshold LOWQUAL_THRESHOLD

                        SCORE for LowQual (PASS for lowqual_threshold <= score

                        < pass_threshold)

> preprocess.py  -h

# preprocess.py  -h

usage: preprocess.py [-h] --mode {train,call} --reference REFERENCE

                     --region_bed REGION_BED --tumor_bam TUMOR_BAM

                     --normal_bam NORMAL_BAM --work WORK

                     [--dbsnp_to_filter DBSNP_TO_FILTER]

                     [--scan_window_size SCAN_WINDOW_SIZE]

                     [--scan_maf SCAN_MAF] [--min_mapq MIN_MAPQ]

                     [--min_dp MIN_DP] [--max_dp MAX_DP] [--good_ao GOOD_AO]

                     [--min_ao MIN_AO] [--snp_min_af SNP_MIN_AF]

                     [--snp_min_bq SNP_MIN_BQ] [--snp_min_ao SNP_MIN_AO]

                     [--ins_min_af INS_MIN_AF] [--del_min_af DEL_MIN_AF]

                     [--del_merge_min_af DEL_MERGE_MIN_AF]

                     [--ins_merge_min_af INS_MERGE_MIN_AF] [--merge_r MERGE_R]

                     [--truth_vcf TRUTH_VCF] [--tsv_batch_size TSV_BATCH_SIZE]

                     [--matrix_width MATRIX_WIDTH]

                     [--matrix_base_pad MATRIX_BASE_PAD]

                     [--min_ev_frac_per_col MIN_EV_FRAC_PER_COL]

                     [--ensemble_tsv ENSEMBLE_TSV] [--long_read] [--restart]

                     [--first_do_without_qual] [--num_threads NUM_THREADS]

                     [--scan_alignments_binary SCAN_ALIGNMENTS_BINARY]

 

Preprocess input alignments for train/call

 

optional arguments:

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

  --mode {train,call}   train/call mode

  --reference REFERENCE

                        reference fasta filename

  --region_bed REGION_BED

                        region bed

  --tumor_bam TUMOR_BAM

                        tumor bam

  --normal_bam NORMAL_BAM

                        normal bam

  --work WORK           work directory

  --dbsnp_to_filter DBSNP_TO_FILTER

                        dbsnp vcf (will be used to filter candidate variants)

  --scan_window_size SCAN_WINDOW_SIZE

                        window size to scan the variants

  --scan_maf SCAN_MAF   minimum allele freq for scanning

  --min_mapq MIN_MAPQ   minimum mapping quality

  --min_dp MIN_DP       min depth

  --max_dp MAX_DP       max depth

  --good_ao GOOD_AO     good alternate count (ignores maf)

  --min_ao MIN_AO       min alternate count

  --snp_min_af SNP_MIN_AF

                        SNP min allele freq

  --snp_min_bq SNP_MIN_BQ

                        SNP min base quality

  --snp_min_ao SNP_MIN_AO

                        SNP min alternate count for low AF candidates

  --ins_min_af INS_MIN_AF

                        INS min allele freq

  --del_min_af DEL_MIN_AF

                        DEL min allele freq

  --del_merge_min_af DEL_MERGE_MIN_AF

                        min allele freq for merging DELs

  --ins_merge_min_af INS_MERGE_MIN_AF

                        min allele freq for merging INSs

  --merge_r MERGE_R     merge af ratio to the max af for merging adjacent

                        variants

  --truth_vcf TRUTH_VCF

                        truth vcf (required for train mode)

  --tsv_batch_size TSV_BATCH_SIZE

                        output files batch size

  --matrix_width MATRIX_WIDTH

                        target window width

  --matrix_base_pad MATRIX_BASE_PAD

                        number of bases to pad around the candidate variant

  --min_ev_frac_per_col MIN_EV_FRAC_PER_COL

                        minimum frac cov per column to keep columm

  --ensemble_tsv ENSEMBLE_TSV

                        Ensemble annotation tsv file (only for short read)

  --long_read           Enable long_read (high error-rate sequence) indel

                        realignment

  --restart             Restart the process. (instead of continuing from where

                        we left)

  --first_do_without_qual

                        Perform initial scan without calculating the quality

                        stats

  --num_threads NUM_THREADS

                        number of threads

  --scan_alignments_binary SCAN_ALIGNMENTS_BINARY

                        binary for scanning alignment bam

root@e8719675e96e:/# 

 

 

実行方法

1、Preprocess step

python preprocess.py \
--mode call \
--reference GRCh38.fa \
--region_bed region.bed \
--tumor_bam tumor.bam \
--normal_bam normal.bam \
--work work_call \
--min_mapq 10 \
--number_threads 10 \
--scan_alignments_binary ../bin/scan_alignments

 

 2、Call variants

python call.py \
--candidates_tsv work_call/dataset/*/candidates*.tsv \
--reference GRCh38.fa \
--out work_call \
--checkpoint work_train/some_checkpoint.pth \
--num_threads 10 \
--batch_size 100

 

 3、Postprocess step

python postprocess.py \
--reference GRCh38.fa \
--tumor_bam tumor.bam \
--pred_vcf work_call/pred.vcf \
--candidates_vcf work_call/work_tumor/filtered_candidates.vcf \
--output_vcf work_call/NeuSomatic.vcf \
--work work_call

 

 

引用

Deep convolutional neural networks for accurate somatic mutation detection
Sayed Mohammad Ebrahim Sahraeian, Ruolin Liu, Bayo Lau, Karl Podesta, Marghoob Mohiyuddin, Hugo Y. K. Lam
Nature Communications 10, Article number: 1041 (2019)