macでインフォマティクス

macでインフォマティクス

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

2005-2017年の各大学のバイオイオンフォマティクス系論文出版数と内容を視覚化した BIOLITMAP

 

 科学および技術のさまざまな分野の科学的貢献に価値を置くことがますます重要になっている。バイオインフォマティクスの急速な発展とその本質的な学際的性質のために、さまざまな応用分野や機関の貢献がどのように発展しているかを視覚化することは特に難しい。組織、助成機関、その他の社会関係者による分析と計画を容易にするために、視覚化がますます重要になっている。

 バイオインフォマティクスのさまざまな研究分野におけるpublicationsの探査を容易にするために、本著者らはBIOLITMAPを発表する。 BIOLITMAPは、年ごとの探査、ジャーナルおよび特定のバイオインフォマティクスのトピックを可能にする、機関によって発行された世界的なバイオインフォマティクス研究をキャプチャする。 BIOLITMAPを使用して、バイオインフォマティクスにおけるさまざまな研究ライン、機関、国の進化のグローバルな概要を提供している。これはこれまで利用できなかった新しい機能である。さらに、それは、期間、ジャーナル、トピックおよび機関に基づいて科学出版物(scientific publications)を見つけ、それらを整頓、視覚化し、検索結果をプリントして見るためにユーザが情報を容易にナビゲートするインタラクティブな方法を提供する。関連する統計このツールを使用して、他の以前のツールでは不可能だった方法で科学論文を調査および取得するための新しいベースラインを提供し、それをさらに分析する可能性を開く。

 このシステムはScopus(https://www.scopus.com/)(wiki)データベースから得られたデータに基づいている。信頼できる引用元として認められているためである(Falagas et al、2008)。 このシステムの最初のバージョンでは、Bioinformatics OUP、BMC Bioinformatics、BMC Genomics、PLoS Computational Biology、Nucleic Acids Researchなど、5つの代表的なジャーナルに掲載された論文が掲載されている。 Scopusには、「ジャーナル」と「レビュー」という2つのカテゴリの文書がある。 Googleの検索では、「Reviews」を除く「Article」タイプのドキュメントのみを保持するようにフィルタを適用した。 その結果、Oxford BioinformaticsのApplication Notesと、ジャーナルの special issues が含まれる。 現在のデータベースは2005年から2017年までをカバーしており、合計46,552のarticlesが含まれている。

 

ソースコード

Github


使い方

http://socialanalytics.bsc.es/biolitmap/にアクセスする。

f:id:kazumaxneo:20190320221414p:plain

 

例えばアイビーリーグあたりを拡大して見てみる。ハーバード大学のscientific publicationsを可視化。MITと隣接しているので、十分に拡大しないとplotが被ってしまう。

f:id:kazumaxneo:20190321101258p:plain

 

右上から年代を選択。Select Allで全て表示した。

f:id:kazumaxneo:20190321101604p:plain

2005年-2017年の推移。

 

日本を見てみる。Bioinformatics、BMC Bioinformatics、BMC Genomics、PLoS Computational Biology、Nucleic Acids Researchのうち、BioinformaticsとBMC Bioinformaticsだけ選択してみた。調べてみると東大が一番多かった。

f:id:kazumaxneo:20190321101845p:plain

日本のplotが非常に少ないが、数が一定数ないとカウント対象にならないと思われる。また、大学以外の研究機関はごく一部以外集計対象になっていない。

 

view articlesをクリックするとカウント対象になった論文一覧を表示/プリントできる。

f:id:kazumaxneo:20190321102147p:plain

 

 

 

右上からは、年代、5 journal、研究カテゴリーを選択できる。

f:id:kazumaxneo:20190321071856p:plain

上はカテゴリー選択時の画面

 

左端のgraphicボタンをクリックすると、地図に重ねてグラフを表示できる。

ジャーナル内訳

f:id:kazumaxneo:20190321072306p:plain

 

年毎のscientific publications件数

f:id:kazumaxneo:20190321102633p:plain

 

カテゴリー内訳

f:id:kazumaxneo:20190321072312p:plain

引用

BIOLITMAP: a web-based geolocated, temporal and thematic visualization of the evolution of bioinformatics publications
Adrián Bazaga Alfonso Valencia María- JoséRementeria
Bioinformatics, Published: 05 December 2018

 

Pacbioのpolishingツール Quiver / ArrowとバリアントコーラーPlurality

 

 Quiverは、Pacbioがテンプレートリードを前提として、最大準尤度テンプレートシーケンスを見つける、より洗練されたアルゴリズムである。 PacBioのリードは、テンプレートシーケンスを指定してリードの準尤度をスコア付けする条件付きランダムフィールドアプローチを使用してモデル化される。各リードの基本シーケンスに加えて、Quiverはbasecaller元が提供するいくつかの追加のQV共変量を使用する。これらの共変量を使用すると、各リードに関する追加情報が提供され、より正確なコンセンサスコールが可能になる。Quiverは、マッパーによって提供されるアライメント(通常はBLASR)を使用しない。ただし、マクロなレベルでリードをまとめどのようにグループ化するか決定する場合を除く。暗黙的に独自のリアライメントを実行するので、indelを含むすべてのバリアント型に非常に敏感である。

 Arrowは近い将来Quiverに取って代わることを意図した新しいモデルである。 Quiverとの主な違いは、CRFの代わりにHMMモデルを使用し、本当の可能性を計算し、さらに少数の共変量を使用することである。 Arrowに関するホワイトペーパーがもうすぐ利用可能になる予定である。

 Pluralityはシンプルなバリアンコールアルゴリズムである。それはアライメントされたリード(BLASRによって生成、または代替マッピングツールで生成されたもの)をpileupし、各ポジションで、リファレンスベースで最も豊富な(つまり複数の)塩基をコンセンサスとしてコールする。

 

 

矢筒(Quiver)と矢(Arrow)のセットです。

  • How to install and use GenomicConsensus 

GenomicConsensus/HowTo.rst at develop · PacificBiosciences/GenomicConsensus · GitHub

  • Quiver FAQ

https://github.com/PacificBiosciences/GenomicConsensus/blob/develop/doc/FAQ.rst

  • Quiver is described in detail in the supplementary material to the HGAP paper(ref.1).

関連ツイート

 

インストール

ubuntu16.04のPython 3.6.8環境でテストした(ホストOS macos10.14)。

 

Github

GenomicConsensus

PacBio Secondary Analysis Tools on Bioconda(Bicondaに対応したPacbioのオフィシャルツール)

#condaで導入
conda install -c bioconda genomicconsensus

> quiver -h

$ quiver -h

usage: variantCaller [-h] [--version] [--emit-tool-contract]

                     [--resolved-tool-contract RESOLVED_TOOL_CONTRACT]

                     [--log-file LOG_FILE]

                     [--log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL} | --debug | --quiet | -v]

                     --referenceFilename REFERENCEFILENAME -o OUTPUTFILENAMES

                     [-j NUMWORKERS] [--minConfidence MINCONFIDENCE]

                     [--minCoverage MINCOVERAGE]

                     [--noEvidenceConsensusCall {nocall,reference,lowercasereference}]

                     [--coverage COVERAGE] [--minMapQV MINMAPQV]

                     [--referenceWindow REFERENCEWINDOWSASSTRING]

                     [--alignmentSetRefWindows]

                     [--referenceWindowsFile REFERENCEWINDOWSASSTRING]

                     [--barcode _BARCODE] [--readStratum READSTRATUM]

                     [--minReadScore MINREADSCORE] [--minSnr MINHQREGIONSNR]

                     [--minZScore MINZSCORE] [--minAccuracy MINACCURACY]

                     [--algorithm {quiver,arrow,plurality,poa,best}]

                     [--parametersFile PARAMETERSFILE]

                     [--parametersSpec PARAMETERSSPEC]

                     [--maskRadius MASKRADIUS] [--maskErrorRate MASKERRORRATE]

                     [--pdb] [--notrace] [--pdbAtStartup] [--profile]

                     [--annotateGFF] [--reportEffectiveCoverage] [--diploid]

                     [--queueSize QUEUESIZE] [--threaded]

                     [--referenceChunkSize REFERENCECHUNKSIZE]

                     [--fancyChunking] [--simpleChunking]

                     [--referenceChunkOverlap REFERENCECHUNKOVERLAP]

                     [--autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE]

                     [--aligner {affine,simple}] [--refineDinucleotideRepeats]

                     [--noRefineDinucleotideRepeats] [--fast]

                     [--skipUnrecognizedContigs]

                     inputFilename

 

Compute genomic consensus and call variants relative to the reference.

 

optional arguments:

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

  --version             show program's version number and exit

  --emit-tool-contract  Emit Tool Contract to stdout (default: False)

  --resolved-tool-contract RESOLVED_TOOL_CONTRACT

                        Run Tool directly from a PacBio Resolved tool contract

                        (default: None)

  --log-file LOG_FILE   Write the log to file. Default(None) will write to

                        stdout. (default: None)

  --log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}

                        Set log level (default: WARN)

  --debug               Alias for setting log level to DEBUG (default: False)

  --quiet               Alias for setting log level to CRITICAL to suppress

                        output. (default: False)

  -v, --verbose         Set the verbosity level. (default: None)

 

Basic required options:

  inputFilename         The input cmp.h5 or BAM alignment file

  --referenceFilename REFERENCEFILENAME, --reference REFERENCEFILENAME, -r REFERENCEFILENAME

                        The filename of the reference FASTA file (default:

                        None)

  -o OUTPUTFILENAMES, --outputFilename OUTPUTFILENAMES

                        The output filename(s), as a comma-separated

                        list.Valid output formats are .fa/.fasta, .fq/.fastq,

                        .gff, .vcf (default: )

 

Parallelism:

  -j NUMWORKERS, --numWorkers NUMWORKERS

                        The number of worker processes to be used (default: 1)

 

Output filtering:

  --minConfidence MINCONFIDENCE, -q MINCONFIDENCE

                        The minimum confidence for a variant call to be output

                        to variants.{gff,vcf} (default: 40)

  --minCoverage MINCOVERAGE, -x MINCOVERAGE

                        The minimum site coverage that must be achieved for

                        variant calls and consensus to be calculated for a

                        site. (default: 5)

  --noEvidenceConsensusCall {nocall,reference,lowercasereference}

                        The consensus base that will be output for sites with

                        no effective coverage. (default: lowercasereference)

 

Read selection/filtering:

  --coverage COVERAGE, -X COVERAGE

                        A designation of the maximum coverage level to be used

                        for analysis. Exact interpretation is algorithm-

                        specific. (default: 100)

  --minMapQV MINMAPQV, -m MINMAPQV

                        The minimum MapQV for reads that will be used for

                        analysis. (default: 10)

  --referenceWindow REFERENCEWINDOWSASSTRING, --referenceWindows REFERENCEWINDOWSASSTRING, -w REFERENCEWINDOWSASSTRING

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd (default: entire reference).

                        (default: None)

  --alignmentSetRefWindows

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd will be pulled from the alignment

                        file. (default: False)

  --referenceWindowsFile REFERENCEWINDOWSASSTRING, -W REFERENCEWINDOWSASSTRING

                        A file containing reference window designations, one

                        per line (default: None)

  --barcode _BARCODE    Only process reads with the given barcode name.

                        (default: None)

  --readStratum READSTRATUM

                        A string of the form 'n/N', where n, and N are

                        integers, 0 <= n < N, designating that the reads are

                        to be deterministically split into N strata of roughly

                        even size, and stratum n is to be used for variant and

                        consensus calling. This is mostly useful for Quiver

                        development. (default: None)

  --minReadScore MINREADSCORE

                        The minimum ReadScore for reads that will be used for

                        analysis (arrow-only). (default: 0.65)

  --minSnr MINHQREGIONSNR

                        The minimum acceptable signal-to-noise over all

                        channels for reads that will be used for analysis

                        (arrow-only). (default: 2.5)

  --minZScore MINZSCORE

                        The minimum acceptable z-score for reads that will be

                        used for analysis (arrow-only). (default: -3.5)

  --minAccuracy MINACCURACY

                        The minimum acceptable window-global alignment

                        accuracy for reads that will be used for the analysis

                        (arrow-only). (default: 0.82)

 

Algorithm and parameter settings:

  --algorithm {quiver,arrow,plurality,poa,best}

  --parametersFile PARAMETERSFILE, -P PARAMETERSFILE

                        Parameter set filename (such as ArrowParameters.json

                        or QuiverParameters.ini), or directory D such that

                        either D/*/GenomicConsensus/QuiverParameters.ini, or

                        D/GenomicConsensus/QuiverParameters.ini, is found. In

                        the former case, the lexically largest path is chosen.

                        (default: None)

  --parametersSpec PARAMETERSSPEC, -p PARAMETERSSPEC

                        Name of parameter set (chemistry.model) to select from

                        the parameters file, or just the name of the

                        chemistry, in which case the best available model is

                        chosen. Default is 'auto', which selects the best

                        parameter set from the alignment data (default: auto)

  --maskRadius MASKRADIUS

                        Radius of window to use when excluding local regions

                        for exceeding maskMinErrorRate, where 0 disables any

                        filtering (arrow-only). (default: 3)

  --maskErrorRate MASKERRORRATE

                        Maximum local error rate before the local region

                        defined by maskRadius is excluded from polishing

                        (arrow-only). (default: 0.7)

 

Verbosity and debugging/profiling:

  --pdb                 Enable Python debugger (default: False)

  --notrace             Suppress stacktrace for exceptions (to simplify

                        testing) (default: False)

  --pdbAtStartup        Drop into Python debugger at startup (requires ipdb)

                        (default: False)

  --profile             Enable Python-level profiling (using cProfile).

                        (default: False)

  --annotateGFF         Augment GFF variant records with additional

                        information (default: False)

  --reportEffectiveCoverage

                        Additionally record the *post-filtering* coverage at

                        variant sites (default: False)

 

Advanced configuration options:

  --diploid             Enable detection of heterozygous variants

                        (experimental) (default: False)

  --queueSize QUEUESIZE, -Q QUEUESIZE

  --threaded, -T        Run threads instead of processes (for debugging

                        purposes only) (default: False)

  --referenceChunkSize REFERENCECHUNKSIZE, -C REFERENCECHUNKSIZE

  --fancyChunking       Adaptive reference chunking designed to handle

                        coverage cutouts better (default: True)

  --simpleChunking      Disable adaptive reference chunking (default: True)

  --referenceChunkOverlap REFERENCECHUNKOVERLAP

  --autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE

                        Disable the HDF5 chunk cache when the number of

                        datasets in the cmp.h5 exceeds the given threshold

                        (default: 500)

  --aligner {affine,simple}, -a {affine,simple}

                        The pairwise alignment algorithm that will be used to

                        produce variant calls from the consensus (Quiver

                        only). (default: affine)

  --refineDinucleotideRepeats

                        Require quiver maximum likelihood search to try one

                        less/more repeat copy in dinucleotide repeats, which

                        seem to be the most frequent cause of suboptimal

                        convergence (getting trapped in local optimum) (Quiver

                        only) (default: True)

  --noRefineDinucleotideRepeats

                        Disable dinucleotide refinement (default: True)

  --fast                Cut some corners to run faster. Unsupported! (default:

                        False)

  --skipUnrecognizedContigs

                        Do not abort when told to process a reference window

                        (via -w/--referenceWindow[s]) that has no aligned

                        coverage. Outputs emptyish files if there are no

                        remaining non-degenerate windows. Only intended for

                        use by smrtpipe scatter/gather. (default: False)

> arrow -h

$ arrow -h

usage: variantCaller [-h] [--version] [--emit-tool-contract]

                     [--resolved-tool-contract RESOLVED_TOOL_CONTRACT]

                     [--log-file LOG_FILE]

                     [--log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL} | --debug | --quiet | -v]

                     --referenceFilename REFERENCEFILENAME -o OUTPUTFILENAMES

                     [-j NUMWORKERS] [--minConfidence MINCONFIDENCE]

                     [--minCoverage MINCOVERAGE]

                     [--noEvidenceConsensusCall {nocall,reference,lowercasereference}]

                     [--coverage COVERAGE] [--minMapQV MINMAPQV]

                     [--referenceWindow REFERENCEWINDOWSASSTRING]

                     [--alignmentSetRefWindows]

                     [--referenceWindowsFile REFERENCEWINDOWSASSTRING]

                     [--barcode _BARCODE] [--readStratum READSTRATUM]

                     [--minReadScore MINREADSCORE] [--minSnr MINHQREGIONSNR]

                     [--minZScore MINZSCORE] [--minAccuracy MINACCURACY]

                     [--algorithm {quiver,arrow,plurality,poa,best}]

                     [--parametersFile PARAMETERSFILE]

                     [--parametersSpec PARAMETERSSPEC]

                     [--maskRadius MASKRADIUS] [--maskErrorRate MASKERRORRATE]

                     [--pdb] [--notrace] [--pdbAtStartup] [--profile]

                     [--annotateGFF] [--reportEffectiveCoverage] [--diploid]

                     [--queueSize QUEUESIZE] [--threaded]

                     [--referenceChunkSize REFERENCECHUNKSIZE]

                     [--fancyChunking] [--simpleChunking]

                     [--referenceChunkOverlap REFERENCECHUNKOVERLAP]

                     [--autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE]

                     [--aligner {affine,simple}] [--refineDinucleotideRepeats]

                     [--noRefineDinucleotideRepeats] [--fast]

                     [--skipUnrecognizedContigs]

                     inputFilename

 

Compute genomic consensus and call variants relative to the reference.

 

optional arguments:

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

  --version             show program's version number and exit

  --emit-tool-contract  Emit Tool Contract to stdout (default: False)

  --resolved-tool-contract RESOLVED_TOOL_CONTRACT

                        Run Tool directly from a PacBio Resolved tool contract

                        (default: None)

  --log-file LOG_FILE   Write the log to file. Default(None) will write to

                        stdout. (default: None)

  --log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}

                        Set log level (default: WARN)

  --debug               Alias for setting log level to DEBUG (default: False)

  --quiet               Alias for setting log level to CRITICAL to suppress

                        output. (default: False)

  -v, --verbose         Set the verbosity level. (default: None)

 

Basic required options:

  inputFilename         The input cmp.h5 or BAM alignment file

  --referenceFilename REFERENCEFILENAME, --reference REFERENCEFILENAME, -r REFERENCEFILENAME

                        The filename of the reference FASTA file (default:

                        None)

  -o OUTPUTFILENAMES, --outputFilename OUTPUTFILENAMES

                        The output filename(s), as a comma-separated

                        list.Valid output formats are .fa/.fasta, .fq/.fastq,

                        .gff, .vcf (default: )

 

Parallelism:

  -j NUMWORKERS, --numWorkers NUMWORKERS

                        The number of worker processes to be used (default: 1)

 

Output filtering:

  --minConfidence MINCONFIDENCE, -q MINCONFIDENCE

                        The minimum confidence for a variant call to be output

                        to variants.{gff,vcf} (default: 40)

  --minCoverage MINCOVERAGE, -x MINCOVERAGE

                        The minimum site coverage that must be achieved for

                        variant calls and consensus to be calculated for a

                        site. (default: 5)

  --noEvidenceConsensusCall {nocall,reference,lowercasereference}

                        The consensus base that will be output for sites with

                        no effective coverage. (default: lowercasereference)

 

Read selection/filtering:

  --coverage COVERAGE, -X COVERAGE

                        A designation of the maximum coverage level to be used

                        for analysis. Exact interpretation is algorithm-

                        specific. (default: 100)

  --minMapQV MINMAPQV, -m MINMAPQV

                        The minimum MapQV for reads that will be used for

                        analysis. (default: 10)

  --referenceWindow REFERENCEWINDOWSASSTRING, --referenceWindows REFERENCEWINDOWSASSTRING, -w REFERENCEWINDOWSASSTRING

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd (default: entire reference).

                        (default: None)

  --alignmentSetRefWindows

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd will be pulled from the alignment

                        file. (default: False)

  --referenceWindowsFile REFERENCEWINDOWSASSTRING, -W REFERENCEWINDOWSASSTRING

                        A file containing reference window designations, one

                        per line (default: None)

  --barcode _BARCODE    Only process reads with the given barcode name.

                        (default: None)

  --readStratum READSTRATUM

                        A string of the form 'n/N', where n, and N are

                        integers, 0 <= n < N, designating that the reads are

                        to be deterministically split into N strata of roughly

                        even size, and stratum n is to be used for variant and

                        consensus calling. This is mostly useful for Quiver

                        development. (default: None)

  --minReadScore MINREADSCORE

                        The minimum ReadScore for reads that will be used for

                        analysis (arrow-only). (default: 0.65)

  --minSnr MINHQREGIONSNR

                        The minimum acceptable signal-to-noise over all

                        channels for reads that will be used for analysis

                        (arrow-only). (default: 2.5)

  --minZScore MINZSCORE

                        The minimum acceptable z-score for reads that will be

                        used for analysis (arrow-only). (default: -3.5)

  --minAccuracy MINACCURACY

                        The minimum acceptable window-global alignment

                        accuracy for reads that will be used for the analysis

                        (arrow-only). (default: 0.82)

 

Algorithm and parameter settings:

  --algorithm {quiver,arrow,plurality,poa,best}

  --parametersFile PARAMETERSFILE, -P PARAMETERSFILE

                        Parameter set filename (such as ArrowParameters.json

                        or QuiverParameters.ini), or directory D such that

                        either D/*/GenomicConsensus/QuiverParameters.ini, or

                        D/GenomicConsensus/QuiverParameters.ini, is found. In

                        the former case, the lexically largest path is chosen.

                        (default: None)

  --parametersSpec PARAMETERSSPEC, -p PARAMETERSSPEC

                        Name of parameter set (chemistry.model) to select from

                        the parameters file, or just the name of the

                        chemistry, in which case the best available model is

                        chosen. Default is 'auto', which selects the best

                        parameter set from the alignment data (default: auto)

  --maskRadius MASKRADIUS

                        Radius of window to use when excluding local regions

                        for exceeding maskMinErrorRate, where 0 disables any

                        filtering (arrow-only). (default: 3)

  --maskErrorRate MASKERRORRATE

                        Maximum local error rate before the local region

                        defined by maskRadius is excluded from polishing

                        (arrow-only). (default: 0.7)

 

Verbosity and debugging/profiling:

  --pdb                 Enable Python debugger (default: False)

  --notrace             Suppress stacktrace for exceptions (to simplify

                        testing) (default: False)

  --pdbAtStartup        Drop into Python debugger at startup (requires ipdb)

                        (default: False)

  --profile             Enable Python-level profiling (using cProfile).

                        (default: False)

  --annotateGFF         Augment GFF variant records with additional

                        information (default: False)

  --reportEffectiveCoverage

                        Additionally record the *post-filtering* coverage at

                        variant sites (default: False)

 

Advanced configuration options:

  --diploid             Enable detection of heterozygous variants

                        (experimental) (default: False)

  --queueSize QUEUESIZE, -Q QUEUESIZE

  --threaded, -T        Run threads instead of processes (for debugging

                        purposes only) (default: False)

  --referenceChunkSize REFERENCECHUNKSIZE, -C REFERENCECHUNKSIZE

  --fancyChunking       Adaptive reference chunking designed to handle

                        coverage cutouts better (default: True)

  --simpleChunking      Disable adaptive reference chunking (default: True)

  --referenceChunkOverlap REFERENCECHUNKOVERLAP

  --autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE

                        Disable the HDF5 chunk cache when the number of

                        datasets in the cmp.h5 exceeds the given threshold

                        (default: 500)

  --aligner {affine,simple}, -a {affine,simple}

                        The pairwise alignment algorithm that will be used to

                        produce variant calls from the consensus (Quiver

                        only). (default: affine)

  --refineDinucleotideRepeats

                        Require quiver maximum likelihood search to try one

                        less/more repeat copy in dinucleotide repeats, which

                        seem to be the most frequent cause of suboptimal

                        convergence (getting trapped in local optimum) (Quiver

                        only) (default: True)

  --noRefineDinucleotideRepeats

                        Disable dinucleotide refinement (default: True)

  --fast                Cut some corners to run faster. Unsupported! (default:

                        False)

  --skipUnrecognizedContigs

                        Do not abort when told to process a reference window

                        (via -w/--referenceWindow[s]) that has no aligned

                        coverage. Outputs emptyish files if there are no

                        remaining non-degenerate windows. Only intended for

                        use by smrtpipe scatter/gather. (default: False)

> plurality -h

$ plurality -h

usage: variantCaller [-h] [--version] [--emit-tool-contract]

                     [--resolved-tool-contract RESOLVED_TOOL_CONTRACT]

                     [--log-file LOG_FILE]

                     [--log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL} | --debug | --quiet | -v]

                     --referenceFilename REFERENCEFILENAME -o OUTPUTFILENAMES

                     [-j NUMWORKERS] [--minConfidence MINCONFIDENCE]

                     [--minCoverage MINCOVERAGE]

                     [--noEvidenceConsensusCall {nocall,reference,lowercasereference}]

                     [--coverage COVERAGE] [--minMapQV MINMAPQV]

                     [--referenceWindow REFERENCEWINDOWSASSTRING]

                     [--alignmentSetRefWindows]

                     [--referenceWindowsFile REFERENCEWINDOWSASSTRING]

                     [--barcode _BARCODE] [--readStratum READSTRATUM]

                     [--minReadScore MINREADSCORE] [--minSnr MINHQREGIONSNR]

                     [--minZScore MINZSCORE] [--minAccuracy MINACCURACY]

                     [--algorithm {quiver,arrow,plurality,poa,best}]

                     [--parametersFile PARAMETERSFILE]

                     [--parametersSpec PARAMETERSSPEC]

                     [--maskRadius MASKRADIUS] [--maskErrorRate MASKERRORRATE]

                     [--pdb] [--notrace] [--pdbAtStartup] [--profile]

                     [--annotateGFF] [--reportEffectiveCoverage] [--diploid]

                     [--queueSize QUEUESIZE] [--threaded]

                     [--referenceChunkSize REFERENCECHUNKSIZE]

                     [--fancyChunking] [--simpleChunking]

                     [--referenceChunkOverlap REFERENCECHUNKOVERLAP]

                     [--autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE]

                     [--aligner {affine,simple}] [--refineDinucleotideRepeats]

                     [--noRefineDinucleotideRepeats] [--fast]

                     [--skipUnrecognizedContigs]

                     inputFilename

 

Compute genomic consensus and call variants relative to the reference.

 

optional arguments:

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

  --version             show program's version number and exit

  --emit-tool-contract  Emit Tool Contract to stdout (default: False)

  --resolved-tool-contract RESOLVED_TOOL_CONTRACT

                        Run Tool directly from a PacBio Resolved tool contract

                        (default: None)

  --log-file LOG_FILE   Write the log to file. Default(None) will write to

                        stdout. (default: None)

  --log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}

                        Set log level (default: WARN)

  --debug               Alias for setting log level to DEBUG (default: False)

  --quiet               Alias for setting log level to CRITICAL to suppress

                        output. (default: False)

  -v, --verbose         Set the verbosity level. (default: None)

 

Basic required options:

  inputFilename         The input cmp.h5 or BAM alignment file

  --referenceFilename REFERENCEFILENAME, --reference REFERENCEFILENAME, -r REFERENCEFILENAME

                        The filename of the reference FASTA file (default:

                        None)

  -o OUTPUTFILENAMES, --outputFilename OUTPUTFILENAMES

                        The output filename(s), as a comma-separated

                        list.Valid output formats are .fa/.fasta, .fq/.fastq,

                        .gff, .vcf (default: )

 

Parallelism:

  -j NUMWORKERS, --numWorkers NUMWORKERS

                        The number of worker processes to be used (default: 1)

 

Output filtering:

  --minConfidence MINCONFIDENCE, -q MINCONFIDENCE

                        The minimum confidence for a variant call to be output

                        to variants.{gff,vcf} (default: 40)

  --minCoverage MINCOVERAGE, -x MINCOVERAGE

                        The minimum site coverage that must be achieved for

                        variant calls and consensus to be calculated for a

                        site. (default: 5)

  --noEvidenceConsensusCall {nocall,reference,lowercasereference}

                        The consensus base that will be output for sites with

                        no effective coverage. (default: lowercasereference)

 

Read selection/filtering:

  --coverage COVERAGE, -X COVERAGE

                        A designation of the maximum coverage level to be used

                        for analysis. Exact interpretation is algorithm-

                        specific. (default: 100)

  --minMapQV MINMAPQV, -m MINMAPQV

                        The minimum MapQV for reads that will be used for

                        analysis. (default: 10)

  --referenceWindow REFERENCEWINDOWSASSTRING, --referenceWindows REFERENCEWINDOWSASSTRING, -w REFERENCEWINDOWSASSTRING

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd (default: entire reference).

                        (default: None)

  --alignmentSetRefWindows

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd will be pulled from the alignment

                        file. (default: False)

  --referenceWindowsFile REFERENCEWINDOWSASSTRING, -W REFERENCEWINDOWSASSTRING

                        A file containing reference window designations, one

                        per line (default: None)

  --barcode _BARCODE    Only process reads with the given barcode name.

                        (default: None)

  --readStratum READSTRATUM

                        A string of the form 'n/N', where n, and N are

                        integers, 0 <= n < N, designating that the reads are

                        to be deterministically split into N strata of roughly

                        even size, and stratum n is to be used for variant and

                        consensus calling. This is mostly useful for Quiver

                        development. (default: None)

  --minReadScore MINREADSCORE

                        The minimum ReadScore for reads that will be used for

                        analysis (arrow-only). (default: 0.65)

  --minSnr MINHQREGIONSNR

                        The minimum acceptable signal-to-noise over all

                        channels for reads that will be used for analysis

                        (arrow-only). (default: 2.5)

  --minZScore MINZSCORE

                        The minimum acceptable z-score for reads that will be

                        used for analysis (arrow-only). (default: -3.5)

  --minAccuracy MINACCURACY

                        The minimum acceptable window-global alignment

                        accuracy for reads that will be used for the analysis

                        (arrow-only). (default: 0.82)

 

Algorithm and parameter settings:

  --algorithm {quiver,arrow,plurality,poa,best}

  --parametersFile PARAMETERSFILE, -P PARAMETERSFILE

                        Parameter set filename (such as ArrowParameters.json

                        or QuiverParameters.ini), or directory D such that

                        either D/*/GenomicConsensus/QuiverParameters.ini, or

                        D/GenomicConsensus/QuiverParameters.ini, is found. In

                        the former case, the lexically largest path is chosen.

                        (default: None)

  --parametersSpec PARAMETERSSPEC, -p PARAMETERSSPEC

                        Name of parameter set (chemistry.model) to select from

                        the parameters file, or just the name of the

                        chemistry, in which case the best available model is

                        chosen. Default is 'auto', which selects the best

                        parameter set from the alignment data (default: auto)

  --maskRadius MASKRADIUS

                        Radius of window to use when excluding local regions

                        for exceeding maskMinErrorRate, where 0 disables any

                        filtering (arrow-only). (default: 3)

  --maskErrorRate MASKERRORRATE

                        Maximum local error rate before the local region

                        defined by maskRadius is excluded from polishing

                        (arrow-only). (default: 0.7)

 

Verbosity and debugging/profiling:

  --pdb                 Enable Python debugger (default: False)

  --notrace             Suppress stacktrace for exceptions (to simplify

                        testing) (default: False)

  --pdbAtStartup        Drop into Python debugger at startup (requires ipdb)

                        (default: False)

  --profile             Enable Python-level profiling (using cProfile).

                        (default: False)

  --annotateGFF         Augment GFF variant records with additional

                        information (default: False)

  --reportEffectiveCoverage

                        Additionally record the *post-filtering* coverage at

                        variant sites (default: False)

 

Advanced configuration options:

  --diploid             Enable detection of heterozygous variants

                        (experimental) (default: False)

  --queueSize QUEUESIZE, -Q QUEUESIZE

  --threaded, -T        Run threads instead of processes (for debugging

                        purposes only) (default: False)

  --referenceChunkSize REFERENCECHUNKSIZE, -C REFERENCECHUNKSIZE

  --fancyChunking       Adaptive reference chunking designed to handle

                        coverage cutouts better (default: True)

  --simpleChunking      Disable adaptive reference chunking (default: True)

  --referenceChunkOverlap REFERENCECHUNKOVERLAP

  --autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE

                        Disable the HDF5 chunk cache when the number of

                        datasets in the cmp.h5 exceeds the given threshold

                        (default: 500)

  --aligner {affine,simple}, -a {affine,simple}

                        The pairwise alignment algorithm that will be used to

                        produce variant calls from the consensus (Quiver

                        only). (default: affine)

  --refineDinucleotideRepeats

                        Require quiver maximum likelihood search to try one

                        less/more repeat copy in dinucleotide repeats, which

                        seem to be the most frequent cause of suboptimal

                        convergence (getting trapped in local optimum) (Quiver

                        only) (default: True)

  --noRefineDinucleotideRepeats

                        Disable dinucleotide refinement (default: True)

  --fast                Cut some corners to run faster. Unsupported! (default:

                        False)

  --skipUnrecognizedContigs

                        Do not abort when told to process a reference window

                        (via -w/--referenceWindow[s]) that has no aligned

                        coverage. Outputs emptyish files if there are no

                        remaining non-degenerate windows. Only intended for

                        use by smrtpipe scatter/gather. (default: False)

> variantCaller -h

$ variantCaller -h

usage: variantCaller [-h] [--version] [--emit-tool-contract]

                     [--resolved-tool-contract RESOLVED_TOOL_CONTRACT]

                     [--log-file LOG_FILE]

                     [--log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL} | --debug | --quiet | -v]

                     --referenceFilename REFERENCEFILENAME -o OUTPUTFILENAMES

                     [-j NUMWORKERS] [--minConfidence MINCONFIDENCE]

                     [--minCoverage MINCOVERAGE]

                     [--noEvidenceConsensusCall {nocall,reference,lowercasereference}]

                     [--coverage COVERAGE] [--minMapQV MINMAPQV]

                     [--referenceWindow REFERENCEWINDOWSASSTRING]

                     [--alignmentSetRefWindows]

                     [--referenceWindowsFile REFERENCEWINDOWSASSTRING]

                     [--barcode _BARCODE] [--readStratum READSTRATUM]

                     [--minReadScore MINREADSCORE] [--minSnr MINHQREGIONSNR]

                     [--minZScore MINZSCORE] [--minAccuracy MINACCURACY]

                     [--algorithm {quiver,arrow,plurality,poa,best}]

                     [--parametersFile PARAMETERSFILE]

                     [--parametersSpec PARAMETERSSPEC]

                     [--maskRadius MASKRADIUS] [--maskErrorRate MASKERRORRATE]

                     [--pdb] [--notrace] [--pdbAtStartup] [--profile]

                     [--annotateGFF] [--reportEffectiveCoverage] [--diploid]

                     [--queueSize QUEUESIZE] [--threaded]

                     [--referenceChunkSize REFERENCECHUNKSIZE]

                     [--fancyChunking] [--simpleChunking]

                     [--referenceChunkOverlap REFERENCECHUNKOVERLAP]

                     [--autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE]

                     [--aligner {affine,simple}] [--refineDinucleotideRepeats]

                     [--noRefineDinucleotideRepeats] [--fast]

                     [--skipUnrecognizedContigs]

                     inputFilename

 

Compute genomic consensus and call variants relative to the reference.

 

optional arguments:

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

  --version             show program's version number and exit

  --emit-tool-contract  Emit Tool Contract to stdout (default: False)

  --resolved-tool-contract RESOLVED_TOOL_CONTRACT

                        Run Tool directly from a PacBio Resolved tool contract

                        (default: None)

  --log-file LOG_FILE   Write the log to file. Default(None) will write to

                        stdout. (default: None)

  --log-level {DEBUG,INFO,WARNING,ERROR,CRITICAL}

                        Set log level (default: WARN)

  --debug               Alias for setting log level to DEBUG (default: False)

  --quiet               Alias for setting log level to CRITICAL to suppress

                        output. (default: False)

  -v, --verbose         Set the verbosity level. (default: None)

 

Basic required options:

  inputFilename         The input cmp.h5 or BAM alignment file

  --referenceFilename REFERENCEFILENAME, --reference REFERENCEFILENAME, -r REFERENCEFILENAME

                        The filename of the reference FASTA file (default:

                        None)

  -o OUTPUTFILENAMES, --outputFilename OUTPUTFILENAMES

                        The output filename(s), as a comma-separated

                        list.Valid output formats are .fa/.fasta, .fq/.fastq,

                        .gff, .vcf (default: )

 

Parallelism:

  -j NUMWORKERS, --numWorkers NUMWORKERS

                        The number of worker processes to be used (default: 1)

 

Output filtering:

  --minConfidence MINCONFIDENCE, -q MINCONFIDENCE

                        The minimum confidence for a variant call to be output

                        to variants.{gff,vcf} (default: 40)

  --minCoverage MINCOVERAGE, -x MINCOVERAGE

                        The minimum site coverage that must be achieved for

                        variant calls and consensus to be calculated for a

                        site. (default: 5)

  --noEvidenceConsensusCall {nocall,reference,lowercasereference}

                        The consensus base that will be output for sites with

                        no effective coverage. (default: lowercasereference)

 

Read selection/filtering:

  --coverage COVERAGE, -X COVERAGE

                        A designation of the maximum coverage level to be used

                        for analysis. Exact interpretation is algorithm-

                        specific. (default: 100)

  --minMapQV MINMAPQV, -m MINMAPQV

                        The minimum MapQV for reads that will be used for

                        analysis. (default: 10)

  --referenceWindow REFERENCEWINDOWSASSTRING, --referenceWindows REFERENCEWINDOWSASSTRING, -w REFERENCEWINDOWSASSTRING

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd (default: entire reference).

                        (default: None)

  --alignmentSetRefWindows

                        The window (or multiple comma-delimited windows) of

                        the reference to be processed, in the format refGroup

                        :refStart-refEnd will be pulled from the alignment

                        file. (default: False)

  --referenceWindowsFile REFERENCEWINDOWSASSTRING, -W REFERENCEWINDOWSASSTRING

                        A file containing reference window designations, one

                        per line (default: None)

  --barcode _BARCODE    Only process reads with the given barcode name.

                        (default: None)

  --readStratum READSTRATUM

                        A string of the form 'n/N', where n, and N are

                        integers, 0 <= n < N, designating that the reads are

                        to be deterministically split into N strata of roughly

                        even size, and stratum n is to be used for variant and

                        consensus calling. This is mostly useful for Quiver

                        development. (default: None)

  --minReadScore MINREADSCORE

                        The minimum ReadScore for reads that will be used for

                        analysis (arrow-only). (default: 0.65)

  --minSnr MINHQREGIONSNR

                        The minimum acceptable signal-to-noise over all

                        channels for reads that will be used for analysis

                        (arrow-only). (default: 2.5)

  --minZScore MINZSCORE

                        The minimum acceptable z-score for reads that will be

                        used for analysis (arrow-only). (default: -3.5)

  --minAccuracy MINACCURACY

                        The minimum acceptable window-global alignment

                        accuracy for reads that will be used for the analysis

                        (arrow-only). (default: 0.82)

 

Algorithm and parameter settings:

  --algorithm {quiver,arrow,plurality,poa,best}

  --parametersFile PARAMETERSFILE, -P PARAMETERSFILE

                        Parameter set filename (such as ArrowParameters.json

                        or QuiverParameters.ini), or directory D such that

                        either D/*/GenomicConsensus/QuiverParameters.ini, or

                        D/GenomicConsensus/QuiverParameters.ini, is found. In

                        the former case, the lexically largest path is chosen.

                        (default: None)

  --parametersSpec PARAMETERSSPEC, -p PARAMETERSSPEC

                        Name of parameter set (chemistry.model) to select from

                        the parameters file, or just the name of the

                        chemistry, in which case the best available model is

                        chosen. Default is 'auto', which selects the best

                        parameter set from the alignment data (default: auto)

  --maskRadius MASKRADIUS

                        Radius of window to use when excluding local regions

                        for exceeding maskMinErrorRate, where 0 disables any

                        filtering (arrow-only). (default: 3)

  --maskErrorRate MASKERRORRATE

                        Maximum local error rate before the local region

                        defined by maskRadius is excluded from polishing

                        (arrow-only). (default: 0.7)

 

Verbosity and debugging/profiling:

  --pdb                 Enable Python debugger (default: False)

  --notrace             Suppress stacktrace for exceptions (to simplify

                        testing) (default: False)

  --pdbAtStartup        Drop into Python debugger at startup (requires ipdb)

                        (default: False)

  --profile             Enable Python-level profiling (using cProfile).

                        (default: False)

  --annotateGFF         Augment GFF variant records with additional

                        information (default: False)

  --reportEffectiveCoverage

                        Additionally record the *post-filtering* coverage at

                        variant sites (default: False)

 

Advanced configuration options:

  --diploid             Enable detection of heterozygous variants

                        (experimental) (default: False)

  --queueSize QUEUESIZE, -Q QUEUESIZE

  --threaded, -T        Run threads instead of processes (for debugging

                        purposes only) (default: False)

  --referenceChunkSize REFERENCECHUNKSIZE, -C REFERENCECHUNKSIZE

  --fancyChunking       Adaptive reference chunking designed to handle

                        coverage cutouts better (default: True)

  --simpleChunking      Disable adaptive reference chunking (default: True)

  --referenceChunkOverlap REFERENCECHUNKOVERLAP

  --autoDisableHdf5ChunkCache AUTODISABLEHDF5CHUNKCACHE

                        Disable the HDF5 chunk cache when the number of

                        datasets in the cmp.h5 exceeds the given threshold

                        (default: 500)

  --aligner {affine,simple}, -a {affine,simple}

                        The pairwise alignment algorithm that will be used to

                        produce variant calls from the consensus (Quiver

                        only). (default: affine)

  --refineDinucleotideRepeats

                        Require quiver maximum likelihood search to try one

                        less/more repeat copy in dinucleotide repeats, which

                        seem to be the most frequent cause of suboptimal

                        convergence (getting trapped in local optimum) (Quiver

                        only) (default: True)

  --noRefineDinucleotideRepeats

                        Disable dinucleotide refinement (default: True)

  --fast                Cut some corners to run faster. Unsupported! (default:

                        False)

  --skipUnrecognizedContigs

                        Do not abort when told to process a reference window

                        (via -w/--referenceWindow[s]) that has no aligned

                        coverage. Outputs emptyish files if there are no

                        remaining non-degenerate windows. Only intended for

                        use by smrtpipe scatter/gather. (default: False)

 

実行方法 

マッピング済みのbamとリファレンスを指定する。ここでは8スレッド指定した。

quiver -j8 aligned_reads{.bam 
> -r reference{.fasta or .xml}
-o variants.gff -o consensus.fasta -o consensus.fastq
  •   -j    The number of worker processes to be used (default: 1)

variantCaller --algorithm=quiver arrow | pluralityでもランできます。

 

出力(Githubより)

  1. A consensus FASTA file containing the consensus sequence
  2. A consensus FASTQ file containing the consensus sequence with quality annotations
  3. A variants GFF file containing a filtered, annotated list of variants identified

引用

ref.1

Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data
Chen-Shan Chin, David H Alexander, Patrick Marks, Aaron A Klammer, James Drake, Cheryl Heiner, Alicia Clum, Alex Copeland, John Huddleston, Evan E Eichler, Stephen W Turner & Jonas Korlach
Nature Methods volume 10, pages 563–569 (2013)

 

参考

Polishing PacBio assemblies with Arrow and Pilon

https://flowersoftheocean.wordpress.com/2018/04/16/polishing-pacbio-assemblies-with-arrow-and-pilon/

 

Quiver と Arrowの話

https://pacbiobrothers.blogspot.com/2016/12/quver-arrow.html

 

関連

 

ロングリードを使ってハプロタイプフェージングを行う HapCol

 

 ヒトなどの二倍体生物は、それぞれの親から1つずつ、2組の染色体を含んでいる。ハプロタイプと呼ばれる、各染色体の2つの異なるコピーを再構築することは、個体のゲノムを特徴付けるために重要である。このプロセスは、フェージングまたはハプロタイピングとして知られており、提供された情報は、遺伝性変異と遺伝子機能との間の関係、または遺伝性変異と疾患感受性との間の関係の分析などの多くの用途にとって根本的に重要であり得る。二倍体種において、ハプロタイピングは、各染色体の2つの親コピーに変異を割り当てることを必要とし、それは一塩基多型(SNP)の点で相違を示す。収集したサンプルからハプロタイプを大規模に直接実験的に再構築することはまだ費用効率が良くないので(Kuleshov et al、2014)、染色体コピーからシーケンシングされたリードのセットを考慮する計算アプローチ(ハプロタイプアセンブリ) が提案されている。可能であれば、予備のマッピング段階でリファレンスゲノムを用いて、リード(フラグメントとも呼ばれる)を未知のハプロタイプに割り当てる必要がある。これは、何らかの方法でシーケンシングおよびマッピングエラーを扱うことを含み、そして一般に最適化問題としてモデル化される計算タスクをもたらす(Lancia et al、2001; Lippert et al、2002)。

 Minimum error correction (MEC) (Lippert et al、2002)は、ハプロタイプ構築のための著名な組み合わせアプローチの一つである。それは、SNPの値に対する最小数の補正で入力データを補正することを目的とし、その結果得られるリードは、それぞれハプロタイプを識別する2つのセットに明確に分割することができる。 wMEC (Greenberg et al、2004) は、問題の重み付き変形であり、各可能な補正は、対応する位置でそのSNP値に割り当てられた信頼度を表す重みと関連付けられる。この信頼度は、そのベースコールについてのシーケンシング中にエラーが発生した確率(フレッドベースのエラー確率)と、そのゲノム位置へのリードマッピングの信頼度との組み合わせである。そのような重みの使用は、精度を向上させる強力な方法として実験的に検証されている(Zhao et al、2005)。

 次世代シーケンシング(NGS)技術の出現は、二倍体生物のゲノムを組み立てる費用効果の高い方法を提供した。しかしながら、正確なハプロタイプアセンブリするためには、いくつかの異なるヘテロ接合位置にまたがるのに十分な長さのリードを有することが必要である(Duitama et al、2012)。この種のデータは、PacBio RS II(http://www.pacificbiosciences.com/products/)やOxford Nanoporeフローセルのような一分子リアルタイムシーケンシング技術のような「未来世代」のシーケンシング技術の出現でますます利用可能になりつつある。

 これらの技術は、10000塩基を超えるシングルエンドリードを生成するそれらの能力のおかげで、ペアエンドデータの必要性を排除し、そしてゲノムFinishingおよびハプロタイプアセンブリのようなタスクに既に使用されている(Smith et al、2012)。リード長に加えて、次世代のシーケンシング技術は、シーケンシングエラーが一様分布する新規な特徴を持つフラグメントを生成する。

 最近、MECおよびwMECアプローチがロングリードの文脈で使用されており、長いフラグメントが伝統的なショートリードよりも正確にハプロタイプアセンブリすることを可能にすることを確認している(Aguiar and Istrail、2012; Duitama et al、2012; Patterson et al、2014、 2015)。 MECはNP困難なので(Cilibrasi et al、2007)、厳密な解は指数関数的に複雑になる。問題の計算困難性に取り組む様々なアプローチが文献で提案されている。Integer linear programming テクニックが最近使用されている(Chen et al、2013)が、このアプローチではいくつかの「難しいブロック」を最適に解くことができなかった。(一部略)

 WHATSHAP(Patterson et al、2014)と呼ばれる2番目のアプローチは、ロングリードを処理できるwMECの最初の正確なアルゴリズムである。 20xまでのカバレッジ範囲でシミュレートされたロングリードのデータに対して良好な精度を得ることができ、そして全ての以前の正確なアプローチよりも優れていることが示された。ただし、20倍を超えるカバレッジを処理することはできず、その限界に近づくとそのパフォーマンスは明らかに低下する。(複数段落省略)

 カバレッジ15〜20倍のシミュレーションデータセットの結果は、HAPCOLはWHATSHAPと同じくらい正確だが(平均誤差約2%を達成)、より高速で大幅にメモリ効率が高いことを示している(約2倍高速、約28倍メモリ使用量が少ない)。 HAPCOLの効率により、精度をさらに向上させることができる。実際、実験結果は、HAPCOLが標準ワークステーション/小型サーバー上で25倍のカバレッジでデータセットを処理できることを示している(一方、WHATSHAPはすべての利用可能なメモリ、256 GBを使い果たした)。カバレッジ25倍でHAPCOLによって再構成されたものは、カバレッジ20倍で再構成されたものよりも約9%正確である。

 

HP

http://hapcol.algolab.eu

 

インストール

依存

  • CMake (>= 2.8)
  • GNU make

本体 Github

#1GBほどあるのでダウンロード時は注意
git clone https://github.com/AlgoLab/HapCol.git
cd HapCol/
mkdir -p build
cd build
cmake ../src
make

 > ./hapcol -h

# ./hapcol -h

* INFO (main        :c/HapCol.cpp:167 ) 10:09:30 | HapCol (master@97d4a5e-clean)

* INFO (main        :c/HapCol.cpp:171 ) 10:09:30 | Starting...

Program options:

  -h [ --help ] produce (this) help message

  -i [ --input ] arg file containing the input reads (in WIF

format)

  -o [ --haplotypes ] arg file where the computed haplotypes will

be written to

  -u [ --discard-weights ] discard weights

  -x [ --no-ambiguous ] do not mark ambiguous positions with Xs

  -A [ --all-heterozygous ] all-heterozygous assumption

  -U [ --unique ] input as unique block

  -e [ --error-rate ] arg (=0.05) read error rate

  -a [ --alpha ] arg (=0.01) significance (smaller is better)

 

* INFO (main        :c/HapCol.cpp:178 ) 10:09:30 | Arguments:

* INFO (main        :c/HapCol.cpp:179 ) 10:09:30 | Initialized? False

* INFO (main        :c/HapCol.cpp:180 ) 10:09:30 | Input filename: ''

* INFO (main        :c/HapCol.cpp:181 ) 10:09:30 | Haplotype filename: ''

* INFO (main        :c/HapCol.cpp:182 ) 10:09:30 | Discard weights? False

* INFO (main        :c/HapCol.cpp:183 ) 10:09:30 | Do not add X's? False

* INFO (main        :c/HapCol.cpp:184 ) 10:09:30 | All-heterozygous assumption? False

* INFO (main        :c/HapCol.cpp:185 ) 10:09:30 | Input as unique block? False

* INFO (main        :c/HapCol.cpp:186 ) 10:09:30 | Error rate: 0.05

* INFO (main        :c/HapCol.cpp:187 ) 10:09:30 | Alpha: 0.01

* FATAL(main        :c/HapCol.cpp:190 ) 10:09:30 | Arguments not correctly initialized! Exiting..


 

テストラン

> head ../docs/sample.wif

# head ../docs/sample.wif

36950998 C 0 10 : 36951355 G 1 7 : 36951968 G 0 3 : 36952090 C 0 11 : 36952098 C 1 28 : 36954126 A 0 15 : 36961286 G 0 20 : 36961747 C 0 23 : 36961983 T 1 24 : 36966400 G 0 21 : 36966495 T 0 26 : 36984266 T 0 31 : 36993308 G 0 10 : # 60 : NA

36950998 C 0 12 : 36951355 G 1 7 : 36951968 G 0 22 : 36952090 C 0 33 : 36952098 C 1 23 : 36954126 A 0 4 : 36961286 G 0 26 : 36961747 C 0 34 : 36961983 T 1 33 : 36966400 G 0 29 : 36966495 T 0 24 : # 60 : NA

36950998 C 0 12 : 36951355 G 1 9 : 36951968 G 0 26 : 36952090 C 0 36 : 36952098 C 1 38 : 36954126 A 0 27 : 36961286 G 0 29 : 36961747 C 0 24 : 36961983 T 1 37 : 36966400 G 0 32 : 36966495 T 0 8 : 36984266 T 0 38 : 36993308 G 0 35 : # 60 : NA

36950998 C 0 13 : 36951355 G 1 37 : 36951968 G 0 19 : 36952090 C 0 31 : 36952098 C 1 29 : 36954126 A 0 17 : 36961286 G 0 33 : 36961747 C 0 24 : 36961983 T 1 35 : 36966400 G 0 34 : 36966495 T 0 11 : # 60 : NA

36950998 C 0 14 : 36951355 G 1 8 : 36951968 G 0 20 : 36952090 C 0 17 : 36952098 C 1 23 : 36954126 A 0 14 : 36961286 G 0 26 : 36961747 C 0 35 : 36961983 T 1 23 : 36966400 G 0 19 : 36966495 T 0 21 : 36984266 T 0 19 : # 60 : NA

36950998 C 0 16 : 36951355 G 1 32 : 36951968 G 0 24 : 36952090 C 0 6 : 36952098 C 1 11 : 36954126 A 0 32 : 36961286 G 0 6 : 36961747 C 0 20 : 36961983 T 1 31 : 36966400 G 0 11 : 36966495 T 0 10 : # 60 : NA

36950998 C 0 17 : 36951355 G 1 31 : 36951968 G 0 24 : 36952090 C 0 11 : 36952098 C 1 36 : 36954126 A 0 24 : 36961286 G 0 24 : 36961747 C 0 22 : 36961983 T 1 27 : 36966400 G 0 26 : 36966495 T 0 26 : # 60 : NA

36950998 C 0 17 : 36951355 G 1 37 : 36951968 G 0 19 : 36952090 C 0 36 : 36952098 C 1 23 : 36954126 A 0 15 : 36961286 G 0 3 : 36961747 C 0 38 : 36961983 T 1 8 : 36966400 G 0 4 : 36966495 T 0 15 : 36984266 T 0 36 : 36993308 G 0 9 : # 60 : NA

36950998 C 0 20 : 36951355 G 1 34 : 36951968 G 0 24 : 36952090 C 0 32 : 36952098 C 1 26 : 36954126 A 0 36 : 36961286 G 0 9 : 36961747 C 0 33 : 36961983 T 1 2 : 36966400 G 0 35 : 36966495 T 0 28 : 36984266 T 0 15 : 36993308 G 0 33 : # 60 : NA

36950998

どのリードがどのSNPアレルをどれだけカバーしているかをwifフォーマット (introduced in (Patterson et al., RECOMB, 2014)) で記載している。このwifフォーマットについてはGithubの"Input format"の項目で説明されている。

 

.wifファイルを指定して実行

./hapcol -i ../docs/sample.wif -o haplotypes.txt
  • --input (or -i), which specifies the file containing the reads in input (in WIF format);
  • --output (or -o), which specifies the file for the computed haplotypes.

出力

>head haplotypes.txt 

# head haplotypes.txt 

101101110111110110111100101111011100011001110110100010010X10100100111001110010011100100111110011010001100010000010100010010011001100110011110011110101101001011011011010111101001001110111101001010101111101111101000110100110101000010101011100111100111001100001111000100110000101110010010111010111011000110110010101000001010011000111010101000001010010110101111100000101001101100111000100001001110100100111111000000111011001100111000000010001011100001010100100010001110010010011100110101111100011010011000100011000111101101010001010101001011111110010000

010010001000001001000011010000100011100110001001011101101001011011000110001101100011011000001100101110011101111101011101101100110011001100001100001010010110100100000101000010110110001000010110101010000010000010111001011001010111101010100011000011000110011110000111011001111010001101101000101000100111001001101010111110101100111000101010111110101101001010000011111010110010011000111011110110001011011000000111111000100110011000111111101110100011110101011011101110001101101100011001010000011100101100111011100111000010010101110101010110100000001101111

 

引用

HAPCOL: accurate and memory-efficient haplotype assembly from long reads
Yuri Pirola, Simone Zaccaria, Riccardo Dondi, Gunnar W. Klau, Nadia Pisanti, Paola Bonizzoni
Bioinformatics, Volume 32, Issue 11, 1 June 2016, Pages 1610–1617,

 

関連

 

 

pacbioのロングリードの構造変異検出ツール pbsv

 

pbsvは、PacBio一分子リアルタイムシークエンシング(SMRT)リードから二倍体ゲノムの構造変異をコールして分析するための一連のツールである。 このツールは、PacBioのSMRT Link GUIのStructural Variant Calling分析ワークフローを強化する。

pbsvは挿入、欠失、逆位、duplications、および転座をコールする。シングルサンプルコールとジョイント(マルチサンプル)コールの両方がサポートされている。

pbsvの効果的な範囲

  • 20 bpから10 kbまでの挿入
  • 20塩基対から100塩基対までの欠失
  • 200 bpから10 kbへの逆位
  • 20 bpから10 kbまでのduplications
  • 異なる染色体間の転座、あるいは単一の染色体上で100kb以上離れている 

 

Structural Variant Detection in SMRT Link 5 with pbsv

f:id:kazumaxneo:20190317193217p:plain

 

 

インストール

ubuntu16.04でテストした(ホストOS macos10.12)。

PacBio Secondary Analysis Tools on Bioconda

本体 Github

#anacondaを使っているならcondaで導入可能(linux only)
conda install -c bioconda pbsv

#pbmm2がないならこちらも入れる
conda install -c bioconda pbmm2

pbsv

# pbsv 

pbsv - PacBio structural variant (SV) calling and analysis tools

 

Usage: pbsv <tool>

 

Basic Options:

  -h, --help   Output this help.

  --version    Output version information.

 

Tools:

  discover     Find structural variant signatures in read alignments (BAM to SVSIG).

  call         Call structural variants from SV signatures and assign genotypes (SVSIG to VCF).

 

Typical workflow:

  1. Align PacBio reads to a reference genome, per movie.

     Subreads BAM input:

     $ pbmm2 align ref.fa movie1.subreads.bam ref.movie1.bam --sort --median-filter --sample sample1

 

     CCS BAM input:

     $ pbmm2 align ref.fa movie1.ccs.bam ref.movie1.bam --sort --preset CCS --sample sample1

 

     CCS FASTQ input:

     $ pbmm2 align ref.fa movie1.Q20.fastq ref.movie1.bam --sort --preset CCS --sample sample1 --rg '@RG\tID:movie1'

 

  2. Discover signatures of structural variation (per movie or per sample):

     $ pbsv discover ref.movie1.bam ref.sample1.svsig.gz

     $ pbsv discover ref.movie2.bam ref.sample2.svsig.gz

 

  3. Call structural variants and assign genotypes (all samples), for CCS input append "--ccs":

     $ pbsv call ref.fa ref.sample1.svsig.gz ref.sample2.svsig.gz ref.var.vcf

root@908d5cac1b6e:/# 

pbsv discover

pbsv discover

Usage: pbsv discover [options] <ref.in.bam|xml> <ref.out.svsig.gz>

pbsv discover - Find structural variant (SV) signatures in read alignments (BAM to SVSIG).

 

Basic Options:

  -h,--help                       Output this help.

  --version                       Output version information.

  --log-file                      Log to a file, instead of stderr.

  --log-level                     Set log level: "TRACE", "DEBUG", "INFO", "WARN", "FATAL". ["WARN"]

 

Sample Options:

  -s,--sample                     Override sample name tag from BAM read group.

 

Alignment Filter Options:

  -q,--min-mapq                   Ignore alignments with mapping quality < N. [20]

  -m,--min-ref-span               Ignore alignments with reference length < N bp. ["100"]

 

Downsample Options:

  -w,--downsample-window-length   Window in which to limit coverage, in basepairs. ["10K"]

  -a,--downsample-max-alignments  Consider up to N alignments in a window; 0 means disabled. [20]

 

Discovery Options:

  -r,--region                     Limit discovery to this reference region: CHR|CHR:START-END.

  -l,--min-svsig-length           Ignore SV signatures with length < N bp. ["7"]

  -b,--tandem-repeats             Tandem repeat intervals for indel clustering.

  -k,--max-skip-split             Ignore alignment pairs separated by > N bp of a read or reference. ["100"]

 

Options:

  --emit-tool-contract            Emit tool contract.

  --resolved-tool-contract        Use args from resolved tool contract.

 

Arguments:

  ref.in.bam|xml                  Coordinate-sorted aligned reads in which to identify SV signatures.

  ref.out.svsig.gz                Structural variant signatures output.

 

pbsv call

# pbsv call    

Usage: pbsv call [options] <ref.fa|xml> <ref.in.svsig.gz|fofn> <ref.out.vcf>

pbsv call - Call structural variants from SV signatures and assign genotypes (SVSIG to VCF).

 

Basic Options:

  -h,--help                                   Output this help.

  --version                                   Output version information.

  --log-file                                  Log to a file, instead of stderr.

  --log-level                                 Set log level: "TRACE", "DEBUG", "INFO", "WARN", "FATAL". ["WARN"]

  -j,--num-threads                            Number of threads to use, 0 means autodetection. [0]

  -z,--chunk-length                           Process in chunks of N reference bp. ["1M"]

 

Variant Options:

  -t,--types                                  Call these SV types: "DEL", "INS", "INV", "DUP", "BND", "CNV". ["DEL,INS,INV,DUP,BND,CNV"]

  -m,--min-sv-length                          Ignore variants with length < N bp. ["20"]

  --min-cnv-length                            Ignore CNVs with length < N bp. ["1K"]

  --max-ins-length                            Ignore insertions with length > N bp. ["10K"]

  --max-dup-length                            Ignore duplications with length > N bp. ["100K"]

 

SV Signature Cluster Options:

  --cluster-max-length-perc-diff              Do not cluster signatures with difference in length > P%. [25]

  --cluster-max-ref-pos-diff                  Do not cluster signatures > N bp apart in reference. ["200"]

  --cluster-min-basepair-perc-id              Do not cluster signatures with basepair identity < P%. [10]

 

Consensus Options:

  -x,--max-consensus-coverage                 Limit to N reads for variant consensus. [20]

  -s,--poa-scores                             Score POA alignment with triplet match,mismatch,gap. ["1,-2,-2"]

  --min-realign-length                        Consider segments with > N length for re-alignment. ["100"]

 

Call Options:

  -A,--call-min-reads-all-samples             Ignore calls supported by < N reads total across samples. [2]

  -O,--call-min-reads-one-sample              Ignore calls supported by < N reads in every sample. [2]

  -S,--call-min-reads-per-strand-all-samples  Ignore calls supported by < N reads per strand total across samples [1]

  -P,--call-min-read-perc-one-sample          Ignore calls supported by < P% of reads in every sample. [20]

  --ccs                                       CCS optimized parameters: -A 1 -O 1 -S 0 -P 10.

 

Genotyping:

  --gt-min-reads                              Minimum supporting reads to assign a sample a non-reference genotype. [1]

 

Annotations:

  --annotations                               Annotate variants by comparing with sequences in fasta. Default annotations are ALU, L1, SVA.

  --annotation-min-perc-sim                   Annotate variant if sequence similarity > P%. [60]

 

Variant Filtering Options:

  --min-N-in-gap                              Consider >= N consecutive "N" bp as a reference gap. ["50"]

  --filter-near-reference-gap                 Flag variants < N bp from a gap as "NearReferenceGap". ["1K"]

  --filter-near-contig-end                    Flag variants < N bp from a contig end as "NearContigEnd". ["1K"]

 

Options:

  --emit-tool-contract                        Emit tool contract.

  --resolved-tool-contract                    Use args from resolved tool contract.

 

Arguments:

  ref.fa|xml                                  Reference genome assembly against which to call variants.

  ref.in.svsig.gz|fofn                        SV signatures from one or more samples.

  ref.out.vcf                                 Variant call format (VCF) output.

 

 

 

実行方法

1、リファレンスにマッピングする。pbmm2(紹介)(インストール参照)の使用が推奨されている。

subreads.bam(説明)を入力とする。

pbmm2 align ref.fa movie1.subreads.bam ref.movie1.bam --sort --median-filter --sample sample1

またはccs.bamを入力とする。

pbmm2 align ref.fa movie1.ccs.bam ref.movie1.bam --sort --preset CCS --sample sample1

またはfastqを入力とする。

pbmm2 align ref.fa movie1.Q20.fastq ref.movie1.bam --sort --preset CCS --sample sample1 --rg '@RG\tID:movie1'

unaligned.bam/fastqからマッピング結果のbamが出力される。 

 

2、pbsv discoverでSVのシグネチャを検出

sample1と2があるなら、それぞれに対して実行。タンデムリピートのbedファイルがあるなら、感度とrecall改善のためにフラグを立てることが推奨されている("--tandem-repeats repeat.bed")。

pbsv discover ref.movie1.bam ref.sample1.svsig.gz
pbsv discover ref.movie2.bam ref.sample2.svsig.gz

ref.sample1.svsig.gz、 ref.sample2.svsig.gzが出力される。

 

3、pbsv callでSV検出(joint callも可能)

sample1と2があるなら両方指定する。

pbsv call ref.fa ref.sample1.svsig.gz ref.sample2.svsig.gz ref.var.vcf

 

 ラージゲノムの場合、2以降のプロセスは並行して実行することでかなり高速化できる(*1)。実行例はpbsvのGithubを参照してください。検出アルゴリズムについてもGithubで説明されています。

 

引用

GitHub - PacificBiosciences/pbsv: pbsv - PacBio structural variant (SV) calling and analysis tools

 

関連


*1

マッピングをchromosomeごとに分けてやらないこと。

植物ゲノムアノテーションwebサービス MEGANTE

 

ハイスループットシーケンシング技術の出現により、植物ゲノムシーケンシングは加速し、そしてデータは作物改良のために利用されてきている(Bevan and Uauy 2013)。大量の植物ゲノム配列の蓄積は、比較ゲノミクスデータベースの構築(Mihara et al、2010、Nagamura et al、2011、Rouard et al、2011、Goodstein et al、2012)および効率的なデータ統合のための植物特異的なvocabulary開発をもたらした(Cooper et al、2013)。しかし、ハイスペックなコンピュータ、膨大な量のデータストレージ、およびコンピュータサイエンス分子生物学の両方における専門知識の必要性から、データ管理と分析のコストは増大している。これらのデータ解析において、ゲノムアノテーションは最も基本的かつ不可欠なステップの1つであり(Yandell and Ence 2012)、分子進化解析、トランスポゾンタギングおよびマイクロアレイ実験などのさらなる研究に直接影響を及ぼす。いくつかの解析プログラムを実施した後、結果を統合して遺伝子構造を予測し、遺伝子機能をアサインする必要があるため、アノテーション手順にはより高いレベルのバイオインフォマティクススキルが必要になる。したがって、研究者がゲノムアノテーションを実行し、結果をグラフィカルビューアで視覚化しアノテーション結果を解釈する、バイオインフォマティクスの専門知識を必要としない使いやすいアノテーションプラットフォームが不可欠である。

 現在、植物ゲノムアノテーションのために数種類の分析ツールがオンラインで利用可能である(論文執筆当時)。例えば、AUGUSTUS(Stanke and Waack 2003)、Fgenesh(Salamov and Solovyev 2000)、GeneMark.hmm(Lukashin and Borodovsky 1998)などのオンラインバージョンのab initio遺伝子予測プログラムを使用してゲノム配列からオープンリーディングフレーム(ORF)を見つけることができる。FPGP(Amano et al、22010)は、双子葉植物および単子葉植物の全長cDNA(FLcDNA)配列をクエリ配列にアライメントさせる。 Gramene(Youens-Clark et al、2011)およびPlantGDB(Duvick et al、2008)は、植物ヌクレオチドまたはタンパク質データベースに対する類似検索のWebサービスを提供している。分析結果のグラフ表示には、WebGBrowse(Podicheti et al。2009)が良い候補である。このようなWebサービスはゲノムアノテーションには役立つが、研究者が複数のWebサイトにアクセスしてその結果を1つずつ解釈するのは時間がかかる。さらに、そのようなアノテーション手順は、バイオインフォマティクススキルのない人々にとっては、入力シーケンスに対して適切なツールおよびパラメータセットを選択することが困難である。したがって、植物遺伝子のポジショナルクローニングなどのゲノム解析をサポートするには、一連の解析プログラムを自動的に実行するための統合解析ツールが必要である(Chen et al、2009、Xu et al、2011)。

 植物ゲノムシーケンスには、いくつかのWebベースのアノテーションパイプラインが利用可能である。それらのいくつかは、特定の植物ゲノムアノテーション用に設計されている。 RiceGAAS(Sakata et al。2002)はライス用で、TriAnnot(Leroy et al、2012)はコムギ用である。植物だけでなく他の種にも適応できる、より用途の広いゲノムアノテーションツールもある。 DNA subway(Goff et al、2011)は動物と植物の両方にパラメータセットを提供している。 MAKER(Cantarel et al、2008)は、分析プログラム用のリファレンスデータベースとパラメータを選択するための高度に設定可能なウェブインタフェースを持っている。ただし、既存のアノテーションパイプラインでサポートされている植物種はほとんどない。

 ここでは、クエリ配列に対していくつかの解析プログラムを実行し、結果を統合し、ゲノムブラウザ上でアノテーション情報を視覚化する、MEGANTEと呼ばれる新しい植物ゲノムアノテーションWebサービスについて説明する。既存のツールと比較して、MEGANTEの注目すべき機能の1つはそのシンプルなインターフェースである。MEGANTEはエキスパートでなくても使いやすい。さらに、このサービスは多種多様な植物種を対象としており、長さが最大10 Mbまでの大きなクエリシーケンスを受け入れる。

 

 アップロードされたシーケンスは最初にキューに入れられ、次にクエリがround-robin方式でアプリケーションサーバー上で処理され、プロセスが公平にスケジュールされる。現在のシステムでは、5つのアノテーションプロセスを並行して実行できる。アノテーションプロセス全体は、1Mbの配列については約150分以内に、そして10Mbの配列については約15時間以内に完了できる。アノテーションプロセスが終了すると、転写物およびタンパク質配列のアラインメント。予測される遺伝子構造既知のタンパク質との類似点機能ドメインGene Ontology(GO)(Ashburner et al。2000)、が報告される。すべての結果は、システムに統合されている広く使用されているゲノムブラウザ、GBrowse(Stein et al。2002)で視覚化される。さらに、システムはアノテーション結果をダウンロード用の単一のZIPファイルにアーカイブする。ファイルには、Microsoft ExcelとGFF3(http://www.sequenceontoloty.org/gff3.html)の両方の形式でアノテーション情報が含まれている。ユーザがクエリを送信するときに電子メール通知のオプションを選択した場合は、注釈の完了時に電子メールが送信されてユーザに通知される。 Webブラウザとサーバー間のデータ転送はSSL暗号化によって保護されている。

 MEGANTEは、ゲノムアノテーションのために、INSDC(Nakamura et al。2013)から得られたFLcDNAおよびEST、Swiss-ProtとTrEMBLとUniProtKBの植物部門(Magrane and UniProt Consortium 2011)、タンパク質ファミリーおよびドメインデータベース、Interpro(Hunter et al、2011)を含む複数のリファレンスデータベースを使用している。データベースは定期的に更新している。配列数などのデータベースの最新の詳細は、MEGANTE Webサイトに記載されている。論文表1に記載されている各種のFLcDNAとESTを取得した後、SeqCleanスクリプト(http://sourceforge.net/projects/seqclean/)を実行して、転写産物からpoly(A)テイル、ベクター、複雑さが少ない短い配列を削除する。

 全体的なアノテーションワークフローを論文図2に示す。アノテーションプロセスは、MIPS Repeat Element データベース(Nussbaumer et al。2013)とRepeatMasker(http://repeatmasker.org)によって検出されたRepeat Element を除外することから始まる。次に、エクソン - イントロン構造を予測するために、システムはBLATを使用して種内FLcDNAをクエリ配列にアラインメントし、カットオフ率≧98%の同一性と網羅率を得る(Kent 2002)。種内FLcDNAは正確な遺伝子予測に有効であるが、多くの場合、配列の数は遺伝子全体を網羅するのに十分ではない。このため、本著者らは、ab initio遺伝子予測にAUGUSTUS(Stanke and Waack 2003)、GeneZilla(Allen et al、2006)、Glimmer HMM(Allen et al、2006)およびSNAP(Korf、2004)、ProSplign(Sayers et al、2006)を使用する。 (

以下略)

 

ヘルプ

https://megante.dna.affrc.go.jp/help

 

f:id:kazumaxneo:20190317165912p:plain

Overview of the genome annotation workflow in MEGANTE. 論文より転載

 

対応ゲノム

f:id:kazumaxneo:20190317173329p:plain

f:id:kazumaxneo:20190317173331p:plain

 

MEGANTEに関するツイート

 

使い方

https://megante.dna.affrc.go.jp/homeにアクセスする。

f:id:kazumaxneo:20190317170838p:plain

https://megante.dna.affrc.go.jp/sigin からサインインする。初めてならばまず登録する。

右上のuploadからマルチfastaをアップロードする。

f:id:kazumaxneo:20190317171055p:plain

1配列最大10Mb、最大100配列まで同時アップロードできる。それ以上の配列があるなら、複数回繰り返す必要がある。

-----------参考-----------

10Mb以上の配列があるならseqkit(紹介) などで除く。10Mb以下の配列を出力

seqkit seq -M 10000000 scaffolds.fa > filtered.fa
  • -M, --max-len int only print sequences shorter than the maximum length (-1 for no limit) (default -1) 

f:id:kazumaxneo:20190317155754p:plain

配列はパラレルに処理される。

f:id:kazumaxneo:20190317172035p:plain

終わるとダウンロードできるようになる(↑上から3つ目が終わってる)。

 

 

表示数は左上から変更する。

f:id:kazumaxneo:20190317172125p:plain

 

ジョブが終わった配列はダウンロード / ゲノムブラウザで可視化ができる(他の配列が終わってなくても可能)。

f:id:kazumaxneo:20190317172527p:plain

 

exampleファイルのアノテーションファイルの中身をチェックする。

f:id:kazumaxneo:20190317155925p:plain

gene.xlsx

f:id:kazumaxneo:20190317155934p:plain

domain.xlsx

f:id:kazumaxneo:20190317160022p:plain

blast.xlsx

f:id:kazumaxneo:20190317160116p:plain

function.xlsx

f:id:kazumaxneo:20190317160146p:plain

 

ゲノムブラウザでのORF可視化結果。

f:id:kazumaxneo:20190317155902p:plain

 

100配列投げたら100回完了のメールが来るのでびっくりしますが、そうゆう仕様のようです。

引用

MEGANTE: A Web-Based System for Integrated Plant Genome Annotation

Hisataka Numa, Takeshi Itoh

Plant Cell Physiol. 2014 Jan; 55(1): e2.

 

ゲノムアノテーションウェブサービスMEGANTEの 果樹への応用(著者らによる使い方の説明)

https://www.naro.affrc.go.jp/publicity_report/publication/files/40f6ec588a1b579067ebd49d3b271113.pdf

 

参考

生物情報工学II  遺伝子の配列解析

https://www.agr.nagoya-u.ac.jp/~bioinfo/Ashikari/7th/2017.7-1.pdf

 

メタゲノムシーケンシングリードからCRISPRスペーサーを検出する MetaCRAST

 

 原核生物のゲノムに見られる clustered regularly interspaced short palindromic repeat (CRISPR) arraysは、我々がより多くの生態系において重要なウイルス - 微生物相互作用をよりよく理解するのを助け得る。ウイルスは溶菌感染を介して細胞内の栄養を生態系に放出し、viral shunt (Weitz & Wilhelm, 2012)と呼ばれる生態学的短絡を形成する可能性がある。このようにして、ウイルスは個々の生態系における栄養循環に貢献するだけでなく、より広い規模で地球生物化学的循環を維持することにも貢献する。 CRISPRアレイに組み込まれたウイルスDNAの短いスペーサーは過去の感染の歴史的記録を形成し、したがってウイルスを宿主に結び付ける(Sorek、Kunin&Hugenholtz、2008; Makarova、Wolf&Koonin、2013)。ウイルスの宿主特異性を決定するCRISPRスペーサーのこの能力は最近、多くの生態系からのメタゲノムを使用して探索されている(Anderson、Brazelton&Baross、2011; Sanguino et al、2015; Edwards et al、2015)。アセンブリされたゲノムにおいてCRISPRを検出するための多くのツールが存在するが(Bland et al、2007; Edgar、2007; Grissa、Vergnaud&Pourcel、2007; Rousseau et al、2009)、メタゲノムリードにおけるCRISPR検出方法はほとんど存在しない (Rho et al., 2012; Skennerton, Imelfort & Tyson, 2013; Skennerton, 2006).。

 CRISPRの反復性はそれらをメタゲノムからアセンブリすることを困難にし、アセンブリされていないリードにおいて検出する特別なツールを必要とする。リードからCRISPRアレイを検出し、そしてアセンブルするいくつかのツールが開発された。 MinCED(Mining CRISPRs in Environmental Datasets)はCRT(Bland et al、2007)の修正版であり、一方、ツールCrassは、CRT (Bland et al., 2007) とCRISPRFinder (Grissa, Vergnaud & Pourcel, 2007b)のハイブリッドアルゴリズムでスペーサーを検出する(一部略)。ロングリード(> 177 bp)では、前述のCRT戦略を使用してリピートを検索する。一方、ショートリード(<177 bp)では、適切な間隔の完全長リピート(20〜50 bp)を検索し、これらのリピートを同一のヌクレオチドでのみ拡張することで、CRTアルゴリズムによって引き起こされる潜在的なエラーを回避する。次にCrassは、単一のリピートを含むリードをさらに検索し、一致するダイレクトリピートを決定し、検出されたスペーサーの最初と最後のk-merを使用してスペーサー配置のグラフを作成し、このグラフを使用してCRISPRアレイをアセンブリする。 MinCEDとCrassはどちらもダイレクトリピート配列に関する以前の知識に頼っていないde novo検出法となっている。代わりに、これらは発見されたリピートが確かにCRISPRであるかどうかを決定するために発見的方法を使う。そのような発見的方法には、短く偽のCRISPRアレイを避けるための閾値アレイ長、およびスペーサーが繰り返しにあまりにも類似するCRISPRよりもむしろマイクロサテライトを示すかもしれないアレイを避けるための閾値リピート - スペーサ類似性が含まれる(Bland et al., 2007; Grissa, Vergnaud & Pourcel, 2007a; Skennerton, Imelfort & Tyson, 2013)。

 この論文では、Metagenomic CRISPR Reference-Aided Search Tool(MetaCRAST)、アセンブルされていないメタゲノムシーケンシングリードにおけるCRISPRスペーサー検出を改善する新しいリファレンスガイドツールを紹介する。著者らの知る限りでは、以前の研究は既知のダイレクトリピートを使用してCRISPR検出を改善したが、MetaCRASTは最初のリファレンスガイド、リード依存メタゲノムCRISPR検出ツールである。ゲノムCRISPR同定アルゴリズムCRISPRDetectは、新たに同定されたダイレクトリピートをリファレンスライブラリーと照合して、リピート境界を洗練しアレイを検証する (Biswas et al., 2016)。(一部略)MinCEDやCrassとは異なり、リファレンスガイド方式として、MetaCRASTはユーザーが指定したダイレクトリピート(DR)をメタゲノムから検索することでスペーサー検出を制限する。これらのツール間の関係および使用におけるそのような違いは、論文図1にさらに示されている。そのような特定のDRは、メタゲノムリードのアセンブリまたは分類学的プロファイリングに基づいて選択することができる。 MetaCRASTは、ユーザーがメタゲノムの分類学的構成を制御できるようにすることで、CRISPRアノテーションを改善する。それはまた、de novo検出法に必要とされる発見的方法のために起こり得る真のCRISPR拒絶を回避できる。さらに、CrassおよびMinCEDとは異なり、MetaCRASTは異なるリード長のIlluminaデータセットに対して一貫したパフォーマンスを示す。 

 

 

f:id:kazumaxneo:20190317120834p:plain

Figure 2: A comparison of per-read CRISPR detection strategies (A) between MetaCRAST and existing de novo detection tools (e.g., Crass, MinCED) and an outline of the MetaCRAST workflow (B).   論文より転載 

 

インストール

ubuntu16.04のPython 3.6.8環境でテストした(ホストOS macos10.14)。

依存

  • perl
  • fasta-splitter.pl
  • fastq-splitter.pl
  • cd-hit
  • (CPAN): Text::Levenshtein::XS, String::Approx, Getopt::Std, Bio::SeqIO, Bio::Perl, MCE, MCE::Loop, and MCE::Shared

Fasta-splitter.pl and fastq-splitter.pl are included in the repository. CD-HIT can be installed by entering sudo apt-get install cd-hit

Github

git clone https://github.com/molleraj/MetaCRAST.git
cd MetaCRAST/
sh local_install.sh

MetaCRAST -h

# MetaCRAST -h

MetaCRAST -piod [-tqhrlcan] 

 -p patterns.fasta/q 

 -i infile.fasta/q 

 -o output_dir 

 [-t] tmp_dir 

 -d dist_allowed 

 [-q] (if FASTQ file input) 

 [-h] (use Hamming Distance) 

 [-r] reverse_complement 

 [-l] max_spacer_length 

 [-c] cd_hit_similarity_threshold 

 [-a] total_spacer_cd_hit_similarity_threshold 

 [-n] num_procs 

 [Optional parameters are in brackets] 

dockerイメージも上げておきます。

docker pull kazumax/metachip 

#ホストのカレントディレクトリとイメージの/dataをシェアして起動(pullを飛ばして以下を実行してもOK)
docker run -itv $PWD:/data/ -w /root/MetaCHIP kazumax/metachip

> source ~/.bash_profile
> MetaCHIP PG -h

 

テストラン

MetaCRAST -p query/AMDquery.fa -i data/simAMDmetagenome-600-454.fasta -o test -d 3 -l 60 -c 0.9 -a 0.9
  • -p   Pattern file containing query DR sequences in FASTA or FASTQ format
  • -i    Input metagenome in FASTA or FASTQ format
  • -d   Allowed edit distance (insertions, deletions, or substitutions) for initial read detection with the Wu-Manber algorithm and subsequent DR detection steps
  • -l    Maximum spacer length in bp
  • -c    CD-HIT similarity threshold for clustering spacers detected for each query direct repeat (value from 0 to 1)
  • -a    CD-HIT similarity threshold for clustering all detected spacers (value from 0 to 1)
  • -n    Number of processors to use for parallel processing (and number of temporary metagenome parts)

出力

> head test/totalSpacersCD90.fa

# head test/totalSpacersCD90.fa

>P0S0

AAAAAAGAGTATTGTTCTGGTAAACTGTTGCACTTGC

>P0S5

AATATCTTATAGGTCTCACTGCAACCGTCAGGGAAT

>P0S6

AGAAAATTCAACGGTTTCATGAAGATGGCGAGAT

>P0S7

CGTGCCTCAATGCCAAGGAACAGATCCCTTGTGCC

>P0S9

CTGAAAAATTAAGGGATTACAAAAACCAGCTTTTAAA

 

検出されたCRISPRスペーサー数の確認

grep -c ">" test/totalSpacersCD90.fa

# grep -c ">" test/totalSpacersCD90.fa

117

 

引用

MetaCRAST: reference-guided extraction of CRISPR spacers from unassembled metagenomes
Abraham G. Moller, Chun Liang

PeerJ. 2017; 5: e3788. Published online 2017 Sep 7

 

関連


参考

CRISPR関連文献メモ_2016/07/17(6件) : crisp_bio

 

メタゲノムbinsからHGTを検出する MetaCHIP

 

 非培養微生物のゲノム再構築(ビニング)は、微生物群集DNA(メタゲノムDNA)の包括的なシーケンシングおよび新規の計算手法により最近になって実現可能になった[ref. 1-3]。再構成されたゲノムビンは、以前には特徴付けられていなかった微生物群の生化学、生理学および適応への新しい洞察を提供した[ref.4-8]。さらに、それらは、培養されていない微生物コミュニティ内の水平遺伝子伝達(HGT)を研究する機会を提供する。

 生物間の遺伝情報伝播であるHGTは、抗生物質耐性と病原性の発達を含む、微生物の進化と適応の重要な推進力であると考えられている[ref.9、10]。いくつかのバイオインフォマティクスツールが、 HGTを同定するための一連のアルゴリズムおよび特徴を用いて開発されてきた。例えば、GIST [ref.11]とIslandViewer [ref.12](紹介)はゲノム配列の構成上の特徴を利用してHGTイベントを予測するが、DarkHorse [ref.13]とHGTector [ref.14]はHGT予測に配列類似性(最良一致)を使用する。明示的な系統学的アプローチはRanger-DTL [ref.15]とAnGST [ref.16]によって採用されている。これらは遺伝子ツリーと対応するspeciesツリーとの調和を通してHGTを予測する。

 しかしながら、現在のHGT検出方法は、コミュニティ全体に適用することができず、またはリファレンスゲノムを必要とする。たとえば、HGTector [reff.14]は定義された遠位グループのメンバーから定義された自己グループメンバーまでのHGTしか検出できないため、微生物コミュニティ内のすべてのメンバーのHGTを予測する用途は限られている。培養されていない微生物には利用できないことが多いHGTを予測する。

そこで本著者らは、ここでMetaCHIP(「メタゲノム」の「Meta」、「コミュニティレベルのHGT識別パイプライン」の「CHIP」)を開発した。これはリファレンスに依存しない、コミュニティレベルのHGT識別のためのパイプラインである。シミュレートされたデータと実際のデータを分析した結果、MetaCHIPはコミュニティからのHGTを高い信頼性で検出し、新しい生物学的および生態学的洞察を得ることができることを示した。

 MetaCHIPのワークフローは論文図1に示されている。MetaCHIPはHGT検出のためにベストマッチと系統発生アプローチの両方を使用する(上記参照)。 その入力は、メタゲノムデータおよびそれらの分類学的分類に由来する一組のゲノムまたはゲノムビンの配列ファイルである。 系統学的に較正されたゲノム分類データベース(GTDB)[ref.18]に基づく最近開発されたGTDB-Tkツール[ref.17]は、入力ゲノムの分類の分類に推奨されている。 入力ゲノムは、最初にユーザー指定のランク(e.g. class, order, family or genus)でそれらの分類に従いMetaCHIPによって分類される。(以下略)

 

f:id:kazumaxneo:20190311203522p:plain

Fig.1 Workflow of MetaCHIP.  論文より転載

 

manual

https://github.com/songweizhi/MetaCHIP/blob/master/manual/MetaCHIP_User_Manual_v1.1.10.pdf

 

インストール

ubuntu16.04のminiconda3-4.3.21環境(python 3.6.1)でテストした(docker使用、ホストOS mac os10.14)。

依存

Python libraries

R packages

  • optparse: command line option parser in R.
  • ape: package for analyses of phylogenetics and evolution in R.
  • circlize: package for circular visualization.

 Third-party software

MetaCHIP makes use of the following 3rd party dependencies and assumes these are on your system path. Specify full path to their executables in the config file if they are not on the system path.

  • Prodigal: protein-coding gene prediction tool for prokaryotic genomes.
  • HMMER: tool for biosequence analysis using profile hidden Markov models.
  • MAFFT: multiple sequences alignment program.
  • FastTree: tool for inferring phylogenies from alignments .
  • BLAST+: you know what it is!
  • Ranger-DTL 2.0: software for inferring gene family evolution.

本体 Github

 

pip install MetaCHIP

進捗log

Successfully built MetaCHIP ete3 subprocess32

Installing collected packages: numpy, biopython, python-dateutil, pytz, backports.functools-lru-cache, cycler, subprocess32, kiwisolver, matplotlib, scipy, pillow, reportlab, ete3, MetaCHIP

Successfully installed MetaCHIP-1.1.10 backports.functools-lru-cache-1.5 biopython-1.73 cycler-0.10.0 ete3-3.1.1 kiwisolver-1.0.1 matplotlib-2.2.4 numpy-1.16.2 pillow-5.4.1 python-dateutil-2.8.0 pytz-2018.9 reportlab-3.5.13 scipy-1.2.1 subprocess32-3.5.3

ヘルプ

> MetaCHIP -h

$ MetaCHIP -h

 

        ...::: MetaCHIP v1.1.10 :::...

        

    HGT detection modules:

       PI      ->      Prepare Input files 

       BM      ->      Best-Match approach 

       PG      ->      PhyloGenetic approach

 

    # for command specific help

    MetaCHIP PI -h

    MetaCHIP BM -h

    MetaCHIP PG -h

 

    

kazu@b1a29c465ac7:/bin$ 

> MetaCHIP PI -h

$ MetaCHIP PI -h

usage: MetaCHIP PI [-h] -i I [-taxon TAXON] -p P [-r R] [-g G] [-x X]

                   [-grouping_only] [-nonmeta] [-noblast] [-t T] [-qsub]

                   [-force] [-quiet]

 

Prepare input files

 

optional arguments:

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

  -i I            input genome folder

  -taxon TAXON    taxonomic classification

  -p P            output prefix

  -r R            grouping rank

  -g G            grouping file

  -x X            file extension

  -grouping_only  run grouping only, deactivate Prodigal and Blastn

  -nonmeta        annotate Non-metagenome-assembled genomes (Non-MAGs)

  -noblast        not run all-vs-all blastn

  -t T            number of threads, default: 1

  -qsub           run blastn with job scripts, only for HPC users

  -force          overwrite previous results

  -quiet          not report progress

 

Example: MetaCHIP PI -h

> MetaCHIP BM -h

$ MetaCHIP BM -h

usage: MetaCHIP BM [-h] -p P [-r R] [-g G] [-cov COV] [-al AL] [-flk FLK]

                   [-ip IP] [-ei EI] [-t T] [-plot_iden] [-NoEbCheck] [-force]

                   [-quiet] [-tmp]

 

Best-match approach

 

optional arguments:

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

  -p P        output prefix

  -r R        grouping rank

  -g G        grouping file

  -cov COV    coverage cutoff, default: 75

  -al AL      alignment length cutoff, default: 200

  -flk FLK    the length of flanking sequences to plot (Kbp), default: 10

  -ip IP      identity percentile cutoff, default: 90

  -ei EI      end match identity cutoff, default: 95

  -t T        number of threads, default: 1

  -plot_iden  plot identity distribution

  -NoEbCheck  disable end break and contig match check for fast processing,

              not recommend for metagenome-assembled genomes (MAGs)

  -force      overwrite previous results

  -quiet      Do not report progress

  -tmp        keep temporary files

 

Example: MetaCHIP BM -h

MetaCHIP PG -h

$ MetaCHIP PG -h

usage: MetaCHIP PG [-h] -p P [-r R] [-g G] [-cov COV] [-al AL] [-flk FLK]

                   [-ip IP] [-ei EI] [-t T] [-force] [-quiet]

 

Phylogenetic approach

 

optional arguments:

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

  -p P        output prefix

  -r R        grouping rank

  -g G        grouping file

  -cov COV    coverage cutoff, default: 75

  -al AL      alignment length cutoff, default: 200

  -flk FLK    the length of flanking sequences to plot (Kbp), default: 10

  -ip IP      identity percentile, default: 90

  -ei EI      end match identity cutoff, default: 95

  -t T        number of threads, default: 1

  -force      overwrite previous results

  -quiet      Do not report progress

 

Example: MetaCHIP PG -h

dockerイメージも上げておきます。

docker pull kazumax/metachip 

#ホストのカレントディレクトリとイメージの/dataをシェアして起動(pullを飛ばして以下を実行してもOK)
docker run -itv $PWD:/data/ -w /root/MetaCHIP kazumax/metachip

> source ~/.profile
> MetaCHIP PG -h

 

実行方法

調べたい全てのゲノムのfastaファイルが入ったディレクトリを指定して実行する。

テストディレクトリの中身。ラン時はfastaファイルのprefixを-xで指定する(e.g., -x fasta)。

f:id:kazumaxneo:20190316214801p:plain

ランにはfastaの他にbinsにtaxonomic情報をアサインしたTSVファイルも必要になる。GTDB-Tk(github)を使うことが推奨されている。テストランのhuman_gut_bins_GTDB.tsvの中身。

f:id:kazumaxneo:20190316221427p:plain

 

 

1、MetaCHIP PI

準備

MetaCHIP PI -i human_gut_bins -x fasta -taxon human_gut_bins_GTDB.tsv -r c -p Soil -t 8

出力解説(githubより)

Output files:
1. Grouping result is exported to [prefix]_grouping_[taxon_rank][group_num].txt.
2. Gene calling results in GenBank and FASTA format.
3. A SCG protein tree of input genomes.
4. A bar plot shows the number of input genomes in each group at provided taxonomic rank.
5. Blast results

1の出力ディレクトリを-p Soilと指定すると、"Soil_MetaCHIP_wd"ディレクトリができる。以降のコマンドは"Soil"だけ指定すれば認識する。

出力ディレクト

f:id:kazumaxneo:20190316214426p:plain

f:id:kazumaxneo:20190316214659p:plain

 

2、MetaCHIP BM

ベストマッチアプローチでHGT検出

MetaCHIP BM -p Soil -r c -t 8

結果は1の出力に追加される。

f:id:kazumaxneo:20190316213529p:plain

> cat Soil_MetaCHIP_wd/Soil_c5_HGTs_ip90_al200bp_c75_ei90bp_f10kbp/Soil_c5_HGTs_BM.txt

$ cat Soil_MetaCHIP_wd/Soil_c5_HGTs_ip90_al200bp_c75_ei90bp_f10kbp/Soil_c5_HGTs_BM.txt

Gene_1 Gene_2 Gene_1_group Gene_2_group Identity end_match full_length_match

bin15_01531 bin96_00912 A C 72.587 no no

bin30_00569 bin829_00466 A E 66.221 no no

> cat Soil_MetaCHIP_wd/Soil_c5_HGTs_ip90_al200bp_c75_ei90bp_f10kbp/Soil_c5_HGTs_BM_aa.fasta

$ cat Soil_MetaCHIP_wd/Soil_c5_HGTs_ip90_al200bp_c75_ei90bp_f10kbp/Soil_c5_HGTs_BM_aa.fasta 

>bin15_01531

MSYVDEVLAYVVAKNPAQPEFHQAVKEVLESLRVVIEANEEEYRKDALLERLITPERVIM

FRVPWVDDKGNVQVNNGFRVQFNSAIGPYKGGLRFHPSVNLGIIKFLGFEQIFKNSLTGL

PIGGGKGGSDFDPKGKSDREVMAFCQSFMTELCKHIGADTDVPAGDIGVGGREIGFLYGQ

YKRIRNLYEGVLTGKGLTYGGSLARTEATGYGLLYLTEEMLKCNGKDIAGKTIAVSGSGN

VAIYAIQKAQQLGAKPVTCSDSTGWVYDPEGIDVALLKEVKEVHRARLTEYAAKRPSAEY

HDKATEGTNQWSVKVDIALPCATQNELNIDDAKALVANGVFAVAEGANMPTTLEATEYFQ

NNGILFCPGKASNAGGVATSALEMSQNSERLSWTFEEVDSKLKNIMVNIFHNLDDASKKY

GMEGNYVAGANIAGFLKVAEAMKAQGIV*

>bin30_00569

VRGFIRSVPKRRLKMSYIDEVLNRTTTRYDYQPEFCQAVTEVLKSIEPAVERNPQYQKAA

LLERLVAPEKATVFRVPWVDDNGTVHVNRGYRVQFNSAIGPYKGGLRFHPSVNMSIIKFL

GFEQTFKNSLTGLPIGGGKGGSDFDPKGKSDYEIMRFCQCFMTELYKVIGPNSDVPAGDI

GVGGREIGYLFGQYKKITGRHEGVLTGKGLSYGGSLARTEATGYGLIYLVEEMLKNHGNS

IEGKTVAVSGSGNVAIYAIEKAQMFGAKVVTASDSSGYVYDKDGIDIALLKQVKEQERAR

IVRYTELKPTAKFVPGKRVWEVPCDVALPCATQNELSLDDAKELIKNGCIAVGEGANMPS

TIDATNAFLQSKVLFAPAKAANAGGVATSALEMSQNSARMIWTFDEVDEKLKDIMESIYG

HMANAAKEYSTPDDFVAGANIAGFLKVADAMMAQGIV*

>bin829_00466

VAVPLGRLGRILVDGLHDLLQLGVDLLEGPGEPCGVLAHLEGGGGDTSGVGGLGGCEEDS

RRLVLGDGLGGGGHVCSLSDRVASVLDQDLGGLLVDLVLGRAGECDVAGDGPDAVAALGV

LGVLSEVVVEVGLYPVPLLLLDELEVPVVDTVVVLDVSVGVGDRDDLGSELGGLLAGVDG

DVSGSGDDDLLSLEGLAVGLQHLVDEVAETVSGGLGPGEGSAGADGLSGQDAGELVSEPL

VLTEHVSDLPSAGTDVSCGDVGVGTDVSEELGHEGLAEAHDLVVGLALGVEVGSSFSSSD

GEGGEAVLEDLLETEEFQDGQVDGGVQPESSLVGSDRRVELDAVSAVDLDLSVVVHPGHA

EHDDPLGFNEPLDDSVLLDLGAGLDDGLQGDEDFLDGLEELGLVCVALFQTVVDGFQVLV

VDCHHKCLKDLLTFVYCR*

>bin96_00912

MSYVDDVIELTVKQNPSEPEFHQAVKEVLESLRVVIEANEEEYKKNALLERLVNPERQLK

FRVPWVDDNGQVQVNTGYRVQFNSAIGPYKGGLRFHPSVNVGIIKFLGFEQIFKNSLTGL

AIGGGKGGSDFDPKGKSDREIMAFCQSFMTELFKYIGADTDVPAGDIGVGGREIGFLYGQ

YKRIRGLSEGVLTGKALSYGGSLARTEATGYGLLYFTDAMLKANDIDIKGKTIAVSGAGN

VAIYAIEKAQQLGGNPVTCSDSTGWIYDPEGIDVELLKEVKEVKRARLTEYAEARPSAEY

HEGKGVWSVKCDIALPCATQNELLLDDAKQLVANGVVAVAEGANMPTSIEATEYLQDNDV

LFGPGKASNAGGVATSALEMAQNSQRLSWDFDKVDKRLKVIMENIFANVDEAAKTYGFEK

NYVVGANIAGFEKVVDAMNAQGIV*

HGT候補のDNA配列も出力される。

f:id:kazumaxneo:20190316222535p:plain

 

3、MetaCHIP PG

PhyloGenetic approachでHGT検出。BMとPGを最後に比較するため2を先に実行しておく必要がある。

MetaCHIP PG -p Soil -r c -t 8

結果は1、2の出力に追加される。

f:id:kazumaxneo:20190317123201p:plain

BMのランと同様にPG候補のDNA配列と予測コード領域のアミノ酸配列も出力される。

まとめ

f:id:kazumaxneo:20190316221023p:plain

f:id:kazumaxneo:20190316221238p:plain

f:id:kazumaxneo:20190316221221p:plain

f:id:kazumaxneo:20190316221302p:plain

 

引用

MetaCHIP: community-level horizontal gene transfer identification through the combination of best-match and explicit phylogenetic tree approaches

Weizhi Song, Bernd Wemheuer, Shan Zhang, Kerrin Steensen, Torsten Thomas

Microbiome 2019 7:36