macでインフォマティクス

macでインフォマティクス

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

メタゲノムデータをサブサンプリングして繰り返し アセンブリする Spherical

 

 過去10年間、研究者らは、ハイスループットシーケンシングを利用して、世界中の多様な環境からの微生物群集の構造と機能を調べてきた[論文より ref.1、2、3]。これらの研究は、微生物の働きについてユニークで斬新な洞察を提供してきたが、入手可能なツールが、データが表す完全な機能的または分類学的多様性を記述していないというコンセンサスが増えている[ref.4]。 「失われた」情報をキャプチャし、サンプリングされたメタゲノム全体を表すアセンブリの生成を可能にするために、アセンブラとワークフローを開発する必要がある。

 数学的には、ゲノムの新しいアセンブリは、効率的なアルゴリズムが知られていない問題のクラス(NP-hard)[ref.5]にあり、さまざまなヒューリスティック解の提案につながっている[ref.6,7,8,9]。これらは、単純なオーバーラップレイアウトコンセンサスアプローチ[ref.5]、de Bruij ngraph[ref.10]などから構成されている。しかし、微生物群全体のゲノム情報(メタゲノミクス)をサンプリングするためのシーケンシング手法が出現し始めているため、新しいアプローチが必要であることは明らかである[re.12 pubmed]。

 メタゲノムシーケンシングデータセットの最大の問題は、天然の微生物群集の種の複雑さと不均一な分布に起因する、微生物群の不均一なカバレッジである[ref.4]。これは、コミュニティで支配的な種の過剰なシーケンス解析と、少数の種のゲノムの断片化したアセンブリに繋がる[13]。

 この問題に対処するための有望なソリューションは、データをより簡単に管理可能な部分に分割するための「divide and conquer」アプローチである[ref.13]。例えば、環境サンプルからのデータは、コミュニティからの異なる分類群を表す「ビン」に分割することができる[ref.14]。シーケンスリードは、kmer頻度やそれらに含まれるGuanineとCytosines(GC)のパーセンテージなどのプロパティに基づいてビンにソートできる[ref.15]。これは、豊富な種の集積率を高める可能性を秘めているが、データの正確な分割に大きく依存する[ref.16]。実際、このようにして作られたメタゲノムデータのビンは、コミュニティの複雑さに応じて単一の種または門全体を表すかもしれない[ref.16]。 SLICEMBLER [ref.8]はカバレッジに基づいてメタゲノムの所定の「スライス」を生成することによって「divide and conquer」アプローチを実装する。各スライスは別々にアセンブリされ、頻繁に発生する文字列が識別され、コンティグをマージするために使用される。これは、カバレッジが分かっている深くシーケンス解析されたゲノムデータセットに対してはうまくいくが、未知の微生物群集からのメタゲノミクスデータセットでは、カバレッジは一般に未知の変数である[ref.17]。アセンブリの反復ビニングを利用する試みも開発されているが、公に利用可能ではない[ref.18]。

 デジタル正規化は、複雑なデータセットに取り組むために一般的に使用されるもう1つのアプローチである。 Khmerパッケージ[ref.19]に実装されているこのメソッドは、データの重複部分の発生を識別して減らすためにkmersを使用する。これは、過度に表現された分類群のカバレッジを減らし、より均一にするためにカバレッジを正規化する効果を有する[19]。ただし、この前処理ステップではデータセットのサイズを縮小できるが、データの複雑さを軽減することはできない[rerf.19]。後続のアセンブリでは、コミュニティの過小評価の問題が発生しやすくなる。

 著者らは、この問題に対処するため、上記のアセンブリ手法のいずれかと組み合わせて使用​​できる簡単で効果的なヒューリスティックを提案する。反復的なアプローチに基づいて、 著者らのワークフロー(Sphericalと呼ばれる)は、アセンブリ内で考慮されていないリードを識別し、それをアセンブリの連続ラウンドの入力として使用する。 Sphericalは、「divide and conquer」方法で全体のランダムなサブセットをアセンブリすることによって必要とされるリソースを最小限に抑えるオプションも提供する。このアプローチの本質的な欠点は、アセンブリを複数回に渡り反復実行することによって、時間要件も増加することである。

 

f:id:kazumaxneo:20180603174415j:plain

図1。フローチャート。論文より転載。

 

論文中では、メモリの削減量とその影響、繰り返しアセンブリの効果等が議論されています。

 

インストール

cent os6、python2.7.12でテストした。

依存

  • Python 2.7
  • Velvet (tested using version 1.2.10)
  • Bowtie2 (tested using version 2.2.3)

Python modules;

  • Numpy
  • HTSeq

 

本体 Github

https://github.com/thh32/Spherical

git clone https://github.com/thh32/Spherical.git 
cd Spherical/

ヘルプはない。

Githubより転載。

f:id:kazumaxneo:20180603172747j:plain

 

 

 

ラン

 

テストデータが用意されている。ランするには、シーケンスデータを解凍して以下のように打つ。

gzip -dv test_data.fa.gz

python Spherical.py -fasta -m -k 21 -R 1 -align 99 -iter 5 -i test_data.fa -velvet -bowtie2 -o test_assembly

アセンブリにはvelvetを指定し(-velvet)、リードのマッピングにはbowtie2を指定している(-bowtie2)。アセンブリのk-merは21にして(-k 21)、リードはサブサンプリングせず100%使用している(-R 1)(メモリが足りなければ減らす)。アセンブリしたcontigにリードの99%がマッピングされない限り(-align 99)、最大5回サイクルは繰り返される(-iter 5)。

 

引用

Spherical: an iterative workflow for assembling metagenomic datasets

Thomas C. A. Hitch and Christopher J. Creevey.

BMC Bioinformatics. 2018; 19: 20.

 

臨床向けアンプリコンシーケンス自動解析パイプライン Canary

 

 臨床診断は、ヌクレオチドレベルで患者DNAを分析することができる技術によって変化している。臨床シーケンシングの精度、処理時間および再現性は、rawシーケンシングデータを有意義なバリアントに変換するバイオインフォマティクスパイプラインに大きく依存している。これらのパイプラインは、複数のソフトウェア依存、ポータビリティの欠如、複雑なパラメータ調整、および並列実行のためのクラスタコンピューティング環境を必要とすることが多い[論文より ref.1]。これらの属性により、現場での展開が困難なパイプラインが発生する。

 ここでは、マルチツール・パイプラインの機能を実行するスタンドアロンJavaユーティリティであるCanaryを紹介する。Canaryはアンプリコン・アッセイから生成された圧縮FASTQファイルから変異遺伝子のアノテーション付きVCFファイルを直接生成することができる。 CanaryはJavaランタイムのみを必要とするため、現在のパイプラインの無数の依存関係とは対照的に、Javaをインストールした任意のコンピュータに展開できる。さらに、パブリックリポジトリからDocker [ref.2]イメージとして入手でき、単一のコマンドでDockerをサポートするどのプラットフォームでもインストールおよび実行できる(%docker run -v /tmp/data:/canary.data dockercanary / canary)。

 ショットガンシーケンス解析パイプラインを用いたアンプリコンシーケンスデータの処理は、スピードと品質の点で最適以下の結果を導く[ref.3,4]。 Illumina BaseSpace [ref.5]のような市販のプラットフォームの外でアンプリコンデータを処理するオプションは比較的少数である。このプラットフォームは、Amplicon TruSeqやTruSightなどの独自のIlluminaアッセイにのみ対応し、他の遺伝子を標的とするカスタムパネルや、アンプリコン解析を社内のパイプラインに組み込むことはできない。非営利のアンプリコンソフトウェアには、アラインメントを行わないMutascope [ref.4]、バリアントコールを行わないAmpliVar [ref.6]、ノーマライズしないUNDR ROVER [ref.3]がある。 Canaryは、圧縮されたFASTQファイルから、臨床的なキュレーションに適したアノテーション付きVCFファイルを、単一のコマンドで実行し必要なパイプライン手順を簡潔にする。 FASTQファイルは、FASTQC [ref.7]などのプログラムで事前にクオリティ管理されているものとする。

 著者らが知る限り、Canaryは、単一の実行可能プログラムで、アライメント、バリアントコール、正規化、トランスクリプトの選抜、およびアノテーションの必要なパイプラインステップすべてを実行できる唯一のツールである。

 

実行可能な処理の詳細(論文リンク)。

 

インストール

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

依存

  • Java JDK 1.7 from Oracle
  • Groovy 2.1.9 from here
  • Gradle for building Canary (we use 1.10 at time of writing) Gradle 1.10
  • Genome Analysis Toolkit (Currently GATK 3.3) and the Sting utility JAR (currently 2.1.8) available from here
  • The PathOS Core library available from PathosCore-all-1.3.jar and maintained here.
  • JNI wrapper to the striped Smith-Waterman alignment library SSW see here

本体 Github

https://github.com/PapenfussLab/Canary

依存が多いので、ここではdockerコンテナでテストする。

docker run -v /tmp/data:/canary.data dockercanary/canary

#ここではコンテナに入って作業を行う。公式の説明のようにしてもOK。
docker run -it --entrypoint 'bash' dockercanary/canary

Canary

# Canary 

error: Missing required options: a, p

usage: Canary [options] read1.fastq.gz read2.fastq.gz

 

Available options (use -h for help):

 -a,--amplicon <arg>        Amplicon FASTA file [required]

 -ano,--annotation <arg>    File of MyVariant annotation fields

 -b,--bam <arg>             Optional BAM file of alignment

 -c,--complex               Coalesce complex events aka MNPs

 -cols,--columns <arg>      File of VCF field names to output to TSV (one

                            per line with optional alias after comma)

 -d,--debug                 Turn on debugging (Note: will generate large

                            file of alignments [debug.out])

 -f,--flank <arg>           Size of flanking region (bp) [5]

 -fastq,--fastq <arg>       Optional FASTQ output files prefix

 -filt,--filter <arg>       List of comma separated amplicon names to use

 -h,--help                  This help message

 -maxmut,--maxmut <arg>     Maximum number of mutations allowed per read

                            pair [10]

 -minpair,--minpair <arg>   Min read pairs for variants [10]

 -mnpgap,--mnpgap <arg>     Maximum size of inter mutation gap for complex

                            mutations [15]

 -mnpmax,--mnpmax <arg>     Maximum size of complex mutations [30]

 -mut,--mutalyzer <arg>     Mutalyzer annotation server host

                            https://mutalyzer.nl

 -n,--nocache               Dont use read cache

 -norm,--normalise <arg>    Generates annotated VCF file from VCF output

 -o,--output <arg>          Output report file

 -p,--primers <arg>         Amplicon Primers file [required]

 -r,--reads <arg>           Percent of reads to process [100]

 -t,--tsv <arg>             TSV (Tab separated variable) file of VCF

                            output

 -ts,--transcript <arg>     File of transcripts mapping genes -> refseq

                            (without version)

 -v,--vcf <arg>             Found variants VCF file [canary.vcf]

 -vaf,--vaf <arg>           Minimum VAF for variants [3.0%]

 -ver,--version             Display Canary version and exit

 

 

 

ラン

テストラン用のシェルスクリプトを使って動作を確認する。

mkdir ~/test_dir
cd ~/test_dir

#テストラン実行
/opt/Canary/bin/runCanary.sh

 このスクリプトを実行すると、Miseq のTruSeqライブラリで実行されたがん原遺伝子48個のアンプリコンシーケンスデータ(/opt/Canary/Fastq)の分析が自動で行われる。

シェルスクリプトの中身は以下のようになっている。

> cat /opt/Canary/bin/runCanary.sh

#!/bin/bash

#

# Run Canary for testing

#

 

Canary --mutalyzer 'https://mutalyzer.nl' \

--amplicon   ${CANARY_HOME}/Amplicon/amplicon.fa \

--primers    ${CANARY_HOME}/Amplicon/amplicon.primers.tsv \

--transcript ${CANARY_HOME}/etc/transcript.tsv \

--columns    ${CANARY_HOME}/etc/cols.txt \

--annotation ${CANARY_HOME}/etc/myvariant.txt \

--reads      100 \

--vaf        3.0 \

--normalise  canary.norm.vcf \

--tsv        outvcf.tsv \

--output     canary.tsv \

--vcf        canary.vcf \

--bam        canary.bam \

$* \

${CANARY_HOME}/Fastq/*R1_001.fastq.gz \

${CANARY_HOME}/Fastq/*R2_001.fastq.gz

 

 

head ${CANARY_HOME}/Amplicon/amplicon.fa

# head ${CANARY_HOME}/Amplicon/amplicon.fa

>1:43814982-43815163

CCGTCCTGGGCCTGCTGCTGCTGAGGTGGCAGTTTCCTGCACACTACAGGTACCGCCCCC

GCCAGGCAGGAGACTGGCGGTGGACCAGGTGGAGCCGAAGGCCTGTAAACAGGCATTCTT

GGTTCGCTCTGTGACCCCAGATCTCCGTCCACCGCCCGTGCGCACCTACGGCTTCGCCAC

TT

>1:115256500-115256680

ATTGGTCTCTCATGGCACTGTACTCTTCTTGTCCAGCTGTATCCAGTATGTCCAACAAAC

AGGTTTCACCATCTATAACCACTTGTTTTCTGTAAGAATCCTGGGGGTGTGGAGGGTAAG

GGGGCAGGGAGGGAGGGAAGTTCAATTTTTATTAAAAACCACAGGGAATGCAATGCTATT

G

head ${CANARY_HOME}/Amplicon/amplicon.primers.tsv

# head ${CANARY_HOME}/Amplicon/amplicon.primers.tsv

1:43814982-43815163 24 26 MPL1_2.chr1.43815008.43815009_tile_1

1:115256500-115256680 26 27 NRAS1_7.chr1.115256528.115256531_tile_1

1:115258702-115258884 26 29 NRAS8_13.chr1.115258730.115258748_tile_1

2:29432636-29432822 27 27 ALK1.chr2.29432664.29432664_tile_1

2:29443667-29443836 25 26 ALK2.chr2.29443695.29443695_tile_1

2:132181238-132181750 30 25 Off_target_7_GNAQ_5.chr9.80409379.80409508_tile_1-GNAQ_7.chr9.80336240.80336429_tile_1

2:132181332-132181839 25 24 Off_target_9_GNAQ_7.chr9.80336240.80336429_tile_2-GNAQ_5.chr9.80409379.80409508_tile_2

2:132181448-132181600 25 28 Off_target_8_GNAQ_6.chr9.80343430.80343583_tile_1-GNAQ_7.chr9.80336240.80336429_tile_3

2:209113084-209113264 27 27 IDH1_1_2.chr2.209113112.209113113_tile_1

2:212288912-212289100 27 27 ERBB4_1_2.chr2.212288942.212288955_tile_1

独自に解析するには、これらのファイルを準備する。primerのフォーマット詳細はGithub参照(リンク)。

 

 

ランが終わると、bamに加え、バリアントをアノテーション付きでまとめたファイル(canary.vcf、outvcf.tsv)が出力される。canary.vcfは変異遺伝子のアノテーション付きVCF。outvcf.tsvはタブ区切りのファイルで、スプレッドシートやデータベースに読み込みやすいように工夫されている。

# head -n 50 outvcf.tsv |tail -n 10

##FORMAT=<ID=RDF,Number=1,Type=Integer,Description="Depth of reference-supporting bases on forward strand (reads1plus)">

##FORMAT=<ID=RDR,Number=1,Type=Integer,Description="Depth of reference-supporting bases on reverse strand (reads1minus)">

##FORMAT=<ID=ADF,Number=1,Type=Integer,Description="Depth of variant-supporting bases on forward strand (reads2plus)">

##FORMAT=<ID=ADR,Number=1,Type=Integer,Description="Depth of variant-supporting bases on reverse strand (reads2minus)">

##unpack="expanded by org.petermac.util.Vcf.unpack()"

#chr position ID REF altbases QUAL FILTER GT GQ HGVSg HGVSc HGVSp gene genename consequence cadd_raw mutdb.mutpred_score mutdb.cosmic_id mutdb.uniprot_id mutdb.strand mutdb.rsid clinvar.allele_id clinvar.gene.id clinvar.gene.symbol clinvar.cytogenic clinvar.variant_id clinvar.rcv

2 212578379 . TA T . PASS 0/1 chr2:g.212578380del

2 212578379 . TAA T . PASS 0/1 chr2:g.212578380_212578381del

2 212578379 . TAAA T . PASS 0/1 chr2:g.212578380_212578382del

2 212578379 . T TA . PASS 0/1 chr2:g.212578379_212578380insA

 

引用

Canary: an atomic pipeline for clinical amplicon assays

Doig KD, Ellul J, Fellowes A, Thompson ER, Ryland G, Blombery P, Papenfuss AT, Fox SB

BMC Bioinformatics. 2017 Dec 15;18(1):555.

 

CGDV

 

 次世代シークエンシング(NGS)技術の進歩により、前例のない量の異なる形式のデータが生成されている。大規模なNGSデータの解釈は複雑で困難である。可視化はNGSデータを解釈する手段の1つであり、データ分析において重要な役割を果たしている。円グラフは、大規模なデータとそれらの相互関係を1つのフレームで表示するのに非常に便利である。円形ビューでデータを視覚化するためのさまざまなWebベースのツールがある(論文 表1)。 Circos [論文より ref.1]に基づくOnline Circos(http://mkweb.bcgsc.ca/tableviewer/)は、円形ビューでデータを視覚化するためのWebツールだが、Circosの使用方法に関する詳細な知識が必要となる。 CiVi [ref.2](紹介)のようなツールは特定のゲノミクスデータのみを扱うことができ、微生物ゲノムからのデータをプロットすることに限定されている。もう1つのWebtoolであるCliCo FS [ref.3]はgene bankファイルのみをサポートしている。他の種類のファイルの場合、自動化されていないため、アップロード前にファイルをフォーマットする必要がある。さらに、ClicO FSはデータ駆動型ではなく視覚化駆動型である。さらに、プロットを生成する前に複数のクリックが必要である(紹介)。 J-Circos [ref.4]などのデスクトップベースのアプリケーションは、実行する前にインストールする必要がある。さらに、J-Circosは、すべてのタイプのゲノミクスおよびトランスクリプトミクスのデータフォーマットをサポートしておらず、モデル生物の限られたセットだけをサポートしている。したがって、これらのツールのいずれもCircroのインストールと使用に関する知識がまったくないか、または最小限の知識を持つ生物学者にとって、循環型の視覚化の形でデータを便利に解釈するための、自動化されたガイド付きの様々なゲノミクスとトランスクリプトミックのraw出力ファイルをサポートしていない。
 Circos [ref.1]を用いて任意のプロットを生成するためには、リファレンスゲノムまたはコンティグの長さ、サイズ、色および各染色体またはコンティグの適切な標識の染色体長さなどの基本情報を定義するカリオタイプ(karyotype)ファイル(karyotype wiki)が必要である。必要な別のファイルは、その内容に基づいてデータを視覚化する方法に関する情報を含むconfiguration(config)ファイルとなる。ユーザーがデータ、コンテンツ、さまざまな視覚化オプションを徹底的に理解する必要があるため、設定ファイルを作成するのは複雑になる。われわれはCircosのラッパーであるCGDVを開発し、大規模なゲノミクスと転写オミクスデータをシームレスに環状可視化する自動化されたガイド付きのツールを提供する。 CGDVは、様々なモデル生物のために予めパッケージ化されたカリオタイプファイルを提供するだけでなく、ユーザーによって提供されるゲノミクスおよびトランスクリプトミックスデータに基づいて設定ファイルを生成する。 CGDVは、SVGおよびPNG形式のデータ特有の環状ビジュアルを生成するために、ゲノミクスおよびトランスクリプトミックスデータの大部分の標準の生の出力ファイルを入力として取り込む(論文 図1)。

 CGDVはApache Webサーバー上で動作する。 CGDVのウェブインタフェースは、ユーザ電子メールID(オプションでゲストユーザとしても実行できる)、モデル生物、円形ダイアグラムを作成したいデータ型などの他のパラメータと共に入力ファイルを必要とする。 入力ファイルから関連情報を抽出し、構成ファイルとデータファイルを作成する。 標準ゲノムのカリオタイプ情報はSQLiteデータベースに保存されている。 モデル生物の選択ごとに、特定のカリオタイプの詳細がデータベースから取り出される。 CGDVは、設定ファイル、データファイル、カリオタイプファイルを使用して、Circos [ref.1]をバックグラウンドで実行し、与えられた入力ファイルから環状の図を作成する。 CGDVは、SVGおよびPNG形式で画像を生成する。 ユーザーが電子メールIDを入力すると、出力は送信日から15日間アーカイブされ、その後に削除される。

 

 

ラン 

https://cgdv-upload.persistent.co.in/cgdv/

f:id:kazumaxneo:20180531215823j:plain

FirefoxChromeに最適化されている。

 

以下のゲノムのkaryotypeに対応している。

f:id:kazumaxneo:20180601153201j:plain

データトラックとして、Exampleの以下のようなデータフォーマットに対応している。

f:id:kazumaxneo:20180603121756j:plain

 

 

ストアカウントでテストデータを何度かアップロードしてみたが、いずれもエラーになる。そもそも動作しているのか不安になってアーカイブを確認してみたが、いくつかのデータはジョブが成功しており、ダウンロードすると可視化できていた。しかしデータトラックの描画までできているデータは見つからなかった(15日しか保存されないので、たまたま初期テストの結果しかなかったのかもしれない)。

改善したら追記します。

 

引用

CGDV: a webtool for circular visualization of genomics and transcriptomics data

Jha V, Singh G, Kumar S, Sonawane A, Jere A, Anamika K

BMC Genomics. 2017 Oct 24;18(1):823.

 

 

トランスクリプトームのblast比較結果を統合し、ベン図を描く VennBLAST

 

 ハイスループットシークエンシングは広範な技術となり、進化的研究を含む様々な研究分野でアクセス可能となっている。ゲノムが利用できない生物の転写産物をシーケンスし、注釈を付ける能力は、分子進化の分野における生物学者、特に非モデル生物を含むルーチンの仕事となっている。 RNA-Seqを介した非モデル生物の集団ゲノム解析の詳細なウォークスルーパイプラインは、De Wit [論文より ref.1 CrossRef]によって記述されている。このプロセスの第1ステップは、多数の配列のde novoアセンブリによる一連の転写セットの構築である。これは、Trinity [ref.2]、Trans-ABySS [ref.3]、Velvet-Oasis [ref.4]、MIRA [ref.5]、Newbler for 454[ref.6]などのさまざまなツールを使用して行うことができる。 Blast2Go [ref.7]、Trinotate [ref.8](紹介)、T-Ace [ref.9 pubmed]、annot8r [ref,10 pubmed]、FastAnnotator [ref.11 pubmed]などのようないくつかのソフトウェアスイートを使用して実行することができる。これらのツールは相同性検索のためのよく参照されてきた方法を使い、、既知配列データベースに対する相同性検索、タンパク質ドメイン同定、GO termおよびpathway解析などを行うが、全てBLASTライクな分析に大きく依存している。

 BLAST(Basic Local Alignment Search Tool)[ref.12]は、未知ゲノム生物のトランスクリプトームシーケンシングからアセンブリされたESTのアノテーションを含む、配列間の局所的な相同性を示す領域を見出すため最も広く使用されるプログラムセットである。BLASTの結果はクエリ配列と個々のマッチ(被験体)の各々との間の類似性のレベルを示しているが、全トランスクリプトーム結果の洞察を達成する能力は乏しい。典型的には、種間の進化的保存の程度を決定し、研究対象の種において利用可能な経路の景観を推定し、あるいは種間発現パターンを比較することができる。

ESTデータベースの2つ以上のセット間のdifferential expressionプロファイルの確立は、アセンブリソフトウエアおよびDESeq [ref.13]およびedgeR [ref.14]のようなデジタル遺伝子発現統計ツールからなるパイプラインを用いて達成することができる。 TrinityとedgeRツールを使用するプロトコルは、Trinityウェブサイト[ref.2 CrossRef]に記載されている。これは、前提としてコンセンサストランスクリプトームを作成する必要がある。 T-Aceはトランスクリプトーム間の統計的な比較を提供する。しかしながら、いくつかの種からのraw発現データを統合してアセンブルされたコンセンサスのトランスクリプトームに依存する。

 発現解析だけでなく、アセンブリソフトウェアやアノテーションソフトウェアの開発にも大きな努力が注がれているが、トランスクリプトーム全体の比較などの下流の調査ははるかに利用できない。ここでは、高速並列化BLASTフィルタリング・ユーティリティと全トランスクリプトのアラインメント比較を組み合わせた、ユーザーフレンドリーな統合ソフトウェア、VennBLASTを紹介する。 VennBLASTは直感的なベン図で結果を示し、BLAST比較での数関係を提供する。 VennBLASTは、全トランスクリプトームの進化的関連性の鳥瞰図を提供するが、遺伝子を有意義なサブグループに解剖し、様々なツールを使用してさらに解析することも可能である。例えばGene set enrichment analysisは、それらの背後にある生物学的意味を理解するために、VennBLAST選択遺伝子リストで実施することができる。これらのタスクを実行することは、バイオインフォマティクスにとっては直接的なことかもしれないが、VennBLASTの対象である非専門家にとっては時間がかかるものである。

 VennBLASTは、対話型ユーザーインターフェイスを備えたデスクトップアプリケーションとして実行するよう実装されている。これはC#プログラミング言語で実装されており、ベン図を描画するために "Venn diagram plotter"コード[ref.15]を使用している。

VennBLASTはhttp://www.ariel.ac.il/research/fbl/softwareで非営利目的で自由に利用できる。

 

ダウンロード

mac os10.13のparallels 12-windows 10のエミュレーション環境でテストした。

公式ページ(プログラムとword形式のマニュアル)

Software

 

ラン

 

下の図のようなワークフローで進める。

f:id:kazumaxneo:20180602173142j:plain

図1。論文より転載。

 

比較したい配列のblast解析を行う。マニュアルの通り、アミノ酸配列データベースに対してblastx解析する。例えば3つのデータベースA、B、 Cに対してコンティグ配列contigg.faをクエリとしてblastx解析を行う。

#Aのindexの作成
makeblastdb -in inputA.aa -dbtype prot
#Aのblastx
blastx -query contig.fa -db inputA.aa -out out.blastxA -outfmt "6 std qcovs" -max_target_seqs 1 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -num_threads 20

#Bのindexの作成
makeblastdb -in inputB.aa -dbtype prot
#Bのblastx
blastx -query contig.fa -db inputB.aa -out out.blastxB -outfmt "6 std qcovs" -max_target_seqs 1 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -num_threads 20

#Cのindexの作成
makeblastdb -in inputA.aa -dbtype prot
#Cのblastx
blastx -query contig.fa -db inputC.aa -out out.blastxC -outfmt "6 std qcovs" -max_target_seqs 1 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -num_threads 20
  • -outfmt "6 std qcovs"     getting the output in a tabular format which includes the default (standard)output parameters as well as the query coverage.
  • -max_target_seqs    control the number of matches recorded in the alignment (1 in this case)
  • best_hit_overhang  and  best_hit_score_edge  are used here in order to search for only the best matches for each query region reporting matches, and to avoid short hits

論文では2つの例で説明している。1つはone to manyの例で、もう1つはmany to oneの例である。one to manyの例では、非モデル生物Stylophora pistillata(ショウガサンゴ)のトランスクプトームデータからアセンブリして得たcontig配列をクエリとして、イソギンチャク目、イシサンゴ目、Hydrozoa目?の3つの目の様々なゲノムに対して、blastx解析をしている。それから、VennBLASTを使って同じ目の結果をマージしてe-valueでフィルタリングし、そのフィルタリング結果を比較してどの目に近いか等を議論している(論文を参照)。

 

 dllとexeファイルが入っている。exeファイルをたたき起動する。

f:id:kazumaxneo:20180602121753j:plain

 

Uploadボタンからblast結果をアップロードし、 必要であれば中央のFilterボタンをクリックしフィルタリングする。例えばE valueが1e-10より高ければ除く。

f:id:kazumaxneo:20180602222741j:plain

save outputでフィルタリング後のblast結果は出力できる。

 

右下のMergeボタンを押すと下の画面が出現する。3つの群のblast結果を入力する。

f:id:kazumaxneo:20180602222936j:plain

 Venn Diagramボタンをクリックするとベン図が出力される(上は動作確認のテスト結果)。

 

引用

VennBLAST—whole transcriptome comparison and visualization tool

Zahavi T, Stelzer G, Strauss L, Salmon AY, Salmon-Divon M.

Genomics. 2015 Mar;105(3):131-6.

 

関心のあるバクテリアゲノムのシグネチャを迅速に検出する Neptune

 

 安価かつ迅速に大量のシーケンスを生成する能力は、生物、特にバクテリアのような比較的小さなゲノムを有する生物全体のゲノムを研究する能力を可能にした。計算生物学者は、歴史的に、少数のバクテリアゲノムを比較し、ヌクレオチド、遺伝子およびゲノムスケールで基本的な特徴付けを行うために、幅広いバイオインフォマティクスソフトウェアツールを使用してきた。しかしながら、バクテリアゲノム全体の効率的な比較分析および特徴付けを行うためのバイオインフォマティクスソフトウェアの必要性が現在存在する。いくつかのツールが登場した。これらのツールのほとんどは、リファレンスマッピングアプローチ(論文より ref.1(CrossRef), 2(pubmed))を用いた1塩基変異(SNV)の同定、または正確なサブストリング(k-mer)(ref.3-5(ref.5 pubmed))に基づく距離推定に焦点を当てていて、スケール変更にも単純な並列化戦略で対応することができる。微生物のゲノム集団を解析してゲノムの特徴を表現型の形質と相関させる微生物genome-wide association studies (GWAS) は、例えば長距離の連鎖不平衡( long range linkage disequilibrium)およびクローン集団構造(ref.6)などに取り組む最近の数学的手法の発達により、表現形質とゲノムの特徴を相関づけるためバクテリアのゲノム集団を分析することを可能にした。 SNVsまたはk-mersと生物学的形質とを関連付けるバクテリアGWAS用のソフトウェアツールがいくつか開発されている(ref.7 pubmed)。しかしながら、バクテリアGWASの場合、特に、新規な生物学的形質を獲得するために水平遺伝子導入に関与する大部分のバクテリアについて、より大きな規模のゲノムの利益および損失を含むバクテリアゲノム変動の全ての様式を同定することが重要である。ある個体群を他の個体群と区別しながらそれらの遺伝子座内の対立遺伝子変異を許容する大規模ゲノム遺伝子座を迅速に抽出することができるスケーラブルなソフトウェアは、バクテリアGWASを達成する上で貴重であり、標的分子診断の開発などの多くの他の用途に有用である。

 この課題に取り組むために、著者らはゲノムシグネチャディスカバリー( genomic signature discovery)の分野を検討した。ここでのシグネチャとは、背景となる配列グループから目的の配列グループを識別することができる配列として定義される。シグネチャーは、ゲノムまたは遺伝子間領域に存在し得、ゲノムアイランド、ファージ領域またはオペロン全体に対応し得る。しかし、シグネチャが機能的に意味のあるコンテンツを含む必要はなく、そのシーケンスによって2つのグループが効果的に区別さえすればよい。効果的なシグネチャ発見アルゴリズムとは、感受性の高さと特異性の高さ両方を持ち、迅速な計算も可能なものである。しかしながら、実際には、これらの3つの属性すべてを有するアルゴリズムを開発することは依然として困難である。病原体検出診断アッセイ(ref.8)を生成する特定の目的で、シグネチャー発見のための初期のアルゴリズムアプローチが開発された。一般に、これらのアプローチは、BLAST(ref.9)のようなアライメントベースの方法を使用して全ての配列を徹底的に比較して、除外群に含まれていない包含群内のシグネチャ領域を見つけることを含む。しかしながら、これらのアプローチは効率的に拡大縮小せず、固定長の分子診断プライマーを生成することに焦点を当てている。

 他の洗練されたアプローチでは、迅速に検索可能なデータ構造のゲノムから固定サイズの部分文字列をコード化し、次にユニークな部分文字列(ref.10)についてこれらのデータ構造を分析する計算的に最適化された文字列処理アプローチを使用することにより、これらの手法は非常に速くスケール性が良くなるが、ターゲットシーケンスの変動性に対応できず、人為的に固定長シグネチャに限定されている。標的の可変性は、複数の配列アラインメント(ref.8)または他のクラスタリング操作(ref.11)を用いて類似の配列をグループ化することによって達成することができる。しかしながら、これらの一般的なクラスタリング技術は、高い計算コストを伴い、うまく拡張できない。いくつかのアルゴリズムは、不要な計算量を減らすために、クラスタリングの前にデータ削減ステップを組み込んでいる。たとえば、Insignia(ref.10 pubmed)、TOFI(ref.12 CrossRef)およびTOPSI(ref.13 CrossRef)は、包含ターゲット内の正確な一致および排除バックグラウンドを事前計算するために効率的なsuffix ツリーを使用する。しかしながら、バックグラウンドデータベースのサイズに依存して、これは計算上高価な演算のままである可​​能性がある。興味深い新規実装の1つはCaSSiS(ref.11 pubmed)である。これは、他のシグネチャ発見パイプラインよりも完全にシグネチャ発見問題に近づいている。このソフトウェアは、系統樹のような階層的にクラスタ化されたデータセット内のすべての位置について同時にシグネチャを生成し、それによってすべての可能なサブグループの候補シグネチャを生成する。しかし、このプロセスでは、入力データを計算上高価な系統樹などの階層的にクラスタ化された形式で提供する必要がある。効率と感度のトレードオフに加えて、これまでに発見されたプログラムのほとんどには、ゲノム集団間の共通の変異を特定するのに不適切な欠点がある。例えば、それらは分析を単一の包含ゲノム(ref.12)に制限するか、標的識別のためのユーザ供給ゲノムを許可しないか(ref.10)、またはエンドユーザにソフトウェアを提供しない可能性がある(ref.8)。

 著者らは、識別可能なバクテリア配列シグネチャを発見し、任意のゲノム配列群の比較分析を効率的で正確な斬新な方法で行うためのシステムとしてNeptuneを設計した。 Neptuneは、ユーザー指定の関心がある配列グループ間でに共有されるが、背景グループには欠けているゲノム遺伝子座を同定する。 事前計算、ターゲットへの制限、および低速クラスタリングアプローチとは無関係に、Neptuneは基準に基づいた並列化された完全一致のk-mer戦略を速度に適用し、不正確な一致には感度を向上させる。 ネプチューンシグネチャ発見は、統計的な信頼の尺度を用いて決定を行う確率論的モデルによって導かれる。 Neptuneはgithub.com/phac-nml/neptuneで無料で入手できるオープンソースのソフトウェアであり、細菌集団の迅速な比較評価に広く適用できる。

 

f:id:kazumaxneo:20180602103943j:plain

 図1。オーバービュー。論文より転載。

 

インストール

mac os 10.13のPython 2.7.12 :: Anaconda 4.2.0でテストした。

依存

本体 Github

https://github.com/phac-nml/neptune

MinicondaとBiocondaをいれてない人は公式マニュアル参照(リンク)。

condaでインストールする。

conda install neptune

neptune -h

$ neptune -h

usage: neptune-conda -i INCLUSION [INCLUSION ...] -e EXCLUSION 

    [EXCLUSION ...] -o OUTPUT

 

Neptune identifies signatures using an exact k-mer matching strategy. Neptune

locates sequence that is sufficiently present in many inclusion targets and

sufficiently absent from exclusion targets.

 

optional arguments:

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

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

 

REQUIRED:

  -i INCLUSION [INCLUSION ...], --inclusion INCLUSION [INCLUSION ...]

                        The inclusion targets in FASTA format.

  -e EXCLUSION [EXCLUSION ...], --exclusion EXCLUSION [EXCLUSION ...]

                        The exclusion targets in FASTA format.

  -o OUTPUT, --output OUTPUT

                        The directory to place all output.

 

KMERS:

  -k KMER, --kmer KMER  The size of the k-mers.

  --organization ORGANIZATION

                        The degree of k-mer organization in the output files.

                        This exploits the four-character alphabet of

                        nucleotides to produce several k-mer output files,

                        with all k-mers in a file beginning with the same

                        short sequence of nucleotides. This parameter

                        determines the number of nucleotides to use and will

                        produce 4^X output files, where X is the number of

                        nucleotides specified by this parameter. The number of

                        output files directly corresponds to the amount of

                        parallelization in the k-mer aggregation process.

 

FILTERING:

  --filter-percent FILTER-PERCENT

                        The maximum percent identity of a candidate signature

                        with an exclusion hit before discarding the signature.

                        When both the filtered percent and filtered length

                        limits are exceed, the signature is discarded.

  --filter-length FILTER-LENGTH

                        The maximum shared fractional length of an exclusion

                        target alignment with a candidate signature before

                        discarding the signature. When both the filtered

                        percent and filtered length limits are exceed, the

                        signature is discarded.

  --seed-size SEED-SIZE

                        The seed size used during alignment.

 

EXTRACTION:

  -r REFERENCE [REFERENCE ...], --reference REFERENCE [REFERENCE ...]

                        The FASTA reference from which to extract signatures.

  --reference-size REFERENCE-SIZE

                        The estimated total size in nucleotides of the

                        reference. This will be calculated if not specified.

  --rate RATE           The probability of a mutation or error at an arbitrary

                        position. The default value is 0.01.

  --inhits INHITS       The minimum number of inclusion targets that must

                        contain a k-mer observed in the reference to begin or

                        continue building candidate signatures. This will be

                        calculated if not specified.

  --exhits EXHITS       The maximum allowable number of exclusion targets that

                        may contain a k-mer observed in the reference before

                        terminating the construction of a candidate signature.

                        This will be calculated if not specified.

  --gap GAP             The maximum number of consecutive k-mers observed in

                        the reference during signature candidate construction

                        that fail to have enough inclusion hits before

                        terminating the construction of a candidate signature.

                        This will be calculated if not specified and is

                        determined from the size of k and the rate.

  --size SIZE           The minimum size of all reported candidate signatures.

                        Identified candidate signatures shorter than this

                        value will be discard.

  --gc-content GC-CONTENT

                        The average GC-content of all inclusion and exclusion

                        targets. This will be calculated from inclusion and

                        exclusion targets if not specified.

  --confidence CONFIDENCE

                        The statistical confidence level in decision making

                        involving probabilities when producing candidate

                        signatures.

 

PARALLELIZATION:

  -p PARALLELIZATION, --parallelization PARALLELIZATION

                        The number of processes to run simultaneously. Note

                        that this is only applicable when running Neptune in

                        non-DRMAA mode (default).

 

DRMAA:

  --drmaa               Whether or not to run Neptune in DRMAA-mode and

                        attempt to schedule jobs through DRMAA. This will

                        require setting up DRMAA in advance.

  --default-specification DEFAULT-SPECIFICATION

                        The default DRMAA parameters.

  --count-specification COUNT-SPECIFICATION

                        The DRMAA parameters specific to k-mer counting.

  --aggregate-specification AGGREGATE-SPECIFICATION

                        The DRMAA specific parameters specific to k-mer

                        aggregation.

  --extract-specification EXTRACT-SPECIFICATION

                        The DRMAA parameters specific to candidate signature

                        extraction.

  --database-specification DATABASE-SPECIFICATION

                        The DRMAA parameters specific to database construction

                        and querying.

  --filter-specification FILTER-SPECIFICATION

                        The DRMAA parameters specific to candidate signature

                        filtering.

  --consolidate-specification CONSOLIDATE-SPECIFICATION

                        The DRMAA parameters specific to filtered signature

                        consolidation.

——

 

ラン

ランには、ターゲットとなるゲノムのFASTA、除外したいゲノムのFASTAディレクトリ単位で与える必要がある。

 Githubレポジトリ(リンク)のテストデータをランする。inc_dir/にあるFASTAにあって、exc_dir/にあるFASTAにない特徴が自動検出され、統合された結果がoutput/のconsolidated.fasta に出力される。

neptune -p 12 --organization 3 --inclusion inclusion_dir/ --exclusion exclusion_dir/ --output output/
  • -i    A list of inclusion targets in FASTA format. You may list multiple file or directory locations following the parameter. Neptune will automatically include all files within directories. However, Neptune will not recurse into additional directories.-i A list of inclusion targets in FASTA format. You may list multiple file or directory locations following the parameter. Neptune will automatically include all files within directories. However, Neptune will not recurse into additional directories.
  • -e    A list of exclusion targets in FASTA format. You may list multiple file or directory locations following the parameter. Neptune will automatically include all files within directories. However, Neptune will not recurse into additional directories. 
  • -o    The location of the output directory. If this directory exists, any files produced with existing names will be overwritten. If this directory does not exist, then it will be created.
  • -p    The number of processes to run simultaneously. Note that this is only applicable when running Neptune in non-DRMAA mode (default).
  • --organization    The degree of k-mer organization in the output files. This exploits the four-character alphabet of nucleotides to produce several k-mer output files, with all k-mers in a file beginning with the same short sequence of nucleotides. This parameter determines the number of nucleotides to use and will produce 4^X output files, where X is the number of nucleotides specified by this parameter. The number of output files directly corresponds to the amount of parallelization in the k-mer aggregation process.

 

出力されるファイル一覧。

f:id:kazumaxneo:20180602095840j:plain

詳細はマニュアル(リンク)で解説されている。

すべての参照のソートされたシグネチャは、 統合され単一の consolidated.fastaファイルに出力される。この出力が、フィルタリングされソートされたNeptuneの最終出力となる。

 > cat output/consolidated/consolidated.fasta 

$ cat output/consolidated/consolidated.fasta 

>1.0 score=1.0000 in=1.0000 ex=0.0000 len=103 ref=inclusion2 pos=99

TAGTCTCCAGGATTCCCGGGGCGGTTCAGATAATCTTAGCATTGACCGCCTTTATATAGAAGCTGTTATTCAAGAAGCATTTTCAAGCAGTGATGTAAGAAAA

>1.1 score=0.9979 in=0.9979 ex=0.0000 len=640 ref=inclusion2 pos=3494

CGCGGGCGATATTTTCACAGCCATTTTCAGGAGTTCAGCCATGAACGCTTATTACATTCAGGATCGTCTTGAGGCTCAGAGCTGGGAGCGTCACTACCAGCAGATCGCCCGTGAAGAGAAAGAGGCAGAACTGGCAGACACATGGAAAAAGGCCTGCCCCAGCACCTGTTTTGAATCGCTATGCATCGATCATTTGCAACGCCACGGGGCCAGCAAAAAAGCCATTACCCGTGCGTTTGATGACGATGTTGAGTTTCAGGAGCGCATGGCAGAACACATCCGGTACTGGTTAAACCATTGCTCACCACCAGGTTGATATTGATTCAGAGGTATAAAACGAATGAGTACAGCACTCGCAACGCTGGCAGGGAAGCTGGCTGAACGTGTCGGCATGGATTCTGTCGACCCACAGGAACTGATCACCACTCTTCGCCAGACGGCATTTAAAGGTGATGCCAGCGATGCGCAGTTCATCGCATTGCTGATCGTCGCCAACCAGTACGGTCTTAATCCGTGGACGAAAGAAATTTACGCCTTTCCTGATAAGCAGAACGGCATCGTTCCGGTGGTGGGCGTTGATGGCTGTCCCGTATCATCAATGAAAACCAGCAGTTTGAGGCATGGTACTTTGAGCAGGACA

>0.2 score=0.9966 in=0.9966 ex=0.0000 len=98 ref=inclusion1 pos=5209

GCGAGTTTTGCGAGATGGTGCCGGAGTTCATCGAAAAAATGGACGAGGCACTGCTGAAATTGGTTTTGTATTTGGGGAGCAATGGCGATGAAGCATCC

——

スコアで降順に出力される(sensitivityとspecificityが高いとスコアが高くなる)。

 

ディレクトリ内の特定のFASTAだけ指定することもできます(マニュアル)。

Walkthrough(出力について)

https://phac-nml.github.io/neptune/walkthrough/

 

引用

Neptune: a bioinformatics tool for rapid discovery of genomic variation in bacterial populations

Eric Marinier, Rahat Zaheer, Chrystal Berry, Kelly A. Weedmark, Michael Domaratzki, Philip Mabon, Natalie C. Knox, Aleisha R. Reimer, Morag R. Graham, Linda Chui, Laura Patterson-Fortin, Jian Zhang, Franco Pagotto, Jeff Farber, Jim Mahony, Karine Seyer, Sadjia Bekal, Cécile Tremblay, Judy Isaac-Renton, Natalie Prystajecky, Jessica Chen, Peter Slade, and Gary Van Domselaar

Nucleic Acids Res. 2017 Oct 13; 45(18): e159.

 

ロングリードを使ってSVを検出する Picky

 

 ゲノム構造変異の獲得(SV)は、ガンゲノムの主要な特徴であるが、ショートリードシーケンシングデータから再構成することは困難である。ここでは、カスタマイズされたパイプライン、Picky(https://github.com/TheJacksonLaboratory/Picky)を使用し、ナノポアプラットフォームのロングリードを使って乳癌モデルにおける多様なアーキテクチャのSVを明らかにした。著者らはショートリード解析に比べて優れた特異性および感度を有するSVの全スペクトルを同定し、主要な変異源としてリピートDNAを明らかにした。 SVに関連する共通の構造的特徴として、ヌクレオチド解明されたマイクロ挿入におけるゲノムワイドなブレークポイントの検討。ゲノムにわたるブレークポイント密度は、染色体間連結性の傾向と関連しており、ゲノムのプロモーターおよび転写領域が豊富であることが判明した。さらに、染色体二重交叉から段階的SVに至る相互転座の過剰発現を観察した。著者らは、Picky分析が、ロングリードデータからガンゲノムのSVを包括的に検出する有効な手段であることを示す。

 

Github Motivationsより

 2016年4月時点では、Nanoporeのロングリードに利用できるエンドツーエンドの構造変異検出パイプラインはなかった。著者らはbwa-memのような既存のツールを使って解析を始めたが、機能はするがいくつかの点では最適ではないことに気づいた。例えば、ショートリードの世界では、リードがtandem duplicationを完全にスパンする可能性はほとんど無視できる。従って、大部分の構造変異同定ソフトウェアはtandem duplicationを同定できない。最も重要なのは、SVが識別できるかどうかは、SVコーラーに提供されたアラインメントに依存することである(アライメントが一番大事)。ほとんどのアライナーは低いシークエンシングエラー率を仮定しており、ロングリードでは最適なアライメントを提供しない可能性がある。そこで、ロングリードアラインメントとSVコールの2つの領域を探索することにした。

 SVをコールするには適切なアライメントが重要であるため、bwa-mem、lastal、damapper、graphmap、NextGenMap-LR、blastnを使ってSV(tandem duplication, inversion, deletion, and insertion)を含むロングリードアライメント能力を探索した。

lastal、blastn、そしてdamapperは、非常に似たアライメントパフォーマンス(計算時間パフォーマンスではない)を持つことがわかった。 Bwa-memも近いが、感度の差が無視できなかった。graphmapは、「グローバル」アラインメントのフィットソートを強制的に行う傾向がある。 NextGenMap-LRは有望だが、PacBioのロングリードに焦点を当てていて、ONTのリードでは不安定だったため、使用できなかった。これは徹底的な比較ではないが、実際にアライナーを選択するには役立つ。最終的に、ゲノムとゲノムとのアライメントの文脈で使用されるlastalを、感度と計算時間両方を考慮して使うことに決定した。

 いろいろなSVコーラーを調べると、ほとんどのツールがSVの1つのサブセットに焦点を絞っていて、そこで非常にうまく機能していることが明らかである。また、ショートリードのシナリオで最適化されているため、ロングリードによるSV検出シナリオで問題になることもある。 .vcf形式のレポートは良いスタンダードだが、主にブレークポイントレベルで提示された場合、生物学者が概念レベルでSVを追跡することを難しくしている。著者はコミュニティがより多くのシナリオを想定してスタンダオードの基準を改善し続けることは疑いがない。今のところ、本ツールは都合が良いタブ区切りのテキストファイル形式でレポートし、近い将来に.vcf形式に対応する。

 

インストール

cent os6でテストした。

依存

  • LAST

LASTはbrewで導入できる(LAST紹介)。

 

本体 Github

GitHub - TheJacksonLaboratory/Picky: Structural Variants Pipeline for Long Reads

git clone https://github.com/TheJacksonLaboratory/Picky.git
cd Picky/src/

> perl picky.pl -h

$ perl picky.pl -h

Please specify the command.

 

picky.pl <command> -h

 

<command> [hashFq, selectRep, callSV]

hashFq    : hash read uuids to friendly ids

lastParam : Last parameters for alignment

selectRep : select representative alignments for read

callSV    : call structural variants

xls2vcf   : convert Picky sv xls file to vcf

preparepbs: chunk last fastq file and write pbs script for cluster submission

script    : write a bash-script for single fastq processing

 

 

      License = JACKSON LABORATORY NON-COMMERCIAL LICENSE AGREEMENT

                https://github.com/TheJacksonLaboratory/Picky/blob/master/LICENSE.md

Documentation = https://github.com/TheJacksonLaboratory/Picky/wiki

   Repository = https://github.com/TheJacksonLaboratory/Picky/

 

 

ラン

faToFastqを使い、フェイクqualityのfastqを作成する。なければダウンロードする。

curl -O http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToFastq 
chmod u+x faToFastq

  ScrappieでベースコールしたONTのfastaをダウンロードし、fastqにする。

curl -O http://s3.amazonaws.com/nanopore-human-wgs/na12878.chr20ScrappieFiltered.fasta

#フェイクfastqに変換
./faToFastq -qual=H na12878.chr20ScrappieFiltered.fasta LongRead.fastq

 hg19のindexを作成する。

lastdb -v -P 10 hg19.lastdb hg19.fa

#sequence dictionaryの作成
samtools dict -H hg19.fa > hg19.seq.dict

以下のコマンドでアライメントとSV検出を実行するためのシェル スクリプトが作成される。

./picky.pl script --fastq LongRead.fastq --thread 20 > LongRead.sh

 実行権をつける。

chmod u+x LongRead.sh

中身を確認。

# general installation

export LASTAL=last-755/src/lastal

export PICKY=./picky.pl

 

# general configuration

export LASTALDB=hg19.lastdb

export LASTALDBFASTA=hg19.fa

export NTHREADS=20

 

# FASTQ file

export RUN=LongRead

 

# reads alignments

time (${LASTAL} -C2 -K2 -r1 -q3 -a2 -b1 -v -v -P${NTHREADS} -Q1 ${LASTALDB} ${RUN}.fastq 2>${RUN}.lastal.log | ${PICKY} selectRep --thread ${NTHREADS} --preload 6 1>${RUN}.align 2>${RUN}.selectRep.log)

 

# call SVs

time (cat ${RUN}.align | ${PICKY} callSV --oprefix ${RUN} --fastq ${RUN}.fastq --genome ${LASTALDBFASTA} --exclude=chrM --sam)

 

# generate VCF format

${PICKY} xls2vcf --xls ${RUN}.profile.DEL.xls --xls ${RUN}.profile.INS.xls --xls ${RUN}.profile.INDEL.xls --xls ${RUN}.profile.INV.xls --xls ${RUN}.profile.TTLC.xls --xls ${RUN}.profile.TDSR.xls --xls ${RUN}.profile.TDC.xls > ${RUN}.allsv.vcf

hg38でindexを作っているなら、hg19のとこも修正する。 

 

LASTのパスはbrewで通していて間違っていたので、2行目の

export LASTAL=last-755/src/lastal

 ↓

export LASTAL=lastal

に修正した。

 

 スクリプトを実行。

./LongRead.sh

 

ランが終わると VCFが出力される(すでにVCFに対応している)。模擬fastq(合計19億bp)からhg19のアライメントには3時間ほどかかった。

$ head -n 30 LongRead.allsv.vcf |tail -n 10

##FILTER=<ID=singleton,Description="Single read">

##FILTER=<ID=CIPOS,Description="The CIPOS is greater or less than 20">

##FILTER=<ID=CIEND,Description="The CIEND is greater or less than 20">

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT LongRead

chr16 33968730 ID1 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33968756;SVTYPE=DEL;RE=3;RNAMES=nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch253_read1117_strand.fast5,DEAMERNANOPORE_20161206_FNFAB49164_MN16450_sequencing_run_MA_821_R9_4_NA12878_12_06_16_71094_ch489_read1323_strand.fast5,nanopore2_20161128_FNFAB49712_MN17633_sequencing_run_20161128_Human_Qiagen_1D_R9_4_64849_ch247_read294_strand.fast5;SVLEN=27;CIPOS=-2,434;CIEND=-3,743;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(3) GT ./.

chr16 33978912 D2 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33979358;SVTYPE=DEL;RE=2;RNAMES=LomanLabz_PC_20161125_FNFAF01132_MN17024_sequencing_run_20161124_Human_Qiagen_1D_R9_4_96293_ch470_read2257_strand.fast5,PLSP61583_20161015_FNFAB42316_MN17048_sequencing_run_Hum94_ss_jt_86199_ch225_read262_strand1.fast5;SVLEN=447;CIPOS=-90,90;CIEND=-211,212;NOTE=CIPOS_CIEND;ISVTYPE=DEL(1),INDEL(1) GT ./.

chr16 33979001 D3 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33980791;SVTYPE=DEL;RE=9;RNAMES=DEAMERNANOPORE_20161206_FNFAB49164_MN16450_sequencing_run_MA_821_R9_4_NA12878_12_06_16_71094_ch74_read4369_strand.fast5,LomanLabz_PC_20161128_FNFAF01253_MN17024_sequencing_run_20161128_Human_Qiagen_1D_R9_4_50330_ch444_read8778_strand.fast5,PLSP61583_20161128_FNFAB49914_MN17048_sequencing_run_Hu_Nott_Bi_FC4_tune_45519_ch481_read148_strand.fast5,MinION2_20161115_FNFAB46664_MN20093_sequencing_run_Chip105_Genomic_R9_4_450bps_tune_16109_ch48_read301_strand.fast5,FMH_15Le080325s_20161103_FNFAB42798_MN17638_sequencing_run_161103_Human5_LSK108R9_4_13493_ch170_read2217_strand.fast5,nanopore2_20161209_FNFAF04090_MN17633_sequencing_run_20161209_Human_Rapidv2_29324_ch79_read58_strand.fast5,nanopore2_20161122_FNFAF01169_MN17633_sequencing_run_20161122_Human_Qiagen_1D_R9_4_65629_ch138_read425_strand.fast5,MinION2_20161013_FNFAB42706_MN16454_sequencing_run_Chip99_Genomic_R9_4_480bps_69747_ch509_read42_strand.fast5,nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch42_read1895_strand.fast5;SVLEN=1791;CIPOS=-5,11;CIEND=-1,5;NOTE=CIPOS_CIEND;ISVTYPE=DEL(9) GT ./.

chr16 33983845 ID4 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33985021;SVTYPE=DEL;RE=2;RNAMES=MinION2_20161013_FNFAB42706_MN16454_sequencing_run_Chip99_Genomic_R9_4_480bps_69747_ch509_read42_strand.fast5,LomanLabz_PC_20161128_FNFAF01253_MN17024_sequencing_run_20161128_Human_Qiagen_1D_R9_4_50330_ch444_read8778_strand.fast5;SVLEN=1177;CIPOS=-120,121;CIEND=-6,7;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(2) GT ./.

chr16 33985139 ID5 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33985230;SVTYPE=DEL;RE=2;RNAMES=nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch121_read1284_strand.fast5,nanopore2_20161125_FNFAF01127_MN17633_sequencing_run_20161124_Human_Qiagen_1D_R9_4_80573_ch42_read1895_strand.fast5;SVLEN=92;CIPOS=-21,21;CIEND=-38,38;NOTE=CIPOS_CIEND;ISVTYPE=INDEL(2) GT ./.

chr16 33985190 D6 N <DEL> . PASS IMPRECISE;SVMETHOD=picky;END=33986368;SVTYPE=DEL;RE=2;RNAMES=nanopore2_20161209_FNFAF04090_MN17633_sequencing_run_20161209_Human_Rapidv2_29324_ch79_read58_strand.fast5,nanopore2_20161122_FNFAF01169_MN17633_sequencing_run_20161122_Human_Qiagen_1D_R9_4_65629_ch138_read425_strand.fast5;SVLEN=1179;CIPOS=-2,3;CIEND=-10,11;NOTE=CIPOS_CIEND;ISVTYPE=DEL(2) GT ./.

(anaconda2-4.2.0) [uesaka@cyano src]$ 

 

 

引用

Picky comprehensively detects high-resolution structural variants in nanopore long reads.

Gong L, Wong CH, Cheng WC, Tjong H, Menghi F, Ngan CY, Liu ET, Wei CL

Nat Methods. 2018 Apr 30. doi: 10.1038/s41592-018-0002-6.

 

CircosをWeb上 で利用できる ClicO FS

 2020 7/24 関連文献追記

 

 Circos(Krzywinski et al、2009)(HP)は、ビジュアルデータを環状形式で表現するPerl言語ベースのツールである。ネイティブのCircosソフトウェアは、コマンドラインインターフェイスCLI)を介して提供されている。ソフトウェアのインストールと設定は、コマンドラインの経験を持つユーザーにとっては難しくないが、CLIに精通していないユーザーにはある程度の難しさがある。

 Circoletto、CIRCUS、RCircos(Darzentas、2010; Naquin et al、2014; Zhang et al、2013)を含む環状画像のプロットを可能にするCircosの他の選択肢がある。 CircolettoとCIRCUSは、より特殊化された機能のために開発された。前者(Circoletto)はBLAST(Altschul et al、1990)の結果とCircosをウェブベースのユーザインターフェース(UI)を介して統合する一方で、後者(CIRCUS)はゲノムの構造変異をプロットすることに特化している。しかし異なるデータ型を持つプロット生成に対する柔軟性。を持っていない。RCircosはCircosのdata trackのプロットをサポートしているが、RパッケージのCLIに基づいている。

 上記のアプリケーションの中には、アプリケーションを十分に活用するためにセットアップ、設定、コーディングのスキルが必要なものがある。 CLI環境での経験がほとんどまたはまったくないユーザーは、厳しい学習曲線に遭遇する可能性がある。この問題に対処するために、ClicO FS(Open Circular Layout Interactive Converter)が開発された。 ClicO FSは、ユーザーがCircosプロットを簡単に生成できる、使いやすいWebサービスとなっている。 Circosのすべての機能がClicO FSによって実行されるわけではない(例えば、染色体軸破損など)。

 ClicO FSは、バックエンド技術としてRuby on Railsを使用し、フロントエンドフレームワークとしてTwitter Bootstrapを使用して構築されている。 ClicO FSは、ユーザーアップロードによって必要なファイルを受け取る。 Circosでは、(i)カリオタイプ(karyotype  wiki)、(ii)データおよび(iii)構成ファイルの3種類のファイルが必要だが、ClicO FSでは最初の2種類のファイルのみが必要となる。カリオタイプファイルは、典型的には染色体、配列コンティグ、または生物学的文脈におけるクローンである軸を定義する固有の識別子、ラベル、サイズ、色など、各軸の情報を含む7列のフォーマットファイルである。データファイルは、さまざまなタイプのデータを描画するための別の入力ファイルとなる。現在、ClicO FSでサポートされているデータタイプには、リンク、ヒストグラム、ラインプロット、散布図、ヒートマップ、タイル、テキスト、コネクタおよびハイライトが含まれる。データファイルフォーマットは単純であり、一般的には描画されるデータトラックに応じて3〜7個の列を含む。ファイル形式の詳細な説明は、ClicO FSインターフェイスで指定する。 ClicO FSで受け入れられるすべてのファイルは、2つのプログラムを切り替えるときの使いやすさを高めるため、Circosと同じフォーマットとなっている。

 ClicO FSでは、ユーザーが設定ファイル(configuration file)をアップロードする必要はない。すべての設定はユーザーインターフェイスで設定できる。 ClicO FSはユーザの利便性を重視して設計されているため、一般的ではない設定にはデフォルトのパラメータが与えられる。しかし、上級ユーザは、ClicO FSの特別なテキストボックスに、背景、ルール、軸、ラベルスニグリング、ハイライトなどの多くの設定を入力することができる。

 いくつかの追加の機能が提供されている。例えば、様々なカラースキームをカリオタイプに割り当てることができる。また、設定を保存しなくても、いつでも出力をプレビューできる。登録されたClicO FSユーザーは、複数の独立したプロジェクトを作成できる。関連するすべてのデータおよび構成は、それぞれのプロジェクト専用のデータベースに格納される。したがって、ユーザーは、次回のログイン時にプロジェクトを検索したり、引き続き作業を続けることができるようになっている。それに加えて、ClicO FSは、いくつかの一般的な形式のプレゼンテーションのスピードアップを可能にする、事前設定されたテンプレートベースのアプローチを提供する。例えば、GenBankデータベースは、配列データの最も一般的なソースの1つである。ユーザーはGenBank形式のファイルをClicO FSにアップロードできる。 5回のクリックで、遺伝子、GC含量、ゲノムの構造RNAなどの情報を視覚化することができる。

Circosの機能性は多様でダイナミックである。 Circoの機能に合わせてClicO FSの機能を拡張する。この問題に対処するために、ClicO FSのユーザーは、作成したプロジェクトのすべての構成ファイルとデータファイルをダウンロードすることができる。これらのファイルはテンプレートとして、ClicO FSでサポートされていない要素を追加したり、ローカル端末でCircosを使って実行できるようになる。

 

 

ラン

利用するには、はじめにユーザ登録が必要となる。

ClicO FS | Circular Layout Interactive Converter Free Services

f:id:kazumaxneo:20180531235108j:plain

数分で登録作業は終了する。パスワードを設定したらログインできるようになる。

 

ヘルプ

http://codoncloud.com:3000/help

 

ログインしたら、Create Projectをクリックし、新しいプロジェクトを作成する。

f:id:kazumaxneo:20180531221737j:plain

 

プロジェクト名を記載する。

f:id:kazumaxneo:20180531221943j:plain

右のopenマークをクリックする。

f:id:kazumaxneo:20180531222104j:plain

 

 

karyotypeファイルをアップロードする。

f:id:kazumaxneo:20180531222112j:plain

 

ここでは、circosのlibraryのhumanのファイルを使った(circos/0.67-7/libexec/data/karyotype/karyotype.human.hg38.txt)。ただし性染色体、オルガネラゲノム、bandingなどの行は除いた。

f:id:kazumaxneo:20180531222539j:plain

boxにレ点チェックをつけて、右の→をクリック。

f:id:kazumaxneo:20180531223434j:plain

 

chr名、色などを設定する。手動でもできるが、ここでは設定一覧からUCSCを選んで、Apply & Saveを押し、自動で色を割り振った。上のPreviewを押せばプレビューできる。最後に左上のもう1つあるsaveボタンを押すのを忘れないようにする(忘れると反映されない)。右端の矢印をクリックし、次の設定に進む。

f:id:kazumaxneo:20180531224736j:plain

 

一般設定。設定を変更して、プレビューで確認できるのは非常に便利。saveを押さないと設定は保持されないことに注意する。

f:id:kazumaxneo:20180531224215j:plain

Label settings。右矢印をクリックして、どんどん右の設定に進んでいく。

f:id:kazumaxneo:20180531225851j:plain

 

Ticks settings。

f:id:kazumaxneo:20180531230026j:plain

Chromosome settings。

f:id:kazumaxneo:20180531230309j:plain

 

Chromosomes Ordering。ドラッグして順番を変更できる。

f:id:kazumaxneo:20180531230352j:plain

 

 

カリオタイプの設定が終わったら、次にデータトラックを追加する。データトラックについては、ユーザーによって何を描くかで最適なタイプが異なってくる(例えば染色体間の組み替えならリンクプロット、オーム解析データならヒートマッププロットやスキャッタープロットなど)。

下の画像を参考に、どんな種類のプロットを描画するか決める(例えば、各染色体間で高い相同性を持つ領域を描画したければ "Links" が使える)。

f:id:kazumaxneo:20180531234300j:plain

 

項目のフォーマットの詳細はcircosのHP(リンク)や、画面右上にあるClicO FSのhelpを押して確認できる。

f:id:kazumaxneo:20180601000201j:plain

 

アップロードするタイプを決めて、アップロードする。ここではcircosのexampleのヒートマップをアップロードした(circos/example/data/heatmap.hs.mm.5e6.txt)。

f:id:kazumaxneo:20180601121046j:plain

 

チェックをつけてセッティングに進む。

f:id:kazumaxneo:20180601120719j:plain

 

最後に右上のOutputボタンをクリックして画像を出力する。

f:id:kazumaxneo:20180601121413j:plain

データはSVGadobeツールと互換性はない)かPNGでダウンロードできる。 

 

リンクも追加してみる。

f:id:kazumaxneo:20180601155836j:plain

 

透明度は0.5、リボンをONにした。

f:id:kazumaxneo:20180601155859j:plain

f:id:kazumaxneo:20180601155931j:plain

 さらに別ファイルのリンクを追加。

f:id:kazumaxneo:20180601161435j:plain

追加リンクはリボンオフにして色も変えた。このように、同じリンクでも、別ファイルで追加することで別の設定を適用できる。

 

 

Genbankからuploadするには、上のメニューからGet startを選択する。

f:id:kazumaxneo:20180601122621j:plain

Plot from Genbankを選択。

f:id:kazumaxneo:20180601122629j:plain

テストした時は、エラーが繰り返し出て進めなかった。

 

 

ログイン方式なので、セーブしておけば結果を別の端末で読み込むこともできます。保存期間については規約を確認してください。複雑なCIrcosを使いやすくする良ツールだと思います。

 

 

 

 

追記

Circaというツールもあります(macwindowslinux対応)。100ドルしますが使いやすそうです。

引用

ClicO FS: an interactive web-based service of Circos

Wei-Hien Cheong, Yung-Chie Tan, Soon-Joo Yap, and Kee-Peng Ng

Bioinformatics. 2015 Nov 15; 31(22): 3685–3687.

 

2020 7/24 

Galaxyで使えるCIrcosが発表されています。また機会をみて紹介します。

Galactic Circos: User-friendly Circos plots within the Galaxy platform
https://academic.oup.com/gigascience/article/9/6/giaa065/5856406