macでインフォマティクス

macでインフォマティクス

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

Small RNA-Seqのfeatureを定量する mmannot

 

 真核生物のスモールRNA(sRNA)は、長さ200bp未満の、通常は翻訳されていないRNAとして定義されている。これらのRNAは、細胞生活の多くの側面に関与していることが示されている[1, 2]。

 一般的には、特定の大きさの範囲、生合成、機能パスウェイによって分類されている。その中でも特にマイクロRNA(miRNA)が最も研究されている。成熟したmiRNAは、長い一次転写産物であり、primary transcriptsによって認識されるPre-miRNAステムループ二次構造へと折り畳まれる。他のsRNAクラスもまた、それらがよく特徴づけられた構造の一部:small nucleolar RNAs(snoRNA)、small nuclear RNAs(snRNAまたはU-RNA)、および tRNA-derived small RNAs(tsRNA)であるため、正確に同定することができる。

 逆に、他のsRNAは構造が知られていない。いくつかのクラスは、タンパク質をコードする遺伝子、または長いノンコーディングRNAとの相対的な位置関係によってのみ定義されている[3]。いくつかのsRNAは、TSSの近く、エクソン/イントロンジャンクション部の近く、遺伝子のプロモーター内、または下流領域に蓄積することが知られている。最後に、いくつかのクラス:small interfering RNAs(siRNAs)、piwi-associated RNAs(piRNAs)、およびリピートassociated siRNAs(rasiRNAs)は、非常に少数の特徴的なパターン(明確な二次構造を持たず、一次転写産物を容易に検出できない)を持っている。

 非常に多様であるが、これらのsRNAは単一の実験で捕捉することができ、このタスクを実行するためにいくつかのプロトコルが考案されている[4]。簡単に言えば、プロトコルには、RNAの精製、サイズ選択、場合によっては末端5'リン酸分子の濃縮、および配列決定が含まれている。使用する技術やマルチプレックスにもよるが、数十から数百万のリードが生成される。その後、各sRNAクラスの豊富さを評価するのが一般的である。これを行う最も簡単な方法は、リードをゲノムにマッピングし(通常は BWA [5] または bowtie [6] を使用)、マッピングされたリードをゲノム上の既知のアノテーションされた sRNA の位置を含むアノテーションファイルと比較することである。これが、各リードに対して実行されるアノテーションステップである。次に、各クラスのメンバーと共局在するリードの数をカウントすることにより、sRNAのクラスを定量化する。シンプルで広く使われている手法だが、この手法はいくつかの曖昧なケースでは機能しない。

  1. ゲノムの2つの異なる領域が同一である場合(通常はgenome duplicationの後)、1つのリードが異なる位置でも同じようにマップされることがある(論文図1A)。
  2. ゲノム内の2つの異なるアノテーション位置が部分的に重なっていて、その両方にヒット(すなわち、リードマッピング)が重なっている場合。この場合、ヒットはどちらかのアノテーションに起因している可能性がある(論文図1B)。
  3. アノテーションが重なっていないがタンデムに並んでいる場合。ヒットは通常、アノテーションの境界線上にある(論文図1C)。

 最初の曖昧さの原因は、おそらく最もよく知られていて、取り組むのが難しいものである。いくつかのsRNAは、重複領域(piRNAやsiRNAなど)と共局在したり、重複遺伝子(miRNAやtsRNA)に含まれたりすることが知られている。これらの問題を回避するために、いくつかの戦略が使用された。

  • マルチマッピングリードを廃棄する。
  • ランダムヒットを使う
  • 各ヒットに重み付けをする(リードがn回マップされた場合、各ヒットは1/nとしてカウント)。

 しかし、すべての戦略には明らかなバイアスがある。最初のストラテジーは通常はかなりの割合のリードを破棄し、2番目のストラテジーはランダムな推測を行い、最後のストラテジーは間違ったクラスを過大評価することがあるが、真のクラスは過小評価する。遺伝子がduplicateされていて、リードがduplicateされた遺伝子の両方に一意にマップされている場合、最後の2つの戦略は機能し、期待された定量化を提供することに注意する。

 原則として、MMR [7]のようなツールによって実装された第4の戦略がある。MMRは、putative lociの発現に基づいて、リードのための最適なマッピング位置を推測する。しかし、MMRはロングRNAシーケンスのために開発されたもので、リードの分布が転写物全体にわたってある程度一定であることを前提としている。これは、ショートRNA-Seqの現実の状況とはかけ離れており、特に小さな領域のみが標的とされるsRNAs silencing genesでは、真実とは程遠いものである。第二に、MMRの本当の利点は、曖昧でない領域から観察されたカバレッジを使用して、曖昧な領域にわたるリードの分布を推論することである。この戦略は、各リードが成熟したトランスクリプトの全長を占めるmiRNAに対しては単純に機能しない。

 ShortStack [8]で実装された最終的な戦略は、以前のアイデアをsRNA-Seqに適用したものである。ここでは、プロファイルの発現の推定が不規則な発現プロファイルをモデル化し、マルチマッピングリードを割り当てるために使用される。

[9]のようないくつかの論文では、根本的に異なるパイプラインを使用している。彼らはまず、異なるヌクレオチドデータベース(通常はクラスごとに1つずつ)を収集し、そのデータベースにリードをマッピングする。あるクラスのデータベースにマッピングされたリードは、そのリードがそのクラスに属していることを示唆している。しかし、このパイプラインでは、1つのリードが複数のデータベースにマッピングされる可能性があるため、複数のリードがマッピングされているという問題を解決することはできない。

 2つ目の曖昧さの原因は、複数のアノテーションが重なっている場合、2つの原因が考えられる。第一に、アノテーションに矛盾があるかもしれない。例えば、miRNAがsnoRNAと共局在することは知られていない。しかし、snoRNAがmiRNAとして検出されることがあることがわかったので、同じアノテーションデータセットの中で2つのアノテーションが競合している可能性がある。第二に、アノテーションのレベルが異なる場合がある。例えば、miRNAはイントロンの内部に局在している場合がある。ユーザーがイントロン上のmiRNAとオーバーラップするリードを見つけた場合、そのリードはおそらくイントロンではなくmiRNAに起因するものと思われる。

 最後の曖昧さの原因は、アノテーションが正しく定義されていない場合に生じる可能性がある。このような曖昧さは(長い)RNA-Seqのコンテキストでは存在するが、通常のRNA-Seqの定量ツール、例えばfeatureCounts [10]のようなものは容易に使用することができない。第一に、これらの定量化ツールは遺伝子ごとのリード数をカウントするものであり、クラスごとのリード数をカウントするものではない。第二に、イントロンや遺伝子のフランキング領域のような小さなRNAを産生する領域は、アノテーションファイルには存在しないため、これらのツールでは無視される。第三に、定量化ツールは、通常、ライブラリがstrandedであること、またはunstrandedであることを期待する。したがって、遺伝子に対してアンチセンスであるリードや、任意の方向(例えば、フランキング領域)を採用できるリードを一回の実行で定量化することは(直接的に)不可能である。最後に、マルチマッピングされたリードは通常破棄されるか(これがデフォルトの動作である)、または遺伝子ごとに一度だけカウントされる。

 これらの曖昧さを解決するために、著者らはmmannotと呼ばれるsRNAクラス定量ツールに実装された新しい戦略を提案する。ここでは紹介した5つの手法と比較した。

 

GIthubより

mmannotはscRNA seqのデータでmiRNA、rRNA、tRNAなどのシークエンスを何回行ったかを調べるツールである。リードの膨大な割合が、実際には複数の場所でマップされている可能性がある。このような複数のマップを持つリードは、通常、類似の定量ツールではうまく処理できない。mmannotは、あるリードが複数の場所にマッピングされている場合、これらのすべての場所を検査する。

 

インストール

macos10.14でテストした。

ビルド依存

Compile everything with make.

  • You will need a C++11 compiler, and zlib.

Github

git clone https://github.com/mzytnicki/mmannot.git
cd mmannot/
make -j

> ./mmannot

$ ./mmannot

Usage: mmannot [options]

Compulsory options:

-a file: annotation file in GTF format

-r file1 [file2 ...]: reads in BAM/SAM format

Main options:

-o output: output file (default: stdout)

-c config_file: configuration file (default: config.txt)

-n name1 name2...: short name for each of the reads files

-s strand: string (U, F, R, FR, RF, FF, defaut: F) (use several strand types if the library strategies differ)

-f format (SAM or BAM): format of the read files (default: guess from file extension)

-l integer: overlap type (<0: read is included, <1: % overlap, otherwise: # nt, default: -1)

-d integer: upstream region size (default: 1000)

-D integer: downstream region size (default: 1000)

-y string: quantification strategy, valid values are: default, unique, random, ratio (default: default)

-e integer: attribute a read to a feature if at least N% of the hits map to the feature (default: 100%)

Output options:

-p: print progress

-m file: print mapping statistics for each read (slow, only work with 1 input file)

-M file: print mapping statistics for each interval (slow, only work with 1 input file)

-t integer: # threads (default: 1)

-h: this help

 

 テストラン

bam、gtf、そしてconfigファイルを指定する。

cd mmannot/
mmannot -a test_dataset.gtf -r test_dataset.bam -c configHS38.txt -o output

> cat output

f:id:kazumaxneo:20200924013153p:plain

 

引用

mmannot: How to improve small–RNA annotation?
Matthias Zytnicki , Christine Gaspin
PLoS One. 2020; 15(5)

 

配列の豊富さを含むsequence indexを作る REINDEER

 

 本研究では、配列の索引付けを行い、データセットのコレクションに渡ってその豊富さを記録する新しい計算手法であるREINDEERを紹介する。これまでのところ、他の方法では、大規模なデータセットに対して効率的なインデックス付けを行うことができなかったが、本研究では、REINDEERを用いて、大規模なデータセットに対して効率的なインデックス付けを行うことができる。REINDEERを用いて、わずか56GBのRAMを使用して、2,585件のヒトRNA-seq実験から得られた配列の豊富さを45時間でインデックス化した。これにより、REINDEERは、2,585件のデータセットに渡って40億のk-mer規模の配列を記録した初めての方法となった。また、REINDEERはk-merの正確なpresence / absenceクエリもサポートしている。簡単に説明すると、REINDEERは各データセットの圧縮de Bruijnグラフ(DBG)を構築し、それらのDBGを一つのグローバルなものに統合する。そして、REINDEERはmonotigsを構築し、インデックスを作成する。

 

REINDEERは、データセット(例えば、生のRNA-seqやメタゲノムリード)のコレクションのk-merとその豊富さをインデックス化するデータ構造を構築する。その後、FASTA配列は、インデックス化された各データセットにおけるその存在と豊富さを検索することができる。他のツール(SBT、BIGSIなど)も大規模なk-merの有無検索用に設計されているが、豊富さの検索はこれまでのところサポートされていなかった(単一のデータセットでKMC、Jellyfishなどのk-merカウンタを使用して調べるケースを除いて)。REINDEERは、高速なクエリ、小さなインデックスサイズ、インデックス作成とクエリ時のメモリフットプリントの低さを兼ね備えてる。

 

インストール

ubntu18.04でテストした。

ビルド依存

  • GCC >= 4.8
  • CMAKE > 3.10.0

本体 Github

git clone --recursive https://github.com/kamimrcht/REINDEER.git
cd REINDEER
sh install.sh

./Reindeer

# ./Reindeer 

############# REINDEER version v1.0.2 #############

Command line was: ./Reindeer 

 

You must choose: either indexing (--index) or querying (--query)

 

******************* REINDEER v1.0.2**********************************

******************* REad Index for abuNDancE quERy ******************

 

                    INDEX BUILDING

 

      * Mandatory parameters

--index                 :     Indexing mode

-f <file>               :     File of file for colors. Either:

                                  i) you've already computed each DBG on your samples, in this case the fof is a list of unitig files

                              OR ii) you need Bcalm to be run to obtain unitigs per sample, in this case use --bcalm option

      * Optional parameters

-k                      :     k-mer size (default 31)

--nocount               :     Only compute presence/absence, not abundance

--bcalm                 :     Launch bcalm on each single read dataset

--paired-end            :     Index using paired-end files (provide pairs of files one after another in the fof). Works only with --bcalm.

--disk-query            :     Index for on-disk query (default: in-memory). To be used for large indexes that won't fit in RAM.

--quantization          :     Quantize the abundances in bins (to use only with --count).

--log-count             :     Record the log of the counts, gives approximate counts that save space (to use only with --count).

      * Output options

-o <file>               :     Directory to write output files (default: output_reindeer)

 

      * Advanced parameters (we recommend not to change these values unless you are very aware of REINDEER's inner components)

--minimizer-size <integer>          :    MPHF option: minimizer size

--buckets <integer>          :    MPHF option: number of buckets (log)

                    QUERY

 

      * Mandatory parameters

--query                 :     Query mode

-l                      :     Reindeer index directory (should be output_reindeer if you've not used -o during indexing)

-q <FASTA>              :     FASTA query file with query sequences

      * Optional parameters

-P                      :     Threshold: at least P% of the positions in the query must be covered by k-mers present in a dataset for the dataset to be reported (default: 40%)

-o <file>               :     Directory to write output files (default: output_reindeer/)

--disk-query            :     On-disk query (default: in-memory). To be used for large indexes that won't fit in RAM, if the index was constructed with the same option.

 

                    General

 

-t <integer>            :     Number of threads (default 1)

--help                  :     Show help

(base) root@336d73c32c08:/data/REINDEER# 

 

thread 'main' panicked at 'McFly error: "/var/folders/dy/5_xmxfts2jbdzlc1pb92j8mw0000gn/T/mcfly.XXXXXXXX.LVBfNWio" file not found: Os { code: 2, kind: NotFound, message: "No such file or directory" }', src/libcore/result.rs:1009:5

note: Run with `RUST_BACKTRACE=1` for a backtrace.

 

 

 

テストラン

1、build index

Reindeer --index -f test/fof_unitigs.txt -o quick_out

#or
sh test.sh

fof_unitigs.txtは2つのFASTAファイルのパスを記載したファイル。

出力されるindex

quick_out

f:id:kazumaxneo:20200924001512p:plain

 

2、query

クエリの配列とindexを指定する。

Reindeer --query -q test/query_test.fa -l quick_out -o quick_query

 

引用

REINDEER: efficient indexing of k-mer presence and abundance in sequencing datasets

Camille Marchet, Zamin Iqbal, Daniel Gautheret, Mikael Salson, Rayan Chikhi

bioRxiv, Posted March 30, 2020

 

long RNA sequencingリードの正確なアラインメントを行う uLTRA

 

 ロングリードRNAシークエンシング技術は、トランスクリプトームのランドスケープを研究するための主要なシークエンシング技術として急速に確立されつつある。このような解析の多くは、ゲノムに対するリードのスプライスアラインメントに依存している。しかし、ロングリード技術のエラー率やシークエンシング長は、これらのリードを正確にアラインメントするための新たな課題を生み出している。本研究では、シミュレーションデータおよび合成データを用いて、小さなexonに対して、最先端技術よりも高い精度を示したアラインメント手法uLTRAを紹介する。生物学的データを用いて、他のアライナーでは検出できないエクソン構造を持つ既知および新規のアイソフォームにuLTRAをアラインメントした例をいくつか示す。

 

インストール

依存

  • parasail
  • pysam (>= v0.11)
  • dill
  • gffutils
  • slaMEM

本体 Github

#依存 slaMEM
git clone https://github.com/fjdf/slaMEM.git
cd slaMEM
make -j
export PATH=$PATH:$PWD

conda create -n ultra python=3 pip
conda activate ultra
pip install ultra-bioinformatics

#mummer
conda install --yes -c bioconda mummer

uLTRA --help

$ uLTRA --help

usage: uLTRA [-h] [--version] {pipeline,prep_splicing,prep_seqs,align} ...

 

uLTRA -- Align and classify long transcriptomic reads based on colinear chaining algorithms to gene regions

 

positional arguments:

  {pipeline,prep_splicing,prep_seqs,align}

                        Subcommands for eaither constructing a graph, or align reads

    pipeline            Perform all in one: prepare splicing database and reference sequences and align reads.

    prep_splicing       Prepare all splicing structures from annotation

    prep_seqs           Prepare reference sequences to align to.

    align               Classify and align reads with colinear chaining to DAGs

 

optional arguments:

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

  --version             show program's version number and exit

 

 

テストラン

#容量が大きいのでclone時は注意
git clone https://github.com/ksahlin/ultra.git
cd ultra/test/

#データベースの準備とアラインの一括実行
uLTRA pipeline SIRV_genes_C_170612a.gtf SIRV_genes.fasta reads.fa outfolder/

 

引用

Accurate spliced alignment of long RNA sequencing reads

Kristoffer Sahlin, Veli Mäkinen

bioRxiv, Posted September 03, 2020

 

真菌のgenomeとtranscriptomeのデータベース Fungi.guru

 

 Fungi kingdomは真核生物の従属栄養生物で構成されており、生態系のバランスを整える役割を担い、分解者として大きな役割を果たしている。また、真菌は、抗生物質や薬理学的な性質を持つ二次代謝物を多種多様に産生している。しかし、真菌の遺伝子機能に関する知識が不足しているため、真菌を私たちのニーズに合わせて調整したり、真菌の代謝多様性を活用したりすることができない。これを解決するために、著者らは最も広く研究されている19の真菌のゲノムと遺伝子発現データを収集し、保存されたパスウェイ、機能的遺伝子モジュール、遺伝子ファミリーの種間識別のためのツールを含むデータベース、fungi.guruを構築した。このデータベースは、Aspergillus fumigatusのグリオトキシンを産生する二次代謝経路、Coprinopsis cinereaセルロースの異化経路、Fusarium graminearumPyricularia oryzaeの保存されたDNA複製経路を同定することで、様々な生物学的プロセスに関与する遺伝子の分子機能、生物学的プロセス、細胞成分を明らかにすることができることを例示している。このデータベースは www.fungi.guru で利用可能である。

 

Tutorials

https://github.molgen.mpg.de/proost/CoNekT/blob/master/docs/tutorials/overview.md

Feature

https://fungi.sbs.ntu.edu.sg/features

 

webサービス

https://fungi.sbs.ntu.edu.sg にアクセスする。

f:id:kazumaxneo:20200818013014p:plain

 

exampleにあるEAL85149(nonribosomal brevianamide peptide synthase FtmA)を検索してみる。

出力

f:id:kazumaxneo:20200818015506p:plain

 利用できるtranscriptome条件でのTPM正規化された発現プロファイルが表示される。

f:id:kazumaxneo:20200920123001p:plain

 

値はその条件での最小、平均、最大TPMで表示される。

f:id:kazumaxneo:20200920123438p:plain

 

Perturbation specificity

特異性測定値(SPM)は0から1までの範囲で、各条件/組織について計算することができる。この値は、発現プロファイル全体に対する1つの条件/組織の寄与がどれだけ大きいかを示す。CoNekTデータベースでは、各条件/組織のSPMが計算され、最も高いものだけが保存される。任意のしきい値と組み合わせてSPMを使用して、条件/組織特異的な遺伝子を見つけることができる。SPMはTiSGeD(Pubmed)で使用されている。(helpより)

f:id:kazumaxneo:20200920183740p:plain

 

Co-expression

f:id:kazumaxneo:20200920123006p:plain


 共発現パートナーは、共発現ネットワークの直接の隣人を第1レベルのneighbor、「隣人の隣人」を第2レベルのneighborとみなしている。グラフでは、デフォルトでは第2レベルのneighborが表示されるが、オプションで第1レベルのみを表示することもできる。

f:id:kazumaxneo:20200920123844p:plain

 

 他にも様々な機能がある。Toolsから選択する。

f:id:kazumaxneo:20200920184507p:plain

 

Compare Profiles

同じ種の遺伝子間で発現プロファイルを比較できる。種を選択して遺伝子を記入する。ウィンドウ下の正規化のチェックはつけることが推奨されている(デフォルトではON)。

f:id:kazumaxneo:20200920191453p:plain

 

出力

デフォルトでは、プロファイルは各プロファイルの中で最も発現している遺伝子に対して正規化される。

f:id:kazumaxneo:20200920191702p:plain

f:id:kazumaxneo:20200920191700p:plain

 

 

ダウンロードボタンを押すと、プロットのデータをタブ区切りのテキストファイルでダウンロードできる。

 

Heatmap

発現プロファイルの比較は50遺伝子に制限されているが、ヒートマップは1つの種内、または種間でより大きな遺伝子のセットを比較できる。種名と遺伝子名を選ぶ。

f:id:kazumaxneo:20200920192146p:plain

 

作成したいヒートマップのタイプに応じて、異なる正規化、または正規化なしオプションが利用できる。

 

出力

青色のセルは、遺伝子の発現が平均以下であるサンプルを示し、赤色のセルは、サンプル中で遺伝子が平均以上に発現していることを示す。濃い灰色のセルは、生の発現がゼロ(対数変換できない)を示す。白のセルは平均値を表す。(helpより)

f:id:kazumaxneo:20200920192412p:plain

 

Comparativeヒートマップでは、複数の種間の比較が可能。ただし、特定の条件セットのみを対象としている。

 

Compare Specificity

どの遺伝子が種や条件で特異的に発現していて、異なる種や条件で共通するホモログを持っているかどうかを調べることができる。

共通した発現プロファイルをもつ遺伝子を探すため、種名、条件を選択する。SPM (Specificity Metric) スライダーを調整することで感度は変更可能。

f:id:kazumaxneo:20200920202917p:plain

種間で発現が保存されている遺伝子をピックアップするだけでなく、重複後に機能が変化した遺伝子を検出することができる。

 

出力

f:id:kazumaxneo:20200920205032p:plain

 

intersectionは1つのみ。

f:id:kazumaxneo:20200920205239p:plain

OrthoFinder

f:id:kazumaxneo:20200920205815p:plain

Enriched Clusters

特定のGOに富む共発現クラスタを探す。

f:id:kazumaxneo:20200920210310p:plain

出力

f:id:kazumaxneo:20200920210332p:plain

f:id:kazumaxneo:20200920210350p:plain

 

f:id:kazumaxneo:20200920210435p:plain

 

他にBLASTサーチ機能がある。クエリの配列と相同性のあるfungi.guruの遺伝子を見つけることができる。

f:id:kazumaxneo:20200920193326p:plain

 

引用

Fungi.guru: comparative genomic and transcriptomic database for the Fungi kingdom

Jolyn Jia Jia Lim, Jace Koh, Jia Rong Moo, Erielle Marie Fajardo Villanueva, Dhira Anindya Putri, Yuen Shan Lim, Wei Song Seetoh, Sriya Mulupuri, Janice Wan Zhen Ng, Nhi Le Uyen Nguyen, Rinta Reji, Margaret Xuan Zhao, Tong Ling Chan, Edbert Edric Rodrigues, Ryan Kairon, Natasha Cassandra Chee, Ann Don Low, Zoe Chen Hui Xin, Shan Chun Lim, Vanessa Lunardi, Fong Tuck Choy, Cherlyn Xin’Er Chua, Kenny Koh Ting Sween, Jonathan Wei Xiong Ng, Marek Mutwil

bioRxiv, Posted June 30, 2020

 

Metalign

 

 サンプル中の微生物の存在と相対的な存在量を予測するメタゲノムプロファイリングは、マイクロバイオーム解析の重要な第一歩である。アラインメントベースのアプローチは、多くの場合、正確ではあるが計算が困難であると考えられている。ここでは、効率的かつ正確なアラインメントベースのメタゲノムプロファイリングを行う新しい手法、Metalignを紹介する。著者らは、新しい封じ込めミニハッシュアプローチを使用して、アライメントの前に参照データベースを事前にフィルタリングし、一意にアライメントされたリードとマルチアライメントされたリードの両方を処理して、正確なアバンダンス推定値を生成する。実際のデータセットとシミュレーションされたデータセットの両方で性能評価を行ったところ、Metalignは、すべてのデータセットにおいて高い性能と競争力のある実行時間を維持した唯一の評価手法であることがわかった。

 

方法

macos10.14でcondaの仮想環境を作ってテストした。

Github

conda create -n Metalign python=3.8 -y
conda activate Metalign
conda install -c bioconda Metalign -y

> metalign.py -h

$ metalign.py -h

usage: metalign.py [-h] [--cutoff CUTOFF] [--db_dir DB_DIR] [--dbinfo_in DBINFO_IN] [--keep_temp_files] [--input_type {fastq,fasta,AUTO}] [--length_normalize]

                   [--low_mem] [--min_abundance MIN_ABUNDANCE] [--no_quantify_unmapped] [--output OUTPUT] [--pct_id PCT_ID] [--precise] [--rank_renormalize]

                   [--read_cutoff READ_CUTOFF] [--sampleID SAMPLEID] [--sensitive] [--strain_level] [--temp_dir TEMP_DIR] [--threads THREADS] [--verbose]

                   reads data

 

Runs full metalign pipeline on input reads file(s).

 

positional arguments:

  reads                 Path to reads file.

  data                  Path to data/ directory with the files from setup_data.sh

 

optional arguments:

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

  --cutoff CUTOFF       CMash cutoff value. Default is 0.01.

  --db_dir DB_DIR       Directory with all organism files in the full database.

  --dbinfo_in DBINFO_IN

                        Location of db_info file. Default: data/db_info.txt

  --keep_temp_files     Retain KMC files after this script finishes.

  --input_type {fastq,fasta,AUTO}

                        Type of input file (fastq/fasta). Default: try to auto-determine

  --length_normalize    Normalize abundances by genome length.

  --low_mem             Run in low memory mode, with inexact multimapped processing.

  --min_abundance MIN_ABUNDANCE

                        Minimum abundance for a taxa to be included in the results. Default: 10^(-4).

  --no_quantify_unmapped

                        Do not factor in unmapped reads in abundance estimation.

  --output OUTPUT       Output abundances file. Default: abundances.tsv

  --pct_id PCT_ID       Minimum percent identity from reference to count a hit.

  --precise             Run in precise mode. Overwrites --read_cutoff and --min_abundance to 100 and 0.1.

  --rank_renormalize    Renormalize abundances to 100 pct. at each rank, e.g if an organism has a species but not genus label.

  --read_cutoff READ_CUTOFF

                        Number of reads to count an organism as present.

  --sampleID SAMPLEID   Sample ID for output. Defaults to input file name(s).

  --sensitive           Run in sensitive mode. Sets --cutoff value to 0.0.

  --strain_level        Profile strains (off by default).

  --temp_dir TEMP_DIR   Directory to write temporary files to.

  --threads THREADS     Number of compute threads for Minimap2/KMC. Default: 4

  --verbose             Print verbose output.

 

 

 テストラン

wget https://ucla.box.com/shared/static/ybz1xgke32kh56p4lqsg41t4g1bq8xy0.gz && mv ybz1xgke32kh56p4lqsg41t4g1bq8xy0.gz reads.fna.gz
metalign.py reads.fna.gz data/ --output metalign_results.tsv

テスト中

 

引用

Metalign: efficient alignment-based metagenomic profiling via containment min hash

Nathan LaPierre, Mohammed Alser, Eleazar Eskin, David Koslicki & Serghei Mangul

Genome Biology volume 21, Article number: 242 (2020)

 

NCBI Genome のBrowse by Organism機能

 

 ハイスループットシークエンシング技術の普及により、NCBIなどの塩基配列データベースに登録されるゲノム数は爆発的に増大している。BLAST検索をやり直したら少し前は無かったゲノム情報が出てきた、という話も度々耳にする。特にバクテリアアーキアはゲノム決定がしやすくなってきているため、病原性バクテリアを中心に、同じ種の様々な株が登録されるようになっている。しかし登録されるゲノムの多くは種の基準となるタイプストレインではなく表現型が不明の株だったりするため(論文の引用がないものも多い)、その情報が何に利用できるのか、という疑問は多くの研究者にあると思う。それはさておき、ここでは、NCBI塩基配列データベースに登録された細菌の莫大なゲノムの一覧を調べるのに便利な、NCBI GenomeのBrowse by Organism機能を紹介する。

 

webサービス

NCBI Genome https://www.ncbi.nlm.nih.gov/genome/に行き、Browse by Organismをクリックする(写真の左上の方のリンク)。

f:id:kazumaxneo:20200918235055p:plain

 

登録されているゲノム一覧が表示される。2020年9月17日現在、Eukaryotes 12838、 Prokaryotes 267477、Viruses 41058、Plasmids 23075、Organelles genome 16977となっている。ゲノムアセンブリのレベルはコンティグから完全長まで様々ある。

f:id:kazumaxneo:20200918235319p:plain

 

興味がある分類群を絞り込んで、どのくらいのゲノム情報が利用できるのか調べてみる。Prokaryotesをクリック。原核生物は現在267477ゲノム利用できる。

f:id:kazumaxneo:20200919002214p:plain

 

絞り込んでいく。右端のFilterをクリック。

f:id:kazumaxneo:20200919002256p:plain

 

 

様々な条件で絞り込める。

f:id:kazumaxneo:20200919002338p:plain

 

ここではまず一番登録数が多いProteobacteriaを選択。

f:id:kazumaxneo:20200919002507p:plain

 

次にガンマプロテオバクテリアを選択。

f:id:kazumaxneo:20200919002643p:plain

 

Assembly levelはCompleteを選択。この時点でまだ7303ゲノムある。

f:id:kazumaxneo:20200919002718p:plain

 

Hostはhumanを選択。

f:id:kazumaxneo:20200919002823p:plain

 

他にもRefseq reference sequence、representative sequenceだけ残すフィルタリングや、異常な配列を除外するフィルターがある。

 

得られたリストはCSVファイルとしてダウンロードできる。

f:id:kazumaxneo:20200919003520p:plain

 

もちろんキーワード検索だけで絞り込むことも可能。また、フィルター機能とキーワード検索を併用して探すこともできる。

f:id:kazumaxneo:20200919003902p:plain

大腸菌の完全長ゲノムは1065登録されていた。

 

assembly levelに注意する。⚫️マークが完全長になる。コンティグで、さらに質が非常に低いゲノムもあったりする。

f:id:kazumaxneo:20200919005454p:plain

 

ダウンロードしたCSVファイルをエクセルで開いた。

f:id:kazumaxneo:20200919004159p:plain

 

表の右端の列にはRefSeqとGenBankFTPアクセス情報が記載されており、カスタムリストの全ゲノムをダウンロードするために利用できる。

f:id:kazumaxneo:20200919004741p:plain

 

引用

Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation

Nuala A O'Leary, Mathew W Wright, J Rodney Brister, Stacy Ciufo, Diana Haddad, Rich McVeigh, Bhanu Rajput, Barbara Robbertse, Brian Smith-White, Danso Ako-Adjei, Alexander Astashyn, Azat Badretdin, Yiming Bao, Olga Blinkova, Vyacheslav Brover, Vyacheslav Chetvernin, Jinna Choi, Eric Cox, Olga Ermolaeva, Catherine M Farrell, Tamara Goldfarb, Tripti Gupta, Daniel Haft, Eneida Hatcher, Wratko Hlavina, Vinita S Joardar, Vamsi K Kodali, Wenjun Li, Donna Maglott, Patrick Masterson, Kelly M McGarvey, Michael R Murphy, Kathleen O'Neill, Shashikant Pujar, Sanjida H Rangwala, Daniel Rausch, Lillian D Riddick, Conrad Schoch, Andrei Shkeda, Susan S Storz, Hanzhen Sun, Francoise Thibaud-Nissen, Igor Tolstoy, Raymond E Tully, Anjana R Vatsan, Craig Wallin, David Webb, Wendy Wu, Melissa J Landrum, Avi Kimchi, Tatiana Tatusova, Michael DiCuccio, Paul Kitts, Terence D Murphy, Kim D Pruitt

Nucleic Acids Res. 2016 Jan 4;44(D1):D733-45

 

関連


 

 

 

 

3990xの計算機でprokka を並列ランする(ベンチマーク)

2020 9/17、9/20 文章修正

 

 3990x(64コア)のような多くのCPUコアが利用できるプロセッサを積んだ計算機を使いこなすには、高度に並列化された計算が欠かせない。しかし並列化は技術的に可能なケースと不可能なケースが存在する。並列化が捗りそうなHigh-Throughput Sequencing(HTS)関連の計算分野でも、CPUコア数増加に比例してほぼ線形の時間短縮が可能なタスクというものは少ない(*1)。しかしながら、1つのジョブでそれなりの計算時間を要するHTSの計算では、1つのジョブを走らせている間に別の処理を行うことも多い(*2)。そのため、利用可能なCPU数に余裕があることは悪いことではない(*3)。

 CPUコア数増加に応じた計算の高速化が難しいとしても、CPUリソースをより簡単かつ効率的に使い切れる方法がある。それは計算の単純並列化である。計算を行いたいファイルがたくさんある際に、計算を並列に行うということである。技術的には面白みがない単純な並列化だが、単純な並列化であるが故、並列数を増すだけほぼリニアに計算時間は短くなる。特に、処理したいファイルが何百も存在する場合はこの単純並列化の効果は大きく、物理64コアある3990xのようなCPUを使っている場合、圧倒的なランタイム短縮が期待できる。

 ここでは3390xの計算機を使って単純並列化の効果を簡単に調べてみる。原核生物のパンゲノム解析を想定し、prokkaを使ってバクテリアのゲノム配列128個にアノテーションをつけ、タンパク質配列情報を取得するまでの処理時間を調べる。この処理を選んだのは、個人的にこの処理を行うことが多く、また、Novaseqで1レーンシーケンスすれば(アセンブルとビニングにより)相当な量のドラフトゲノム情報が得られるようになっており、この処理の高速化が重要であることも背景にある。128ものゲノムもあると、計算の並列化に対応したprokkaであっても、全ファイルの逐次アノテーションには一晩近くかかる。このアノテーション付加処理を単純並列で実行し、どのくらい早くなるのかを調べる。まず手順を説明してから結果をまとめる。

 

方法

 128バクテリアゲノム(平均ゲノムサイズはおよそ4Mbほど)のアノテーションを行う。並列化にはGNU parallelを使用した。GNU parallelを使い、prokkaを並列ランする。prokkaの使用スレッドを1とし、10個を単純並列してランするなら、以下のコマンドをテキスト保存し、時間を図りながら実行する(time ./code.sh)。

#!/bin/sh 
#prooka annotation
ls *.fasta | parallel -j 10 'prokka {} --outdir prokka_{.} --addgenes --cpus 1 --quiet'
同じディレクトリに128ゲノムのfastaファイルを置いて実行。実行内容を記録する目的も兼ね、ファイルを叩いて実行している。実際に使う際はターミナルで直接叩いてしまって問題ない。

 

計測には以下の構成の計算機を使用した。

  • OS: ubuntu18.04LTS 
  • CPU: Ryzen Threadripper 3990X  64コア (合計128スレッド)
  • CPUクーラー  ROG RYUJIN 360
  • memory: DDR4 3000MHz 16GB x 8 (合計128 GB)
  • GPU: GTX 1600
  • motherboard: TRX40 Taichi
  • system: drive: Crucial M2 SSD:1TB

ゲノムのfastaファイルはSATA3接続の内臓6TB HDDに置き、そのディスク上で読み書きを行なった。

 

結果 

並列処理数を1から128まで変えて、128ゲノム全てのアノテーションがつくまでの時間を調べた結果が下の表になる。prokkaの使用スレッドは1、2、4の3つを検討した。

f:id:kazumaxneo:20200917235335p:plain

 まず左端の列がprokkaの使用スレッドを1として単純並列の数を変えた時のランタイムになる。1スレッドで128のゲノムを順番に処理するには835分(13.9時間)かかったが、並列処理を取り入れると、並列処理数が増えるにしたがってトータルのランタイムは短くなっていき、128並列処理時では21.2分で全ての計算を終えた。単純並列化の時間短縮効果は圧倒的と言える。次に、使用スレッドと並列処理数のバランスについても調べるため、prokkaの使用スレッドを2または4に増やして単純並列のタイムを調べた。それが中央と右端の列となる。ランタイムが一番短かったのは、2、4スレッドいずれも64個並列処理した時で、それぞれ19.1分と19.3分だった。もう1つ注目すべき結果は、1スレッドで4並列した時の231.2分と比較して、4スレッドで1つずつ処理した時は264.3分と、1スレッドで4並列した時の方がよりトータルのランタイムは短かったことである。プログラム内部の並列化より単純並列化した時の方が短い傾向が出ている(*4)。

 次に、CPUコアの多い3990xと他の計算機で単純並列化の効果を比較する。別のゲノム128個のアノテーションにかかる時間を、3990xの計算機とxeon E5 v4 dual(28C56T)の計算機で調べた。以下のコマンドを実行した。

#3990x
ls *.fasta | parallel -j 128 'prokka {} --outdir prokka_{.} --addgenes --cpus 1 --quiet'

#xeon E5 v4 2680 dual
ls *.fasta | parallel -j 64 'prokka {} --outdir prokka_{.} --addgenes --cpus 1 --quiet'

結果は3990xの計算機が20m.59.8sxeon E5 v4 2680 dual の計算機が43m38.9sで、3990xの計算機の方が半分以下の時間で処理を終えた。

 

まとめ

 単純並列化するだけでバクテリアゲノムのアノテーションにかかる時間を大きく短縮できることがわかった。コア数の多い3990xの計算機では特に効果が絶大で、4スレッドの逐次処理では4時間かかる計算を20分まで短縮できた。これはあくまでケーススタディだが、似た計算内容(粒度)の計算では同様の傾向が期待できると思われる。1つ気をつけたいのはメモリの使用量で(*8)、メモリ使用量の多いタスクを単純並列化すると、メモリ不足で最悪計算機が落ちてしまう危険性がある。特にTR3990xなどは256GBしか物理メモリを詰めないため、コアあたりのメモリ利用可能量は一般的なシングルノードの計算機の最大量より少ない。このこともあり、共有計算機で無闇矢鱈に単純並列化する事は勧められない(他のユーザが実行中の計算を阻害してしまうリスクも高い)。あくまで個人用のワークステーションクラス計算機で、様子を見ながら行うものだということを強調しておきます。

 

参考

 *1

64コアのCPUを単一のジョブで使い切る処理も当然存在する。この半年間3990xを使ってきたが、大きなアプリの全ビルド(Boostのビルド)、IQtree、MAFFT、CONSENTなどのエラーコレクションツールのエラー修正、Nextflowを使って実装した一部のツール、そして特にDIamond blastxのラン等は3990XのCPUリソースをかなり効率的に使う。また、Trinityのように非常に効率的にコアを使用するが、計算時間は短くならないものもあった(理由は不明)。

*2

例えば、あるユーザーが16スレッド指定でスモールゲノムのDe novoアセンブリを行なっている間に(数時間かかる)、別のユーザーが同じ計算機でRNA seqの発現解析用にfastqの16スレッド指定マッピングを3つ並行して行う、という処理が急ぎで必要になったとする。物理16コアのryzen最上位3950xでは厳しいことが予想されるが、物理64コアでメモリ帯域も多いTR3990Xなら余裕を持って行うことができる。 

*3

TRやEPYCはメモリのアクセス速度も速くなっており、日常的な複数タスクの処理を前提にCPUとその足回りを設計している節がある。

 *4

単純にprokkaをループ処理しても、128ゲノム程度なら一晩待つだけでアノテーションファイルを得ることができる。しかし、アノテーションをつけるだけで研究が終わるわけではない。20分で終わるなら当日中にかなりの事ができるようになる。

prokkaのランを使用コア数を増やしてランしても、内部で高速化が起きるのはhmmsearchによる相同タンパク質の検索などの一部のステップだけであり、リニアに処理が速くなるわけではない。

 *5

xeon E5 v4 2680は、2016年に発売されたインテルxeon v4 CPUシリーズになる。上で使ったのは、このCPUを1つのボードに2つ積んだ計算機になる。2つ積んでいるので、対等なCPUの性能比較ではないが、2680v4の値段は1つ23万円程度であり、2つあたりの値段は3990xに近い。

*6 

1桁表記は曖昧すぎるので、そろそろ単位のピコメートル表記への切り替えまたは別の指標が必要。面積になるので、6.6nmと7.4nmでは集積率に大きな差が出る。

 

*7

今後より半導体製造が微細化し、製造プロセスが1nm以下のプロセッサ(*6)や積層プロセッサのようなものが出てくれば、2030年ごろには、1ノードあたりで512-1024物理コア利用できるハイエンド環境も出てくるかもしれない(現在はintelではxeonの4-way、AMDではdual quad EPYCが市販品最大)。あるいは現在のwebサーバのように仮想化技術によって混雑度により利用可能コア数を自在に指定できるクラウドサービスがより一般的になることで、必要に応じてメニーコア環境が利用できる環境がより一般的になるかもしれない(*7)。

*8

prokkaのランは、ゲノムサイズにもよるが4スレッドで最大1GB前後のメモリーを使用する。128並列だとピーク時に100GB以上使う可能性があることになるが、実際は128並列でもメモリ使用量は30GB程度にとどまっており、100GBになることはなかった。処理するゲノムが異なるため、メモリのピーク使用量のタイミングが同調しないためと思われる。