macでインフォマティクス

macでインフォマティクス

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

blast結果をインタラクティブなヒートマップで視覚化する BLASTmap

 

 植物と病原体の相互作用の結果を決定する多数の遺伝子が現在発見されている。たとえば、免疫受容体、感受性因子、病原体エフェクター、およびそれらの宿主標的など。ターゲットエンリッチメントシーケンスは、遺伝子型固有のゲノムアセンブリを最初に生成する必要なく、これらの目的の遺伝子を優先的に再シーケンスする手段を提供する。 Basic Local Alignment Search Tool(BLAST)は、ここで開発されたBLASTmapと組み合わせて、ターゲット種または最も近い関連属をリファレンスとして使用することにより、そのような遺伝子を特異的にターゲットとするプローブを設計するために使用できる。

 BLASTアラインメントの視覚化のために多数のプログラムが利用可能である。ただし、現在、ビットスコアなどの大規模なBLAST出力属性を視覚的に比較する専用のプログラムはない。数千のBLAST結果を迅速かつ効率的に比較する必要性により、BLASTmapが開発された。これは、Blast結果をインタラクティブなヒートマップとしてクラスタリングおよび表示するためにカスタマイズされたShiny Rパッケージを使用して作成されたインタラクティブなWebアプリケーションである。論文では、カスタムDNA / RNAプローブシーケンスを分析し、4つのプローブでジャガイモR2耐病性遺伝子ファミリーの特定の包括的エンリッチメントに十分であることを視覚的に判断するため、BLASTmapが正常に適用された例を示す。  

 

推奨インターネットブラウザはChromeVivaldiOperaとされている。

 

使い方

1、BLAST実行

ローカルBLASTを行う。出力はタブ出力を指定する(Query name, hit name, percentage identity, alignment length, mismatches, gaps, query start, query end, hit start, hit end, e-value, bitscore、が記載されている事)。

#blastデータベース(protein)
makeblastdb -dbtype prot -in database_eq.fasta -out blastn_database
#blastデータベース(nucletoide)
makeblastdb -dbtype nucl -in database_eq.faa -out blastp_database

#blastp
blastp -db blastp_database -query query.faa -outfmt 6 -out out.txt -num_threads 8 -evalue 1e-5
#blastn
blastn -db blastn_database -query query.fasta -outfmt 6 -out out.txt -num_threads 8 -evalue 1e-5

#diamond blastx (紹介)
diamond blastx --query input.fa \
--db uniprot_ref_proteomes.diamond.dmnd \
--outfmt 6\
--sensitive \
--evalue 1e-5 \
> blast.out

 

2、BLASTmapへのアクセス

https://ics.hutton.ac.uk/blastmap/ にアクセスする。タブを切り替える方式になっている。

f:id:kazumaxneo:20200303004449p:plain

 

まずimportにアクセスしてデータを読み込む。ここではテストデータを使う。チェックボックスにチェックを入れる。

f:id:kazumaxneo:20200304111411p:plain

 

interactive heatmapのタブに切り替え、ヒートマップに描画する要素を選択する。

f:id:kazumaxneo:20200304110330p:plain

決めたら左上のPlot heat mapボタンをクリックする。

 

インタラクティブなヒートマップで視覚化された。

f:id:kazumaxneo:20200303223143p:plain

ヒートマップのmatrixはクエリシーケンス数 x ヒットした配列数で表現されている。クエリ配列数やヒット配列数が少なすぎると視覚化されないので注意する。

 

左のメニューから様々な条件でフィルタリングができる。

f:id:kazumaxneo:20200303223746p:plain

項目はかなり多く、スクロールしないと見えないので注意。

 

アラインメント長を2000-bp以上にした。

f:id:kazumaxneo:20200303223930p:plain

 

入力ファイルが巨大だと作図でエラーになる。その時は、ファイルを読み込んでそのままimnportのタブでフィルタリングする。

f:id:kazumaxneo:20200304110611p:plain

赤字の部分がエラーからOKに変わった。

 

描画された。

f:id:kazumaxneo:20200304110605p:plain

 

Exportから作図したマップは出力できる。

f:id:kazumaxneo:20200303223556p:plain

こちらのタブのヒートマップはインタラクティブではない。あくまで出力の最終確認に使う。

 

 

コメント

uniprotの配列をデータベースにするとそのままIDで出力されます。これではヒートマップでIDしか表示されず、視覚化する価値が半減してしまいます。blastデータベース作成時に ”--taxonmap <proteomes.taxids>” を 指定してtaxids付きでblastデータベースを作っておくと良いかと思います。手順については、以前紹介したblobtoolsのデータベース作成パートを確認してください。taxidを使わずデータベースを作ってしまったなら、UniprotのID変換webサービスを使って後からtaxonomyに変換することもできますが、やや面倒です。

引用
BLASTmap: A Shiny-Based Application to Visualize BLAST Results as Interactive Heat Maps and a Tool to Design Gene-Specific Baits for Bespoke Target Enrichment Sequencing.

Baker K, Stephen G, Strachan S, Armstrong M, Hein I

Methods Mol Biol. 2018;1848:199-206

 

関連


 

 

アセンブリ結果を評価するwebサービス gVolante

2021 5/12 ツイート追記

 

 全ゲノムやトランスクリプトームなどの包括的な配列情報へのアクセスが増加するとともに、それらの品質を評価する必要性が高まっている。N50などのシーケンス長に基づくメトリックが標準になったが、これはアセンブリ品質の1つの側面のみを評価する。逆に、事前に選択されたリファレンスタンパク質をコードする遺伝子のカバレッジを分析すると、重要なコンテンツベースの品質評価が提供されるが、この目的で現在利用可能なパイプラインであるCEGMAおよびBUSCOには、ユーザーにとって使いやすいインターフェイスを持たない。ここでは、(i)以前に開発されたパイプラインCEGMAおよびBUSCOによるシーケンスセットのオンデマンドの完全性評価、および(ii)事前計算された完全性スコアの閲覧のためのオンラインツールを提供する、新しいWebサーバーgVolanteを紹介する。参照遺伝子のカバレッジだけでなく、配列の長さ(N50長さなど)にもとづいてgVolanteレポートのスコアで実行される完全性評価により、複数の側面での品質管理が可能になる。 gVolanteを使用すると、元のアセンブリの品質を複数のバージョン間(プログラムの選択やパラメーターの調整などで取得)で比較し、データベースセクションにあるパブリックリソースのスコアと比較して評価できる。gVoalteはhttps://gvolante.riken.jp/で無料で利用できる。

 

 

 

HP

https://gvolante.riken.jp

Tutorial

 

webサービス

https://gvolante.riken.jp/analysis.html にアクセスする。

f:id:kazumaxneo:20200229225206p:plain

 

1、multi-fastaファイルを指定する。各種圧縮形式にも対応している (.gz, .tgz, .bz2, .tbz, .tar or .zip). 

f:id:kazumaxneo:20200229224950p:plainUPLOADをクリック。

 

2、ジョブタイトルとメールアドレスを記載する。メールアドレスは任意だが、記載しておくとラン後にメールが届く。Cut-off lengthは、すべての配列を使用する場合は"1"を入力。テストデータはコード領域のアミノ酸配列のためPeptideを選択。

f:id:kazumaxneo:20200229225705p:plain

他のフィールドについては上のリンク先にあるチュートリアルを参照。

 

出力

ランが終わるとメールが届く。

BUSCOのスコア、配列長などがまとめられる。

f:id:kazumaxneo:20200301104844p:plain

出力内容についてはFAQ も参照。

 

パイチャート右のリンクからオロソログの詳細を確認できる。完全長検出されたオロソログ、partialにしか検出されなかったオルソログ、そして検出されなかったオルソログに分けられているので分かりやすい。

f:id:kazumaxneo:20200301104847p:plain

 

Misssingのオルソログを1つ見てみる。chorNOG00418のリンクをクリック。

f:id:kazumaxneo:20200301110436p:plain



eggNOG(オロソログデータベース。系統樹の様々な位置に存在するオルソログがどのような関係にああるか調べることができる)

f:id:kazumaxneo:20200301110841p:plain

 

aLeaves (目的の配列のホモログを検索 => ツリーを構築)

f:id:kazumaxneo:20200301111401p:plain

f:id:kazumaxneo:20200301113041p:plain

aLeavesの使い方は統合TVを確認して下さい。丁寧に解説されています。

 

 

gVolanteの結果はLINK TO RAW RESULTS〜からダウンロードできます。注意点ですが、ランタイムはデータのサイズ、ツールやデータベースの選択、サーバーの混雑度によって変化します。何度も同じジョブを投げないようにして下さい。

引用
gVolante for standardizing completeness assessment of genome and transcriptome assemblies.

Nishimura O, Hara Y, Kuraku S

Bioinformatics. 2017 Nov 15;33(22):3635-3637

 
aLeaves facilitates on-demand exploration of metazoan gene family trees on MAFFT sequence alignment server with enhanced interactivity

Kuraku S1, Zmasek CM, Nishimura O, Katoh K

Nucleic Acids Res. 2013 Jul;41(Web Server issue):W22-8

 

eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses
Jaime Huerta-Cepas, Damian Szklarczyk, Davide Heller, Ana Hernández-Plaza, Sofia K Forslund, Helen Cook, Daniel R Mende, Ivica Letunic, Thomas Rattei, Lars J Jensen, Christian von Mering, Peer Bork

Nucleic Acids Res. 2019 Jan 8; 47

 

参考

 

関連


パンゲノム解析のためbacteria populationsをシミュレートする SimPan

 

 細菌ゲノムは、広範な相同組換え、水平遺伝子導入、遺伝子損失、遺伝子重複などの複雑な進化の歴史によって形作られている。細菌ゲノムの定義されたセット内のすべての遺伝子で構成されるパンゲノムは、系統学的推論および集団研究の基礎を提供できる。ここでは、遺伝的に多様化された数千の細菌ゲノムからパンゲノムを構築できるパイプラインであるPEPPAを紹介する。 PEPPAは、ツリーベースおよびシンテニーベースのアプローチを組み合わせて実装し、パラロガス遺伝子を特定および除外する。これにより、全ゲノムおよびコアゲノムMLSTタイピングスキームの構築が容易になる。 PEPPAは、個々のゲノム内の遺伝子および偽遺伝子の一貫したアノテーションをサポートする類似性ベースの遺伝子予測も実装する。これは、専門の手動キュレーションとほぼ同じくらい正確である。 PEPPAパッケージには、PEPPA_parserが含まれている。これは、アクセサリー遺伝子の内容とMLSTコアゲノムの対立遺伝子の違いに基づいてツリーを計算する。また、細菌のパンゲノムの進化をシミュレートするための新しいパイプラインであるSimPanも紹介する。 PEPPAは、経験的およびシミュレートされたデータセットの両方で、4つの最先端のパンゲノムパイプラインと比較された。それは他のどのパイプラインよりも高い精度と特異性を示し、パンゲノムの計算ではそれらとほぼ同じ速さであった。その能力の実証として、PEPPAを使用して、少なくとも80種にわたる3,170の代表的なゲノムから40,000を超える遺伝子の連鎖球菌パンゲノムを構築した。結果として生じる遺伝子と対立遺伝子のツリーは、この属のゲノムの多様性の前例のない概要を提供する。

 

PEPPAは以前紹介しました。ここではSimPanを紹介します。

http://kazumaxneo.hatenablog.com/entry/2020/02/16/073000

 

インストール

pytho3.7の仮想環境でテストした(ホストOS; ubuntu18.04LTS)。

依存

SimPan runs in Python with versions >= 3.5 and requires two libraries:

  • numpy
  • ete3 

本体 Github

git clone https://github.com/zheminzhou/SimPan.git
cd SimPan/

python SimPan.py -h

# python SimPan.py -h

usage: SimPan.py [-h] [-p PREFIX] [--genomeNum GENOMENUM] [--geneLen GENELEN]

                 [--igrLen IGRLEN] [--backboneBlock BACKBONEBLOCK]

                 [--mobileBlock MOBILEBLOCK] [--operonBlock OPERONBLOCK]

                 [--aveSize AVESIZE] [--nBackbone NBACKBONE] [--nCore NCORE]

                 [--nMobile NMOBILE] [--pBackbone PBACKBONE]

                 [--pMobile PMOBILE] [--tipAccelerate TIPACCELERATE]

                 [--rec REC] [--recLen RECLEN] [--seqRec SEQREC]

                 [--insRec INSREC] [--delRec DELREC] [--noSeq]

                 [--idenOrtholog IDENORTHOLOG] [--idenParalog IDENPARALOG]

                 [--idenDuplication IDENDUPLICATION] [--indelRate INDELRATE]

                 [--indelLen INDELLEN] [--freqStart FREQSTART]

                 [--freqStop FREQSTOP]

 

SimPan is a simulator for bacterial pan-genome. 

Global phylogeny and tree distortions are derived from SimBac and the gene and intergenic sequences are simulated using INDELile.

 

optional arguments:

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

  -p PREFIX, --prefix PREFIX

                        prefix for all intermediate files and outputs. {DEFAULT: SimPan]

  --genomeNum GENOMENUM

                        No of genome in population [DEFAULT: 20]

  --geneLen GENELEN     [negative bionomial with r=2] mean,min,max sizes of genes [DEFAULT: 900,150,6000]

  --igrLen IGRLEN       [negative bionomial] mean,min,max sizes of intergenic regions [DEFAULT: 50,0,300]

  --backboneBlock BACKBONEBLOCK

                        [geometric] mean,min,max number of backbone genes per block [DEFAULT: 3,0,30]

  --mobileBlock MOBILEBLOCK

                        [geometric] mean,min,max number of mobile genes per block [DEFAULT: 10,0,100]

  --operonBlock OPERONBLOCK

                        [geometric] mean,min,max number of continuous genes that share the same coding strand [DEFAULT: 3,0,15]

  --aveSize AVESIZE     average gene number per genome (greater than nBackbone). [DEFAULT: 4500]

  --nBackbone NBACKBONE

                        number of backbone genes (present in common ancestor) per genome. [DEFAULT: 4000]

  --nCore NCORE         sizea of core gene (smaller than the size of backbone genes). [DEFAULT: 3500]

  --nMobile NMOBILE     size of mobile gene pool for accessory genome. [DEFAULT: 20000]

  --pBackbone PBACKBONE

                        propotion of paralogs in backbone (core) genes. [DEFAULT: 0.05]

  --pMobile PMOBILE     propotion of paralogs in mobile (accessory) genes. [DEFAULT: 0.4]

  --tipAccelerate TIPACCELERATE

                        grandient increasing of gene indels in recent times. [DEFAULT: 100]

  --rec REC             expected coverage of homoplastic events in pairwise comparisons. [DEFAULT: 0.05]

  --recLen RECLEN       expected size of homoplastic events. [DEFAULT: 1000]

  --seqRec SEQREC       Use homoplastic events to infer sequences. Use 0 to disable [DEFAULT: 1]

  --insRec INSREC       Use homoplastic events to infer gene insertions. Use 0 to disable [DEFAULT: 1]

  --delRec DELREC       Use homoplastic events to infer gene deletions. Use 0 to disable [DEFAULT: 1]

  --noSeq               Do not infer sequence but only the gene presence/absence. [DEFAULT: False]

  --idenOrtholog IDENORTHOLOG

                        average nucleotide identities for orthologous genes. [DEFAULT: 0.98]

  --idenParalog IDENPARALOG

                        average nucleotide identities for paralogous genes. [DEFAULT: 0.6]

  --idenDuplication IDENDUPLICATION

                        average nucleotide identities for recent gene duplications. [DEFAULT: 0.995]

  --indelRate INDELRATE

                        average frequency of indel events relative to mutation rates. [DEFAULT: 0.01]

  --indelLen INDELLEN   average size of short indel events within each gene. [DEFAULT: 10]

  --freqStart FREQSTART

                        frequencies of start codons of ATG,GTG,TTG. DEFAULT: 0.83,0.14,0.03

  --freqStop FREQSTOP   frequencies of stop codons of TAA,TAG,TGA. DEFAULT: 0.63,0.08,0.29

 

 

実行方法

それぞれ50遺伝子からなる10のゲノムをシミュレートする。

python SimPan.py --aveSize 50 --nBackbone 30 --nMobile 1000 -p test --genomeNum 10

 

はじめに、SimBacを使用して細菌ゲノムのグローバルな系統発生および組換えイベントがシミュレートされる。次に、ランダムなindelイベントにより、コアゲノムとアクセサリーゲノムの両方の遺伝子コンテンツが生成される。http://abacus.gene.ucl.ac.uk/software/indelible/を使い、これらのpan遺伝子の配列がfill inされる。

出力についてはGithubで確認して下さい。

引用

Accurate reconstruction of the pan- and core- genomes of bacteria with PEPPA

Zhemin Zhou* and Mark Achtman

bioRxiv preprint, Posted January 03, 2020

 

 

scientific dataのインタラクティブな視覚化ツール BPG interactive visualization of scientific data

 

  生物学的実験はますます大きく、多面的なデータセットを生み出している。そのようなデータを探索して観察結果を伝達することはますます困難になっており、堅牢な科学的データ視覚化の必要性は加速している[ref.1、2、3、4]。特にWebベースのインタフェースやローカルソフトウェアパッケージとして、無数のデータ視覚化ツールが存在する。残念ながら、これらは広く使用されているBioconductor [ref.5]のように、Rベースの統計パイプラインには簡単に統合できない。 R内には、ベースグラフィックス[ref.6]、ggplot2 [ref.7]、ラティス[ref.8]、Sushi [ref.9]、circlize [ref.10]、multiDimBio [ref.11]、NetBioV [ref.12]、GenomeGraphs [ref.13]、そしてそれらを含む多くの視覚化パッケージが存在する。 特定のタスクや分析の種類に焦点を当てた固有の視覚化パッケージもある[ref.15、16、17、18、19、20、21、22、23、24]。これらの中には、高解像度、適切なラベルサイズ、グレースケール用途に適し、赤 - 緑の色覚異常を持つ人々に見えるようなデフォルトのカラーパレットなどの、publication品質のデフォルトがないものもある。多くは重要なパラメータ化を必要とする。他のものは限られたプロットタイプを含んでいるか、マルチパネルの図を自動生成する機能は限定的であるか、または特定のデータタイプに制限されている。特定のプロット要素をハイライト表示し、それらをカスタマイズするために利用可能な一連のパラメータを自動的に識別し、プロットの変化をリアルタイムで視覚化するGUIインタフェースを介してインタラクティブにRコードを生成できるインタラクティブな視覚化ツールはほとんどない。したがって、これらの可視化パッケージのそれぞれには大きな価値があるが、それぞれ計算生物学者やデータ科学者にとって有益な機能がいくつか欠けている。

優れたビジュアライゼーションソフトウェアは、利用可能なデータタイプの多様性に合わせるために、さまざまなチャートタイプを作成する必要がある。高度にカスタマイズされた図に対して柔軟なパラメータ化を提供し、高解像度の出力を生成するなど、適切でpublication に適したデフォルト設定を採用しながら複数の出力フォーマットを可能にするべきである。さらに、それは既存の計算パイプラインとシームレスに統合しながら、直感的でインタラクティブなモードを提供する必要がある。周期的な開発を可能にして、パイプラインとインタラクティブなモード間で移行する能力があるべきである。最後に、適切な色の選択や特定のユースケースに合わせたレイアウトの提案など、優れたデザイン原則を奨励する必要がある。ユーザーがすぐに習熟するためには、詳細な例、チュートリアル、リアルタイムのインタラクティブなプロット調整機能、およびアプリケーションプログラミングインターフェイスAPI)が必要になる。今日まで、既存のビジュアライゼーションスイートでこれらのニーズを完全に満たしたものはない。

   このギャップを解消するために、本著者らはBPG(BoutrosLab.plotting.general)ライブラリーを作成した。これはグリッドグラフィックスシステムとラティスフレームワークを使ってRで実装されている。棒グラフや箱ひげ図などの一般的なプロットからマンハッタンプロットなどの特殊なプロットまで、幅広い種類のチャートタイプを生成できる(論文図1;コードは追加ファイル1にあり)。これらには、ドットマップを含む、いくつかの新しいプロットタイプが含まれる。マトリックスの内側に挿入された円のグリッドで、4次元データの表現を可能にする(論文図1n)。各プロット関数は高度にパラメータ化されており、プロットの美観を正確に制御できる。BPGのデフォルトパラメータは、公開に適した高解像度(1600 dpi)のTIFFファイルを生成する。ファイルの種類は、単にファイル拡張子を指定することによって指定される。目盛りの追加、フォントの選択、およびプロジェクト全体で一貫した印刷スタイルを作成するために連携するデフォルトの色など、その他のデフォルト値はグラフィックの一貫性に役立つ。デフォルト値は、高品質の数値を生成するように最適化されているため、手動で調整する必要はない。しかしながら、良いデフォルトでさえすべてのユースケースに適しているわけではない[ref.15]。追加ファイル2:図S1は、デフォルト設定または最適化設定(BPG、ベースRグラフィック、ggplot2、およびラティス)の4つの別々のグラフィックフレームワークを使用して作成された単一の散布図を示している。BPGは、デフォルトプロットと最適化プロットの両方に対して、他のフレームワークの半分のコードのみ必要とし、少なくとも同等の品質のプロットを作成した(追加ファイル3)。

  迅速なグラフィカルプロトタイピングを容易にするために、オンラインのインタラクティブなプロットインタフェースも作成した(http://bpg.oicr.on.ca)。このインタフェースにより、ユーザはパラメータ値を調整した結果を簡単かつ迅速に確認でき、それによってプロットの美観を正確に向上させることができる。このインタフェースによって生成されたRコードもダウンロードできる。

  BPGの重要な機能の1つは、複数のプロットを1つの図にまとめる機能である。これは、publicationsで広く使用されている手法である。これは、自動的にプロットを整列させ、最終的なFigure内のすべてのプロット要素にわたって線幅やフォントサイズなどのパラメータを標準化するcreate.multiplot関数によって実現されている。これは、PowerPointLaTeX、または他の同様のソフトウェアを使用した、遅くてエラーが発生しやすい数字の手動の組み合わせ、またはlayout()のような関数を使用して手動でプロット位置を直接整列するという時間のかかるパラメーター化を置き換えられる。複数のプロットを組み合わせる必要があるのは、データセットが複雑なためである。高次元のデータでは、関連するすべての情報を1つのチャートタイプにまとめることは困難である。複数のチャートタイプを組み合わせると、データをより詳細に視覚化できる。例えば、1つ目のプロットは異なるサンプルに存在する突然変異の数を伝える。2番目のプロットでは異なる変異タイプの割合を追加し、3番目のプロットではサンプルレベルの情報を追加する(論文図2)。論文図2のビジュアライゼーションの作成に使用したものを含む一連のサンプルデータセットをBPGに直接含めている。これらのデータセットからこのプロットを作成するためのソースコードは、追加ファイル4にある。

 

HP

https://labs.oicr.on.ca/boutros-lab/software/bpg

Guide

https://cran.r-project.org/web/packages/BoutrosLab.plotting.general/vignettes/PlottingGuide.pdf

Video tutorial

help

http://bpg.oicr.on.ca/API/BoutrosLab.plotting.general/5.3.4/index.html

 

 

webサービス

ここではオンラインのshiny アプリ版について、使い方を簡単に説明する。

http://bpg.oicr.on.ca にアクセスする。

f:id:kazumaxneo:20190410004831p:plain

 

サイトの構成を見るに、まず最初からロードされているデモデータを使って出力図を決定し、それから自分のデータを使うことを前提に設計されているように見える。そこで、この流れに従って使い方を見ていく事にする。

 

1、BPGにアクセスすると、defaultで以下のグラフが描画される(数秒~10秒程度かかる)。

f:id:kazumaxneo:20190607224020p:plain

 

最初に自動で描かれたファイルは以下の1-4の工程が自動実行されて描画されている。

f:id:kazumaxneo:20190607231513p:plain

他のパラメータはdefaultになっている。

 

2、パラメータは素早く変更できるように工夫されており、図の変更したい箇所をダブルクリックすると、左のメニューが指定パラメータのウィンドウに変化する。

f:id:kazumaxneo:20190607232456p:plain

Font Optionsタブからx軸の目盛りのフォントサイズを1.25=>2、カラーをblack => redに変更。

 

3、自分のデータを使う前に、デモデータの段階でどのようなグラフが適切か考える。左上からグラフの種類を切り替える。

f:id:kazumaxneo:20190607232833p:plain

 

barplotからboxplotに切り替えた。データはcars。他に定番のirisなども選べる。

f:id:kazumaxneo:20190607233225p:plain

(同じcarsのデータです)

 
他の例を見ていく。

hexbinplot

f:id:kazumaxneo:20190607233228p:plain

 

qqplot

f:id:kazumaxneo:20190607233548p:plain

 

histgram

f:id:kazumaxneo:20190607233421p:plain

 

manhattanplot

f:id:kazumaxneo:20190607233752p:plain

 

seqplot

f:id:kazumaxneo:20190607233826p:plain

 

densityplot

f:id:kazumaxneo:20190607233916p:plain

 

stripplot

f:id:kazumaxneo:20190607234021p:plain

 

violinplot

f:id:kazumaxneo:20190607234112p:plain

 

heatmap

f:id:kazumaxneo:20190607234144p:plain

 

自分のデータを使うにはBrowseからアップロードする。

f:id:kazumaxneo:20190608000623p:plain

Preview Dataで先頭を確認。

f:id:kazumaxneo:20190608000650p:plain

セパレータの種類とヘッダーの有無を選択する。セパレータのdefaultはコンマになっている。先頭にHeader行があるなら、Headerにチェックをつける。

f:id:kazumaxneo:20190607214045p:plain

 

それぞれのグラフでオススメの例を表示できる。下の画像のExampleのところでColourを選択、Loadをクリック。

f:id:kazumaxneo:20190608011555p:plain

 

右下を見ると使用データがexample1に切り替わっている(boxplotの例1(Colour))。

f:id:kazumaxneo:20190608011809p:plainこのパラメータ設定を真似れば、同じグラフが書ける。

 

右下のタブからR CODEタブを選択するとコードが表示される。再現性のある解析のため、図の作成工程を記録として残すのは重要である。また、GUIを使わずコンソールでランする時にも役に立ってくる。

f:id:kazumaxneo:20190608012048p:plain

f:id:kazumaxneo:20200302201553p:plain

f:id:kazumaxneo:20200302201751p:plain

f:id:kazumaxneo:20200302201755p:plain

f:id:kazumaxneo:20200302201800p:plain

f:id:kazumaxneo:20200302201811p:plain

f:id:kazumaxneo:20200302201905p:plain

f:id:kazumaxneo:20200302201903p:plain




 

機能が豊富すぎて説明しきれません。詳細は上に張ったリンクから確認して下さい。リンクのhelpには豊富な作成があります。

ローカルでの使い方についても、機会があれば説明したいと思います。

引用

BPG: Seamless, automated and interactive visualization of scientific data

Christine P’ng, Jeffrey Green, Lauren C. Chong, Daryl Waggott, Stephenie D. Prokopec, Mehrdad Shamsi, Francis Nguyen, Denise Y. F. Mak, Felix Lam, Marco A. Albuquerque, Ying Wu, Esther H. Jung, Maud H. W. Starmans, Michelle A. Chan-Seng-Yue, Cindy Q. Yao, Bianca Liang, Emilie Lalonde, Syed Haider, Nicole A. Simone, Dorota Sendorek, Kenneth C. Chu, Nathalie C. Moon, Natalie S. Fox, Michal R. Grzadkowski, Nicholas J. Harding, Clement Fung, Amanda R. Murdoch, Kathleen E. Houlahan, Jianxin Wang, David R. Garcia, Richard de Borja, Ren X. Sun, Xihui Lin, Gregory M. Chen, Aileen Lu,
Yu-Jia Shiah, Amin Zia, Ryan Kearns and Paul C. 

Boutros BMC Bioinformatics 2019 20:42

 

関連 

 

 

 

 

 

 

(モデル生物 )MNase-seqやchip-seeのアラインメントを2D plotで視覚化する plot2DO

 

 ヌクレオソーム、これは147 bpのDNAがA約1.7ターンでヒストンオクタマーに包まれる真核生物のDNAパッケージングの基本単位である(Luger、et al、1997)。標的部位へのDNA結合因子のアクセスは、これらの部位がヌクレオソームフリー領域(NFR)にある場合、約10〜20倍高くなる(Liu、et al、2006)。したがって、ヌクレオソームとNFRの正確な位置を知ることは、DNA結合と遺伝子調節を理解するために非常に重要である。
 現在、ヌクレオソームマッピングに最もよく使用される方法はMNase-seqである。これは、Micrococcal nuclease消化とそれに続くディープシーケンス(MNase-seq)ある。クロマチンがMicrococcal nuclease(MNase)で消化され、残りの未消化DNAフラグメントがハイスループットシーケンスにかけられる。残念ながら、MNaseには強い配列優先性があり(Dingwall e tal、1981;Hörzand Altenburger、1981)、MNase-seq実験から生じるヌクレオソーム断片はMNase消化の程度の影響を受ける(Chereji et al、2016; Chereji、et al、2017)。さらに、穏やかな消化の後、ゲノムの大部分はまだモノヌクレオソームDNA断片(〜150bpの長さ)に分割されておらず、さらなる分析から破棄されるが、大規模な消化の後、A / Tリッチ配列を占める多くのヌクレオソームモノヌクレオソーム断片のサンプルから過剰に消化され、失われる(Chereji、et al、2017)。そのため、MNase-seq実験では、消化レベルを慎重に制御する必要があり、特に複数のサンプルを比較する場合は、さまざまな程度の消化を常に考慮する必要がある。
 ここでは、ゲノムデータの2D占有率(2DO)をプロットするツールであるplot2DOを紹介する。これは、MNase-seqデータの初期品質チェックとして消化の程度を評価するだけでなく、MNase消化キネティクスからゲノムの機能領域の近くでのヌクレオソーム組織の洞察を得るためにも役立つ。

 Plot2DOは、Rで記述されたオープンソースパッケージであり、ターミナルのコマンドラインから起動できる。ユーザーは、プロットする分布のタイプ(未消化DNAフラグメントの占有率/カバレッジ、フラグメント中心の分布(ヌクレオソームダイアド)、またはフラグメントの5 '/ 3'末端の分布)、アラインメント(転写開始サイト(TSS)、転写終結サイト(TTS)、+ 1ヌクレオソーム、または特定のユーザー提供サイトのリスト)を選択する。ユーザーは、ヌクレオソーム集団の代表として使用されるフラグメントのサイズ制限を指定することにより、プロットされるウィンドウの幅を選択でき、消化されていないDNAのin silicoサイズ選択も実行できる。Plot2DOを使用すると、さまざまな生物(酵母、ハエ、mouse、線虫、マウス、およびヒト)から生成され、以下のゲノムバージョンのいずれかにマッピングされたペアエンドシーケンスデータを調査できる:sacCer3、dm3、dm6、ce10、ce11、 mm9、mm10、hg18、hg19。

 

デモデータ(bam)

https://onedrive.live.com/?authkey=%21ADHdTJiN14nIAFw&id=7713D31A609D5B2F%21102271&cid=7713D31A609D5B2F

 

インストール

CRANとBioconducterからたくさんのパッケージが導入されるため、サブのmacosマシン(mac pro2012, 10.13)でテストした。(*1)

本体 Github

git clone https://github.com/rchereji/plot2DO.git
cd plot2DO/
Rscript plot2DO_setup.R

Rscript plot2DO.R --help 

$ Rscript plot2DO.R --help  

Usage: plot2DO.R [options]

 

 

Options:

-f FILE, --file=FILE

Name of the file containing aligned sequencing data [options: BAM or BED file]

 

-t TYPE, --type=TYPE

Type of distribution to plot [options: occ, dyads, fivePrime_ends, threePrime_ends; default = occ]

 

-g GENOME, --genome=GENOME

Genome version

[options: sacCer3 (default) (S. cerevisiae); EF2 (S. pombe); dm3, dm6 (D. melanogaster);

ce10, ce11 (C. elegans); mm9, mm10 (M. musculus);

hg18, hg19, hg38 (H. sapiens); tair10 (A. thaliana)]

 

-r REFERENCE, --reference=REFERENCE

Reference points to be aligned [options: TSS (default), TTS, Plus1]

 

-s SITES, --sites=SITES

User-provided sites to be aligned (BED file)

 

-a ALIGN, --align=ALIGN

Points of the provided intervals to be aligned? [options: center (default), fivePrime, threePrime]

 

--siteLabel=SITELABEL

Label for the aligned sites [default = Sites]

 

-l MINLENGTH, --minLength=MINLENGTH

The smallest DNA fragment to be considered [default = 50]

 

-L MAXLENGTH, --maxLength=MAXLENGTH

The largest DNA fragment to be considered [default = 200]

 

-u UPSTREAM, --upstream=UPSTREAM

Length of the upstream region to be plotted [default = 1000]

 

-d DOWNSTREAM, --downstream=DOWNSTREAM

Length of the downstream region to be plotted [default = 1000]

 

-m COLORSCALEMAX, --colorScaleMax=COLORSCALEMAX

Maximum value on the color scale (e.g. 0.02)

 

--simplifyPlot=SIMPLIFYPLOT

Simplify the plot (show only the 2D heat map) [options: on, off (default)]

 

--squeezePlot=SQUEEZEPLOT

Simplify the plot and squeeze the heat map [options: on, off (default)]

 

-h, --help

Show this help message and exit

 

 

 

 

実行方法

plot2DOのほとんどすべての引数は任意。必須の引数は"--file=" (or "-f") で指定するアライメントのBAMファイルのみになる。リファレンスは、Saccharomyces cerevisiae S288Cがデフォルト対応している。 "--genome="で変更可能。

Rscript plot2DO.R --file=input.bam --genome=sacCer3
  • --genome=<Genome version>      

 [options: sacCer3 (default) (S. cerevisiae); EF2 (S. pombe); dm3, dm6 (D. melanogaster);ce10, ce11 (C. elegans); mm9, mm10 (M. musculus); hg18, hg19, hg38 (H. sapiens); tair10 (A. thaliana)]

出力

f:id:kazumaxneo:20200229170255p:plain

 

出力を見てみる(macos)。

open output/2D_occ_TSS/OCC_matrix.TSS.50_200.yeast_50U_MNase_SRR3649301.5M_reads.pdf \
-a /Applications/Preview.app/ 

f:id:kazumaxneo:20200229161219p:plain

3つのパネル が生成される。

(1)2D占有率(2DO)プロット(中央のヒートマップ)。指定された長さのDNAフラグメントの相対カバレッジを表す。ここでは転写開始部位(TSS)が基準点で、TSS周辺1 kbの範囲が含まれている。 赤色はカバレッジが高いことを示し、青色はカバレッジがゼロであることを示す。

(2)ヒートマップに示されているすべての長さのDNAフラグメントを積み重ねることによって生成される1次元の占有率(上部パネル)。

(3)シーケンスリードのサンプル全体からの各DNAフラグメントサイズに対応するパーセンテージを示すフラグメント長ヒストグラム(右パネル)。 

 

 


 複数描画して並べるなら"squeezePlot=on"フラグを立てて1パネル出力する。 上限は統一する。

#No.1
Rscript plot2DO.R --file=yeast_50U_MNase_SRR3649301.5M_reads.bam --squeezePlot=on -g sacCer3 -r Plus1 -m 0.03 -u 100 -d 100

#No.2
Rscript plot2DO.R --file=yeast_50U_MNase_SRR3649301.5M_reads.bam --squeezePlot=on -g sacCer3 -r Plus1 -m 0.03 -u 100 -d 100
  • -u     Length of the upstream region to be plotted [default = 1000]
  • -m    Maximum value on the color scale (e.g. 0.02)
  • --squeezePlot   Simplify the plot and squeeze the heat map [default: off ]

 

f:id:kazumaxneo:20200229180944p:plain

 

引用

plot2DO: a tool to assess the quality and distribution of genomic data

Răzvan V. Chereji

bioRxiv, Posted September 15, 2017.

 

*1

環境を汚したくないなら、rockerプロジェクトのrstudioなどを使う手もあります。また記事にします。

 

関連


 

インタラクティブなヒートマップを描く heatmaply

2020 2/29 誤字修正

 

 クラスターヒートマップは、高次元のデータを視覚化するための一般的なグラフィカルな方法である。その中で、数値のテーブルは、色付きセルのタイル状のマトリックスとしてスケーリングおよびエンコードされる。マトリックスの行と列は、パターンを強調するように順序付けられており、よく樹形図とカテゴリー注釈の追加の列を伴う。 1世紀以上にわたるこの象徴的な視覚化の継続的な開発は、すべてのバイオインフォマティクスディスプレイの中で最も広く使用されているディスプレイ方法の1つの基盤を提供している(Wilkinson and Friendly、2009)。統計計算にR言語を使用する場合(R Core Team、2016)、静的なヒートマップを作成するための多くのパッケージが存在する。たとえば、stats、gplots、heatmap3、fheatmap、pheatmapなどである。最近リリースされたパッケージでは、より複雑なレイアウトも可能である。これらには、gapma、 superheat、およびComplexHeatmapが含まれる(Gu et al、2016)。次の進化のステップは、インタラクティブクラスターヒートマップの作成であり、いくつかのソリューションが既に利用可能である。ただし、idendro Rパッケージ(Sieger et al、2017)などのこれらのソリューションは、多くの場合、研究者のパーソナルコンピューターでのみ探索できるインタラクティブな出力を提供することに焦点を当てている。共有可能なインタラクティブヒートマップを作成するためのソリューションがいくつかある。ただし、これらはXCMS Onlineなどの特定のオンラインプロバイダーに依存しているか、InCHlibなどのJavaScriptの知識が必要である。実際には、学術雑誌に掲載する場合、読者には静的な数字のみが残される(多くの場合、pngまたはpdf形式)。

 このギャップを埋めるために、インタラクティブクラスターヒートマップを含む共有可能なHTMLファイルを簡単に作成するためのheatmaply Rパッケージを開発した。対話機能は、次のコマンドを実行した後、ユーザーのデータに基づいて生成されるクライアント側のJavaScriptコードに基づいている。
 HTMLファイルには、ユーザーがセルの上にマウスを移動したときに値を表示できるだけでなく、ズームインできるようにする、出版用のインタラクティブな図が含まれている。この自己完結型のHTMLファイルは、研究者のホームページにアップロードするか、雑誌のサーバーの補足資料として興味のある読者が利用できるようにすることができる。同時に、このインタラクティブな図は、RStudioのビューアペインに表示したり、Shinyアプリケーションに含めたり、knitr / RMarkdown HTMLドキュメントに埋め込むことができる。

 このペーパーの残りの部分では、効果的なクラスターヒートマップ視覚化を作成するためのガイドラインを提供する。論文図1は、プロジェクトTychoのデータに関するこのセクションの提案を示している(van Panhuis et al、2013)。オンラインの補足資料には、インタラクティブなバージョンと、実世界の生物学データでのパッケージの使用例が含まれている。

  

 

CRAN manual

https://cran.r-project.org/web/packages/heatmaply/heatmaply.pdf

 

インストール

macos10.14のRstudioを使ってテストした。

依存

  • Depends R (>= 3.0.0), plotly (>= 4.7.1), viridis

Imports ggplot2 (>= 2.2.0), dendextend (>= 1.12.0), magrittr (>= 1.0.1), reshape2, scales, seriation, utils, stats, grDevices, methods, colorspace, RColorBrewer, htmlwidgets, webshot, assertthat, egg

本体 Github

#CRAN
> install.packages('heatmaply')

 >help(heatmaply)

heatmaply              package:heatmaply               R Documentation

 

Cluster heatmap based on plotly

 

Description:

 

     An object of class heatmapr includes all the needed information

     for producing a heatmap. The goal is to separate the

     pre-processing of the heatmap elements from the graphical

     rendering of the object, which could be done

 

     (Please submit an issue on github if you have a feature that you

     wish to have added)

 

     heatmaply_na is a wrapper for `heatmaply` which comes with

     defaults that are better for exploring missing value (NA)

     patterns. Specifically, the grid_gap is set to 1, and the colors

     include two shades of grey. It also calculates the is.na10

     automatically.

 

     heatmaply_cor is a wrapper for `heatmaply` which comes with

     defaults that are better for correlation matrixes. Specifically,

     the limits are set from -1 to 1, and the color palette is RdBu.

 

Usage:

 

     heatmaply(x, ...)

     

     heatmaply_na(x, grid_gap = 1, colors = c("grey80", "grey20"), ...)

     

     heatmaply_cor(x, limits = c(-1, 1), colors = cool_warm, ...)

     

     ## Default S3 method:

     heatmaply(

       x,

       colors = viridis(n = 256, alpha = 1, begin = 0, end = 1, option = "viridis"),

       limits = NULL,

       na.value = "grey50",

       row_text_angle = 0,

       column_text_angle = 45,

       subplot_margin = 0,

       cellnote = NULL,

       draw_cellnote = !is.null(cellnote),

       cellnote_color = "auto",

       cellnote_textposition = "middle right",

       cellnote_size = 12,

       Rowv,

       Colv,

       distfun = dist,

       hclustfun = hclust,

       dist_method = NULL,

       hclust_method = NULL,

       distfun_row,

       hclustfun_row,

       distfun_col,

       hclustfun_col,

       dendrogram = c("both", "row", "column", "none"),

       show_dendrogram = c(TRUE, TRUE),

       reorderfun = function(d, w) reorder(d, w),

       k_row = 1,

       k_col = 1,

       symm = FALSE,

       revC,

       scale = c("none", "row", "column"),

       na.rm = TRUE,

       row_dend_left = FALSE,

       margins = c(NA, NA, NA, NA),

       ...,

       scale_fill_gradient_fun = NULL,

       grid_color = NA,

       grid_gap = 0,

       srtRow,

       srtCol,

       xlab = "",

       ylab = "",

       main = "",

       titleX = TRUE,

       titleY = TRUE,

       hide_colorbar = FALSE,

       key.title = NULL,

       return_ppxpy = FALSE,

       row_side_colors,

       row_side_palette = NULL,

       col_side_colors,

       col_side_palette = NULL,

       ColSideColors = NULL,

       RowSideColors = NULL,

       seriate = c("OLO", "mean", "none", "GW"),

       heatmap_layers = NULL,

       side_color_layers = NULL,

       branches_lwd = 0.6,

       file,

       width = NULL,

       height = NULL,

       long_data,

       plot_method = c("ggplot", "plotly"),

       label_names = NULL,

       fontsize_row = 10,

       fontsize_col = 10,

       cexRow,

       cexCol,

       subplot_widths = NULL,

       subplot_heights = NULL,

       colorbar_len = 0.3,

       colorbar_xanchor = if (row_dend_left) "right" else "left",

       colorbar_yanchor = "bottom",

       colorbar_xpos = if (row_dend_left) -0.1 else 1.1,

       colorbar_ypos = 0,

       showticklabels = c(TRUE, TRUE),

       dynamicTicks = FALSE,

       grid_size = 0.1,

       node_type = "heatmap",

       point_size_mat = NULL,

       point_size_name = "Point size",

       label_format_fun = function(...) format(..., digits = 4),

       labRow = NULL,

       labCol = NULL,

       custom_hovertext = NULL,

       col = NULL,

       dend_hoverinfo = TRUE,

       side_color_colorbar_len = 0.3

     )

     

     ## S3 method for class 'heatmapr'

     heatmaply(

       x,

       colors = viridis(n = 256, alpha = 1, begin = 0, end = 1, option = "viridis"),

       limits = NULL,

       na.value = "grey50",

       row_text_angle = 0,

       column_text_angle = 45,

       subplot_margin = 0,

       row_dend_left = FALSE,

       margins = c(NA, NA, NA, NA),

       ...,

       scale_fill_gradient_fun = scale_fill_gradientn(colors = if (is.function(colors))

         colors(256) else colors, na.value = na.value, limits = limits),

       grid_color = NA,

       grid_gap = 0,

       srtRow,

       srtCol,

       xlab = "",

       ylab = "",

       main = "",

       titleX = TRUE,

       titleY = TRUE,

       hide_colorbar = FALSE,

       key.title = NULL,

       return_ppxpy = FALSE,

       draw_cellnote = FALSE,

       cellnote_color = "auto",

       cellnote_textposition = "middle right",

       cellnote_size = 12,

       row_side_colors = x"row_side_colors",

       row_side_palette = NULL,

       col_side_colors = x"col_side_colors",

       col_side_palette = NULL,

       plot_method = c("ggplot", "plotly"),

       ColSideColors,

       RowSideColors,

       heatmap_layers = NULL,

       side_color_layers = NULL,

       branches_lwd = 0.6,

       label_names = c("row", "column", "value"),

       fontsize_row = 10,

       fontsize_col = 10,

       subplot_widths = NULL,

       subplot_heights = NULL,

       colorbar_xanchor = if (row_dend_left) "right" else "left",

       colorbar_yanchor = "bottom",

       colorbar_xpos = if (row_dend_left) -0.1 else 1.1,

       colorbar_ypos = 0,

       colorbar_len = 0.3,

       showticklabels = c(TRUE, TRUE),

       dynamicTicks = FALSE,

       node_type = c("scatter", "heatmap"),

       grid_size = 0.1,

       point_size_mat = x"matrix""point_size_mat",

       point_size_name = "Point size",

       label_format_fun = function(...) format(..., digits = 4),

       custom_hovertext = x"matrix""custom_hovertext",

       dend_hoverinfo = TRUE,

       side_color_colorbar_len = 0.3

     )

     

長くなるのでargumentsは省略。CRANのヘルプを参照して下さい。

 

 

テストラン

library("heatmaply")
heatmaply(mtcars)

webページとして保存できる。

f:id:kazumaxneo:20200227003307p:plain

 

 

実行方法 

1、heatmaplyのロード

library("heatmaply")

 

2、データの読み込み

 ここでは以下のANIの総当たり比較結果のデータを読み込む。

  Genome1 Genome2 Genome3
Genome1 100 95 98
Genome2 95 100 92
Genome3 98 92 100

コピペ読み込み。上の表を選択してコピーしておく。以下のコマンドを実行(コードをコピペするなら実行直前に上の表をコピーすること)。

#mac
x <- read.table(pipe("pbpaste"), header = TRUE)

#windows
x <- read.table("clipboard"), header = TRUE)

またはファイルinput.tsvに保存し、データフレームに読み込む。  

x <- read.table("input.tsv", header=T, sep="\t")

#ファイルがカレントにないならフルパス指定. セパレータがコンマならsep=","
x <- read.table("/home/kazu/data/input.csv", header=T, sep=",")

 

3、ヒートマップで視覚化

heatmaply(x)

#font sizeのみ調整(default=10)
heatmaply(x, fontsize_row = 20, fontsize_col = 20)

Rstudioやjupyter notebookで実行していない場合、デフォルトブラウザに直接出力される。 

f:id:kazumaxneo:20200229160304p:plain

 

 

heatmaply_corはheatmaplyのラッパーで、相関マトリックスに適したデフォルト設定、 具体的には、値が-1〜1に制限され、カラーパレットはRdBu(RColorBreweのカラー設定)になる。

heatmaply_cor(
cor(mtcars),
xlab = "Features",
ylab = "Features",
k_col = 2,
k_row = 2
)
  • k_col    an integer scalar with the desired number of groups by which to color the den- drogram’s branches in the columns (uses color_branches) If NA then find_k is used to deduce the optimal number of clusters. 
  • k_row    an integer scalar with the desired number of groups by which to color the den- drogram’s branches in the rows (uses color_branches) If NA then find_k is used to deduce the optimal number of clusters.
  • xlab    A character title for the x axis.
  • ylab    A character title for the y axis.
  • main    A character title for the heatmap.

 

f:id:kazumaxneo:20200227014224p:plain

他にも様々な例があります。上に載せたリンクから確認して下さい。

 

 

補足

shinyHeatmaply(ファーストオーサーの友人が開発と書かれています)も公開されています。簡単に使えるので興味があれば試してみて下さい。ただし、まだ開発中であり、環境によっては不安定な事もあるようです。

 

メモ

2000 x 2000 程度なら数分待てば動作する。

f:id:kazumaxneo:20200227104149p:plain

引用

heatmaply: an R package for creating interactive cluster heatmaps for online publishing
Tal Galili, Alan O’Callaghan, Jonathan Sidi, Carson Sievert
Bioinformatics, Volume 34, Issue 9, 01 May 2018, Pages 1600–1602

 

関連

 

(ヒトゲノム)個人のサンプルが汚染または交換されている可能性があるかどうかを調べる verifybamid

 

 DNAサンプル汚染の検出と推定は、高品質の遺伝子型コールと信頼性の高いダウンストリーム分析を確保するための重要なステップである。既存の方法は、汚染率の正確な推定のために母集団対立遺伝子頻度情報に依存している。シーケンス解析の初期段階で各個体の集団対立遺伝子頻度を正しく指定することは、多様な集団にわたる複数の研究からのサンプルを同時に処理する大規模なシーケンシングセンターにとって非現実的または不可能ですらある。一方、対立遺伝子の頻度を誤って指定すると、推定汚染率にかなりの偏りが生じる可能性がある。たとえば、既存の方法では、遺伝的祖先が誤って指定されている場合、一般的な3%の汚染除外しきい値で10%の汚染サンプルを特定できないことがよくある。汚染されたサンプルのこのような不完全なスクリーニングは、深くシーケンシングされたゲノムおよびエクソームにおいてさえ、遺伝子型決定エラー率を大幅に増大させる。

 DNA汚染を正確に推定し、目的のサンプルまたは汚染サンプルの遺伝的祖先にとらわれない堅牢な統計的手法を提案する。この方法は、参照遺伝子型から主成分座標に投影された個体固有の対立遺伝子頻度を活用することにより、統一された尤度フレームワークで遺伝的祖先とDNA汚染の推定を統合する。この方法が、さまざまな集団および汚染率における汚染率を堅牢かつ正確に推定することを実証する。さらに、汚染が存在する場合、汚染が無視されると遺伝的祖先の定量的推定値(主成分座標など)が大幅に偏ることがあり、提案された方法がこの偏りを補正することを実証する。このメソッドはhttp://github.com/Griffan/verifyBamIDで公開されている。

 

wiki

https://genome.sph.umich.edu/wiki/VerifyBamID

 

インストール

ubuntu18.04のdocker環境でテストした(ホストOSはmacos10.14)。

本体 Github

#bioconda (link)
conda install -c bioconda -y verifybamid

verifyBamID

# verifyBamID 

verifyBamID 1.1.3 -- verify identity and purity of sequence data

(c) 2010-2014 Hyun Min Kang, Goo Jun, and Goncalo Abecasis

 

 

Available Options

                             Input Files : --vcf , --bam , --bai ,

                                           --subset , --smID

                    VCF analysis options : --genoError [1.0e-03],

                                           --minAF [0.01],

                                           --minCallRate [0.50]

   Individuals to compare with chip data : --site, --self, --best

          Chip-free optimization options : --free-none, --free-mix [ON],

                                           --free-refBias, --free-full

          With-chip optimization options : --chip-none, --chip-mix [ON],

                                           --chip-refBias, --chip-full

                    BAM analysis options : --ignoreRG, --ignoreOverlapPair,

                                           --noEOF, --precise, --minMapQ [10],

                                           --maxDepth [20], --minQ [13],

                                           --maxQ [40], --grid [0.05]

                 Modeling Reference Bias : --refRef [1.00], --refHet [0.50],

                                           --refAlt [0.00]

                          Output options : --out , --verbose

                               PhoneHome : --noPhoneHome,

                                           --phoneHomeThinning [50]

 

 

FATAL ERROR - 

--vcf [vcf file] required

 

 

 

実行方法

ランには外部遺伝子型または対立遺伝子頻度情報を含むVCFファイルとBAMファイルが必要。VCFには常染色体のSNP情報のみ記載されている必要がある(MNPも不可)。VCFのIDとBAMのIDは一致していなければならない。

verifyBamID --vcf input.vcf --bam input.bam --out output --verbose --ignoreRG

 

 

引用

Ancestry-agnostic estimation of DNA sample contamination from sequence reads

Fan Zhang, Matthew Flickinger, Sarah A. Gagliano Taliun, InPSYght Psychiatric Genetics Consortium, Gonçalo R. Abecasis, Laura J. Scott, Steven A. McCaroll, Carlos N. Pato, Michael Boehnke, Hyun Min Kang

Genome Research. Published in Advance January 24, 2020