macでインフォマティクス

macでインフォマティクス

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

エラーの多いロングリードをタンデムリピートにマッピングする tandemmapper

 

 タンデムリピートは、不均等なクロスオーバーによってしばしば生成される複数の連続するほぼ同一のシーケンスによって形成される(Smith、1976)。初期のDNAシーケンスプロジェクトから、タンデムリピートが真核生物ゲノムに豊富にあることが明らかになった(Yunis and Yasmineh、1971; Bacolla et al、2008)。タンデムリピートの最近の研究は、さまざまな細胞プロセスにおける役割を明らかにし、タンデムリピートの突然変異が遺伝的障害につながる可能性があることを実証した(McFarland et al、2015; Giunta and Funabiki、2017; Song et al、2018; Black et al、 2018)。
 広範囲に研究されたショートタンデムリピート(Willems et al。2014; Gymrek et al。、2016; Saini et al。、2018)と長さが数万から数百万のウルトラロングタンデムリピート(ETR)を区別する。 ETRはアセンブリが難しいため、それらの大部分は、他の種はもちろんのこと、ヒトゲノムでもアセンブリされていないままである。セントロメアおよびペリセントロメアには、ヒトゲノムの約3%を占め、メガベース長の領域に及ぶ最長のETRが含まれている(Miga、2019)。それらは、これまでに配列を決定するすべての試みを回避したヒトゲノムの「ダークマターであり、リファレンスヒトゲノムの最大のギャップである(Hayden et al、2013; Miga et al、2019)。
パシフィックバイオサイエンス(PacBio)やオックスフォードナノポアテクノロジー(ONT)などのロングシーケンシング技術の出現により、全ゲノムシーケンスの状況が大きく変わった。ロングリードアセンブラ(Chin et al、2016; Koren et al、2017; Kolmogorov et al、2019; Ruan and Li、2019)およびロングリードとショートリードを組み合わせたハイブリッドアセンブラ(Antipov et al、 2016; Zimin et al、2017)は、ショートリードアセンブリと比較して、アセンブリされたゲノムの連続性を大幅に向上させた。さらに、ロングリードは、ヒト動原体を再構築するための半手動アプローチの成功に貢献した(Jain et al、2018a; Miga et al、2019)。 Flyeアセンブラは、ロングリードにまたがるブリッジされたタンデムリピート、およびロングリードにまたがらないいくつかの非ブリッジタンデムリピートを正常に解決する(Kolmogorov et al、2019)。 centroFlyeアセンブラ(BzikadzeおよびPevzner、2019)は、セントロメアなどの非ブリッジETRを自動的にアセンブルするように設計された。
ETRアセンブリのさまざまな代替戦略と、これらのアセンブリベンチマークのグラウンドトゥルースがないため、品質評価の問題が生じる。

(一部略)

 既存のアセンブリ品質評価ツールは、リードをアセンブリに正確にマッピングするアライナーに依存している(Li and Durbin、2009; Langmead et al。、2009; Li、2016; Li、2018)。ただし、ベンチマークでは、これらのツールはETRで頻繁に失敗することが明らかになった。たとえば、minimap2(Li、2018)は、特にアセンブリエラーのある領域で、ETRへの一部のリードのアライメントが正しくないことを示す。そのため、エラーが発生しやすいロングリードをETRに効率的にマッピングするtandemMapperツールを開発した。 TandemMapperは、tandemQUAST開発を可能にしただけでなく、より正確なリードマッピングとその後のポリッシュによりETRアセンブリの改善にもつながる。
(一部略)
 霊長類の動原体は、レトロトランスポゾンリピートとATリッチアルファサテライト、171 bpモノマーに基づくDNAリピートで構成されている(Manuelidis and Wu、1978)。人間と多くの霊長類では、連続したモノマーがタンデムに並んで高次リピート(HOR)ユニットに配置される(Willard and Waye、1987a)。 HOR内のモノマーの数とその順序は、染色体に固有である。たとえば、DXZ1と呼ばれる染色体X HORは、12個のモノマーで構成されている(Willard and Waye、1987b)。これらの12個のモノマーは、先祖の五量体サテライトリピートABCDEから進化したものであり、「C1D1E1」として表すことができる。ここでは、リードをETRにマッピングするためのツールであるtandemMapperと、ETRアセンブリを評価および改善するためのツールであるtandemQUASTを紹介する。

 

 

インストール

ubuntu18.04LTSでテストした。

本体 Github

git clone https://github.com/ablab/tandemQUAST.git
cd tandemQUAST/

#ここでは仮想環境に導入
conda create -n tandemQUAST -y
conda activate tandemQUAST
conda install --file requirements.txt

 > python tandemquast.py --help

$ python tandemquast.py --help

Usage: tandemquast.py [OPTIONS] [ASSEMBLY_FNAMES]...

 

Options:

  -l TEXT         Comma separated list of assembly labels

  -r PATH         File with reads  [required]

  --hi-fi PATH    File with HiFi reads (optional)

  -m PATH         Monomer sequence

  -o PATH         Output folder

  -t INTEGER      Threads

  --only-polish   Run polishing only

  -f, --no-reuse  Do not reuse old files

  --help          Show this message and exit.

python tandemmapper.py --help

$ python tandemmapper.py --help

Usage: tandemmapper.py [OPTIONS] [ASSEMBLY_FNAMES]...

 

Options:

  -l TEXT         Comma separated list of assembly labels

  -r PATH         File with reads  [required]

  --hi-fi PATH    File with PacBio HiFi reads (optional)

  -o PATH         Output folder

  -t INTEGER      Threads

  -f, --no-reuse  Do not reuse old files

  --help          Show this message and exit.

 

 

実行方法

tandemmapper.py -r reads_file -o output_dir assembly_file1 assembly_file2

 

引用

TandemMapper and TandemQUAST: mapping long reads and assessing/improving assembly quality in extra-long tandem repeats

Alla Mikheenko, Andrey V. Bzikadze, Alexey Gurevich, Karen H. Miga, Pavel A. Pevzner

bioRxiv preprint first posted online Dec. 23, 2019

 

UGENE その2

今回は、1回目で説明できなかった機能について説明する。

 

一旦入力した配列は直接編集できないようになっている。編集するには左端のeditボタンをクリックする。

f:id:kazumaxneo:20200107040909p:plain

選択した配列を消したり、追加できるようになっている。

f:id:kazumaxneo:20200107040955p:plain

編集が終わったらもう一度editボタンをクリックしてロックする。

 

複数ライン表示が見にくければ、1行表示に切り替える。

f:id:kazumaxneo:20200107041732p:plain

 

ORFを表示。

f:id:kazumaxneo:20200107041912p:plain

 

さらにTranslatee selectionにチェックを付け、

f:id:kazumaxneo:20200107041949p:plain

 

配列を囲むと、その領域のアミノ酸配列が素早く確認できる。

f:id:kazumaxneo:20200107042218p:plain

Show all frame

f:id:kazumaxneo:20200107042146p:plain

 

 

Find pattern - Smith–Watermanで相同な領域をサーチ。

f:id:kazumaxneo:20200107033015p:plain

配列を入力して検索する。

f:id:kazumaxneo:20200107034013p:plain

パラメータ、アラインメント感度、ヒット後のアノテーションのカテゴリなどを選択する。


ヒットした領域がアノテーション付けられ、ハイライト表示される。複数ヒットした場合、下のウィンドウのフィーチャ一覧を展開して確認できる。

f:id:kazumaxneo:20200107034105p:plain



build dotplot - 2配列間の相同な領域の視覚化

f:id:kazumaxneo:20200107034709p:plain

 

比較する配列を指定する。multi-fastaを読み込んでいるなら改めてロードせずとも各々の配列を選択可能。

f:id:kazumaxneo:20200107034803p:plain

 

視覚化された。dot plotの図を囲むとその領域の配列が選択される。

f:id:kazumaxneo:20200107034810p:plain

そのままcommand + Cすることで選択した領域の配列をコピーできる。

 

Query NCBI BLAST database - ネット越しのblast検索

f:id:kazumaxneo:20200107035159p:plain

データベースとパラメータを指定する。

f:id:kazumaxneo:20200107035204p:plain

 

primer3 - プライマーデザイン

f:id:kazumaxneo:20200108012336p:plain

パラメータを指定する。

f:id:kazumaxneo:20200108012312p:plain

 

条件を満たすプライマーが自動作成される。

f:id:kazumaxneo:20200108012443p:plain

近いうちにNGSのデータを扱う流れについてもまとめます。

引用

Unipro UGENE: a unified bioinformatics toolkit
Okonechnikov K, Golosova O, Fursov M; UGENE team
Bioinformatics, Volume 28, Issue 8, 15 April 2012, Pages 1166–1167

 


DNA解析ソフト4 次世代シークエンシングデータも扱える Unipro UGENE その1

2020 1/6 タイトル修正

2020 3/2 わかりにくい説明を修正

 

明けましておめでとうございます。今年もよろしくお願い致します。

2020年初回はDNA解析ソフトUGENEを紹介します。発表はかなり前ですが、今でもアップデートが続いており、塩基配列の編集のみならず、サンガーの波形データの表示やNGSのデータの操作も可能になっています。

 

 *1

 Unipro UGENEは、バイオインフォマティクスの専門知識があまりない分子生物学者がデータを管理、分析、視覚化できるよう支援することを主な目的とするマルチプラットフォームオープンソースソフトウェアである。 UGENEは、広く使用されているバイオインフォマティクスツールを共通のユーザーインターフェイスに統合する。このツールキットは複数の生物学的データ形式をサポートし、リモートデータソースからデータを取得できる。アノテーション付きゲノムシーケンス、次世代シーケンス(NGS)アセンブリデータ、複数のシーケンスアラインメント、系統樹、3D構造などの生物学的オブジェクトの視覚化モジュールを提供する。統合されたアルゴリズムのほとんどは、マルチスレッドと特別なプロセッサ命令の使用により、最大のパフォーマンスが得られるように調整されている。 UGENEには、ローカルリソースまたはHigh Performance Computing(HPC)環境で起動できる再利用可能なワークフローを作成するためのビジュアル環境が含まれている。 UGENEはQtフレームワークを使用してC ++で記述されている。組み込みのプラグインシステムと構造化されたUGENE APIにより、ツールキットの新しい機能を拡張できる。UGENEバイナリは、MS WindowsLinux、およびMac OS Xhttp://ugene.unipro.ru/download.htmlから無料で入手できる。 UGENEコードはGPLv2の下でライセンスされている。統合ツールのコードライセンスと著作権に関する情報は、ソースバンドルで提供されるLICENSE.3rd_partyファイルに記載されている。

 

*2

 次世代シーケンス(NGS)テクノロジーの出現により、研究者に新しい可能性が開かれた。しかし、より生物学がよりデータ集約型フィールドになると、より多くの生物学者が複雑な計算ツールを使用してNGSデータを処理および分析する方法を学ぶ必要がある。一般的なパイプラインが利用できる場合でも、多くの場合、ベンチサイエンティストがパイプラインツールをインストールして構成するのは時間のかかる面倒な作業が必要になる。 NGSデータ分析ツールの統合された、デスクトップおよび生物学者に優しいフロントエンドは、この分野の生産性を大幅に向上させると考えられる。ここでは、Unipro UGENEデスクトップツールキットに統合されたNGSパイプライン「SAMtoolsを使用したバリアントコール」、「RNAシーケンスデータ分析用のTuxedoパイプライン」、「ChIPシーケンスデータ分析用のシストロームパイプライン」を紹介する。研究者がこれらのパイプラインを異なるデータセットで実行し、結果を保存および調査し、同じパラメーターでパイプラインを再実行するのに役立つ利用可能なUGENEインフラストラクチャについて説明する。これらのパイプラインツールは、UGENE NGSパッケージに含まれている。これらのパイプラインの個々のブロックは、専門ユーザーが独自の高度なワークフローを作成するためにも利用できる。

 

http://ugene.net HP

f:id:kazumaxneo:20200104205557p:plain

 

Quick Start Guide

https://ugene.net/wiki/display/QSG/Quick+Start+Guide

Documentation

http://ugene.net/documentation.html

YouTube channel

http://youtube.com/uniprougene

 


NGSデータ解析ワークフロー

Metagenomic workflows

DNA editing

 

インストール

All downloadsよりダウンロードする。ここではmacos10.15の.dmgイメージをダウンロードした。

 

ダウンロードした.dmgを実行する。

メタゲノム解析ツールも導入する場合、224GBほど追加のディスクスペースが必要になる。

f:id:kazumaxneo:20200101205007p:plain

CHIP-seqは23.7GB。

f:id:kazumaxneo:20200101205012p:plain

指定したツールがダウンロードされる。

f:id:kazumaxneo:20200101205643p:plain

メタゲノムモジュールを指定しているとダウンロードにかなりの時間を要する。

 

使い方

起動したところ

f:id:kazumaxneo:20191231132849p:plain

 

 

Open FilesでDNA配列を開く。ウィンドウ内にFASTA配列をペーストする。

f:id:kazumaxneo:20200101213507p:plain

直接FASTAファイルを叩いて開いてもOK。

f:id:kazumaxneo:20200101213722p:plain

 

開いたところ。左が配列名。multi-fastaならクロモソームごとに表示される。右下は配列のウィンドウ、右上は配列全体の俯瞰表示になる。

f:id:kazumaxneo:20200101214201p:plain

別のファイルの配列を同時に開いている場合、左のウィンドウや上のタブから切り替えする。

f:id:kazumaxneo:20200101214924p:plain

 

右上の俯瞰表示ウィンドウの領域を囲むと、その領域の配列がハイライト表示される。

f:id:kazumaxneo:20200101215210p:plain

 

俯瞰表示ウィンドウの領域をダブルクリックすると、右下の配列ウィンドウがその領域の配列に切り替わる。

f:id:kazumaxneo:20200101215517p:plain

 

右クリック => Go to positionで数値をして飛ぶこともできる。

f:id:kazumaxneo:20200101215658p:plain

 

ゲノム全体から目的の領域に飛ぶには、上のウィンドウでポジションをタイプするのが早い。120000にジャンプ。

f:id:kazumaxneo:20200106021520p:plain

 

また、一番上のバーの部分をクリックしてドラッグすることで配列全体の高速スクロールが可能。

f:id:kazumaxneo:20200102221606p:plain

 

2020 3/2 補足

一番上のアイコンを使うとうまく移動できる。をスクロールすると、下半分の配列のATGC表示部分がスクロールする。四角の半透明灰色ボックス◽️を動かすと上半分のフィーチャがスクロールする。下半分ウィンドウと上半分ウィンドウは独立しているので、両方のウィンドウで同じ領域を示すにはこの両方のアイコンを移動させる必要がある。

f:id:kazumaxneo:20200302113647p:plain

 

ちなみに、アイコンの配列の領域は、下半分ウィンドウと上半分ウィンドウが同じ領域にあるときは上のウィンドウでも表示される。を半透明灰色ボックス◽️の中に収まるように移動してやると、上のウィンドウ内に破線でボックスが表示される(下の画像のGCプロット中央付近の縦長破線ボックスがの配列表示部分)。

f:id:kazumaxneo:20200302114224p:plain

 

半透明灰色ボックス◽️を正確に動かすには左端のアイコンを使う。左端のアイコンの一番下にあるボタンをクリックすると半透明灰色ボックス◽️ の表示領域を指定できる。

f:id:kazumaxneo:20200302114548p:plain

他は半透明灰色ボックス◽️の拡大/縮小などになる。

 

このように、染色体クラスの大きな配列もできるだけ直感的に扱えるように工夫されている。

 

 

次にボタンについて説明していく。

ORFをクリックすると、3frameでORFが表示される。

f:id:kazumaxneo:20200101220032p:plain

 

ORFの下のRestriction sitesをクリックすると制限酵素部位が表示される。

f:id:kazumaxneo:20200101220154p:plain

 

上のカラフルな円グラフのようなマークを押すと環状表示になる。

f:id:kazumaxneo:20200101221333p:plain

 環状表示にして制限酵素サイトを表示。

f:id:kazumaxneo:20200101221623p:plain

アイコンが見えない場合、右上の>>マークをクリックして隠れたボタンを表示する。

 

線形表示に戻すには、右上の>>マークから隠れたボタンを表示、円マークをクリック。

f:id:kazumaxneo:20200101221357p:plain

 

カラフルな円マークボタンの右隣のボタンをクリックすると、GCとやATプロットが表示される。

f:id:kazumaxneo:20200106021007p:plain



 

 右端に5つあるボタンにも便利な機能がアサインされている。ボタンの一番上は配列検索ボタンになる。

f:id:kazumaxneo:20200101222607p:plain

 

パーフェクトマッチのほか、不完全マッチの配列も検索できる。

f:id:kazumaxneo:20200101225907p:plain

 

上から2つ目のボタンは制限酵素認識部位やメチル化部位などアノテーションの色指定。

f:id:kazumaxneo:20200101222447p:plain

 

上から3つ目のボタンは選択した配列のGCやTmなどの要約統計。

f:id:kazumaxneo:20200102193757p:plain

 

上から4つ目のボタンは選択したプライマー配列を使ったinsilico PCR

f:id:kazumaxneo:20200102194502p:plain

 

配列ウィンドウからプライマーにする配列を選択してcommnad + Cでコピー、Forwardとreverseにそれぞれ貼り付ける。

f:id:kazumaxneo:20200102194646p:plain

 

Revserseは貼り付けてからウィンドウ右下のreverse compボタン(C⇄)を押して逆相補鎖にする。

f:id:kazumaxneo:20200102194712p:plain

 

Forwardについて赤字でlow Tmの警告が出ている。必要に応じて配列を修正する。また最大増幅サイズなどのパラメータを指定する。それから下のFind productsボタンをクリック。ここではForward primerのTmが低い警告が出ており、下のボタンもFind products anywayとなっているが、気にせず実行。

f:id:kazumaxneo:20200102195116p:plain

 

増幅部位のサイズとTaが表示され、配列がハイライト表示される。サイズとTa部分をダブルクリックすると別ウィンドウにジャンプする。

f:id:kazumaxneo:20200102195728p:plain

 

上に並んでいるそのほかのボタンは、メニューの表示/非表示、領域の配列の画像出力、reverse comp、copyなどになります。明日は残りのボタンについて簡単に説明します。  

追記


引用

1   Unipro UGENE: a unified bioinformatics toolkit
Okonechnikov K, Golosova O, Fursov M; UGENE team
Bioinformatics, Volume 28, Issue 8, 15 April 2012, Pages 1166–1167

 

2  Unipro UGENE NGS pipelines and components for variant calling, RNA-seq and ChIP-seq data analyses
Golosova O, Henderson R, Vaskin Y, Gabrielian A, Grekhov G, Nagarajan V, Oler AJ, Quiñones M, Hurt D, Fursov M1, Huyen Y

PeerJ. 2014 Nov 4;2:e644. doi: 10.7717/peerj.644. eCollection 2014

 

インタラクティブなGOエンリッチメント解析webサービス ShinyGO

2021 6/5 バージョンアップのツイート

 

2019年最後の投稿になります。今年もお世話になりました。来年もどうぞよろしくお願い致します。

 

 ゲノムワイドな研究で特定された一連の遺伝子について、遺伝子オントロジー(GO)によって定義されたものなど、特定のパスウェイまたは機能カテゴリの遺伝子がエンリッチされているかどうかを調べるためエンリッチメント解析を行うことができる(Ashburner et al 、2000)。エンリッチメント解析用に多数のツールが開発された(Khatri et al、2012)。これらのツールの小さなサブセットがSupp table.1にリストされている。一部のツールは生物医学研究用に設計されているため、主にヒトとマウスに焦点を当てている。たとえば、Enrichr(Kuleshov et al、2016)には、GO、共発現、組織特異的遺伝子、転写因子(TF)またはmicroRNA(miRNA)標的遺伝子からさまざまなパスウェイデータベースに至るまでの遺伝子セットが含まれている。同様に、PlantGSEAなどのツールは15の植物種に焦点を当てている(Yi、et al、2013)。
g:Profilerは、300を超える植物および動物種のEnsemblの遺伝子アノテーション(Aken et al、2017)に基づいている。 STRING(Szklarczyk et al、2015)は、タンパク質間相互作用(PPI)の大規模なデータベースである。また、5090種(バージョン11)のGOおよびタンパク質ドメインのエンリッチメント解析機能も提供する。もう一方の極端な例は、DAVID(Huang da et al、2009)であり、NCBI、UniProt、KEGG、GO、Biocarta、REACTOMEなどの多くのソースから情報を取得するため、数万の生物の大規模なデータベースを含めることができる。これらのツールは、生物学者が遺伝子リストから洞察を得るのに役立つ。
 視覚化と統計分析のための多くの強力なRパッケージへのアクセスを可能にするShinyフレームワークを利用して、Ensemblとパスウェイのアノテーションデータベースに基づく新しいツールを開発した。 ShinyGOのユニークな機能には、1)KEGGおよびSTRINGへのAPIアクセスに基づいて、パスウェイダイアグラムおよびタンパク質-タンパク質相互作用(PPI)ネットワークにクエリ遺伝子を表示する、2)階層的クラスタリングおよび対話型ネットワークを使用してエンリッチされたパスウェイ間のオーバーラップを視覚化する、および3)遺伝子の種類、長さ、GC含有量、クエリ遺伝子とバックグラウンドの間のクロモソーム分布の統計的に有意な違いを特定する。いくつかのモデル生物について、他の遺伝子セット、特にTFおよびmiRNA標的遺伝子を組み込むことにより、アノテーションデータベースも強化した。
 詳細については、補足ファイルを参照。ソースコードidep/shinyapps/go at master · iDEP-SDSU/idep · GitHubから入手できる。

 

2023/05/03

2022/12/14

 

 

webサービス

http://bioinformatics.sdstate.edu/go/にアクセスする。

f:id:kazumaxneo:20191228231604p:plain

iDEPとほぼ同じインターフェイスなので、iDEPに慣れているとすぐ使い始めることができる。

================参考================

iDEP webサーバーのDEG2タブからもリンクしている。

f:id:kazumaxneo:20191229113725p:plain

左下のShinyGOからリンクしている(直接データが読み込まれるわけではない)。

=====================================================

 

Speciesを選択し、遺伝子リストをペーストする。ここではDemo genesをペーストした。Speciesを選択しないと、ベストマッチの生物種が画面左下に表示される。

f:id:kazumaxneo:20191228231706p:plain

p-value cutoff値と表示する遺伝子数を指定してSubmitをクリック。

 

iDEPで発現変動遺伝子群を絞り込んでいるなら、DEG1タブのGene listsをダウンロードする。

f:id:kazumaxneo:20200109200425p:plain

(iDEP webサービスより)

ダウンロードしたファイルは1行目コメント、2行目誘導遺伝子群、3行目抑制遺伝子群とそれぞれ遺伝子名が横一列にずらっと並んでいるが、横1列全てコピーし、そのままShinyGOのウィンドウにペーストすればOK。

 

その場ですぐにエンリッチされた機能カテゴリが表示される。

f:id:kazumaxneo:20191228233116p:plain

 

左のタブからカテゴリを変更可能。変更するとすぐさま更新される。

f:id:kazumaxneo:20191229111554p:plain

 

REACTOME、pfam、interpro、KEGGなどGOのデータベース以外も選べる。素晴らしい!

f:id:kazumaxneo:20200410183529p:plain

 

hierarchical clustering tree

f:id:kazumaxneo:20191228233121p:plain

Network

f:id:kazumaxneo:20191228233129p:plain

 

クロモソーム毎の遺伝子数、遺伝子タイプ(IncRNA, miRNA, snRNA...)などのグラフf:id:kazumaxneo:20191229104532p:plain

 

locus 

f:id:kazumaxneo:20191229104818p:plain

プロモーター配列

f:id:kazumaxneo:20191229104824p:plain


STRINGとの連携によるprotein-protein network分析。

f:id:kazumaxneo:20191229104830p:plain


PPI network of DEGs

f:id:kazumaxneo:20191229111110p:plain

Click here~からSTRINGにジャンプできる。

f:id:kazumaxneo:20191229111127p:plain

f:id:kazumaxneo:20191229111202p:plain

 

右端の?タブがhelpになっており、データソースなど確認できる。

f:id:kazumaxneo:20191229105747p:plain

 

引用

ShinyGO: a graphical enrichment tool for animals and plants
Steven Xijin Ge, Dongmin Jung, Runan Yao
Bioinformatics, btz931, Published: 27 December 2019

 

関連

 

ラージゲノムにもスケールする高速且つ精度の高いドラフトゲノムポリッシャー hypo

  

 DNAシーケンサーによって生成されたフラグメント(リード)からゲノムを再構築するゲノムアセンブリと、種間または種内の遺伝的変異を調べるためのその解析は、ゲノミクスの中心である。 Pacific Biosciences(PacBio)やOxford Nanopore Technologies(ONT)などの第3世代(DNA)シーケンサーTGS)は、ゲノムの高品質なアセンブリと分析を可能にすることで、ゲノミクスに新しい刺激を与えた。 TGSは、リード長が数百塩基対に制限されていることから生じる第2世代または次世代シーケンサー(NGS; Illuminaなど)の主な制限を克服し、平均長が数万塩基対のリードを生成する。連続的なアセンブリ、よりrepetitiveなエレメントの解決、およびより大きな構造変異の啓示など[Roberts et al、2013、Lee et al、2016]。ごく最近まで、ヒトゲノムのde novoアセンブリは、ナノポアからのウルトラロングリードを使用して行われ(品質を改善する他の補完的な技術とともに)、ヒトリファレンスゲノム(GRCH38)の連続性を超えるだけでなく、テロメアからテロメアへの完全な染色体Xを再構成した [Miga et al、2019]。ただし、99%を超える精度を持つNGSのショートリードとは対照的に、ロングリードにはエラー率が高い(> 10%)という欠点がある[Weirather et al、2017、Jain et al、2018]。さらに、ロングリードのエラープロファイルは、置換よりも挿入-欠失( indel)に偏っている。それらのうちホモポリマーのindelがより顕著である[Weirather et al、2017]。
 ノイズの多いロングリードに対処するために、アセンブラは通常、アセンブリの前にエラーを修正することに頼る。これは、特に大きなゲノムの場合、計算コストがかかる[Sovićet al、2016、Fu et al、2019、Zhang et al、2019]。さらに、アセンブリされたコンティグは通常、最後にコンセンサス生成を使用することでさらにアセンブリの品質を向上させる。最近、いくつかの高速アセンブラ(miniasm [Li、2016]、Ra [Vaser andŠikić、2019]、wtdbg2 [Ruan and Li、2019])が利用可能になった。ベースレベルのエラーは多くなるが(修正されたリードを使用する他のアセンブラに比べての約10倍のエラー)1桁早く動作する。これらの高速アセンブラは、エラー修正をpolishのみに依存している。重要なことに、ロングリードアセンブリでは、エラー率が高く indelがドミナントであるため、エラーを修正してタンパク質予測に重大な影響を与えないようにすることが重要である[Watson and Warr、2019]。したがって、ポリッシングツールは、正確で長いリードのアセンブリ、特に高速でエラー修正のないアセンブラによって生成されるアセンブリを作成する上で重要な役割を果たす。
 ポリッシャーは大まかに「Sequencer-bound」と「General」に分類できる。シーケンサーにバインドされたポリッシャーは、特定のシーケンサーによって生成された生の信号レベル情報を必要とするため、特定のシーケンサーからのリードのみをポリッシュできる。 ONTの場合、NanoPolish [Loman et al、2015]およびその後継のMedaka [Nanopore Technologies、2019]はこのカテゴリに分類される。同様に、PacBioについては、Quiver [Chin et al、2013]およびその後継のArrow [Laird Smith et al、2016]が利用可能である。

 一方、Generalなポリッシャーは、どのようなシーケンサーによって生成されたリードでも処理できる十分な堅牢性を持つ。以前、Pilon [Walker et al、2014]は、細菌および小さな真核生物ゲノムに広く使用されている一般的なポリッシャーであった。しかしながら、Pilonは、Raconに取って代わられつついる(または組み合わせて使用​​されている)[Vaser et al。、2017]。Raconは超高速であるため、大規模なゲノム上でリソース的にうまくスケーリングできる。最近、いくつかの新しいポリッシャーが登場した:wtpoa-cns(wtdbg2のスタンドアロンコンセンサスモジュール)、ntEdit [Warren et al、2019]、Apollo [Firtina et al、2019]。 ntEdit以外の各ポリッシャーは、ドラフト(未修正)アセンブリのリードのアライメント情報に依存する。 Pilonは、ドラフトコンティグの各ベース位置でのリードからのベースのパイルアップに基づいている。 Raconおよびwtpoa-cnsはコンティグをより小さな断片に分割し、より高速且つより実用的にするため、POAのsingle instruction multiple data (SIMD)実装[Lee et al、2002、Lee、2003]を使う。 Apolloは、機械学習アプローチを展開して、ドラフトアセンブリのプロファイルHidden Markov Model(pHMM)[Firtina et al、2018]を構築し、それを使用してエラーを修正する。リードからアセンブリへのアライメントの活用とは対照的に、ntEditはドラフト内のkmers(長さkのシーケンス)のスキャンに基づいてエラーを修正し、リードにkmerを格納するBloom Filterを使用してそれらの有無をチェックする。
 一般的なポリッシャーにはそれぞれ制限がある。 PilonとntEditは、主に非常に正確なショートリードで動作するように設計されている。さらに、前述のように、Pilonはリソースの観点から大きなゲノムではうまくスケーリングできないが、Raconとwtpoa-cnsは大きなゲノムでもスケーラブルし高速であり、ショートリードとノイズの多いロングリードを処理できる。ただし、1回の実行で、両方のポリッシングに使用できるのは、ロングリードのみ、またはショートリードのみである。精度を高めるには、ロングリードポリッシングとショートリードポリッシングを勧める。 Apolloは1回の実行で両方のタイプのリードを使用できるが、非常に時間がかかる。たとえば、ApolloはPac​​Bioリードを使用してE.Coliデータセットを洗練するのに約2時間半かかったが、Raconは約2分しかかからなかった[Firtina et al、2019]。現在、Raconはその速度と比較的優れた精度を考えると、広く使用されているポリッシャーである(また、Raconは全体的に、他のポリッシャーよりも正確な結果を生成することを確認している)。
 ここでは、1回の実行でショートリードとロングリードを利用して、大小のゲノムのロングリードアセンブリをポリッシュするHyPo–a Hybrid Polisherを紹介する。HyPoはユニークなゲノムkmerを利用して、選択的な読み取りセグメントのPOAを開拓して、コンティグのセグメントを選択的にポリッシュする。 Hypoは、Raconに比べて約3分の1の時間で、メモリ要件が約半分で、非常に正確なポリッシュされたアセンブリを生成することを示している。
 

 

インストール

ubuntu18.04LTSでテストした。

依存

Either Mac OS X or Linux are currently supported.

  • Zlib
  • OpenMP
  • GCC (>=8) to support filesystem
  • Following are the commands to update GCC on an Ubuntu machine (from say GCC 5):
sudo apt-get update; sudo apt-get install build-essential software-properties-common -y;
sudo add-apt-repository ppa:ubuntu-toolchain-r/test -y; sudo apt update;
sudo apt install gcc-snapshot -y; sudo apt update
sudo apt install gcc-8 g++-8 -y;
sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-8 60 --slave /usr/bin/g++ g++ /usr/bin/g++-8
sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-5 60 --slave /usr/bin/g++ g++ /usr/bin/g++-5
  •  KMC3
#KMC3
conda install -c bioconda kmc
  • HTSLIB (version >=1.10)

本体 Github

git clone --recursive https://github.com/kensung-lab/hypo hypo 
cd hypo
#Htslibも導入
chmod +x install_deps.sh
./install_deps.sh
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
make -j 8

> ./hypo

# ./hypo 

[Hypo::] Error: Invalid command: Too few arguments!

 

 Usage: hypo <args>

 

 ** Mandatory args:

-r, --reads-short <str>

Input file name containing reads (in fasta/fastq format; can be compressed). A list of files containing file names in each line can be passed with @ prefix.

 

-d, --draft <str>

Input file name containing the draft contigs (in fasta/fastq format; can be compressed). 

 

-b, --bam-sr <str>

Input file name containing the alignments of short reads against the draft (in bam/sam format; must have CIGAR information). 

 

-c, --coverage-short <int>

Approximate mean coverage of the short reads. 

 

-s, --size-ref <str>

Approximate size of the genome (a number; could be followed by units k/m/g; e.g. 10m, 2.3g). 

 

 

 ** Optional args:

-B, --bam-lr <str>

Input file name containing the alignments of long reads against the draft (in bam/sam format; must have CIGAR information). 

[Only Short reads polishing will be performed if this argument is not given]

 

-o, --output <str>

Output file name. 

[Default] hypo_<draft_file_name>.fasta in the working directory.

 

  -t, --threads <int>

Number of threads. 

[Default] 1.

 

  -p, --processing-size <int>

Number of contigs to be processed in one batch. Lower value means less memory usage but slower speed. 

[Default] All the contigs in the draft.

 

  -m, --match-sr <int>

Score for matching bases for short reads. 

[Default] 5.

 

  -x, --mismatch-sr <int>

Score for mismatching bases for short reads. 

[Default] -4.

 

  -g, --gap-sr <int>

Gap penalty for short reads (must be negative). 

[Default] -8.

 

  -M, --match-lr <int>

Score for matching bases for long reads. 

[Default] 3.

 

  -X, --mismatch-lr <int>

Score for mismatching bases for long reads. 

[Default] -5.

 

  -G, --gap-lr <int>

Gap penalty for long reads (must be negative). 

[Default] -4.

 

  -n, --ned-th <int>

Threshold for Normalised Edit Distance of long arms allowed in a window (in %). Higher number means more arms allowed which may slow down the execution.

[Default] 20.

 

  -q, --qual-map-th <int>

Threshold for mapping quality of reads. The reads with mapping quality below this threshold will not be taken into consideration. 

[Default] 2.

 

  -i, --intermed

Store or use (if already exist) the intermediate files. 

[Currently, only Solid kmers are stored as an intermediate file.].

 

  -h, --help

Print the usage. 

 

 

 

実行方法

1、mapping

#short read
minimap2 -t 16 -ax sr draft_genome.fasta pair_*.fq.gz |samtools sort -@ 12 -O BAM - > short.bam

#long(ONT)
minimap2 -t 16 -ax map-ont draft_genome.fasta nanopore.fq.gz |samtools sort -@ 12 -O BAM - > long.bam

 

2、fastq名のテキストファイル作成

echo -e "pair_1.fq.gz\npair_2.fq.gz" > names.txt

 

3、polishing

short

hypo -d draft_genome.fasta -r names.txt -s 3g -c 55 -b short.bam -p 96 -t 8 -o output.fa

short & long

hypo -d draft_genome.fasta -r names.txt -s 3g -c 55 \
-b short.bam -B mapped-lg.sorted.bam -p 96 -t 8 \
-o output.fa
  • -t     Number of threads [Default] 1.
  • -p    Number of contigs to be processed in one batch. Lower value means less memory usage but slower speed.  [Default] All the contigs in the draft.
  • -m   Score for matching bases for short reads [Default] 5.
  • -d    Input file name containing the draft contigs (in fasta/fastq format; can be compressed).
  • -r    Input file name containing reads (in fasta/fastq format; can be compressed). A list of files containing file names in each line can be passed with @ prefix.
  • -b    Input file name containing the alignments of short reads against the draft (in bam/sam format; must have CIGAR information).
  • -B   Input file name containing the alignments of long reads against the draft (in bam/sam format; must have CIGAR information).  [Only Short reads polishing will be performed if this argument is not given.

 

 

引用

HYPO: SUPER FAST & ACCURATE POLISHER FOR LONG READ GENOME ASSEMBLIES

Ritu Kundu, Joshua Casey, Wing-Kin Sung

bioRxiv preprint first posted online Dec. 20, 2019

 


 

 

 

Whisper 2

 

第3世代のシーケンシングの開発にもかかわらず、高いスループットと低いエラーレートのショートリードプラットフォームは多くの生物学的分析に不可欠なままである。 これらは、とりわけ、スモール(Kim et al、2018)および構造(Cameron et al。、2019)変異コールだが、RNA配列(Stark et al。、2019)またはゲノムアセンブリ(Bertrand et al。、2019)にも当てはまる。 バリアントコーラーの大部分はリファレンスゲノムへのリードのマッピングを必要とするため、マッピングの信頼性は、バリアントコールの精度にとって重要である。ここではショートリードマッピングアルゴリズムであるWhisper 2を紹介する。 Whisper 2は Indel処理のための新しい手順を備えており、優れた精度を他と競合力のある実行時間で提供するバリアントコールパイプラインである。

 

 

インストール

ubuntu18.04LTSでテストした。

本体 Github

リリースからバージョンv2.0をダウンロードする。

./whisper

# ./whisper

Whisper v. 2.0 (2019-12-15)

Usage:

   whisper [options] <index_name> @<files> 

   whisper [options] <index_name> file_se 

   whisper [options] <index_name> file_pe1 file_pe2

Parameters:

  index_name   - name of the index (as created by asm_pp)

  files        - name of the file containing list of FASTQ files with seq. reads

  file_se      - FASTQ file (single-end)

  file_pe[1|2] - FASTQ files (paired-end)

Options:

  -b <value> - no. of temporary files (minimum: 100, default: 384)

  -clipping-distance <value> - no. of sigmas for max. additional distance in clipping (default: 14)

  -d[fr/ff/rf] - mapping orientation (default: -dfr (forward - reverse)

  -dist_paired <value> - max. distance for paired read (default: 1000)

  -e <value> - max. no of errors (default: auto)

  -e-paired <value> - max. fraction of errors in paired read (default: 0.09)

  -enable-boundary-clipping <value> - enable clipping at boundaries when a lot of mismatches appears (default: 1)

  -enable-mapping_indels <value> - enable looking for long indels during mapping stages (default: 1)

  -enable-short-indel-refinement <value> - enable short indel refinement after mapping (default: 1)

  -enable-short-reads <value> - enable reads shorter than 90% of the longest reads (default: 0)

  -filter <value> - store only mappings for given chromosome (default: )

  -gap-del-open <value> - score for gap (del) open (default: -5)

  -gap-del-extend <value> - score for gap (del) extend (default: -0.4)

  -gap-ins-open <value> - score for gap (ins) open (default: -5)

  -gap-ins-extend <value> - score for gap (ins) extend (default: -0.4)

  -gzipped-SAM-level <value> - gzip compression level of SAM/BAM, 0 - no compression (default: 0)

  -high-confidence-sigmas <value> - (default: 4)

  -hit-merging-threshold <value> - minimal distance between different mappings (default: 12)

  -hit-merging-wrt-first <value> - calculate distance in marged group w.r.t. first (default: 1)

  -m[f/s/a] - mode: first stratum/second stratum/all strata (default: first stratum)

  -mask-lqb <value> - mask bases of quality lower than value (default: 0)

  -max-indel-len <value> - max. indel length (default: 50)

  -min-clipped-factor <value> - mask bases of quality lower than value (default: 1)

  -out <name> - name of the output file (default: whisper)

  -penalty-saturation <value> - no. of sigmas for max. penalty in matching pairs (default: 7)

  -rg <read_group> - complete read group header line, ? ? character will be converted to a TAB in the output SAM while the read group ID will be attached to every read (example line: ?@RG id:foo SM:bar?)

  -r[s|p] - single or paired-end reads (default: single)

  -score-discretization-threshold (default: 0.5)

  -score-clipping <value> score for clipping (default: -6)

  -score-match <value> - score for matching symbol (default: 1)

  -score-mismatch <value> - score for mismatching symbol (default: -5)

  -sens <value> - turn on/off sensitive mode (default: 1)

  -sens-factor <value> - sensitivity factor (default: 2.5)

  -stdout - use stdout to store the output (default: 0)

  -store-BAM - turn on saving in BAM (default: 0)

  -t <value> - no. of threads (0-adjust to hardware) (default: 0)

  -temp <name> - prefix for temporary files (default: ./whisper_temp_)

  -x - load complete suffix arrays in main memory (default: 0)

Examples:

  whisper human @files

  whisper -temp temp/ human reads1.fq reads2.fq

  whisper -out result.sam -temp temp/ -t 12 human reads1.fq reads2.fq

./whisper-index

# ./whisper-index 

Whisper index construction v. 2.0 (2019-12-15)

Usage: 

   whisper-index <index_name> <ref_seq_file_name> <dest_dir> <temp_dir>

   whisper-index <index_name> <@ref_seq_files_name> <dest_dir> <temp_dir>

Hints:

   * vcf_name can be . if not used

 

 

実行方法

1、indexing

whisper-index hg38-chr20 chr20.fa index-dir temp-dir

  

2、mapping

whisper hg38-chr20 pair_1.fq pair_2.fq

  

引用

Whisper 2: indel-sensitive short read mapping

Sebastian Deorowicz, Adam Gudys

bioRxiv preprint first posted online Dec. 19, 2019

 

Whisper: read sorting allows robust mapping of DNA sequencing data
Sebastian Deorowicz, Agnieszka Debudaj-Grabysz, Adam Gudyś, Szymon Grabowski
Bioinformatics, Volume 35, Issue 12, June 2019, Pages 2043–2050

 

体細胞コピー数変化イベントを調べるFACETSをワンライナーで実行するcnv_facets

 2019 12/27 誤字修正

 

 Cancer Genome Atlas(TCGA)およびInternational Cancer Genome Consortium(ICGC)プロジェクトを含む大規模なシーケンス研究により、腫瘍と正常なサンプルペアの何万もの全ゲノム(WGS)および全エキソーム(WES)が生成された。対立遺伝子特異的コピー数解析は、変異の同定を超えてシーケンシングデータの有用性を大幅に拡張できる。ここでは次世代シーケンス(NGS)データ用のllele-specific copy number analysis (ASCN) パイプラインおよびオープンソースソフトウェアであるFACETS(Fraction and Allele-Specific Copy Number Estimates from Tumor Sequencing)を紹介する。

 ASCN分析には、従来の総コピー数分析と比較していくつかの大きな利点がある。第一に、総コピー数のみを分析することでは検出できないコピー中立的なヘテロ接合性喪失(LOH)イベントを含む、コピー数異常のより包括的な識別を提供する。したがって、ゲノム全体のLOHパターンを体系的に評価できる。さらに、従来の分析では通常、総コピー数比を定性的なコピー数状態(high versus low levelのゲイン、shallow versus deep levelのロス、normal)に変換するが、ASCN分析を使用して、ホモ接合型欠失、ヘテロ接合型欠失、コピー中立LOH、対立遺伝子特異的ゲインおよび対応する整数コピー数増幅の検出のためにゲノムに正確にアノテーションを付けることができる。さらに、ASCN分析は、腫瘍の純度と倍数性のより正確な推定値を提供する。出力は、体細胞点突然変異の強化されたクローン不均一性分析に使用できる。

 初期のASCNメソッドは、主にコピー数アレイプラットフォーム用に開発された(ref.1〜4)。最近では、さまざまな分析戦略に基づいて、次世代のシーケンスデータ用に多くのASCNメソッドが開発された。Patchwork(ref.5)は、総リードカウントに基づいてゲノムをセグメント化し、各セグメント内の対立遺伝子固有のコピー数を推定する。制限は、合計リードカウントをセグメント化するだけでは完全なイメージを提供せず、コピーニュートラルLOHなどの特定のイベントを見逃すことが避けられないことである。 Falcon(ref.6)は、ヘテロ接合SNP遺伝子座からの対立遺伝子リードカウントのために、二項プロセスを使用した共同セグメンテーション手順を提供する。 TITEN(ref.7)を含むいくつかの他の方法は、ベイズ混合モデル(ref.8)、隠れマルコフモデル(ref.7)または他の最尤法(ref.9)を含むさまざまな確率論的モデリングアプローチを使用して、コピー数分析の精度を高めるために腫瘍の純度とクローンの不均一性をさらに考慮した(ref.10)。

 FACETSは、既存のメソッドに対していくつかのユニークな貢献を提供する。 1つは、モデルの仮定に依存せず、ゲノムの変化点を検索するための高速実装を提供する、合計およびアレル固有のリードカウントを直接組み合わせることで、Hotelling T2統計に基づくノンパラメトリックな共同セグメンテーションアプローチを採用している。

 対立遺伝子の不均衡はヘテロ接合部位でしか測定できないため、ASCN分析では通常、SNPベースのアプローチを使用する。データをシーケンスするためのほぼすべてのASCNメソッドは、ヘテロ接合サイトからのリードカウント情報のみを使用する。ただし、ヘテロ接合サイトはサブジェクト固有であり、まばらであるため、コピー総数の情報が失われる。したがって、ヘテロ接合であろうとホモ接合であろうと、すべてのSNPからの対立遺伝子特異的リードカウントの体系的な列挙は、合計および対立遺伝子特異的コピー数の両方に関する完全な情報を提供する。さらに、ターゲットインターバルに沿った一連の疑似SNP(非多型遺伝子座)からのリードカウントも使用して、連続SNP間に大きなギャップがある領域が総コピー数分析で表されるようにする。

(以下省略)

 

 

Githubより

The advantage of cnv_facets over the original facets package is the convenience of executing all the necessary steps, from BAM input to VCF output, in a single command line call.

 

インストール

ubuntu18.04LTSのpython3.7環境でテストした。

 

本体 Github

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

cnv_facets.R -h

$ cnv_facets.R -h

usage: /home/kazu/anaconda3/envs/cnv_facets/bin/cnv_facets.R

       [-h] --out OUT [--snp-tumour SNP_TUMOUR] [--snp-normal SNP_NORMAL]

       [--snp-vcf SNP_VCF] [--snp-mapq SNP_MAPQ] [--snp-baq SNP_BAQ]

       [--snp-count-orphans] [--snp-nprocs SNP_NPROCS] [--pileup PILEUP]

       [--depth DEPTH DEPTH] [--targets TARGETS] [--cval CVAL CVAL]

       [--nbhd-snp NBHD_SNP] [--annotation ANNOTATION]

       [--gbuild {hg19,hg38,mm9,mm10}] [--unmatched] [--rnd-seed RND_SEED]

       [-v]

 

DESCRIPTION 

Detect somatic copy number variants (CNVs) and estimate purity and ploidy in a

tumour sample given a matched normal sample. The core of the analysis is based

on the package FACETS (https://github.com/mskcc/facets).

 

INPUT

 

* A bam file for the tumour sample, a bam file for the matched normal sample, and a vcf

  file of SNP positions. Bam and vcf files must be sorted and indexed

 

Or:

 

* A csv file of counts (pileup) generated by a previous run of cnv_facets.R

  or by the accompanying program snp-pileup

 

EXAMPLE

 

    cnv_facets.R -n normal.bam -t tumour.bam -vcf snps.vcf.gz -o my_cnv

 

See the online documentation for more details and usage.

Version 0.15.0; facets=0.5.14

 

optional arguments:

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

  --out OUT, -o OUT     Output prefix for the output files

  --snp-tumour SNP_TUMOUR, -t SNP_TUMOUR

                        BAM file for tumour sample

  --snp-normal SNP_NORMAL, -n SNP_NORMAL

                        BAM file for normal sample

  --snp-vcf SNP_VCF, -vcf SNP_VCF

                        VCF file of SNPs where pileup is to be computed

  --snp-mapq SNP_MAPQ, -mq SNP_MAPQ

                        Sets the minimum threshold for mapping quality. Default 5

  --snp-baq SNP_BAQ, -bq SNP_BAQ

                        Sets the minimum threshold for base quality. Default 10

  --snp-count-orphans, -A

                        Do not discard anomalous read pairs

  --snp-nprocs SNP_NPROCS, -N SNP_NPROCS

                        Number of parallel processes to run to prepare the read pileup file.

                        Each chromsome is assigned to a process. Default 1

  --pileup PILEUP, -p PILEUP

                        Pileup for matched normal (first sample) and tumour

                        (second sample). Not needed if using BAM input. This file

                        is the <prefix>.cvs.gz output of of a previous run of cnv_facet.R

  --depth DEPTH DEPTH, -d DEPTH DEPTH

                        Minimum and maximum depth in normal sample for a position

                        to be considered. Default 25 4000

  --targets TARGETS, -T TARGETS

                        BED file of target regions to scan. It may be the target regions

                        from WEX or panel sequencing protocols. It is not required, even

                        for targeted sequencing, but it may improve the results

  --cval CVAL CVAL, -cv CVAL CVAL

                        Critical values for segmentation in pre-processing and

                        processing. Larger values reduce segmentation. [25 150] is

                        facets default based on exome data. For whole genome consider

                        increasing to [25 400] and for targeted sequencing consider

                        reducing them. Default 25 150

  --nbhd-snp NBHD_SNP, -snp NBHD_SNP

                        If an interval of size nbhd-snp contains more than one SNP,

                        sample a random one.  This sampling reduces the SNP serial

                        correlation. This value should be similar to the median insert

                        size of the libraries. If "auto" and if using paired end BAM

                        input, use the estimated insert size from the normal bam file.

                        Otherwise use 250. Default auto

  --annotation ANNOTATION, -a ANNOTATION

                        Optional annotation file in BED format where the 4th column

                        contains the feature name (e.g. gene name). CNVs will be

                        annotated with an additional INFO/TAG reporting all the

                        overalapping features

  --gbuild {hg19,hg38,mm9,mm10}, -g {hg19,hg38,mm9,mm10}

                        String indicating the reference genome build. Default hg38.

  --unmatched, -u       Normal sample is unmatched. If set, heterozygote SNPs are

                        called using tumor reads only and logOR calculations are different

  --rnd-seed RND_SEED, -s RND_SEED

                        Seed for random number generator. Default: The name of the input file

  -v, --version         show program's version number and exit

 

Options --snp-* are used only to generate the pileup file. They are ignored with option --pileup

 

 

実行方法

tumorのbam、コントロールのnormal.bam (typically, a blood sample from the same patient: Githubより) を指定する。またcommon SNPsのVCFを指定する。ここではdbSNP150のHuman variations without clinical assertionsを指定した。

cnv_facets.R -t tumour.bam -n normal.bam -vcf common_all.vcf.gz -o output

 

出力

f:id:kazumaxneo:20191226180311p:plain

output.cnv.png

f:id:kazumaxneo:20191226180038p:plain

 

output.cov.pdf

f:id:kazumaxneo:20191226180132p:plain

引用

FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing
Ronglai Shen* and Venkatraman E. Seshan

Nucleic Acids Res. 2016 Sep 19; 44(16): e131