macでインフォマティクス

macでインフォマティクス

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

タンパク質のホモリピートを分析するwebサーバー dAPE

 

 Low Complexity(LC)は、タンパク質中のタンデムリピートおよびcompositionally biased regions(CBR)のようなアミノ酸組成にほとんど多様性がない領域を説明するために使用される一般用語である。ホモリピート、またはpolyX領域は、単一のアミノ酸残基の連続からなるCBRである。ヒトのホモリピートを有するタンパク質の割合はわずかな変動があり、5.7%(Jorda and Kajava、2010)から7.7%(51 778のヒト配列から約4000タンパク質)となっている。この多様性の理由は、偏った領域をホモリピートと見なすための標準化されたカットオフの欠如、およびそれがそのような領域内にミスマッチを含むかどうかということである。ホモリピートを定義するためのいくつかのカットオフが文献に記載されている:少なくとも5つから少なくとも7つまで(Jorda and Kajava、2010)、少なくとも5つの同一残基(Alba andGuigó、2004)、少なくとも6(Lobanov and Galzitskaya、2012)、 10残基ウィンドウで8以上(Schaefer et al、2012)、5つのうち4つ(UniProt、http://www.uniprot.org/help/compbias)。実際、ホモリピート検出のカットオフは、アミノ酸タイプだけでなく種にも依存している可能性がある(Lobanov et al、2016)。

 (一部略)

原核生物よりも真核生物プロテオームにおいてより豊富であると記載されているが(Faux et al、2005; Jorda and Kajava、2010)、ホモリピートがどのように進化するかは明らかではない。それらの経時的な発展はまだ完全には記載されていないので、それらが進化に沿ってそれらの長さおよび機能を増加、減少またはランダムに変化させるかどうかは知られていない。

 本著者らは、ホモリピートの進化とそのタンパク質の状況を評価するのに役立つウェブサーバーとデータベースであるdAPEを開発することによって、これらの問題に取り組む。これは、デフォルトでは弱いカットオフを使用して、出現または消失するホモリピートの識別に役立ち、ユーザーが使いやすい設定で独自のシーケンスを研究することを可能にする。

 

dAPEに関するツイート 


使い方

http://cbdm-01.zdv.uni-mainz.de/~munoz/polyx/ にアクセスする。

f:id:kazumaxneo:20190402173449j:plain

 

f:id:kazumaxneo:20190402230906p:plain


作成中

 

 

 

 

引用

dAPE: a web server to detect homorepeats and follow their evolution
Pablo Mie, Miguel A Andrade-Navarro

Bioinformatics. 2017 Apr 15; 33(8): 1221–1223

初めてコマンドを使う人向けの解説:その1、指定した領域から配列を抽出する

2019 9/20追記 ゲノムのダウンロード

 

この記事では、初めてコマンドで動作するツールを使う方向けにゲノムの指定した領域から配列を抽出する方法について説明します。コンピュータはmacを想定しています。普通はpython3やanacondaを入れ、condaのコマンドを使ってBicoondaのプリビルドバイナリを導入するのが簡単ですが、$PATHの設定などでつまづくかもしれません。そこで、計算速度が遅くなることは気にせず、dokcerを使って仮想環境でBEDToolsを動かします。目的は、指定した領域の塩基配列を抽出し、エンリッチされた配列があるか調べることです。

 

1、macosのバージョンチェック

dockerをインストールするためにはmac10.12以上が必要です(間も無く10.13以上が最低ラインになるかもしれません)。また、メモリも8GBは必要です。できるだけ新しいmacを使ってください。

f:id:kazumaxneo:20190919161537j:plain

上のマシンはmacos10.12なのでインストール可能。メモリは96GBになっている。ストレージ(HDDやSSD)とメモリを混同しないように注意。

 

 

、dockerのインストール

Qiitaの記事を参考に、dockerをインストールします。


dockerインストールに成功したら、右上にクジラのアイコンが表示されているはずです。

f:id:kazumaxneo:20190919163339j:plain


クジラのアイコンをクリックし、dockerの仮想環境で利用するメモリやCPU数を指定します。

f:id:kazumaxneo:20190919162213j:plain

mac proはCPUやメモリのリソースが多いため、CPU18、メモリ84GBまで増やしています。 下のボタンが緑色になったら使用できます。

 

 

、ターミナルを起動します。アプリケーション(commnad + shift + A)=>ユーティリティにあります。

f:id:kazumaxneo:20190919151618p:plain

起動したら、dockerと打ちます。

docker

 ヘルプが出ればインストールに成功しています。

f:id:kazumaxneo:20190919151713p:plain

 

4A 、ヒトゲノムアセンブリのダウンロード

Homo sapiens - Ensembl genome browser 97

f:id:kazumaxneo:20190920133334p:plain

右の方のDownload FASTA を選択。ウィンドウが出てくるのでゲストを選択。

f:id:kazumaxneo:20190920133424p:plain

dnaを選択。全ゲノムならばHomo_sapiens.GRCh38.dna.toplevel.fa.gzをダウンロードする。rmとsmはリピートがハードマスクかソフトマスクされたアセンブリになる。

f:id:kazumaxneo:20190920133713p:plain

解凍しておく。
 

またはNCBIその他からもダウンロードできる。

NCBI

https://www.ncbi.nlm.nih.gov/genome/51

f:id:kazumaxneo:20190920135434p:plain

genomeをクリックする。

 

 

4B 、データを作業フォルダに集める。

ここではデスクトップにhumanというフォルダを作り、そこにヒトのGRCh38アセンブリをコピーしました。

f:id:kazumaxneo:20190919164145j:plain

 

ターミナルで移動します。ターミナルでcd (change directory)と打ち、それからスペースキーを一度打って1スペース空け、humanのフォルダをターミナルにドラッグ&ドロップします。

f:id:kazumaxneo:20190919164803j:plain

 

こんな感じになるので、そのままエンターキーを押したら移動します。

f:id:kazumaxneo:20190919164940j:plain

 

ls(list segments)+Enterと打って用意したファイルがあるか確認します。

ls

 

 

、bedtoolsを含むイメージのpull

 dokcerhubというイメージ共有サイトからbedtoolsがプリインストールされたイメージをローカルにダウンロードします。ここではlindsayliangさんのイメージを使います。以下のコマンドをターミナルにコピペしてEnterして下さい。

docker pull lindsayliang/bedtools

f:id:kazumaxneo:20190919153942p:plain

 

 

、イメージの動作確認

全コンテナのpullが終わったら helpが表示されるか確認します。以下のコマンドをコピペして下さい。

docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools bedtools

これはローカルと仮想環境の共有ディレクトリをカレントディレクトリ(現在いるフォルダ)、run後の初期パスをカレントディレクトリ、exit後に消す設定でイメージをrunし、イメージにインストールされたbedtoolsのヘルプを表示しています。

うまく動けば以下のようなbedtoolsのhelpメッセージが表示されるはずです。

$ docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools bedtools

bedtools: flexible tools for genome arithmetic and DNA sequence analysis.

usage:    bedtools <subcommand> [options]

 

The bedtools sub-commands include:

 

[ Genome arithmetic ]

    intersect     Find overlapping intervals in various ways.

    window        Find overlapping intervals within a window around an interval.

    closest       Find the closest, potentially non-overlapping interval.

    coverage      Compute the coverage over defined intervals.

    map           Apply a function to a column for each overlapping interval.

    genomecov     Compute the coverage over an entire genome.

    merge         Combine overlapping/nearby intervals into a single interval.

    cluster       Cluster (but don't merge) overlapping/nearby intervals.

    complement    Extract intervals _not_ represented by an interval file.

    shift         Adjust the position of intervals.

    subtract      Remove intervals based on overlaps b/w two files.

    slop          Adjust the size of intervals.

    flank         Create new intervals from the flanks of existing intervals.

    sort          Order the intervals in a file.

    random        Generate random intervals in a genome.

    shuffle       Randomly redistrubute intervals in a genome.

    sample        Sample random records from file using reservoir sampling.

    spacing       Report the gap lengths between intervals in a file.

    annotate      Annotate coverage of features from multiple files.

(以下略)

 

 

、指定した領域のBEDファイルの準備

ここからが本番です。取ってきたい領域を指定したファイルを用意します。chr1の三箇所、chr3の一箇所から取ってくるには以下のようにします。excelを使っています。

f:id:kazumaxneo:20190919183006j:plain

 左端に染色体番号、2列目にスタート(上流側)ポジション、3列目にエンド(下流側)ポジションを記入。

 

Miなどのテキストエディタに配列をコピペ。region.bedというファイル名で、作成したフォルダに保存する(commnad + S)。

f:id:kazumaxneo:20190919183425j:plain

右上の改行コードボタンの表示がCR(キャリッジリターン)になっていたら、LF(ラインフィード)に変更。

 

 

配列の抽出

作業前にindexファイルを作成します。これは本の索引に相当するファイルで、プログラムが高速動作するためにゲノムのファイルと同じ場所にある必要があります。

docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
samtools faidx Homo_sapiens.GRCh38.dna.toplevel.fa

ゲノムが大きいので10分以上かかります。 取ってきたい領域が特定の染色体に限定されるなら"-r xxx"という風にindex領域を指定することで時間短縮できますが、ここでは全領域のindexファイルを作成します(コマンドに慣れているならrefgenieというツールも使えます)。 Homo_sapiens.GRCh38.dna.toplevel.fa.faiが出力されます。

 

bedtools(紹介)を動かします。途中にある\はコマンドの途中に改行する時に使います。

docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
bedtools getfasta -fi Homo_sapiens.GRCh38.dna.toplevel.fa \
-bed region.bed > reads.fasta

ファイル名が違うところは変えて下さい。

 

 reads.fastaが出力されます。lessコマンドで中身を見てみます。

less reads.fasta

f:id:kazumaxneo:20190919185437j:plain

出力されてますね。スペースキーで1ページ下にスクロールします。bで1ページ戻ります。qを押すとlessコマンドは終了します。 

 exitでターミナルは終了しますが、commnad + Qでターミナルをアプリごと閉じてもらっても構いません。 

 追記

取り出す領域が1つならsamtools faidxが早いです。

docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
samtools faidx Homo_sapiens.GRCh38.dna.toplevel.fa chr1:10000-10010 \
> reads.fasta

 

 

、モチーフ検索

MEME CHIPにアクセスして、エンリッチされているモチーフがあるか調べます。

http://meme-suite.org/tools/meme-chip

f:id:kazumaxneo:20190919201134j:plain

 

先ほど出力された.fastaファイルをアップロードします。

f:id:kazumaxneo:20190919201937j:plain

 

MEME CHIPについてはNGS_handson2015の説明を参考にして下さい。

https://github.com/suimye/NGS_handson2015/wiki/PeakCallAndMDA

 

関連


 

 

 

ONTのロングリードから抗生物質耐性遺伝子の分布を調べるwebサーバー NanoARG

 

 薬剤耐性(AMR)は、感染症を予防および治療する能力を損ない、世界的な公衆衛生の脅威になる[ref.1]。現在、抗生物質耐性による世界中の年間死亡者数は、2050年までに1,000万人を超えると推定されている[ref.2]。これに対応して、多くの国内および国際機関は、クリニックおよび環境設定の両方で監視を拡大することを求めている。特に、環境モニタリングは、抗生物質耐性細菌と抗生物質耐性遺伝子(ARG)の人間と農業のインプットだけでなく、耐性病原体の進化と拡散に寄与する要因についての洞察を提供する。(一部略)

 分子ベースのモニタリングは、環境内の抗生物質耐性を追跡するための培養ベースの技術よりも多くの利点を提供する。これは、複雑な微生物群集内のARGの運搬と移動に関する豊富な情報を回収できる可能性に関して特に当てはまる。培養ベースの手法は時間がかかり、一度に1つの対象種に関する情報のみを提供するため、AMRの拡散に寄与する重要な微生物生態学的プロセスを見落とす可能性がある。したがって、細菌宿主を超越することを懸念する「汚染物質」としてARGを直接ターゲットにすることが人気を集めている。特に、水平遺伝子導入(HGT)[ref.7]は、新しい耐性株の出現と微生物生態系におけるAMRの普及[ref.8]で重要な役割を果たす。細菌間でのARGの細胞間移動は、トランスポゾン、プラスミド、インテグロンなどの可動性遺伝因子(MGE)を介して促進される[ref.9]。インテグロンは、複数のARGの捕捉を促進するため、重要な遺伝的要素であり、したがって多剤耐性の普及の媒体として効果的に機能する[ref.10]。 HGTに関与するメカニズムには、接合、形質転換、形質導入、および相同組換えが含まれ、DNAは転位、複製、および統合によって組み込まれる[ref.9]。

 多剤耐性は、主要な臨床的課題として浮上している。たとえば、メチシリン耐性黄色ブドウ球菌MRSA)は、特にバンコマイシンに耐性がある場合、主要な病院感染症の原因であり、治療の選択肢はほとんどない[ref.11]。最近、ニューデリーメタロベータラクタマーゼ(blaNDM-1)は、強力なラストリゾートカルバペネム抗生物質に対する耐性をコードし、多剤耐性に関連する非常に可動性の高い遺伝エレメントに搭載されているため、主要な関心事として浮上している。この例は、理想的には、モニタリング技術により、ARGの迅速かつ堅牢な特性評価、およびMGEとの関連性、多剤耐性、病原体宿主による保菌を提供する必要があることを強調している。この点で、ショットガンメタゲノムシーケンス手法は、さまざまな環境で見つかった多様なARGアレイの特性評価の有望なツールとして登場した[ref.4、15、16、17]。特に、イルミナプラットフォーム[ref.18]や454パイロシーケンシング[ref.19、20]などのハイスループット次世代DNAシーケンス技術により、環境でのARGモニタリングの新しい次元展開が可能になった。

 

(複数段落省略)

 ナノポアのロングリードは、共起とモビリティの可能性の観点からARGのコンテキストを調査するユニークな機会を提供する。キメラ配列を生成する可能性のあるショートリードの新規アセンブリ[ref.48]とは異なり、ナノポアシーケンスは本質的に長い配列を生成するため、キメラの可能性が低下する。したがって、ナノポアシーケンスは、ARG、MGE、およびMRGの共存を特定するための強力なツールになる可能性がある。このようなアプローチは、環境モニタリングのアプローチを大幅に進歩させ、ARGおよび他の関連遺伝子と遺伝要素の同時発生と共選択を通じてAMRの潜在的な普及に関する洞察を提供する[ref.49,50,51]。また、ARGとMGEの同時発生により、HGTなどの目的の遺伝的事象の証拠を追跡できる[ref.46]。

 ここでは、ナノポアシーケンスデータを使用して環境サンプル中のARGの包括的なプロファイリングを可能にするユーザーフレンドリーなオンラインプラットフォームであるNanoARGを紹介する。包括的なARGプロファイリングに加えて、NanoARGは、MRG、MGE、分類マーカー、既知の病原体との類似性が高い配列の識別、および同じDNA鎖上のこれらのさまざまな要素間の結合のインタラクティブな視覚化も提供する。環境ARGプロファイリングに対するNanoARGの可能性を実証するために、環境および臨床サンプルを含むいくつかのナノポアシーケンスライブラリを分析した。 Webサービスhttps://bench.cs.vt.edu/nanoargから無料で入手できる。ナノポアシーケンスデータをアップロードして処理するには、ユーザーログインとサブスクリプションが必要である。

 

 

使い方

https://bench.cs.vt.edu/nanoarg/#/にアクセスする。

f:id:kazumaxneo:20190907222330p:plain

 

Loginする。初回はユーザー登録が必要。

f:id:kazumaxneo:20190919010256p:plain

 

Loginしたらプロジェクトを作成する。

続いてview projectをクリック。少しわかりにくいが、この緑のAdd Sampleをクリックする。

f:id:kazumaxneo:20190919221617p:plain

 

nanoporeの1D、または2D / 1D2のリードのFASTAをアップロードする。

f:id:kazumaxneo:20190919221623p:plain

nextを押して次に進めていく。

 

しばらくたつと結果が表示される。

リード長分布

f:id:kazumaxneo:20190920134305p:plain

 

Antibiotic resistance hits distribution

f:id:kazumaxneo:20190920134329p:plain

 

Distribution of Metal Resistance Genes (MRGs)

f:id:kazumaxneo:20190920134422p:plain

 

f:id:kazumaxneo:20190920134738p:plain

 

 

f:id:kazumaxneo:20190920134812p:plain

 

f:id:kazumaxneo:20190920134816p:plain

 

1プロジェクトに複数のデータをアップして、比較解析できるようになっています。マニュアルを読んでみてください。
引用

NanoARG: a web service for detecting and contextualizing antimicrobial resistance genes from nanopore-derived metagenomes
Arango-Argoty GA, Dai D, Pruden A, Vikesland P, Heath LS, Zhang L

Microbiome. 2019 Jun 7;7(1):88

pacbioのbamをfastqに変換する BAM2fastx

 

 

PacificBiosciences/bam2fastx Converting and demultiplexing of PacBio BAM files into gzipped fasta and fastq files. by @PacificBiosciences - Repository | DevHub.io

 

BAM format specification for PacBio(5.1.0)

https://pacbiofileformats.readthedocs.io/en/5.1/BAM.html

 

インストール

依存

Github

#bioconda (link)
conda install -c bioconda bam2fastx

> bam2fasta -h

# bam2fasta -h

Usage: bam2fasta [options] INPUT

Converts multiple BAM and/or DataSet files into into gzipped FASTA file(s).

 

Options:

  -h,--help          Output this help.

  --version          Output version info.

  -o,--output        Prefix of output filenames

  -c                 Gzip compression level [1-9] [1]

  -u                 Do not compress. In this case, we will not add .gz, and we ignore any -c setting.

  --split-barcodes   Split output into multiple FASTA files, by barcode pairs.

  -p,--seqid-prefix  Prefix for sequence IDs in headers

 

Arguments:

  input              Input file.

 

> bam2fastq -h

# bam2fastq -h

Usage: bam2fastq [options] INPUT

Converts multiple BAM and/or DataSet files into into gzipped FASTQ file(s).

 

Options:

  -h,--help          Output this help.

  --version          Output version info.

  -o,--output        Prefix of output filenames

  -c                 Gzip compression level [1-9] [1]

  -u                 Do not compress. In this case, we will not add .gz, and we ignore any -c setting.

  --split-barcodes   Split output into multiple FASTQ files, by barcode pairs.

  -p,--seqid-prefix  Prefix for sequence IDs in headers

 

Arguments:

  input              Input file.

 

 

実行方法

subreads.bamを指定する。

bam2fasta -o projectName m54008_160330_053509.subreads.bam

 

引用

GitHub - PacificBiosciences/bam2fastx: Converting and demultiplexing of PacBio BAM files into gzipped fasta and fastq files.

 

参考

Where is the documentation for PacBio tool bam2fastx?

https://www.biostars.org/p/387385/

rRNAのアンプリコンシーケンスのトリミングを行う FIGARO

 

 マイクロバイオーム研究は、巨視的世界にとっての微生物コミュニティの重要性についての途方もない洞察を提供し続けている。ハイスループットDNAシーケンシング技術(すなわち、次世代シーケンス)は、微生物分類群を同定し、生物学的および環境試料の多様性および組成を計算することができるバイオインフォマティクスツールと組み合わせると、微生物集団の費用効果の高い迅速評価を可能にした。原核生物種および真核生物種をそれぞれ同定するために16Sおよび18S rRNA遺伝子配列がそれぞれ使用されるリボソームRNA遺伝子シーケンシングは、現在、マイクロバイオーム分析において使用されている最も広く使用されている技術の1つである。これらの配列の生物情報学的分析の前に、配列自体の予想されるエラーが最小化される一方で、トリミング後の配列情報が最大化されるようにトリミングパラメータを設定しなければならない。このアプリケーションノートでは、FIGARO:クオリティトリミングおよびフィルタリングした後のリード保持を最大化するように設計されたPythonベースのアプリケーションについて説明する。 FIGAROは、再現性を高め、DADA2ベースのパイプラインのトリミングパラメータ選択における試行錯誤を最小限に抑えるように設計されており、ペアエンドオーバーラップが必要な場合にも、トリミングパラメータの最適化および他のパイプラインのシーケンスエラーの最小化に役立つ。

  ユーザーは、アンプリコンの全長と、トリムされていないペアエンドデータリードを含むフォルダーをFASTQ形式で提供する必要がある。オプションの入力には、出力ファイル名指定、入力ディレクトリと出力ディレクトリ指定、FASTQファイルのサブサンプリングレート、ペアエンドリード間の最小オーバーラップ、およびテストされた各位置でのフィルタリングに使用する予想エラー率が含まれる。

 

 インストール

ubuntu16.04のpython3.6.3環境でテストした(docker使用、ホストOS ubuntu18.04)。

本体 Github

docker イメージをビルドする。

git clone https://github.com/Zymo-Research/figaro.git
cd figaro
docker build -t figaro .

> python3 figaro.py -h

# python3 figaro.py -h

usage: figaro.py [-h] [-o OUTPUTDIRECTORY] -a AMPLICONLENGTH -f

                 FORWARDPRIMERLENGTH -r REVERSEPRIMERLENGTH

                 [-i INPUTDIRECTORY] [-n OUTPUTFILENAME] [-m MINIMUMOVERLAP]

                 [-s SUBSAMPLE] [-p PERCENTILE]

 

optional arguments:

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

  -o OUTPUTDIRECTORY, --outputDirectory OUTPUTDIRECTORY

                        Directory for outputs

  -a AMPLICONLENGTH, --ampliconLength AMPLICONLENGTH

                        Length of amplicon (not including primers)

  -f FORWARDPRIMERLENGTH, --forwardPrimerLength FORWARDPRIMERLENGTH

                        Length of forward primer

  -r REVERSEPRIMERLENGTH, --reversePrimerLength REVERSEPRIMERLENGTH

                        Length of reverse primer

  -i INPUTDIRECTORY, --inputDirectory INPUTDIRECTORY

                        Directory with Fastq files to analyze

  -n OUTPUTFILENAME, --outputFileName OUTPUTFILENAME

                        Output file for trim site JSON

  -m MINIMUMOVERLAP, --minimumOverlap MINIMUMOVERLAP

                        Minimum overlap between the paired-end reads

  -s SUBSAMPLE, --subsample SUBSAMPLE

                        Subsampling level (will analyze approximately 1/x

                        reads

  -p PERCENTILE, --percentile PERCENTILE

                        Percentile to use for expected error model

 

 

実行方法

docker container run --rm -e AMPLICONLENGTH=[amplicon length] -e FORWARDPRIMERLENGTH=[forward primer length] \
-e REVERSEPRIMERLENGTH=[reverse primer length] -v /path/to/fastqs:/data/input \
-v /path/to/output:/data/output figaro

エラーになる。改善できたら追記します。

 

引用

FIGARO: An efficient and objective tool for optimizing microbiome rRNA gene trimming parameters
Michael M. Weinstein, Aishani Prem, Mingda Jin, Shuiquan Tang, Jeffrey M. Bhasin

bioRxiv preprint first posted online Apr. 16, 2019

 

ロングリードRNA seqのアライナー Graphmap2

 

 オックスフォードナノポアテクノロジー(ONT)[ref.1]やパシフィックバイオサイエンス(PacBio)[ref.2]などの企業が達成したシーケンシングテクノロジーの進歩により、長さが10 kbpを超えるロングリードが生成される。当初、このようなロングリードのエラー率は非常に高かったが、着実に低下している。PacBioデバイスの最後の世代は、イルミナのショートリードと同等のリードを生成する[ref.3 link]。さらに、ONTは最近、新しいR10ポアを発表し、この社内での結果は有望である。 RNA-seqの分野ではショートリードが依然として主に使用されているが、アイソフォームの検出と定量化に役立つ長いリードが必要である。アルゴリズム的に、RNA-seqリードを既知の転写産物にマッピングすることは、DNAリードをマッピングすることと同等である。しかし、これらのリードを真核生物のゲノムにマッピングすることは、RNAスプライシングのためにより複雑になる。第三世代のシーケンシング技術によって生成されたロングリード用に開発されたいくつかのRNA-seqスプライス対応マッピングツールがある。これらのツールの評価は[ref.4]および[ref.5]で行われ、GMAP [ref.6]およびMinimap2 [ref.7]が最高のパフォーマンスを発揮するツールである。ただし、特にエクソンのエッジを正しく配置し、転写産物内のすべてのエクソン、特に短いエクソンを見つけるには、まだ改善の余地がある。
 ここでは、Graphmapツールの拡張バージョンであるGraphmap2を紹介する。 Graphmapは、DNAリードをリファレンスゲノムにマッピングするための非常に感度の高いツールであり、RNAリードをマッピングする機能を拡張した。この拡張バージョンは、初期バージョン[ref.8]と同じ5段階の「リードファネル」アプローチを使用し、RNAリードのマッピングに固有のアップグレードを追加する。 Graphmap2のパフォーマンスと精度は、シミュレートされたデータセットと実際のデータセットで評価され、結果は最先端のRNAマッピングツールと比較された。さらに、ロングリードを使用して、潜在的に新しいアイソフォームと遺伝子の同定について分析した。
 図1は、RNAリードをリファレンスゲノムにアラインメントする完全なプロセスを示している。元のGraphmapマッピング方法は2つの段階に分割できる。最初の段階では、短いシードを使用してリファレンスゲノム上の候補位置を見つける。2番目の段階では、Edlibを使用してリードの正確なアライメントを計算する。これは、Myersのビットベクトルアルゴリズムの高速実装である[ref.9]。ステップを含む最初のステージ:(i)領域選択、(ii)グラフマッピング、(iii)k merサブストリング(LCSk)の最長共通サブシーケンス[ref.10]はGraphmap2に保持される。入力データセットのすべてのリードについて、このステージはリードの部分とリファレンスの部分との間の近似一致のセットを生成する。これらの一致はアンカーによって表される。すべてのアンカーは、リードの開始位置と終了位置、およびその一致のリファレンスの開始位置と終了位置で構成される。 Graphmap2は、次の段階でワークフローを拡張することにより続行する:(iv)アンカーフィルタリング(ナップサック)、(v)アンカーアライメント、(vi)アライメントチューニング、(vii)エクソンのグループ化と調整。アンカーフィルター処理ステップでは、Graphmap2はナップザックアルゴリズムのバリエーションを使用して最適なアンカーセットを見つけ、次にKSW2 [ref.11]アライナーを使用してこれらのアンカー間で区分的アフィンギャップアライメントを実行し、第1フェーズのリグメントを生成する。その後、これらのアライメントは、エクソン拡張およびエクソン境界調整ステップで処理され、第一段階のアライメントの品質が向上する。改善されたアライメントは、エクソンに分割され、エクソンのグループ化と調整フェーズでさらにグループ化および変更される。エキソンは、同じ転写産物から転写されたエキソンをグループ化するために、リファレンスゲノム上の位置に基づいてグループ化される。エクソングループは、「エクソンオフセット計算」およびエクソン調整ステップで分析および変更される。エクソンオフセット計算ステップは、エクソンのすべてのグループについて、同じグループ内のすべてのエクソンの正しい開始位置と終了位置を計算する。エクソン調整ステップでは、誤ったエクソンの開始位置または終了位置が、グループ全体の以前に計算された開始位置と終了位置に一致するように調整される。エクソン調整フェーズでは、Graphmap2ツールの出力である最終アライメントが生成される。

(以下略)

 


 

インストール

ubuntu18.0.4でテストした(docker使用、ホストOS macos10.14)。

git clone https://github.com/lbcb-sci/graphmap2 
cd graphmap2
make modules
make -j 8
cd /bin/Linux-x64/

> ./graphmap2 -h

# ./graphmap2 -h

ERROR: Unknown value of 'tool' parameter. Exiting.

 

Usage:

  ./graphmap2 tool

 

Options

    tool       STR   Specifies the tool to run:

                       align - the entire GraphMap pipeline.

                       owler - Overlapping With Long Erroneous Reads.

 

> ./graphmap2 align -h

# ./graphmap2 align -h

GraphMap - A very accurate and sensitive long-read, high error-rate sequence mapper

GraphMap Version: v0.6.1

Build date: Sep 16 2019 at 01:41:10

 

GraphMap (c) by Ivan Sovic, Mile Sikic and Niranjan Nagarajan

GraphMap is licensed under The MIT License.

 

Affiliations: Ivan Sovic (1, 3), Mile Sikic (2), Niranjan Nagarajan (3)

  (1) Ruder Boskovic Institute, Zagreb, Croatia

  (2) University of Zagreb, Faculty of Electrical Engineering and Computing

  (3) Genome Institute of Singapore, A*STAR, Singapore

 

 

Usage:

graphmap [options] -r <reference_file> -d <reads_file> -o <output_sam_path>

 

 

Usage:

  ./graphmap2 [options]

 

Options

  Input/Output options:

    -r, --ref                   STR   Path to the reference sequence (fastq or fasta).

    -i, --index                 STR   Path to the index of the reference sequence. If not specified, index is generated

                                      in the same folder as the reference file, with .gmidx extension. For

                                      non-parsimonious mode, secondary index .gmidxsec is also generated.

    -d, --reads                 STR   Path to the reads file.

    -o, --out                   STR   Path to the output file that will be generated.

        --gtf                   STR   Path to a General Transfer Format file. If specified, a transcriptome will be

                                      built from the reference sequence and used for mapping. Output SAM alignments will

                                      be in genome space (not transcriptome).

    -K, --in-fmt                STR   Format in which to input reads. Options are:

                                       auto  - Determines the format automatically from file extension.

                                       fastq - Loads FASTQ or FASTA files.

                                       fasta - Loads FASTQ or FASTA files.

                                       gfa   - Graphical Fragment Assembly format.

                                       sam   - Sequence Alignment/Mapping format. [auto]

    -L, --out-fmt               STR   Format in which to output results. Options are:

                                       sam  - Standard SAM output (in normal and '-w overlap' modes).

                                       m5   - BLASR M5 format. [sam]

    -I, --index-only             -    Build only the index from the given reference and exit. If not specified, index

                                      will automatically be built if it does not exist, or loaded from file otherwise. [false]

        --rebuild-index          -    Always rebuild index even if it already exists in given path. [false]

        --auto-rebuild-index     -    Rebuild index only if an existing index is of an older version or corrupt. [false]

    -u, --ordered                -    SAM alignments will be output after the processing has finished, in the order of

                                      input reads. [false]

    -B, --batch-mb              INT   Reads will be loaded in batches of the size specified in megabytes. Value <= 0

                                      loads the entire file. [1024]

 

  General-purpose pre-set options:

    -x, --preset                STR   Pre-set parameters to increase sensitivity for different sequencing technologies.

                                      Valid options are:

                                       illumina  - Equivalent to: '-a gotoh -w normal -M 5 -X 4 -G 8 -E 6'

                                       overlap   - Equivalent to: '-a anchor -w normal --overlapper --evalue 1e0

                                      --ambiguity 0.50 --secondary'

                                       sensitive - Equivalent to: '--freq-percentile 1.0 --minimizer-window 1'

 

  Alignment options:

    -a, --alg                   STR   Specifies which algorithm should be used for alignment. Options are:

                                       sg       - Myers' bit-vector approach. Semiglobal. Edit dist. alignment.

                                       sggotoh       - Gotoh alignment with affine gaps. Semiglobal.

                                       anchor      - anchored alignment with end-to-end extension.

                                                     Uses Myers' global alignment to align between anchors.

                                       anchorgotoh - anchored alignment with Gotoh.

                                                     Uses Gotoh global alignment to align between anchors. [anchor]

    -w, --approach              STR   Additional alignment approaches. Changes the way alignment algorithm is applied.

                                      Options are:

                                       normal         - Normal alignment of reads to the reference.

                                       (Currently no other options are provided. This is a placeholder for future

                                      features, such as cDNA mapping) [normal]

        --overlapper             -    Perform overlapping instead of mapping. Skips self-hits if reads and reference

                                      files contain same sequences, and outputs lenient secondary alignments. [false]

        --no-self-hits           -    Similar to overlapper, but skips mapping of sequences with same headers. Same

                                      sequences can be located on different paths, and their overlap still skipped. [false]

    -M, --match                 INT   Match score for the DP alignment. Ignored for Myers alignment. [5]

    -X, --mismatch              INT   Mismatch penalty for the DP alignment. Ignored for Myers alignment. [4]

    -G, --gapopen               INT   Gap open penalty for the DP alignment. Ignored for Myers alignment. [8]

    -E, --gapext                INT   Gap extend penalty for the DP alignment. Ignored for Myers alignment. [6]

    -z, --evalue                FLT   Threshold for E-value. If E-value > FLT, read will be called unmapped. If FLT <

                                      0.0, thredhold not applied. [1e0]

    -c, --mapq                  INT   Threshold for mapping quality. If mapq < INT, read will be called unmapped. [1]

        --extcigar               -    Use the extended CIGAR format for output alignments. [false]

        --spliced                -    Align clusters of anchors independently and report them as separate alignments.

                                      Does not align in-between clusters. Works only for anchored alignment modes. [false]

        --split                  -    Align clusters of anchors independently and report them as separate alignments.

                                      Does not align in-between clusters. Works only for anchored alignment modes. [false]

        --no-end2end             -    Disables extending of the alignments to the ends of the read. Works only for

                                      anchored modes. [false]

        --max-error             FLT   If an alignment has error rate (X+I+D) larger than this, it won't be taken into

                                      account. If >= 1.0, this filter is disabled. [1.0]

        --max-indel-error       FLT   If an alignment has indel error rate (I+D) larger than this, it won't be taken

                                      into account. If >= 1.0, this filter is disabled. [1.0]

 

  Algorithmic options:

    -k                          INT   Graph construction kmer size. [6]

    -l                          INT   Number of edges per vertex. [9]

    -A, --minbases              INT   Minimum number of match bases in an anchor. [12]

    -e, --error-rate            FLT   Approximate error rate of the input read sequences. [0.45]

    -g, --max-regions           INT   If the final number of regions exceeds this amount, the read will be called

                                      unmapped. If 0, value will be dynamically determined. If < 0, no limit is set. [0]

    -q, --reg-reduce            INT   Attempt to heuristically reduce the number of regions if it exceeds this amount.

                                      Value <= 0 disables reduction but only if param -g is not 0. If -g is 0, the value

                                      of this parameter is set to 1/5 of maximum number of regions. [0]

    -C, --circular               -    Reference sequence is a circular genome. [false]

    -F, --ambiguity             FLT   All mapping positions within the given fraction of the top score will be counted

                                      for ambiguity (mapping quality). Value of 0.0 counts only identical mappings. [0.02]

    -Z, --secondary              -    If specified, all (secondary) alignments within (-F FLT) will be output to a

                                      file. Otherwise, only one alignment will be output. [false]

    -P, --double-index           -    If false, only one gapped spaced index will be used in region selection. If true,

                                      two such indexes (with different shapes) will be used (2x memory-hungry but more

                                      powerful for very high error rates). [false]

        --min-bin-perc          FLT   Consider only bins with counts above FLT * max_bin, where max_bin is the count of

                                      the top scoring bin. [0.75]

        --bin-step              FLT   After a chunk of bins with values above FLT * max_bin is processed, check if

                                      there is one extremely dominant region, and stop the search. [0.25]

        --min-read-len          INT   If a read is shorter than this, it will be marked as unmapped. This value can be

                                      lowered if the reads are known to be accurate. [80]

        --minimizer-window      INT   Length of the window to select a minimizer from. If equal to 1, minimizers will

                                      be turned off. [5]

        --freq-percentile       FLT   Filer the (1.0 - value) percent of most frequent seeds in the lookup process. [0.99]

        --fly-index              -    Index will be constructed on the fly, without storing it to disk. If it already

                                      exists on disk, it will be loaded unless --rebuild-index is specified. [false]

        --chain-indel-bw        FLT   Diagonal between two anchors needs to be within this band to chain them. [0.23]

 

  Anchor chaining options:

        --chain-max-dist        INT   Maximum distance between two anchors to chain them. [200]

        --chain-min-cov         INT   Minimum number of bases covered by a chain to keep it. Only checked if number of

                                      anchors in the cluster is <= '--chain-min-num-anchors'. [50]

        --chain-min-num-anchors INT   If number of anchors in chain is <= to this value and the chain covers <

                                      '--chain-min-cov' bases, the chain is discarded. [2]

 

  Other options:

    -t, --threads               INT   Number of threads to use. If '-1', number of threads will be equal to min(24, num_cores/2). [-1]

    -v, --verbose               INT   Verbose level. If equal to 0 nothing except strict output will be placed on stdout. [5]

    -s, --start                 INT   Ordinal number of the read from which to start processing data. [0]

    -n, --numreads              INT   Number of reads to process per batch. Value of '-1' processes all reads. [-1]

    -h, --help                   -    View this help. [false]

 

  Debug options:

    -y, --debug-read            INT   ID of the read to give the detailed verbose output. [-1]

    -Y, --debug-qname           STR   QNAME of the read to give the detailed verbose output. Has precedence over -y.

                                      Use quotes to specify.

    -b, --verbose-sam           INT   Helpful debug comments can be placed in SAM output lines (at the end). Comments

                                      can be turned off by setting this parameter to 0. Different values

                                      increase/decrease verbosity level.

                                      0 - verbose off

                                      1 - server mode, command line will be omitted to obfuscate paths.

                                      2 - umm this one was skipped by accident. The same as 0.

                                      >=3 - detailed verbose is added for each alignment, including timing measurements

                                      and other.

                                      4 - qnames and rnames will not be trimmed to the first space.

                                      5 - QVs will be omitted (if available). [0]

 

 

./graphmap2 owler -h

# ./graphmap2 owler -h

GraphMap - A very accurate and sensitive long-read, high error-rate sequence mapper

GraphMap Version: v0.6.1

Build date: Sep 16 2019 at 01:41:10

 

GraphMap (c) by Ivan Sovic, Mile Sikic and Niranjan Nagarajan

GraphMap is licensed under The MIT License.

 

Affiliations: Ivan Sovic (1, 3), Mile Sikic (2), Niranjan Nagarajan (3)

  (1) Ruder Boskovic Institute, Zagreb, Croatia

  (2) University of Zagreb, Faculty of Electrical Engineering and Computing

  (3) Genome Institute of Singapore, A*STAR, Singapore

 

 

Usage:

graphmap [options] -r <reference_file> -d <reads_file> -o <output_sam_path>

 

 

Usage:

  ./graphmap2 [options]

 

Options

  Input/Output options:

    -r, --ref               STR   Path to the reference sequence (fastq or fasta).

    -d, --reads             STR   Path to the reads file.

    -o, --out               STR   Path to the output file that will be generated.

    -K, --in-fmt            STR   Format in which to input reads. Options are:

                                   auto  - Determines the format automatically from file extension.

                                   fastq - Loads FASTQ or FASTA files.

                                   fasta - Loads FASTQ or FASTA files.

                                   gfa   - Graphical Fragment Assembly format.

                                   sam   - Sequence Alignment/Mapping format. [auto]

    -L, --out-fmt           STR   Format in which to output results. Options are:

                                   mhap - MHAP overlap format.

                                   paf  - PAF (Minimap) overlap format. [mhap]

    -I, --index-only         -    Build only the index from the given reference and exit. If not specified, index will

                                  automatically be built if it does not exist, or loaded from file otherwise. [false]

        --load               -    Loads the generated index from a file specified with the --index parameter. [false]

        --store              -    Stores the generated index to a file specified with the --index parameter. If the

                                  index was loaded from a file, it will not be rewritten. [false]

        --rebuild-index      -    If the index file is not compatible with currently set parameters, automatically

                                  rebuild the index. Otherwise, use the disk index. This parameter only works if --load

                                  is specified. [false]

    -i, --index             STR   Path to the index of the reference sequence. If not specified, path will be set to

                                  the same folder as the reference file, with .owlidx extension.

 

  Algorithmic options:

    -e, --error-rate        FLT   Approximate error rate of the input read sequences. [0.45]

        --shape             STR   Seed shape which will be used for indexing and querying the sequences. [111111111111111]

        --minimizer-window  INT   Length of the window to select a minimizer from. If equal to 1, minimizers will be

                                  turned off. [5]

        --freq-percentile   FLT   Filer the (1.0 - value) percent of most frequent seeds in the lookup process. [0.999]

 

  Overlap filtering:

        --min-read-len      INT   Overlaps between sequences, where one sequence is shorter than min-read-len, will be discarded. [500]

        --min-cov-bases     INT   Minimum number of bases covered by seeds in an overlap. [50]

        --min-overlap-len   INT   Minimum span of overlap region on either of the two sequences. [100]

        --overhang          FLT   Value in range [0.0, 1.0] which specifies an allowed distance of an overlap from the

                                  beginning/ending of a sequence (e.g. (overhang * read_len)). [ 0.20]

        --max-overhang      INT   Allowed overhang length will be upper-bound by this parameter. Essentially, overhangs

                                  of min(max_overhang, overhang * read_len) will be allowed unless max_overhang is < 0,

                                  in which case this threshold is turned off. [1000]

        --min-perc-coverage FLT   Discard overlaps which have a fraction of bases covered by seeds then specified by

                                  this parameter. [ 0.01]

        --min-seeds         INT   Discard overlaps with less seed hits (after filtering) than the number specified by

                                  this parameter. [4]

 

  Other options:

    -t, --threads           INT   Number of threads to use. If '-1', number of threads will be equal to min(24, num_cores/2). [-1]

    -v, --verbose           INT   Verbose level. If equal to 0 nothing except strict output will be placed on stdout. [5]

    -s, --start             INT   Ordinal number of the read from which to start processing data. [0]

    -n, --numreads          INT   Number of reads to process per batch. Value of '-1' processes all reads. [-1]

    -h, --help               -    View this help. [false]

 

  General-purpose pre-set options:

    -x, --preset            STR   Pre-set composite parameters. Valid options are:

                                   disk       - Equivalent to: '--load --store'

                                   no-filters - Disables all filters except minimum number of seeds and

                                                overhang thresholds. Equivalent to:

                                             '--min-read-len 0 --min-cov-bases 0 --min-overlap-len 0 --min-perc-coverage 0.0'

 

  Debug options:

    -b, --verbose-out       INT   Helpful debug comments can be placed in output lines (at the end). Comments can be

                                  turned off by setting this parameter to 0. [0]

    -y, --debug-read        INT   ID of the read to give the detailed verbose output. [-1]

    -Y, --debug-qname       STR   QNAME of the read to give the detailed verbose output. Has precedence over -y. Use

                                  quotes to specify.

 

 

 

 

実行方法

アラインメント

graphmap2 align -r reference.fa -d reads.fasta -t 8 -o output.sam 
  • -t    Number of threads to use. If '-1', number of threads will be equal to min(24, num_cores/2). [-1]

 

ロングリードのオーバーラップ

graphmap2 owler -r reads.fasta -d reads.fasta -o output.mhap 

 

引用

Graphmap2 - splice-aware RNA-seq mapper for long reads

Josip Marić​​, Ivan Sović​​, Krešimir Križanović, Niranjan Nagarajan, Mile Šikić

bioRxiv preprint first posted online Jul. 30, 2019

 

関連


バクテリア/アーキアのゲノム距離を計算するwebツール GGDC

  

 DNA-DNAハイブリダイゼーション(DDH)は、古細菌および細菌種の描写のための分類学的ゴールドスタンダードとして現在も使用されているウェットラボ法である。 2つのそれぞれの生物のゲノムDNAがDDHの類似性が70%未満であることが明らかになった場合、これはそれらを異なる種と見なすための主な論点であり、逆もまた同様である。 DDHは、退屈で面倒で、エラーが発生しやすいと広く考えられている[ref.3、4]。さらに、ゲノムシーケンスとは対照的に、DDH値自体より多くの情報を返さないため、結果として、データを再利用して段階的に作業することはできない。

 DDH技術は現在、少数の専門ラボ(主に微生物サービスコレクション)でのみ確立されており、実験的偏差を起こしやすいため、その実験の統計的信頼性を決定するためにいくつかの実験を繰り返す必要がある。たとえば、微生物学における種の境界に関して、関連するクエスチョンは、DDH値が70%を大幅に下回るか上回るかである。これは、表現型の測定などの他の基準に対してDDHからの証拠をトレードオフする必要がある多相アプローチのコンテキストでは特に重要である[ref.5]。 16S rRNA配列の類似性が特定のしきい値を下回っている場合にのみ、DDH実験を新種の説明で省略できる。これは、DDH値が70%を超えることは期待できないことを示す[ref.2]。

(一段落省略)

技術的な問題と進歩を考慮すると、ウェットラボのDDH手順とDDHのデジタル推定の関係は、16S rRNA配列のDNA:rRNAクロスハイブリダイゼーション融解曲線[ref.9、10]が置き換えられた約30年前に起こったことを思い起こさせる。これは、微生物系統の大幅な進歩をサポートした[ref.11]。

 Genome Blast Distance Phylogenyアプローチ(GBDP)は、完全に(または不完全な)シーケンスされたゲノムの特定のセットから系統樹またはネットワークを推論するためのアプローチとして当初考案され[ref.12]、その後再検討および強化された[ref.8、13 -16]。基本的な原理は次のとおりである。最初のステップでは、BLAST [ref.17]などのツールを使用して2つのゲノムAとBをローカルにアラインし、スコアの高いセグメントペア(HSP;これらはゲノム間一致)のセットを生成する。 2番目のステップでは、特定の距離式を使用して、これらのHSPに含まれる情報(たとえば、同一の塩基対の総数)を単一のゲノム間距離値に変換する。系統樹は、neighbour joiningなどの標準的な手法を使用して、そのような距離行列から推測できる[ref.18]。これらの方法は、相当量のパラロガス遺伝子、大きなリピート、減少したゲノム[ref.12]、および配列内の複雑さの低い領域[ref.16]が存在する場合でも堅牢である。 GBDPはプロテオームデータにも適用でき[ref.13]、単一遺伝子にも適用できる[ref.19]。

 GBDPのさらなる使用が最近評価された。つまり、DDH値のデジタルでの同等な物を推測するためである[ref.8、16]。これらはウェットラボハイブリダイゼーションの結果をうまく模倣し、先行するゲノムシーケンスに基づく方法[ref.6]よりも経験的なDDH値のセットとの高い相関を提供し、かなり不完全なゲノム[ref.8]に対処できることが判明した。微生物学者は、http://ggdc.dsmz.de [ref.16]にある無料のWebサービスを使用してGBDPを使用して、ゲノムペアを送信し、DDHアナログおよびモデルベースのDDH推定値を受信できる。

(以下略)

 

GGDC Scientific Background

http://ggdc.dsmz.de/ggdc_background.php#

FAQ

https://ggdc.dsmz.de/faq.php#tabGGDC

 

使い方

http://ggdc.dsmz.de/ggdc.php#  にアクセスする。

f:id:kazumaxneo:20190913024250p:plain


アラインメントのプログラムを選択する。標準ではBLAST+になっている。

f:id:kazumaxneo:20190913024320p:plain

 

クエリのゲノムと比較するリファレンスゲノムをそれぞれ選択する。ゲノムのFASTAをアップロードするか、Genebank accesion IDsを指定する。

クエリのゲノム

f:id:kazumaxneo:20190914014016p:plain

 

リファレンスゲノム。ローカルのファイルを複数指定する場合は、shiftキーやCtrlキーで複数同時選択してアップロードする。acceson IDを指定する場合は、1行に1生物ずつタイプしていく。

f:id:kazumaxneo:20190914014057p:plain

 

メールアドレスを指定してSubmitする。

f:id:kazumaxneo:20190914014637p:plain

 

DDHの結果はメール中に直接記載される。またCSVの添付ファイルにまとめられる。 

 

引用

Genome sequence-based species delimitation with confidence intervals and improved distance functions
Jan P Meier-Kolthoff, Alexander F Auch, Hans-Peter Klenk, Markus Göker
BMC Bioinformatic svolume 14, Article number: 60 (2013)