macでインフォマティクス

macでインフォマティクス

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

バイオインフォマティクスデータの検索と分析サイト Datasets2Tools

 

 ウェブの導入により、研究成果の伝統的な印刷出版物のソフトウェアベースの拡張が可能になった。a)研究論文をより簡単にコピー及び配布できるよりソフトウエアベースでpublishできる。 b)研究によって収集されたデータは、再利用および統合的で遡った(retrospective )分析のために、他の人がより容易に利用できるようにすることができる。 c)研究データを分析し視覚化するための方法とソフトウェアツールは、より強力かつアクセスしやすくなった。ウェブの導入によって促進されたこれらの発展は、生物医学研究の進歩にとって重要であったが、研究結果がどのようにデジタルで伝達されるかを改善するための多くの作業が残されている。ウェブの導入によりバイオインフォマティクス分野も促進されており、毎年何百もの新しいツールやデータベースが公開されている。同時に、biological sampleから大量の情報を抽出するhigh-content のバイオテクノロジーの進歩によるデータ収集の急速な成長は、この豊富なデータを管理するための集約型生物医学レポジトリの開発につながった。例えば、遺伝子発現トランスクリプトミクスとエピゲノミクスデータのためのOmnibus(GEO)(論文より ref.1)がある。実際、生物医学および生物学研究のすべての分野をカバーする多くのリポジトリが現在存在している。これらのリポジトリは、DataMed(ref.2)やFAIRSharing(ref.3)などのサイトで検索を行うよう索引付けされている。同様に、OMICtools(ref.4)などのバイオインフォマティクスツールを集約して索引付けするためのいくつかの取り組みが行われている。

 集約型の生物医学データリポジトリのデータセットランディングページから開始してpublication-ready な図を生成したり、インタラクティブな視覚化グラフィックを生成する方法があるが、現在ほとんどの集約型生物医学リポジトリは、その中に含まれるデータセットを操作できるツールから切り離されている。バイオインフォマティクス分析のプロセスを合理化する努力がなされているが、最近まで、データセットを処理するために適用されたほとんどのバイオインフォマティクスツールは、分析を実行したユーザーには引き続き他のユーザーと共有できる永続URLを提供しない。しかし、特定の集約型データリポジトリに含まれるデータセットから生成され、特定のバイオインフォマティクスツールによって分析されるような、事前に生成された ‘缶詰型’ の分析結果を表示するウェブページは徐々に蓄積されている。 GeneMANIA(ref.5)、Enrichr(ref.6,7)、Clustergrammer(ref.8)、L1000CDS2(ref.9)などのソフトウェアツールは、分析のために遺伝子シグネチャーを提出するユーザー、または相互作用ネットワークを構築するユーザーに、これらのツール、他者と共有できる永続URLを提供する。これに加えて、Galaxy(ref.10)のようなプラットフォームは、Web上で共有できる複雑なバイオインフォマティクスワークフローやパイプラインを簡単に作成できるインターフェイスを提供する。同様に、JupyterやR Markdown(ref.11)のようなインタラクティブなノートブックは、バイオインフォマティクスの分析をウェブ上で体系的に共有するためのプラットフォームである。再現性のある生物医学的および生物学的研究結果を作成し、伝達し、共有するこれらの新しい永続的方法は、我々が「缶詰め分析( ‘canned analysis’)」と呼ぶ新しい形の出版物と考えることができる。

 ここでは、初期状態で31,473のcanned analysisを揃えたcanned analysisのバイオインフォマティクス分析プラットフォームDatasets2Toolsを紹介する。31,473のcanned analysisは、選択されたツールを使用して6,431のデータセットに適用され、永続的なcanned analysisのURLが提供される。Datasets2Toolsは、4つのリーディングバイオインフォマティクスジャーナルBMC BioinformaticsBioinformaticsDatabase、そして Nucleic Acids Researchに掲載されている4901のソフトウェアツールのインデックスも作成している。Datasets2Toolsを使用すると、関連するデータセット、ツール、およびcanned analysisのデジタルオブジェクトをすばやく見つけることができる。著者らは、Datasets2Toolsがコミュニティへの貢献を通じて拡大することを期待している。 Datasets2Toolsは、Google Chromeブラウザの拡張機能Application Programming InterfaceAPI)としても提供されている(リンク)。

 Datasets2Toolsのユニークな機能の1つは、検索でき、アクセスでき、相互運用でき、再利用可能な(FAIR)原則に準拠する、データセット、ツール、およびcanned analysisを評価する機能である(ref.12)。 DataSet2Toolsに含まれる3種類のデジタルオブジェクト(データセット、ツール、およびcanned analysis)の検出可能性、アクセシビリティ、相互運用性、および再利用性に関するさまざまな側面に関する9つのYesまたはNoの質問に答えることができる。評価はデータベースに保存され、各データセット、ツール、または缶詰分析の近くに記章として表示される。 Datasets2Tools内の索引付けされたすべてのバイオインフォマティクスツールは出版物に由来しているため、Altmetric Attention ScoresおよびPlumXアセスメントならびに各ツールのPubMed引用はツールのカードの近くに表示される。このツール・メトリックの可視化は、ユーザーがバイオインフォマティクス・ツールをよりよく識別し、ランク付けするのに役立つ。 Datasets2Toolsのコンポーネントと機能の概要を論文図1に示す(pubmed)。

 

ヘルプページ

http://amp.pharm.mssm.edu/datasets2tools/help

 

 

概要

 Datasets2Toolsは 30,000以上のバイオインフォマティクス解析と、6,000以上のデータセット、4,000以上の計算ツールの索引付けを行っているデジタルオブジェクトの発見と評価のためのプラットフォーム。 これらには、豊富なエンリッチメント分析、遺伝子相互作用ネットワーク、インタラクティブなデータの可視化、 マイクロアレイ、RNA-seq、プロテオミクスなどのデータセットの計算ツールを提供している。 ユーザーは、webページのインタラクティブな検索インターフェイスを使用してこれらのオブジェクトを検索できる。缶詰分析というのは、バイオインフォマティクス分析の結果を、見つけやすく、アクセス可能で、相互運用可能で、再利用可能な方法で保存するように設計された新しいタイプのデジタルオブジェクトと定義されている(詳細な定義は論文とhelp参照)。

 

使用方法 

Datasets2Tools

f:id:kazumaxneo:20180727084219p:plain

 

GEO accesion number(GEOはNGSやアレイのレポジトリで分析もできる)による検索。 

f:id:kazumaxneo:20180731211906j:plain

ツール検索

http://amp.pharm.mssm.edu/datasets2tools/search?object_type=tool

f:id:kazumaxneo:20180727084109p:plain

 

例えば、NCBI GEOのこのデータ(ID: GSE39509)を調べたいとする(マウスのtumorのRNA seq)。

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39509

f:id:kazumaxneo:20180731212056j:plain

IDで検索。

f:id:kazumaxneo:20180731212229j:plain

hitした。クリックする。

f:id:kazumaxneo:20180731212249j:plain

詳細が表示される。

f:id:kazumaxneo:20180731212331j:plain

Toolsタブではこのデータの分析に使われたツールが表示される。そのツールに飛んで、詳細を調べたり、他にこのツールが使われた解析にさらにリンクすることもできる。

f:id:kazumaxneo:20180731212412j:plain

Viewを押すと、GSE39509のデータから作られたヒートマップに飛ぶ。

f:id:kazumaxneo:20180731212758j:plain

マウスのホイールで拡大縮小できる。これだけでも論文よりずっと便利(論文のヒートマップは格好だけはいいが、拡大しても字が潰れて見えなかったり、データの価値が損なわれていることがある)。

f:id:kazumaxneo:20180731212826j:plain

左のメニューから、縦方向と列方向のソート方法を変えたり、ヒートマップの濃さを調整できる。また、ブラウザで表示できるSVG形式、または行列のテキストファイルとしてデータをダウンロードできる。

まずは興味あるキーワードでラフに検索してみて下さい。そこからさらに絞り込めます。

f:id:kazumaxneo:20180731214615j:plain

 

 

ユーザーがデータをアップロードすることも可能。

http://amp.pharm.mssm.edu/datasets2tools/contribute

f:id:kazumaxneo:20180727083836p:plain

 

chrome extensionとしても供給されている(ダイレクトリンクできる)。

f:id:kazumaxneo:20180727083752p:plain

 

おまけ

Ma'ayan LabのGithubレポジトリには他にもいくつかのツールが登録されている。そのうちの1つBioToolBayは、Datasets2Toolsと同じく、4つの主要バイオインフォマティクスジャーナルにアクセスされた論文のツールを検索できるHPらしい。

http://amp.pharm.mssm.edu/biotoolbay

f:id:kazumaxneo:20180731143431j:plain

ユニークリソース、可視化ツール、データベースなど、などで分類されている。

例えば可視化ツール一覧。まだ下にスクロールできる。

f:id:kazumaxneo:20180731143644j:plain

ツールをクリックすれば詳細を確認できる。

f:id:kazumaxneo:20180731144004j:plain

 

 

 感じたこと

 バイオインフォマティクスに限らないが、論文の図表で提示されたデータの多くは、通常、パースしたりフィルタリングされ、その上、場合によっては他のデータと比較可能なように正規化された、編集済みデータである(ただし生のデータが出される研究もたくさんある)。このようなデータセットは、著者らが適切に解釈できると判断した手順で処理されており、読み手が論文内で主張されていることを理解するには必要十分になっている。しかし、角度を変えて分析したい時には少し困る。このブログで書くと釈迦に説教になってしまうが、一般に、競合する研究者やその研究成果を基礎にしてより新しい研究を行いたい研究者がデータドリブンで新しい現象を見つけ出そうとする時、往往にして、false positiveのリスクを取ってでも最大限の情報を抽出したいと思うものである(特にomic関係でよくある。統計のcut-off valueは重要な役目を持つが、値自体に意味は薄いので、cut-offの向こうには大事な遺伝子が潜んでいる可能性がある)。具体的には、論文とともに公共レポジトリにdepositされた生データをダウンロードし(しばしば無い場合すらある)、寸暇を惜しんで自ら同じようなセットセットを作り上げたりする。これには、意外にバカにならないくらいの時間がかかったりする。しかし、たとえその論文のmethodsを注意深く読みながらデータの再解析を実行したとしても、環境の違いで再現できないことがしばしばある(理由は無数に考えられる)。バージョンを統一してもできない場合があり、それが意図的なものではないかと疑義が立つことすらあるが、それは別の問題なので、触れるだけにとどめる。

 このような問題は、オリジナルデータは一切編集できないが、見かけ上だけフィルタリングやパース条件を変えられるインタラクティブなデータベースサイトに図表データがアップロードされていたり、depositされたデータを瞬時に可視化できるデータベースがあれば低減できる(NCBI GEO自体もそうです)。Datasets2Toolsは、公開されたデータを素早く、ユーザーが理解できる形で分析、再利用できることを目標に構築されている(FAIR Principles )。

 

紛らわしい文章を書きましたが、Datasets2ToolsはNCBI GEOにアクセッションされたデータから図を作っています。公開されているデータの中に、論文の著者らが直接アップしたものがどれだけあるのかは不明です。

 

引用

Datasets2Tools, repository and search engine for bioinformatics datasets, tools and canned analyses

Torre D, Krawczuk P, Jagodnik KM, Lachmann A, Wang Z, Wang L, Kuleshov MV, Ma'ayan A

Sci Data. 2018 Feb 27;5:180023.

高速なgermlineとsomaticのSV検出ツール Manta

 

 ゲノムシーケンシングおよびゲノムエンリッチメントシーケンシングは、臨床での遺伝性および体細胞突然変異発見のためにますます使用されてきているが、このシナリオにおける構造変異(SV)およびindelsの迅速な発見のためのツールは限られている。著者らは、SVと中サイズのindelsと大きなサイズの挿入を統一された迅速なプロセスで正確に検出し、評価するための新しい方法であるMantaを使いこのギャップに取り組む。 Mantaは、シークエンシングアッセイのペアエンドとスプリットマッピングから効率的なパラレルワークフローでバリアントを検出する。現在、研究とpopulation genomicsに焦点を当てた多くの高度な構造変異検出ツールが利用可能である(Layer et al、2014; Rausch et al、2012; Sindi et al、2012; Ye et al、2009)。しかし、著者らの知る限り、関連するサンプルの個々のセットまたは小さなセットに焦点を当てた迅速なワークフローに、多くのバリアントタイプを組み合わせることはできない。臨床パイプラインに重点を置くMantaは、リファレンスゲノムと任意の標準リードマッパーからのアラインメント(bam)のみを使用して、検出、アセンブリ、スコアリングのための完全なソリューションを提供する。これは、二倍体個体の生殖系列分析および腫瘍 - ノーマルサンプルペアの体細胞分析のためのスコアリングモデルを提供し、現在開発中のRNA-Seq、de novo変異、および他に類を見ない腫瘍のためのさらなる応用を提供する。(一部略)。

 Mantaのワークフローは、個々のサンプルまたは小さなサンプルセットで高い並列性を実現するように設計されている。それは2つのフェーズで動作する。(一部略)すべてのワークフローコンポーネントの詳細は補足資料(ダイレクトリンク)のメソッドで説明されている。

 

bamからSVが予測された領域についてgraphを作成し、その graphからバリアントを予測する。正確な予測ができなかった場合、ペアエンド情報だけから予測され、IMPRECISEのタグがつく。

f:id:kazumaxneo:20180730132404j:plain

Mantaのワークフロー。論文補足資料より転載(ダイレクトリンク)。

 

著者らの環境では、よく使われるゴールデンスタンダードindelのplatinum genomes NA12878(家系図 CEPE pedigree)(これ以外にもGenome in a bottle (GIAB) のデータがある)の50xシーケンシングデータ(illumina)のbamを、Mantaは20分未満で処理できるとされる(環境: 20 physical cores using a dual Xeon E5-2680 v2 server with the BAM accessed from a conventional local drive, peak total memory (RSS) for this run was 2.35 Gb)。 

ハードウエア必要条件

  • Memory Typical memory requirements are <1Gb/core for germline analysis and <2Gb/core for cancer/FFPE/highly rearranged samples. The exact requirement depends on many factors including sequencing depth, read length, fragment size and sample quality.
  • CPU Manta does not require or benefit from any specific modern CPU feature (e.g. NUMA, AVX..), but in general faster clock and larger caches will improve performance.
  • I/O I/O can be roughly approximated as 1.1 reads of the input alignment file per analysis, with no writes that are significant relative to the alignment file size.

MantaはWGSとWES向けのツールで、以下のような解析に対応している。

  • Joint analysis of small sets of diploid individuals (where 'small' means family-scale -- roughly 10 or fewer samples)
  • Subtractive analysis of a matched tumor/normal sample pair
  • Analysis of an individual tumor sample 

Mantaは以下のようなSVサブクラス検出に対応している。

  • Deletions
  • Insertions
  • Fully-assembled insertions
  • Partially-assembled (ie. inferred) insertions
  • Inversions
  • Tandem Duplications
  • Interchromosomal Translocations

Mantaは以下のようなSVサブクラスは検出できない。

  • Dispersed duplications
  • Most expansion/contraction variants of a reference tandem repeat
  • Small inversions
  • The limiting size is not tested, but in theory detection falls off below ~200bases. So-called micro-inversions might be detected indirectly as combined insertion/deletion variants.
  • Fully-assembled large insertions
  • The maximum fully-assembled insertion size should correspond to approximately twice the read-pair fragment size, but note that power to fully assemble the insertion should fall off to impractical levels before this size
  • Note that manta does detect and report very large insertions when the breakend signature of such an event is found, even though the inserted sequence cannot be fully assembled.

マニュアル

manta/docs/userGuide at master · Illumina/manta · GitHub

 Mantaに関するツイート。


8/3 dockerコマンド修正

 

インストール

mac os10.13でdokcerを使いテストした。

本体 Github

リリースからソースコードのほか、cent os向けバイナリをダウンロードできる。

https://github.com/Illumina/manta/releases

linuxならcondaでも導入できる(リンク)。

ここではdockerコンテナを使う。

docker pull eitanbanks/manta-1.0.3

 > docker run --rm -it eitanbanks/manta-1.0.3 /bin/manta/bin/configManta.py -h

$  docker run --rm -i -t -v /Users/user/data/:/data eitanbanks/manta-1.0.3 /bin/manta/bin/configManta.py -h

Usage: configManta.py [options]

 

Version: 1.0.3

 

This script configures the Manta SV analysis pipeline.

You must specify a BAM or CRAM file for at least one sample.

 

Configuration will produce a workflow run script which

can execute the workflow on a single node or through

sge and resume any interrupted execution.

 

Options:

  --version             show program's version number and exit

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

  --config=FILE         provide a configuration file to override defaults in

                        global config file

                        (/usr/bin/manta/bin/configManta.py.ini)

  --allHelp             show all extended/hidden options

 

  Workflow options:

    --bam=FILE, --normalBam=FILE

                        Normal sample BAM or CRAM file. May be specified more

                        than once, multiple inputs will be treated as each BAM

                        file representing a different sample. [optional] (no

                        default)

    --tumorBam=FILE, --tumourBam=FILE

                        Tumor sample BAM or CRAM file. Only up to one tumor

                        bam file accepted. [optional] (no default)

    --exome             Set options for WES input: turn off depth filters

    --rna               Set options for RNA-Seq input: turn off depth filters

                        and don't treat anomalous reads as SV evidence when

                        the proper-pair bit is set.

    --unstrandedRNA     Set if RNA-Seq input is unstranded: Allows splice-

                        junctions on either strand

    --referenceFasta=FILE

                        samtools-indexed reference fasta file [required]

    --runDir=DIR        Run script and run output will be written to this

                        directory [required] (default: MantaWorkflow)

 

  Extended options:

    These options are either unlikely to be reset after initial site

    configuration or only of interest for workflow development/debugging.

    They will not be printed here if a default exists unless --allHelp is

    specified

 

    --scanSizeMb=INT    Maximum sequence region size (in megabases) scanned by

                        each task during SV Locus graph generation. (default:

                        12)

    --region=REGION     Limit the analysis to a region of the genome for

                        debugging purposes. If this argument is provided

                        multiple times all specified regions will be analyzed

                        together. All regions must be non-overlapping to get a

                        meaningful result. Examples: '--region chr20' (whole

                        chromosome), '--region chr2:100-2000 --region

                        chr3:2500-3000' (two regions)'

——

us

dokcerイメージから表示するなら

docker run --rm -it  eitanbanks/manta-1.0.3 /bin/manta/bin/configManta.py -h

 

>  runWorkflow.py -h #上記のコマンドを打つと、このpyflowのスクリプトができる。

Usage: runWorkflow.py [options]

 

Version: 1.4.0

 

Options:

  --version             show program's version number and exit

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

  -m MODE, --mode=MODE  select run mode (local|sge)

  -q QUEUE, --queue=QUEUE

                        specify scheduler queue name

  -j JOBS, --jobs=JOBS  number of jobs, must be an integer or 'unlimited'

                        (default: Estimate total cores on this node for local

                        mode, 128 for sge mode)

  -g MEMGB, --memGb=MEMGB

                        gigabytes of memory available to run workflow -- only

                        meaningful in local mode, must be an integer (default:

                        Estimate the total memory for this node for local

                        mode, 'unlimited' for sge mode)

  -d, --dryRun          dryRun workflow code without actually running command-

                        tasks

  --quiet               Don't write any log output to stderr (but still write

                        to workspace/pyflow.data/logs/pyflow_log.txt)

 

  development debug options:

    --rescore           Reset task list to re-run hypothesis generation and

                        scoring without resetting graph generation.

 

  extended portability options (should not be needed by most users):

    --maxTaskRuntime=hh:mm:ss

                        Specify scheduler max runtime per task, argument is

                        provided to the 'h_rt' resource limit if using SGE (no

                        default)

 

Note this script can be re-run to continue the workflow run in case of

interruption. Also note that dryRun option has limited utility when

task definition depends on upstream task results -- in this case the

dry run will not cover the full 'live' run task set.

 

 

ラン

1、はじめにマッピングしてbamを作成する。著者らはマッパーにBWA-MEM version 0.7.5a を使っている。

 

2、Configurationファイルの作成。bin/configManta.pyを使う。入力はマッピングして得たbam(cram)とそのリファレンスfasta。bam(cram)、fasta共にindexファイルも必要。

Single Diploid Sample Analysis

#~/documnets/input_file/とイメージのdata/を共有ディレクトリとする。ただしdockerを使わないなら1行目は不要。またファイルパスの/data/部分も不要。
docker run --rm -itv /Users/uesaka/Documents/input_file:/data eitanbanks/manta-1.0.3 \
/bin/manta/bin/configManta.py \
--bam=/data/NA12878_S1.bam \
--referenceFasta=/data/hg19.fa \
--runDir=/data/output

 Joint Diploid Sample Analysis(例えばCEPE familyのコホート(ここではTrio))

docker run --rm -itv /Users/uesaka/Documents/input_file:/data eitanbanks/manta-1.0.3 \
/bin/manta/bin/configManta.py \
--bam=/data/NA12878_S1.cram \
--bam=/data/NA12891_S1.cram \
--bam=/data/NA12892_S1.cram \
--referenceFasta /data/hg19.fa \
--runDir=/data/output

 Tumor Normal Analysis(case-contorol)

docker run --rm -itv /Users/uesaka/Documents/input_file:/data eitanbanks/manta-1.0.3 \
/bin/manta/bin/configManta.py \
--normalBam=/data/HCC1187BL.cram \
--tumorBam=/data/HCC1187C.cram \
--referenceFasta=/data/hg19.fa \
--runDir=/data/output

 Tumor-Only Analysis

docker run --rm -itv /Users/uesaka/Documents/input_file:/data eitanbanks/manta-1.0.3 \
/bin/manta/bin/configManta.py \
--tumorBam=/data/HCC1187C.cram \
--referenceFasta=/data/hg19.fa \
--runDir=/data/output

 

3、 Configurationファイルの実行(シングルノードでの実行。SGE clusterのランはGithub参照)(*1)。

#dockerを使わないなら1行目は不要
docker run --rm -itv /Users/uesaka/Documents/input_file:/data eitanbanks/manta-1.0.3 \
data/output/runWorkflow.py -m local -j 8 --memGb=20

ジョブが終わると、germline解析では3つのVCFファイルが出力される(results/variants)。turmor/nomal解析ではさらにもう1つVCFができ(tumorSV.vcf)、合計4つのVCFが出力される。VCFのフォーマットはversion4.1に則っている。ポジションはSVのleft-shifted breakend coordinateが報告される。

1、diploidSV.vcf.gz

2、somaticSV.vcf.gz

3、candidateSV.vcf.gzとそのサブセットのcandidateSmallIndels.vcf.gz

4、tumorSV.vcf.gz

他にもStatisticsファイルが出力される。

1、alignmentStatsSummary.txt

2、svLocusGraphStats.tsv

3、svCandidateGenerationStats.tsv

4、svCandidateGenerationStats.xml

複数サンプルがある場合、bamの@RGのサンプル名(@RGのSM:のところ) をspecificな名前にしてbamを作ってください。bamを作った後に修正するならまとめ のsamtools view部分参照。

 

GIhubだけでもかなり丁寧に説明されています。まずGithubのマークダウンのマニュアル(リンク)に目を通し、方法論についての詳細は論文(publishされたのはBioinfomaticsジャーナルの"SEQUENCE ANALYSIS")のsupllementary(ダイレクトリンク)を読んでください。

引用

Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications
Chen X, Schulz-Trieglaff O, Shaw R, Barnes B, Schlesinger F, Källberg M, Cox AJ, Kruglyak S, Saunders CT

Bioinformatics. 2016 Apr 15;32(8):1220-2

 

*1

"docker run"時、必要に応じてメモリlimitのオプションも使って下さい。

60G以上になるとkillするなら -m "60g" 

circulating cell free DNAから超高感度な変異検出を行う SiNVICT

 

 精密腫瘍学(precision oncology)の最も有望な分野の1つは、患者に合わせたカスタムターゲット療法の開発である。このような療法の開発および効率的な適用を成功させるには、患者の腫瘍DNAの治療誘発性変化の効率的かつ安価な同定手段とモニタリング方法を必要とする。残念なことに、特に進行性のガンでは、ガンの罹患率および死亡率の主な原因は、組織サンプリングに容易にアクセスできない複数の転移性病変の発生である。例えば、前立腺ガンでは> 90%以上の転移が骨および/または深部リンパ節で起こる(Bubendorf et al、2000)。このような部位を生検は患者の死亡率と関連してしまうため、したがって一般的には行われない。

 1948年以来、哺乳類の血液中に循環している無細胞DNA(cfDNA)の存在が知られている(Mandel、1948)。 cfDNAは細胞死に至る(壊死/アポトーシス)細胞(正常および腫瘍の両方)から放出されると考えられている - 1994年、癌患者の血液中に変異RAS遺伝子断片が検出された(see Schwarzenbach et al 、2011年)。 cfDNAを生成する非特異的メカニズムは、サンプリング変動の対象となる患者の全ての腫瘍DNAを、そしておそらく腫瘍の血流への完全な表現となる。以前の研究では、例えば、著者たちは、去勢抵抗性前立腺癌(CRPC)患者のcfDNAから様々な変異型AR(アンドロゲン受容体)遺伝子の存在を観察したが(Azad et al、2015)、これをもっとも説明できるのはそれぞれの患者の体内に癌細胞に複数のsubpopulationsが存在していたということである。複数の腫瘍foci / subclonesのこの完全な表現は、腫瘍由来のDNA源としての血漿の使用に重要な利点を提供する。残念なことに、患者の血液中に正常および腫瘍DNA両方が存在していることは、cfDNAシーケンシング分析に重大な問題を提起する。事態を悪化させているのは、腫瘍DNAは複数のサブクローンに由来することが多く、従って高度にheterogeneousであることである(一部略)。

 特定の集団内の1塩基変異(SNV)およびindelsを見出すために、体細胞および生殖細胞系や同時に複数サンプリングしてのWhole Genome Shotgun Sequencingを使用して突然変異コーラーが存在する。例えば、GATK(McKenna et al、2010)、VarScan2(Koboldt et al、2012)、Freebayes(Garrison and Marth、2012)、Strelka(Saunders et al、2012)(Strelka2の簡単な紹介)、MuTect(Cibulskis et al、2013)その他。これらのツールのほとんどは、頻度論またはベイズのアプローチを使用して、遺伝子座のバリアントコールがノイズ由来の偽陽性コールなのか(シーケンシングまたはマッピングエラー)、真の突然変異なのか推定する。その中でも、VarScan2はいくつかの発見的手法を使用して候補セットのサイズを縮小し、フィッシャーの正確検定のようないくつかの統計的検定を腫瘍/良性腫瘍ペアに適用して体細胞突然変異をコールする。さらに、strand biasなどの追加要因に基づきフィルタリングを可能にするポストプロセッシング機能も備えている。 Freebayes、MuTect、Strelkaなどの他のツールは、変異を呼び出すためにベイジアンの文脈で突然変異している場所の前後の確率を利用する。残念なことに、これらのツールは次のようなデータに対応していない、(i)患者の複数の時点からのシーケンシングデータ、(ii)非常に高いリードデプス(平均20〜30 k、90 kまで、おそらく将来的にはより高い)、(iii)極端に低い希釈率(約0.01%の変異の対立遺伝子パーセンテージ(Lipson et al、2014)、(iv)腫瘍内のheterogeneityが高いサンプル、または(v)システマティックなノイズの多いサンプル。さらに、ctDNAレベルは、限局性疾患を有する患者および治療を受けた患者における既存のctDNA検出アプローチの分析感度(Bettegowda et al、2014)よりも低くなり得る。

 上記の問題に対処するために、非常に高いリードデプスと非常に低い希釈率を扱える計算ツールSiNVICTを紹介する。 SiNVICTは、単一の腫瘍サンプル、複数の腫瘍サンプルバッチ、または単一の患者の複数時点でシーケンスされた複数サンプルで実行できる。この機能により、SiNVICTは1人の患者の複数のガンステージの同時シーケンシング解析や、異なる患者グループサンプルの同時シーケンシング解析ができる。これらのサンプルが類似の疾患進行および希釈レベルを有する場合、SiNVICTはその系統的なノイズを特徴付けるためにバッチの信号対雑音比(SNR)を利用し、シーケンシングされた領域のノイズの不均一性による偽陽性コールを減らそうと努力する。

 2つのシークエンシングプラットフォームで得られた同じ腫瘍サンプルのデータを使いSiNVICTの堅牢性を評価した(Illuminaの0.1%置換; IonTorrentの1%indel (Glenn、2011)]。著者らの実験は、SiNVICTが両方のシーケンシングプラットフォームで生成されたデータのコールに非常に敏感であることを示している。(一部略)。

 重要なことであるが、SiNVICTはユニークな問題に対処するものであり、既存の一般的なSNVおよびindelコーラー(例えばゲノム解析ツールキット)には匹敵しない。それらのツールは、シーケンシングデプスが高いデータセットの一部を処理するのが一般的で、 同じリードはPCR duplicatesとしてマークされる。しかしながら、ディープアンプリコンシーケンシングでは、同一のリードは必ずしもPCRアーティファクトとはならない、

 

gene target amplicon sequencing向けのツールです。論文では、SiNVICTの評価に、イルミナのIllumina DesignStudio (リンク)(*1)で作った14遺伝子のプローブ配列mixtureを使ったアンプリコンシーケンシングを行っています。

 

SiNVICTに関するツイート。


インストール

cent os7に導入した。

本体 Github

git clone --recursive https://github.com/sfu-compbio/sinvict.git
cd sinvict/
make

./sinvict -h

$ ./sinvict -h

 

SiNVICT: Ultra Sensitive Detection of Single Nucleotide Variants and Indels in Circulating Tumour DNA.

 

Allowed arguments:

--help or -h : Print help message.

--error-rate or -e : Error rate for the sequencing platform used.

--min-depth or -m : Minimum Read Depth required for high confidence in a call.

--left-strand-bias or -l : Lower limit for the strand bias value interval to be used in assessing the confidence of a call.

--right-strand-bias or -r : Upper limit for the strand bias value interval to be used in assessing the confidence of a call.

--read-end-fraction or -f : Average position of the called base on the reads supporting the call as a fraction.End values such as 0.01 as useful for filtering read end artifacts.

--qscore-cutoff or -q : Cutoff value for the qScore assigned to each call by the Poisson model used.

--tumor-directory-path or -t : Specifies directory for the input files.

--output-directory-path or -o : Specifies directory for the output files.

--use-poisson-germline or -s : Use a more robust poisson model to guess somatic/germline status.

 

--tumor-directory-path and --output-directory-path must be specified

 

Usage: ./sinvict -e=[error-rate] -m=[min-depth] -l=[left-strand-bias] -r=[right-strand-bias] -f=[read-end-fraction] -q=[qscore-cutoff] -t=<tumor-directory-path> -o=<output-directory-path> -s=<use-poisson-germline>

 

ラン

SiNVICTを使うにはbamを作成し、リードカウントを得る必要がある。著者は以下の手順で進めることを推奨している。

  1. Trimming FASTQ files (optional, fastq->fastq).
  2. Mapping (fastq -> bam).
  3. Recalibration and error correction (optional, bam -> bam).
  4. Obtaining mapping statistics per location (bam -> readcount).

(1) クオリティトリミングしアダプターが残っているなら除く。(2) bwaやmrFASTでマッピングしてbamを作成し、(3) ABRAなどで不正確なアライメントを改善しノイズを減らす(紹介)。(4) bam-readcountを使ってリードカウントを得る(紹介)。

 

出力されたテキストファイルをディレクトリ(ここではread-count_dir/ とする)に収納し、SiNVICTのコマンドを実行する。複数サンプルあるなら、全てディレクトリに収納しておく。

mkdir output
sinvict -t read-count_dir -o output
  • --min-depth    Minimum required read depth for a call to be considered reliable. (Default: 100)
  • --error-rate    Error Rate for the sequencing technology used (e.g. Illumina, Ion Torrent, ...) (Default: 0.01)
  • --left-strand-bias and --right-strand-bias    The strand bias values in the range [leftStrandBias, rightStrandBias] will be considered reliable. This [lsb, rsb] interval has to be between [0,1]. (Defaults: 0.3 and 0.7)
  • --read-end-fraction    Despite the trimming step, calls can be marked as low confidence according to the average position of the base on the reads that support a call. This value should be within range 0-1. (Default: 0.01)
  • --qscore-cutoff    The poisson model used by SiNVICT assigns a QScore to every call. Calls with a QScore below the user defined threshold will be considered low confidence. This value should be in range 0-99. (Default: 95)
  • --tumor-directory-path    The path to the directory where the readcount files are located.
  • --output-directory-path    The path to an empty directory where the output files will be generated.

結果はoutput/に出力される。 

$ ls -lth output/

合計 16K

-rw-rw-r--. 1 parallels parallels 7.8K  7月 28 22:55 calls_level1.sinvict

-rw-rw-r--. 1 parallels parallels 7.5K  7月 28 22:55 calls_level2.sinvict

-rw-rw-r--. 1 parallels parallels 7.3K  7月 28 22:55 calls_level3.sinvict

-rw-rw-r--. 1 parallels parallels 7.1K  7月 28 22:55 calls_level4.sinvict

-rw-rw-r--. 1 parallels parallels    0  7月 28 22:55 calls_level5.sinvict

-rw-rw-r--. 1 parallels parallels    0  7月 28 22:55 calls_level6.sinvict

フィルタリングレベルで6つのバリアントコールが出力される。6のcalls_level6.sinvictが全てのフィルタリングをクリアしたファイルとなる。テスト時は5でサイズがゼロになった。

  1. Poisson model: calls_level1.sinvict
  2. Minimum Read Depth filter: calls_level2.sinvict
  3. Strand-bias filter: calls_level3.sinvict
  4. Average position of called location among all reads supporting the call: calls_level4.sinvict
  5. Signal-to-Noise ratio filter: calls_level5.sinvict
  6. Homopolymer Regions filter: calls_level6.sinvict

出力はVCFに似た形式です。各フィールドはGIthubのページで解説されています。

 

引用

SiNVICT: ultra-sensitive detection of single nucleotide variants and indels in circulating tumour DNA.
Kockan C, Hach F, Sarrafi I, Bell RH, McConeghy B, Beja K, Haegert A, Wyatt AW, Volik SV, Chi KN, Collins CC, Sahinalp SC
Bioinformatics. 2017 Jan 1;33(1):26-34

 

参考

*1

DesignStudioを用いたプローブデザインの方法と最適化のヒント

https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2016_ilmn_webinar_designstudio_session1.pdf

バリアントをランク付けする Variant Ranker

 

 変異を特定することは、病気の病因を理解する上で重要である。ハイスループットな次世代ゲノム技術の進歩により、ゲノムシーケンシング、エクソンシークエンシング、RNA-SeqおよびChIP-Seqは、複雑なメンデル症の感受性遺伝子座を同定するための標準となっている。課題は、これらの手法が因果関係の変異を識別するために生成する膨大なデータをふるいにかけることにある。これに加えて、有害性の予測(例えばPolyPhen [論文より ref.1]、SIFT [ref.2]、MutationTaster [ref.3])や保全(例えば ' (PhyloP [ref.4]、SiPhy [ref.5]、GERP [ref.6])、異なるツールからの予測にかなりのばらつきが存在することが問題になる。さらに、バリアント機能の注釈はデータベースごとに異なる傾向がある。 SnpEff [ref.7](簡単な紹介)、SeattleSeq [ref.8]、ANNOVAR [ref.9]のようなバリアントのアノテーションには、非常に有用なツールがいくつか存在するが、バリアントをランク付けする能力はない。 eXtasy [ref.10]やSPRING [ref.11]のようなツールは、同義置換以外をランク付けする。他のケースでは、VAAST [ref.12]やKGGSeq [ref.13]のようなツールは病気の原因となるバリアントの優先順位を決める便利なコマンドラインツールだが、通常はツールをダウンロードして実行するためにはある程度のプログラミング知識が必要となる。

 著者らは、様々なアルゴリズムやデータベースからの変異型の予測とアノテーションをそれぞれ組み合わせる簡単な方法を提供することにより、ゲノムデータを解釈する課題に取り組むWebベースのバイオインフォマティクスツール、Variant Rankerを開発した。最終結果は、機能的研究または実験的検証のために引き継がれる変異のランク付けされたリストである。このツールを使用すると、いくつかのデータベースのバリアントに存在する既存の情報と利用可能な情報を結合した単一のスコアを計算することによって、優先順位付きバリアントのランク付けされたリストが生成される。 Variant Rankerは、de factoVCF [ref.14]とANNOVAR [ref.9]の形式を使用して、すべてのタイプのシーケンシングデータに適用できる。このツールの利点は、使いやすさ、すべてのバリアント(コーディングおよびノンコーディング)スコアリング、およびユーザーに提供されるフィルタリングの柔軟性である。ユーザーは、データベースを介して迅速に結果を照会することができ、バイオインフォマティクスのスキルが限られている人を含む、容易にアクセス可能で解釈可能な出力を提供する。ランク付けされた変異/遺伝子リストから重要な生物学的接続を発見するための下流の機能的濃縮分析の目的で、ネットワークアナライザーが統合されている。ネットワークアプローチを介してDAVID(アノテーション、ビジュアライゼーション、統合ディスカバリ用のデータベース、https://david.ncifcrf.gov)[ref.15、16]の表形式の結果を調べるネットワーク視覚化ツールである。

 バリアントランキングアルゴリズム

 使用可能なアノテーションを使用して、すべてのバリアントが0と1の間のウェイトを割り当てることによってエンコードされる。たとえば、ANNOVARアノテーションでのルールに従ってエクステリアにウェイトが与えられる。対応する重みは、それぞれ、1,5,5,6、4 / 6,3 / 6,2 / 6,1 / 6となる。保存アルゴリズムと予測アルゴリズムからのスコアは、各アルゴリズムを使用して対応する重みに変換される。例えば、変異がGERP [ref.6] score> 2(高度に保存されている)である場合、それに対応する重み1が与えられる。同様に、予測アルゴリズムPolyphen2の場合、重みは1(損傷)、0.5(おそらく損傷) MetaSVM [ref.23]、MutationTaster [ref.3]、MutationAssessor [ref.24]、およびFATHMM [ref.25]は、重み1(有害)および0(耐性)に続いて、SIFT [ref.2]、LRT [ref.22] 。 ENCODE [ref.27]エレメント、転写因子結合部位または保存部位を持つ領域の変異、およびdbSNPに存在しない場合、またはGWASカタログ[ref.28]またはクリニカルバリア[ref.29]データベースに存在する場合には、バイナリー加重(1または0)を適用する。集団頻度データベースでは、希少対立遺伝子に重み付けするために重み付け(1 - 対立遺伝子頻度)が割り当てられる。

 このように、新規であり、いくつかの予測アルゴリズム(異なるアルゴリズムは異なる予測を有する傾向がある)によって有害で​​あると予測される、機能的に重要な変異に対してより高いスコアが与えられる。各バリアントの合計得点は、バリアント当たりの符号化ウェイトの合計を取ることによって得られ、その後、すべてのバリアントは、それらの合計得点でソートされ、ランク付けされる。この方法には、バリアントごとに利用可能な情報に基づいて、単一のスコアが与えられる利点もある。

 

チュートリアル

http://paschou-lab.mbg.duth.gr/html5up/Tutorial33.html 

 

実行方法

Variant Ranker以外に複数のモジュールがある。

1、Variant Ranker: Single sample VCF/List of Variants(リンク

データセットの変異のランク付けとアノテーションを実行し、複数ソースの情報から各変異の有害性、新規性および既存の情報などの優先度を統合する。

VCFを指定し、パラメータを設定してsubmitする。

f:id:kazumaxneo:20180728135114j:plain

Sample Identifier:はなるべくspecificな名前にする。これは、Sample Identifierで検索できるが、その時他のデータもヒットしないようにするため。メールアドレスも記載する。ジョブが終わればしゅつ力リンクつきメールが届く。

ランタイムはバリアントコールの数などで変わってくるが、およそ30000コールあるテストvcfでは数十分で結果が得られた(チュートリアルに説明あり)。

チュートリアルでは、Variant Rankerの使用例として、論文で副産物的に報告(pubmed)された突発性の溶血性貧血(wiki)のWESデータ解析由来VCFを使っている。論文ではPKLRが原因遺伝子として報告されているが、Variant RNakerのVCFランキングでは、PKLRが4位に検出されている。下のリンク先にjob ID"nonregistered-2016-01-13_13:24:42"を貼り付ければ結果が見れる。

Welcome - : Spectral Ranking Home

f:id:kazumaxneo:20180728140732j:plain

出力はここには載せません。自分で実行してください。出力内容についての説明(リンク)。

チュートリアルには、他に、ファイファー症候群(wiki)とミラー症候群(wiki)の患者の合成データ由来VCF解析例がある(3万以上のbiallelic variantsからのランキング)。

 

2、Filtering Multi-sample VCF/Case-Control Filtering 

Cases / Controls間(ケースコントロール:wiki)のバリアントをフィルタリングし、ランク付けされたバリアントリストを取得する。ssnpshiftを"casecontrol"フラグつきで走らせて拡張VCFを作っておく必要がある(e.g., java -jar SnpSift.jar caseControl "++++0-----" cc.vcf )。

http://paschou-lab.mbg.duth.gr/index/login/CC.php

f:id:kazumaxneo:20180728142426j:plain

ラン時に簡単なフィルタリング条件を設定できるが、その時、以下のパラメータをつけることが推奨されている。

  1. Cases and not Controls: NCase>0 and NControl=0
  2. Controls and not Cases: NControl>0 and NCase=0

f:id:kazumaxneo:20180728143241j:plain

 

 

3、Filtering Result Explorer
Variant Rankerの結果をフィルタリングする。保存したバリアントランキング結果をアップロードする以外に、 Variant Rankerの結果のページから直接利用することもできる。ユーザー指定の条件でフィルタリングして、機能的に重要なバリアントを抽出するために使う。

f:id:kazumaxneo:20180728150425j:plain

 

 

4、Network Analyser

ネットワークアナライザーは、関連する遺伝子を調べるため biological connections を視覚化するためのツール。アノテーション、ビジュアライゼーション、DAVID(https://david.ncifcrf.gov)データベースを利用する。 例えばVariant Rankerのトップランク遺伝子をここに入力して関連する遺伝子を分析する。

http://paschou-lab.mbg.duth.gr/index/login/PathwayNR.php

f:id:kazumaxneo:20180728151913j:plain

ネットワークの表示にはFlash playerのプラグインが必要になります。

 

5、SNPtoGene

バリアントのポジションリストから遺伝子名に変換する。入力はBEDファイルに対応している。リストが多すぎると機能しない。

f:id:kazumaxneo:20180728145922j:plain

 

全モジュール。

http://paschou-lab.mbg.duth.gr/html5up/index.html

f:id:kazumaxneo:20180728134652j:plain

 

論文での使用例もある。

https://www.frontiersin.org/articles/10.3389/fnins.2016.00428/full

 

*他にも同じ名前のツールがあります。注意してください。

引用

Variant Ranker: a web-tool to rank genomic data according to functional significance.

Alexander J, Mantzaris D, Georgitsi M, Drineas P, Paschou P

BMC Bioinformatics. 2017 Jul 17;18(1):341.

 

SVの正確なGenotypingを行う SV2

 

 構造変化(SV)は、50bpよりも大きい染色体の構造の変化として定義される。 SVはヒトの遺伝子変異の主要な寄与因子であり、様々なヒト疾患にも関与している(Conrad、et al、2010; Sudmant、et al、2015)。特に、de novo構造変異(子孫固有で親にはない)は、特発性自閉症や知的障害のリスクにつながる(Brandler et al、2016)。新生児SVの正確な検出は、子どもの偽陽性や両親の偽陰性などの遺伝子型判定ミスのために困難である。さらに、ジェノタイピングエラーは、親から子への変異の伝達を強め、家系図の交絡分析にバイアスを与える。
 次世代シーケンシングからSVを特徴付けることは、SVの広範なサイズ(50bp〜50Mb)および広範なタイプのために困難な作業である。したがって、SVの多様性を完全に捉えるためには、スタンドアロンソリューション(GoNL Consortium、2014; Sudmant、et al、2015)として設計された多数のツールが必要になる。複数のSVコールメソッドからのバリアンコールと信頼スコアを調和させ統一されたジェノタイプを報告する方法が欠けている。

 ここでは、ペアエンドシーケンシングデータからの欠失およびduplicationのジェノタイピングのための機械学習アルゴリズムであるSV2について説明する。 SV2は、複数のSV発見アルゴリズムからのバリアントコールを、高度な遺伝子タイピング精度とde novo変異を検出する能力を備えた統一されたコールセットに迅速に統合することができる。

 SV2は、前処理、フィーチャの抽出、その後の遺伝子型分類の3段階で機能する(図1:下の図参照)。 ジェノタイピングは、4つの有益な特徴、すなわち、カバレッジデプス、不一致のペアエンド、スプリットリードおよびヘテロ接合性対立遺伝子比(HAR)を利用する。 ジェノタイピングは、1000ゲノムプロジェクト(1KGP)からの全ゲノム配列(WGS)で訓練された教師ありサポートベクターマシン分類器を用いて行われる。 得られた遺伝子型のVCFファイルには、1KGP(Sudmant et al、2015)によって記録された共通SVの遺伝子、反復要素および変異体識別のアノテーションが含まれている。

 

チュートリアル

https://github.com/dantaki/SV2/wiki/tutorial 

SV2に関するツイート。

f:id:kazumaxneo:20180728111325j:plain

ワークフロー。GIthubより転載。

 

インストール

mac os10.12のPython2.7.14環境でテストした。

依存

  • samtools (pileupで使う)
  • Bedtools 2.25 or later 

python2.7とライブラリ

  • cython
  • numpy
  • pandas
  • pybedtools
  • pysam 0.9+
  • scikit-learn v0.19+

本体

Github

github.com

pip install sv2

 > sv2 -h

$ sv2

                       ____

  _____________   ___ |___

 /   _____/   /   // ___/

 _____      Y   //_____)

 /              /

/_________/   ___/

 

Support Vector Structural Variation Genotyper

Version 1.4.3.4    Author: Danny Antaki <dantaki at ucsd dot edu>    github.com/dantaki/SV2

 

  sv2 -i <bam ...> -v <vcf ...> -b <bed ...> -snv <snv vcf ...> -p <ped ...> [OPTIONS]

 

input arguments: github.com/dantaki/SV2/wiki/Options#input-arguments

 

  Note: input arguments can take multiple files, separated by space <-i BAM1 BAM2 ...>

 

  -i, -bam    ...     bam file(s)

  -v, -vcf    ...     vcf files(s) of SVs

  -b, -bed    ...     bed files(s) of SVs

  -snv        ...     snv vcf files(s), must be bgzipped and tabixed

  -p, -ped    ...     ped files(s)

 

genotype arguments: github.com/dantaki/SV2/wiki/Options#genotype-arguments

 

  -g, -genome         reference genome build <hg19, hg38, mm10> [default: hg19]

  -pcrfree            GC content normalization for pcr free sequences

  -M                  bwa mem -M compatibility, split-reads flagged as secondary instead of supplementary

  -merge              merge overlapping SV breakpoints after genotyping

  -min-ovr            minimum reciprocal overlap for merging [default: 0.8]

  -no-anno            genotype without annotating variants   

 

  -pre                preprocessing output directory, skips preprocessing

  -feats              feature extraction output directory, skips feature extraction

 

classifier arguments: github.com/dantaki/SV2/wiki/Options#classifier-arguments

 

  -load-clf           add custom classifiers (-load-clf <clf.json>)

  -clf                define classifier for genotyping [default:default]

 

config arguments: github.com/dantaki/SV2/wiki/Options#config-arguments

 

  -download           download required data files

  -hg19               hg19 fasta file

  -hg38               hg38 fasta file

  -mm10               mm10 fasta file

  -ini                configuration INI file [default: $SV2_INSTALL_PATH/config/sv2.ini]

  

optional arguments:

 

  -L, -log            log file for standard error messages [default: STDOUT]

  -T, -tmp-dir        directory for temporary files [default: working directory]

  -s, -seed           random seed for preprocessing shuffling [default: 42]

  -o, -out            output prefix [default: sv2_genotypes]

  -O, -odir           output path, location for sv2 output directories [default: working directory]

 

  -h, -help           show this message and exit

 

——

 

ラン

テストランの流れを確認する(リンク)。

準備1: テストリファレンスゲノムデータをダウンロードし、パスを記載する。

#リファレンスのダウンロード(持ってなければ)
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz
gunzip human_g1k_v37.fasta.gz
samtools faidx human_g1k_v37.fasta

#付属データのダウンロード
sv2 -download

#fastaのパスを定義しておく。特にメッセージは出ないので、成功してるかわかりにくい(*1)。 
sv2 -hg19 /full/path/to/human_g1k_v37.fasta

 

準備2: テストデータをダウンロードする。

wget https://github.com/dantaki/SV2/releases/download/sv2v1.4.0/sv2_tutorial.zip 
unzip sv2_tutorial.zip
cd sv2_tutorial/

データは、1000Genomesのlow coverageサンプルをダウンサンプリングしたもので、bwa memでhg19にアライメントされbamになっている。VCFファイル1kgp_phase3_tutorial.vcfは 1000 Genomesの成果であるphase3 VCFをランダムサンプリングしたもの。HG00096_sub.vcf.gzはダウンサンプリングしたbamからFreebayesでバリアントコールしたもの。bgzipで圧縮されtabixでindexが付いている。bamとこれらをsv2の入力として使う。他にGWASで使うPLINKのバイナリファイル(HG00096.pedと1kgp.ped)などが入る。

 

テスト1。Freebayes1つのバリアントコール結果からgenotypingを行う。 

sv2 -i HG00096_sub.bam -v 1kgp_phase3_tutorial.vcf -snv HG00096_sub.vcf.gz -p HG00096.ped -o HG00096_sv2
  • -v  vcf files(s) of SVs
  • -snv vcf files(s), must be bgzipped and tabixed

  • -p  ped files(s)

  • -o  output prefix [default: sv2_genotypes]
  • -min-ovr minimum reciprocal overlap for merging [default: 0.8]

出力。図1の3段階のフローに相当するディレクトリができる。最終出力はsv2_genotypes。genotyeはvcfとbedで出力される(*2)。 f:id:kazumaxneo:20180728111056j:plain

 

スト2。複数のバリアントコール結果を統合してgenotypingを行う。 Manta(複数の間接的なSVシグナルを使い、コントロールケース実験の精度の高いSVコールを行う)のコールと、Lumpy(mantaと同様、複数手法を使う高感度なSVコーラー。論文発表はlumpyが先)のコールを使う。

sv2 -i HG00096_sub.bam -v 1kgp_phase3_tutorial.vcf -snv HG00096_sub.vcf.gz -p HG00096.ped -o HG00096_sv2

 

 # merge breakpoints with 95% reciprocal overlap

sv2 -i HG00096_sub.bam -b lumpy.bed manta.bed -snv HG00096_sub.vcf.gz -p HG00096.ped -o HG00096_95merged -min-ovr 0.95 -pre sv2_preprocessing/ -feats sv2_features/

 

SV2では、各サンプルごとにgenotypingを行い、それから統合し、genotpeの行列を出力するフローを推奨している(コマンド例はリンク先一番下のSV2 Guidelinesを参照)。

 

 

修正

mantaの説明をmetasvと勘違いして記載しておりました。修正します。

引用

SV2: accurate structural variation genotyping and de novo mutation detection from whole genomes.

Danny Antaki William M Brandler Jonathan Sebat

Bioinformatics, Volume 34, Issue 10, 15 May 2018, Pages 1774–1777,

 

*1 /sv2/config/sv2.ini を直接開きFASTAのパスを確認する。直接編集でもOK。

Running without install step (aka "config on the command line") · Issue #10 · dantaki/SV2 · GitHub

 

*2 出力について

https://github.com/dantaki/SV2/wiki/Output#vcf-output

 

解析内容などから適したインフォマティクスツールを探せる Bio-TDS

 

 バイオインフォマティクスと計算生物学は、バイオサイエンスとバイオメディカル研究において重要な役割を果たす。研究者が実験プロジェクトを設計する際には、データから新しい知識発見に至る最も関連性の高いバイオインフォマティクスツールキットを見つけることが大きな課題となる。 Bio-TDS(Bioscience Query Tool Discovery Systems、http://biotds.org/)は、研究者が最も適用可能な分析ツールをフリーテキスト検索して探せるように開発された。 Bio-TDSは、複数のバイオサイエンス領域(例えばゲノム、プロテオーム、バイオイメージング)の統合されたコミュニティリポジトリから12000個の分析ツールを照会する能力をユーザに提供する、柔軟な検索システムである。 Bio-TDSの主なコンポーネントの1つは、アノテーション、キュレーション、クエリ処理、評価のためのオントロジー自然言語処理ワークフローである。 Bio-TDSの科学的影響は、biological data analysisに焦点を当てたサイトであるBiostarsにおいて研究者が質問したサンプルを用いて評価された。 Bio-TDSは、5つの同様のバイオサイエンス分析ツール検索システムと比較され、関連性および完全性に関して他のものよりも優れていた。 Bio-TDSは、知識発見プロセスにおいて、データ分析に必要な、研究者の質問と最も関連性の高い計算ツールセットを結びつける能力を提供する。

 

Getting started

http://biotds.org/help/gettingstarted.xhtml;jsessionid=65897f4ccf9e2f234ddcd3a7883b

Overview

他のチュートリアル動画

http://biotds.org/help/videos.xhtml;jsessionid=672bf8c8cc3e16f417347fd7b17a

Bio-TDSに関するツイート。

 

ラン

http://biotds.orgにアクセスする。

f:id:kazumaxneo:20180726202446p:plain

 

1、キーワード検索

例えばクエリとしてRNAとタイプ。インクリメンタルサーチに対応しており、候補が表示される。RNA secondary structureを選択。Searchボタンをクリック。

f:id:kazumaxneo:20180727004136p:plain

候補のツールが表示される。

f:id:kazumaxneo:20180727004434p:plain

クリックすると、詳細とツールへのリンクが表示される。

f:id:kazumaxneo:20180727004428p:plain

ツール名を直接打って検索することもできる。

 

2、分野から検索。

ツールが研究分野毎に分類されている。画像はBioinformaticsを選択したところ。

f:id:kazumaxneo:20180727005128p:plain

 

3、Methodから検索。

画像はpathwayを選択したところ。

f:id:kazumaxneo:20180727005450p:plain

 

4、Data formatから検索。

Alignmentから検索。

f:id:kazumaxneo:20180727005928p:plain

 

引用
Bio-TDS: bioscience query tool discovery system

Gnimpieba EZ, VanDiermen MS, Gustafson SM, Conn B, Lushbough CM

Nucleic Acids Res. 2017 Jan 4;45(D1):D1117-D1122

 

シーケンスクオリティを絵文字で表す fastqe

そのままのツールです。ツイッターで教えてもらいました。

 

インストール

mac os 10.13のpython3.4でテストした(pyton2環境でも多分動きます)。

依存

  • pyemojify
  • BioPython
  • NumPy

pipで必要なライブラリを入れておく。

pip install pyemojify numpy Biopython

本体 Github

 

ラン

実行するにはfastqファイルを引数で指定する(iIllumina 1.8+/Sanger format)。

python input.fq

 

InSilicoSeq(リンク)を使い発生させたリードを分析してみる。

#ランダムにバクテリアゲノムを1つダウンロードしてリードを発生
iss generate --ncbi bacteria -u 1 --model NovaSeq --output ncbi_reads -p 8

分析実行

python fastqe.py ncbi_reads_R1.fastq 

python fastqe.py ncbi_reads_R1.fastq 

🙀🙀🙀🙀🙀😬😬🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈🙈😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😬😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿😿🙅

対応絵文字はコードから確認して下さい(リンク)。

同じデータをfastqcでも分析。

> fastqc #GUI起動

f:id:kazumaxneo:20180726190026p:plain

NovaSeqのクオリティはかなり良好ですね。

 

ちょっとしたジョークアプリだと思いますが、真剣に導入したくなりました😬😬😬

引用

https://github.com/lonsbio/fastqe