macでインフォマティクス

macでインフォマティクス

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

変異のlandscape visualizationを提供するwebサービス LandScape

 

 シーケンシング技術の莫大な発展はゲノムデータの蓄積を加速させ、指数関数的蓄積を引き起こし、ヒトゲノム研究を著しく加速させた。一方、生物学的研究では、増加するサンプルからのトランスオミクスデータを分析している。多くの適切に設計された視覚化は、結果の結果を明確かつ完全に報告するためにpublicationで広く利用されている。最近のガンゲノム研究では、一般に数百の腫瘍ゲノムデータを分析し、mutation burden、突然変異遺伝子、生物学的パスウェイなどの一連の属性間でそれらを比較すが、これらはよく知られている‘Landscape’ figuresで一般的に示される。ランドスケープ図とその変形は、すべてのサンプルの関心のある多層特性パターンを一度に要約するのに最適であるため、サンプルの大規模コホートの研究で広く適用されている。これらの数値の実装は、十分なプログラミングスキルを必要とし、publicationのために微調整する多大な努力を必要とするため、生物学研究者にとって課題のままである。さらに、既存のソフトウェアパッケージによって生成される数値は静的であるため、インタラクティブなデータ探索が不可能になる。 
 ここでは、LandScapeovizを紹介する。これは、バッチサンプルの統合マルチレイヤーデータをインタラクティブかつリアルタイムに視覚化するためのWebベースのアプリケーションである。 LandScapeは、適切に設計されたフォーマットを使用して複数のデータタイプを処理し、表示をカスタマイズし、結果を調べ、高品質の図を生成する一連の組み込み関数を提供する。

 

 「LandScape」ビジュアライゼーションは、バッチサンプルの複数のレイヤーからの統合データの体系的な図を提供するために頻繁に利用される。これらのデータは、ガンで変異した遺伝子や生物学的パスウェイなどの特定の属性で互いに比較される。 このオンラインの「LandScape」ビジュアライゼーションは、追加のパネル(年齢、性別、組織学など)を備えた固定部分(ヒストグラムと遺伝子パネル)で設計されている。

 

Oviz-Bio Github


 

使い方

Oviz-Bio

f:id:kazumaxneo:20200116001034p:plain

Oviz-Bioは、ランドスケープビュー、SNVシグネチャ、突然変異対遺伝子、CNV、融合遺伝子、SV、および染色体3Dコンフォメーションなどの分析をサポートしている。

 

https://bio.oviz.org/demo-project にアクセスする。

f:id:kazumaxneo:20191209195353p:plain

 

 データを視覚化するには、必要な形式でCSVファイルをアップロードし、サイドバーオプションを使用して表示をカスタマイズする。 

demoのCSVをアップロードした。 

f:id:kazumaxneo:20200115235642p:plain

ダークモードOFF

f:id:kazumaxneo:20200115235646p:plain

 

入力ファイルCSVは1行目がヘッダー、2行目以降に1サンプル1行ずつ記載していく。de movo CSVexcelで開いた。

f:id:kazumaxneo:20200116013048p:plain

拡大した。1列目はsampleID。2列目以降に表示したい遺伝子を記載する。変異なしは"-"で表記する。not availableはN/Aで表記する。

f:id:kazumaxneo:20200116013400p:plain

 

この変異情報の遺伝子パネル情報だけアップロードした場合、変異の垂直トラック図と変異の種類のヒストグラムだけが出力される。

f:id:kazumaxneo:20200116014608p:plain

横方向がサンプル、ここではT119まである。縦方向が遺伝子。

 

サブタイプ、性別、などの属性情報もアップすると、上に示したように追加のヒストグラムが表示される。詳細はmanualと、マニュアルの上からダウンロードできるテストデータを確認してください。

引用

LandScape: a web application for interactive genomic summary visualization
Wenlong Jia, Hechen Li, Shiying Li, Shuaicheng Li

bioRxiv preprint first posted online Dec. 6, 2019

 

関連


 

ゲノムの指定した領域をNでマスクする

 

bedtoolsを使う。

 

Document

 

bedtoolsのインストール

本体 Github

 

#bioconda(link)
condaw install -c bioconda -y bedtools

bedtools maskfasta

$ bedtools maskfasta

 

Tool:    bedtools maskfasta (aka maskFastaFromBed)

Version: v2.29.0

Summary: Mask a fasta file based on feature coordinates.

 

Usage:   bedtools maskfasta [OPTIONS] -fi <fasta> -fo <fasta> -bed <bed/gff/vcf>

 

Options:

-fi Input FASTA file

-bed BED/GFF/VCF file of ranges to mask in -fi

-fo Output FASTA file

-soft Enforce "soft" masking.

Mask with lower-case bases, instead of masking with Ns.

-mc Replace masking character.

Use another character, instead of masking with Ns.

-fullHeader Use full fasta header.

By default, only the word before the first space or tab

is used.

 

 

実行方法

入力の FASTAとマスクしたい領域のBED/GFF/VCFを指定して実行する。

bedtools maskfasta -fi input.fa -fo mask.fa -bed regions.bed

出力

f:id:kazumaxneo:20200114200153p:plain

 

 

ソフトマスク

bedtools maskfasta -fi input.fa -fo mask.fa -bed regions.bed -soft

出力

f:id:kazumaxneo:20200114200655p:plain

 

引用

BEDTools: a flexible suite of utilities for comparing genomic features.

Quinlan AR1, Hall IM.

Bioinformatics. 2010 Mar 15;26(6):841-2. 

https://www.ncbi.nlm.nih.gov/pubmed/20110278

 

関連


ベンチマークその2(2019)

 2020 9/12,9/13 誤字修正

 

 1の続きになります。

 

ピークメモリ値

 以下のグラフは前回の投稿のxeon E5 v4 dualのランログから取った、各ツールのピークメモリ使用量(GB)になる (n=5)。flyeのピークメモリが突出していた。特にエラー修正していないraw ONT リードをアセンブリに使うと60GBを超えていた(*1)。前もってLordecでエラー修正してから使うと、ピークメモリは25GB前後まで減った。

f:id:kazumaxneo:20200113194302p:plain

前回の投稿で、32GBメモリのryzen7計算機は、raw ONT リードのアセンブリのみセグメンテーションフォルト を出して終了したと書いたが、ピークメモリの値からメモリが足りなかった事が主要な原因だとわかった。

 

補足  

ところでLoRDEC(紹介)でエラー修正することでONTのアセンブリは変化するのだろうか?アセンブリ配列をseqkitでチェックする。

$ seqkit stats -T flye_correct/assembly.fasta flye_raw/assembly.fasta 

出力を表にした。

f:id:kazumaxneo:20200113200920p:plain

de novo アセンブリ前にエラー修正することで、contig数は144から49まで減り、連続性も大きく改善していた。最大長のcontigサイズも14Mbから16Mbに増加した。QUAST-LG(紹介)の結果も一番下にまとめた(*5)。アセンブリ前にショートリードでエラー修正する価値は多いにあると言える。

 

 

まとめに入る。今回使ったマシンの長所・短所を簡単にまとめた。

1、3700x自作機

不可の低いジョブで高い性能を示したものの、多くのコアを使う計算集約的なジョブでは遅れが目立った。ECCメモリ(参考)を使えないため、長時間のジョブを走らせるには不安が残る。しかし8万円前後で組めることを考えれば費用効果は抜群に高い。2019夏に購入して連続運用しているが、今のところシャットダウンやパーツの破損などのトラブルもない。Ubuntuがかなり成熟してきたことも使い勝手の良さを後押ししている。目立った欠点は128GBしかメモリが積めないことか(ボードによっては64GB)。またメモリがデュアルチャネル、PCIの帯域なども気になるところ。

2、MousePro

簡単なジョブ、計算集約的なジョブを問わず高い性能を示した。

3、xeon platinum自作ワークステーション

本体価格が一番高い割には計算時間がかかることが多く、今回の比較ではいいところがなかった。flyeのランタイムが長くなった原因は不明。AVX-512対応は長所だが(e.g., BWA MEM2)、AVX-512はクロックが下がる問題があるので手放しで喜べない。

4、mac mini 2018上位

3と同様、特にいいところがなかった。計算集約的なジョブ、特にのlordecのランで遅れが目立ったのは、mac miniの冷却性能が低く、長時間負荷による熱の蓄積でCPUのサーマルスロットリングが起き、上限クロックまで上がらなかったためかもしれない(*4)。osxが使え、各種安定ドライバ、完成度の高い有料アプリが揃っているのが一番のメリットかもしれない。

5、Apple mac pro 2012 early (リンク

すでに引退させた機種だが、予想以上に高いパフォーマンスを示した。osのアップデートができないので今後セキュリティ面が難しくなるが、用途を限定すればまだまだ使えそうな気がした。しかし電気代の面で高コストである。主力運用は控えたい。

 

次に、今回実行したジョブの各計算機のランタイムの平均値と、計算機の価格の図を示した。計算機の発売時期は全く異なり、点の数も全然足りないのだが、参考程度に見て欲しい。図では、安くて性能が高い(費用効果が良い)計算機ほど左下に位置する。ryzen7の計算機がその位置に最も近く、2番目に近いのはxeon E5 2680 dualの計算機である。では、ryzen7 3700xが最も優秀な計算機と結論して良いのだろうか。

f:id:kazumaxneo:20200114003900p:plain

 ここで忘れてはいけないのが、計算機の堅牢性や冗長性、環境負荷(消費電力)など数値に変えにくい因子である。上の図はこれらの因子を抜きにして議論している。いくら安くても、ラン中にメモリエラーを起こしたり、データが頻繁に飛ぶような計算機では使い物にならない。よって、左下に行くほど良い計算機と簡単に判断するのは危険である。とすれば、実は散布図中央の一番下付近がスイートスポットで、この位置する計算機が好ましいのかもしれない。高い計算能力を持ち、堅牢性も備えた計算機を組み立てれば、コストが上がって自然にその辺りにプロットされるからである。今回使った計算機の中では、xeon E5 2680 dualが最もその位置に近いだろうか。

 では結局のところどのような計算機を持つのが理想的だろうか?とりあえず3990xを買って計算集約的なジョブにはこちらを使い、ルーチンの軽いジョブはryzenやcoreiシリーズ、という風に使い分けするのがバランスが良いろうか?この手の話はケースバイケースであり、十把一絡げな議論は難しいが、ここで1つ注意したいのは、ZEN2世代のRyzen Threadripperのメモリ上限が256GBしかないことである。256GBメモリでは、真核生物のアセンブリ(例えばmasurcaのメモリ使用量)、メタゲノムのアセンブリmetaFlyeの質問,  metaspades)では不安が残る。幸い現在ではたくさんのクラウドリソースが利用できるようになっていて、お金を払えば24TBメモリなどリッチクラウドはすぐに利用できる。

この文章自体5年、いや10年遅いが、今や研究室単位で巨大なクラスタを持つ時代ではなくなって来ている。研究室単位ではせいぜいRyzen Threadripper程度のマシンにとどめ大きなデータはクラウドスパコン含む)、という風に使い分けるのが賢い選択だと思う。(分野によります。主観性の強い意見なので参考程度に)。

2020 9/13 追記

OEM向けに2TBまで対応したThreadripper Proが発売・出荷されている。こちらのCPUが搭載されたマシンを使えば、メモリ容量問題はほぼ解決する(3990xよりクロックがやや低いので、TRを使うメリットがやや損なわれているのは問題だと思う)。

 

まとめ

xeon E5 v4 2680 dual > ryzen7 3700x   mac mini i7 > xeon platinum P-8136

 

 

 

 

 

 

 

 

補足

*1

flyeのメモリ使用量はGithubにまとめられている。

https://github.com/fenderglass/Flye

 *2

実はxeon platinumのマシンはその位置を狙って組んだのだが、ジョブによってかなり時間がかかる計算機になってしまった(flye参照)。 

 *3

このmac mini2018のcorei7はオプションであり、デフォルトではよりクロックの低いCPUになっている。appleは時々インテル系CPUの発熱を甘く見てやらかすので(macbook pro2018のcorei9コアモデル)、corei7の熱をうまくさばけていないのかもしれない。

 *4

学内のスパコンサービスや計算機研究者にお願いしてリッチな環境も利用させてもらえるなら解決する問題かもしれない。

 *5

連続性は向上していても、その分エラーが増えていれば意味がない。QUEST-LGで評価する。

quast flye_correct/assembly.fasta flye_raw/assembly.fasta \
-R Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz \
-1 clean_1.fq.gz -2 clean_2.fq.gz \
--nanopore ont_trimmed.fq.gz \
--labels flye_correct,flye_raw \
--eukaryote \
-o quast_output \
-t 40 --large

以下に結果の表を乗せた。エラー修正することで、genome fraction6%近く上がった。ミスマッチの塩基数は増えているものの、indelエラーが大きく減った。また、リファレンスゲノムに全くアラインされないcontigの数が41から1に減った。まだindelエラーは多いが、raconやpilonでエラー修正することでさらにrefineできると思われる。

Genome statistics flye_correct flye_raw
Genome fraction (%) 79.563 73.857
Duplication ratio 1.017 1.002
Largest alignment 3832552 3540314
Total aligned length 96645076 88413363
NG50 14486416 12278303
NG75 11422790 9738359
NA50 488799 274714
NA75 58141 -
NGA50 488799 336387
NGA75 56885 -
LG50 4 5
LG75 7 7
LA50 44 64
LA75 222 -
LGA50 44 54
LGA75 224 -
Misassemblies    
# misassemblies 856 565
   # relocations 428 354
   # translocations 415 199
   # inversions 13 12
# misassembled contigs 20 16
Misassembled contigs length 116668822 112170040
# local misassemblies 7465 4680
# scaffold gap ext. mis. 0 0
# scaffold gap loc. mis. 0 0
# possible TEs 304 166
# unaligned mis. contigs 15 56
Unaligned    
# fully unaligned contigs 1 41
Fully unaligned length 66442 1098364
# partially unaligned contigs 36 92
Partially unaligned length 22758306 36369393
Mismatches    
# mismatches 1264899 973593
# indels 381552 860169
Indels length 1335715 1752268
# mismatches per 100 kbp 1330.59 1103.27
# indels per 100 kbp 401.37 974.74
   # indels (<= 5 bp) 348521 827147
   # indels (> 5 bp) 33031 33022
# N's 200 900
# N's per 100 kbp 0.17 0.71
Statistics without reference    
# contigs 47 139
# contigs (>= 0 bp) 47 144
# contigs (>= 1000 bp) 47 141
# contigs (>= 5000 bp) 46 127
# contigs (>= 10000 bp) 43 112
# contigs (>= 25000 bp) 41 102
# contigs (>= 50000 bp) 38 72
Largest contig 16188110 14781025
Total length 119509464 125898504
Total length (>= 0 bp) 119509464 125904440
Total length (>= 1000 bp) 119509464 125902454
Total length (>= 5000 bp) 119505257 125848957
Total length (>= 10000 bp) 119481288 125734397
Total length (>= 25000 bp) 119445279 125544799
Total length (>= 50000 bp) 119331684 124385925
N50 14486416 12278303
N75 11422790 8802673
L50 4 5
L75 7 8
GC (%) 36.1 36.41
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ベンチマーク(2019)

2020 2/10 追記

2020 3/15 文章修正

2020 9/13  誤字修正

 

最近はZEN世代のCPUが高いパフォーマンスを出している。しかし、公開されている情報の多くはCPUの基礎的な特性を示すベンチマークだったり、ゲーミングや映像編集など需要が高い分野に限定されており、世間的にはニッチなバイオインフォマティクスのツールのパフォーマンスについての情報は乏しい。そこで、このブログ記事では実際にryzen7 3700xの計算機やインテル系CPUの計算機を使ってHTS関連のツールを動かし、そのランタイムから、ZEN2世代CPUがHTSのデータ処理においても有用なのか考えていく。結果が良好であれば、同じ規格のCore Complex (CCX)の組み合わせから成るEPYCやRyzen Threadripperも似た傾向が期待できる。

 

 

今回使用したハードウエア

1、3700x自作機

  •  CPU: Ryzen3700x 8コア (合計16スレッド)
  •  memory: DDR4 2667MHz 16B x 2 (合計32 GB)
  •  GPU: GTX1080
  •  system: drive: Crucial M2 SSD:1TB x1 (link)メモリキャッシュ機能は無効
  •  OS: ubuntu18.04LTS

2、MousePro(リンク

  •  CPU: xeon E5 2680 v4  2.4 GHz/14コア x 2 (合計56スレッド)
  •  memory: DDR4 ECC Registered 8GB x 16(合計128 GB)
  •  GPU: quadroK620
  •  system: drive: SATA3 SSD:1TB x2 (RAID0)
  •  OS: ubuntu18.04LTS

3、自作ワークステーション

  •  CPU: xeon E5 platinum  2.1 GHz/ 28コア x 2 (合計56スレッド)
  •  memory: DDR4 ECC Registered 64GB x 8 (合計512 GB)
  •  GPU: quadro K4000
  •  system: drive: SATA3 SSD:1TB
  •  OS: ubuntu18.04LTS

4、mac mini 2018上位

  •  CPU: core i7 8700B 6コア  (合計12スレッド)
  •  memory: DDR4 2667MHz 16GB x 2 (合計32 GB)
  •  GPU: CPU内蔵グラフィック
  •  system: drive: M.2 SSD 512 GB
  •  OS:  OSX mojave

5、Apple mac pro 2012 early (リンク

  •  CPU: xeon (Westmere EP) 3.46 GHz/6コア x 2 (合計24スレッド)
  •  memory: DDR3 1,333 MHz 8GB x 12 (合計96 GB)
  •  GPU: GTX 680
  •  system: drive: SATA3 SSD:512 GB x 3 (RAID0のストライピング)
  •  OS: OSX sierra

 

2020 2/10 追記

6、3990x自作機

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

後から3990xをベンチマークすると、再現性が疑われる結果になってしまいました。何かバージョン管理できていないライブラリがある可能性があります。追記はしないことにします。申し訳ありません。

 

 

簡単な解説

  • 1、Ryzen3700xの自作機。マザーはASUSmicro-ATXボードB450M-K。メモリはDDR4 2666MHz 16GBx2。3700xは定格で3200MHzのメモリまで対応しているが、安価な2667MHzのメモリと組み合わせている(*)。GuppyのGPU版やニューラルネットワークベースのツールを使うため、グラフィックボードはGTX1080 を選択。スモールゲノム解析用のため、メモリは32GBに留めた(ボードは最大128GBメモリまで認識)。CPUクーラーは純正を使用。
  • 2、マウスコンピュータが販売しているxeon E5 v4のワークステーション。CPUの開発コードネームはBroadwell-EP。この世代のxeon E5 2680v4を2つ搭載している。メモリは8GBのDDR4  2400MHz16枚でトータル128GB。システムディスクは1TB SATA3 SSDのストライピング(最近のQLCではなく書込み耐性の高いSSDを使用)。cinebenchR15のスコアは4000前後。電源は900W。GPUを無視すれば、今回の比較中では絶対性能が高く、そして最もバランスが取れた計算機
  • 3、Intel Xeon P-8136インテルOEM向けCPUだが、中国の知人経由で入手)をシングルで積んだワークステーション。CPUの開発コードネームはSkylake-SP。自分でパーツを選定し組み立てた。電源はgoldの1000W。メモリは64GBx8の512GBで、de novo assemblyなどのタスクで多く要求されるメモリは今回のマシンの中で圧倒的に多い。システムディスクはubuntuが512GBのSSD、windows10 proが1TB SSDでdual bootにしている。本機は大容量メモリを積んでいるにしてはシステムの容量が少ない(*2)。ボードは10Gb Ethernetポート搭載。インテルの新しめのサーバに切り替えるとどのくらい処理が早くなるのかの指標になる。(*3)
  • 4、mac mini2018のcore i7 CPUカスタム。インテルの8700Bが積まれていて6コア12スレッドで動作する。メモリは32GB。オプションの10Gb Ethernetポートに変更している。フットプリントは小さいが、背面端子の種類が豊富で、他の機器との接続がしやすい。現在主力機として使用(*4)。
  • 5、2世代前、つまり2009-2012世代のMac proで、本機は2012年モデル。以後、mac pro2012と記載する。デフォルトのCPUからX5690(3.44GHz)に換装し、メモリ容量は96GBまで増設、システムディスクをPCIスロットに繋ぎ3台でストライピングするなどの改良を施している。2009-2012世代のmac proではほぼ最速構成。しかし、10年前の設計であり、メモリ、バス(用語)などの足回りも含めてあらゆる点で現在の世代の計算機より見劣りする。ロートルな計算機と新しい計算機でどのくらい差がつくのかの参考になる。

 

 他にIvy Bridge-EPやHaswell-EPのxeonxeon platinum8180M dualなどの計算機が利用できたのですが、オーバーホール中だったり、他のジョブのラン中だっため今回のベンチマークには使えませんでした。また、私物でmac pro2009 2.22GHz x2、macbook pro2015、2017等がありましたが、今回のシロイヌナズナゲノムのデータ処理にはスペックが不十分と考えて使いませんでした。

 

ベンチマーク時のハードウエア設定と環境

部屋の温度20-25度。他のジョブは実行しない。再起動して10分後に測定開始。mac proのファンは最大回転速度に設定。自作機のファンはマザーの自動回転に設定。

 

使用したデータ

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5803254/のデータを使った。モデル植物シロイヌナズナのシーケンシングデータになる。 リード数はショートリードが16,841,951x2、ロングリードが300,071。トータル長はショートリードが84億bp、ロングリードが34億bp。変異解析やde novoアセンブリなどのタスクであれば、十分なリード数が確保できていると言える。

ファイルサイズ

f:id:kazumaxneo:20200112221734p:plain

 

始めに、mac mini2018とxeon E5 2680v4 dualの2台だけを使って、unixのコマンドやツール実行時の使用スレッド数とランタイムの関係を調べた。pigzコマンドによる圧縮と、3つのツールの実行時間を調べた。

 

 1、pigzによるfastqの圧縮にかかる時間

 fastqの圧縮/解凍はHTS関連でも頻繁に行う作業になる。並列処理に対応したpigzを使って、fastqをgzippingする時のランタイムを調べた。

 下の図が結果。横軸が使用スレッド数、縦軸がランタイム(timeでrealを測定)。使用スレッドが1-4の少ない時は、xeon E5 2680 v4 dual(以後、xeon E5 dual)よりmac mini2018(以後、mac mini)の方が早くジョブが終わっているが、使用コア数が6コア以上に増えると逆転し、xeon E5  dualの方が早くジョブが終わった。xeon E5 2680 v4 dualはスレッド数増加に従いランタイムが短くなり、1スレッドと比較して40スレッドは27.0倍早くジョブが終わった。mac miniは、1スレッドと比べ12スレッド使用時は7.5倍早くジョブが終わった。mac miniの12スレッドとxeon E5dualの40スレッドを比較すると、xeon E5 dualはmac miniの4.2倍早くジョブが終わった。

 xeon E5 2680 v4 dual計算機の利用可能な物理コア数は28で、論理コア数は56になる。興味深いことに、32スレッドから40スレッドにスレッド数を増やしたときでも87%の時間までランタイムが短縮した。pigzはハイパースレッドを上手く使えるらしい。mac miniはあえて16スレッドまで増やした時のタイムも計測した。16スレッド使用時は12スレッド使用時より遅く、12スレッドと比較して106%のランタイムになった。

f:id:kazumaxneo:20200111102840p:plain

ここには載せていないが、unpigz での解凍は、どれだけ並列化してもランタイムはあまり変わらなかった(2コア以上で80-90%前後の時間に短縮)。

 

 

2、GNU parallelによる並列化

 GNU parallelやxargsを使うと簡単に"なんちゃって"並列化ができる。そこで、ONTのsquiggle plotデータをbasecallして得た大量のfastqファイルを入力とし、GNU parallelで並列しながらnanofiltのクオリティフィルタリングラン(使用スレッド-t=1)を実行した時のランタイムを計測した。

 下の図が結果。傾向はpigzによく似ていた。並列数が4以下だとmac mini2018の方がxeon E5 dualより早く全ジョブを終え、並列数が6以上に増えると逆転、xeon E5 dualの方が早く全ジョブが終わるようになった。mac miniの12スレッドとxeon E5 dualの40スレッドを比較すると、xeon E5 dualはmac miniの1.6倍早く全ジョブが終わった(pigzの4.2倍と比べると並列化効率は低かった)。GNU parallelによる並列化は、様子をみて、10以下、できれば4程度に留めておくのが無難と思われる。リミットが低いのは"なんちゃって"並列化=、真面目な並列化をしているのが原因と思われる。増やしすぎて他のジョブと I/Oリソースを奪い合えば、全部遅くなる(同じくらいのサイズのファイルを同じタイミングで同時に処理するので、ピーク負荷のタイミングが常にほぼ同時になる)。

f:id:kazumaxneo:20200111102842p:plain

 

3、SPAdesランの使用スレッド数

pigzやGNU parallelと同様の結果になった。SPAdesのように複数の段階からなる複雑なツールだと、内部で並列化できないプロセスが多くなる。そのためアムダール則(*5)の影響が強く出るのだろうと思われる。

f:id:kazumaxneo:20200111102844p:plain

 

4、mac mini 2018を使って、minimap2によるマッピング、samtoolsのcoordinate sort、minimap2のマッピングとパイプしたsamtoolsのcoordinate sort、の3つのコマンドランタイムを計測した。データは3つとも同じデータを使用した。n=1

f:id:kazumaxneo:20200112203441p:plain

(追記)samtools sortの並列化による恩恵は少ない(オレンジ)。samtools sortは2ー4スレッド使えば飽和量かもしれない。黒線のminimap2とsamtools のパイプ処理は8スレッド付近まで早くなる傾向が見られるが、minimap2がコア数に応じて早く処理しているのを合わせて考えると(黄緑線)、黒線がコア数に応じて高速化しているのは、minimap2のマッピングだけの高速化だけに由来すると思われる。

 

ここまでのまとめ

pigzのランは論理コア数上限近くまでランタイムが短縮したが、他のジョブは物理コア数を超えた付近でランタイムは縮まらなくなった。たくさんのCPUコアを利用できる環境やツールであっても、使用スレッド数は物理コア数+1かそれ以下にとどめておくのが無難そうである(物理コア+1)。




ここからは、1−5の計算機全てを使ってベンチマークを取っていく。以下のスクリプトを1−5の計算機それぞれ5ループ回し、ランタイム(real)とピークメモリをまとめた。

#!/bin/sh 
#benchmark20191230
ref='Arabidopsis_lyrata.v.1.0.dna.toplevel.fa.gz'
shortf='ERR2173372_1.fastq.gz'
shortr='ERR2173372_2.fastq.gz'
ont='ERR2173373.fastq.gz'
thread='12'

for i in 1 2 3 4 5
do
echo fastp${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v fastp -i $shortf -I $shortr -o clean_1.fq.gz -O clean_2.fq.gz -q 30 -w $thread) 2>> run_log

echo Nanoflint${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v gunzip -c $ont |NanoFilt -q 8 --headcrop 50 -l 500 | gzip > ont_trimmed.fq.gz) 2>> run_log

echo NanoPlot${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v NanoPlot --fastq ont_trimmed.fq.gz --loglength -t $thread) 2>> run_log

echo cat${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v cat clean_1.fq.gz clean_2.fq.gz > merged.fq.gz) 2>> run_log

echo lordec${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v lordec-correct -2 merged.fq.gz -k 19 -s 3 -T $thread -i ont_trimmed.fq.gz -o ont-corrected.fasta) 2>> run_log

echo flye_raw${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v flye --nano-raw ont_trimmed.fq.gz --out-dir flye_raw --genome-size 120M --threads $thread) 2>> run_log

echo flye_correct${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v flye --nano-corr ont-corrected.fasta --out-dir flye_correct --genome-size 120M --threads $thread) 2>> run_log

echo minimap2${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v minimap2 -t $thread -ax sr flye_correct/assembly.fasta clean_1.fq.gz clean_2.fq.gz > short.sam) 2>> run_log

echo samtools_sort${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v samtools sort -@ $thread -O BAM short.sam > short.bam) 2>> run_log

echo samtools_index${i} >> run_log
time(/usr/local/Cellar/gnu-time/1.9/bin/gtime -v samtools index -@ $thread short.bam) 2>> run_log

done



 

1、fastpによる前処理

mac mini2018とryzen7 3700x自作機が一番早くジョブが終わった。mac pro2012だとこれらの計算機の2倍近い時間がかかった。xeon platinumはCPU世代が新しいにしては遅かった(xeon platinumの遅さや大きなばらつきは計測を通して最後まで続いた)。

f:id:kazumaxneo:20200111092638p:plain

                                                                                     12 threads  (n=5) (shorter is better)

 

 2、nanofiltによるONTリードのクオリティフィルタリング

fastpと傾向は同じで、mac mini2018とryzen7 3700x自作機が一番早くジョブが終わった。fastpとnanofiltの結果からは、共通してxeon CPUの計算機が遅いという傾向があるように見える。

f:id:kazumaxneo:20200112223520p:plain

                                                                                                                      12 threads  (n=5)

 

3、nanoplotによるONTリードのクオリティ分析

ryzen7 3700xが最も早くジョブが終わった。mini2018はryzen7より時間がかかっている傾向があるが、ばらつきが大き過ぎて差があるのか不透明だった。mac pro2012はエラーになった。

f:id:kazumaxneo:20200112223636p:plain
                                                                                                                          1 thread  (n=5)

これまでのまとめ

1ー3の結果をまとめると、ryzen7はfastp、nanofilt、nanoplotのような軽い計算では十分高速と言える。次に、より計算集約的なツールを走らせた時の結果を見ていく。

 

 

4、lordecによるONTリードのエラーコレクション

lordecによるONTリードのエラーコレクションを、12スレッド指定で実行した。1−3までの傾向とは異なり、最も早く計算が終わったのはxeon E5 v4 dualだった。            

f:id:kazumaxneo:20200111092647p:plain

 

 次に、同じジョブをそれぞれの計算機のCPU論理コア数の上限か、40コアに指定して実行した。mac miniは12スレッドが上限になるため、上の図のデータをそのまま使用している。

 結果を見ると、ryzen7も12スレッドから16スレッドに増やすとランタイムは縮まったが、それ以上にxeon E5 v4 dualとmac pro2012のランタイム短縮が著しく、303min ≒ 5時間でランを終えた(12スレッド時は690min)。mac mini2018は16時間以上かかるため、結果が出るまで1日待たないといけない。一方、5時間で終わるならパラメータを変えて1日に2回は回せる。計算機の性能によって大きな差が出ていると言える。興味深いことに、mac pro2012はryzen7より100分以上早く計算が終えた。

f:id:kazumaxneo:20200112142727p:plain

                                                                        

 

 5、flyeによるONTリードのde novo assembly

 flyeによるONTリードのde novo assemblyを、12スレッド指定で実行した。lordecでエラー修正したリード(下図の上半分)と、エラー修正していないリード(下図の下半分)を使った。エラー修正すると計算量が減るため、ランタイムは大きく短縮される。

 結果を見ると、最も早く計算が終わったのはxeon E5 v4 dualだった。ryzen7の計算機は、raw ONTリードのアセンブリはメモリが足らず強制終了したが(同じメモリ容量のmac miniは正常にランした)、先にエラー修正したリードを使うと、r正常にランは終了した。

f:id:kazumaxneo:20200113145013p:plain

xeon platinumは信じられないほど遅い。

 

 5のジョブをそれぞれの計算機の論理コア数の上限か、40スレッドに指定して実行した。12スレッド時と同じく、エラー修正したリード(下図の上半分)と、lordecでエラー修正したリードを使っている(下図の下半分)。xeon E5 v4 dualの計算機が最もランタイムを縮めた。これはlordecのランでスレッドを増やした時と同様の傾向になる。"エラー修正したONTリードを使うと、120-Mbゲノムのde novo アセンブリがわずか1時間で終わる"、ということを強調しておく。100ドルゲノムならぬ、植物の1-h whole genome de novo assemblyと言える。

f:id:kazumaxneo:20200113144416p:plain

  

これまでのまとめ 

 ここまでの結果から、ryzen7は軽いジョブでは十分高速だが、エラーコレクションやde novoアセンブリのような計算集約的なジョブではワークステーションクラスの計算機に大きく劣ることがわかった。これは誰もが予想する結果ではあるが、実際にその傾向が出たのは検証した甲斐があったと思う。ryzen7は一般向けの性能と消費電力のバランスが取れたCPUだが、NGSの計算機として十分使えると言ってよい(de novo assemblyなどに使うのでなければ)。

 

 次に、flyeでde novoアセンブリしたコンティグにショートリードをマッピングし、bamファイルを出力するまでのランタイムを計測する。初期はおかしなデータが出たが、まとめる時に間違ってしまっていた。やり直した結果が以下になる。

6、minimap2によるマッピング

 minimap2を使い、flyeで作成したcontig配列にショートリードをマッピングした。スレッドはmac minimac pro、ryzen7がCPUの論理コア数の上限、xeonマシンは40コアに指定している。これまでと同じように、xeon E5 v4 dualの計算機が最もランタイムが短くなった。注目すべきはxeon platinumで、これまでのジョブではスレッドを増やしても遅かったに対し、今回は相応にランタイムが短くなった(*6)。

f:id:kazumaxneo:20200113182122p:plain

                                                                         

 

7、samtoolsによるsam => bam変換とsort(samtools sort -O BAM -@ <thread>)

 6のステップで得たsamファイルを入力として、samtoolsによるsam=>bam変換とcoordinate sortに要する時間を計測した。スレッド数は、各マシンともステップ6と同じにした。マッピングと似た傾向が見られるが、xeon platinumとxeon E5 dualのタイムは逆転し、samtools sortはxeon platinumの方が早く計算を終えた。

f:id:kazumaxneo:20200113171901p:plain


2に続く。


 


関連

 

*1

CPU内部のCCX間のインターコネクトであるinfinity fabric(マーケティング的な名称で信号は2種類)のクロックがメモリのクロック(メモリコントローラのクロック)と同期するため(Infinity Fabric Dividerオンで1:2、すなわち半分にもできる)、低クロックメモリと組み合わせた場合、ZEN CPUの性能は伸びにくくなる。

https://www.corsair.com/corsairmedia/sys_master/productcontent/Ryzen3000_MemoryOverclockingGuide.pdf 

 

*2

初期はもっと大容量だったが、色々あってこうなってしまった。

 

*3

蛇足だが、cinebench R20の全CPUスコアは7650前後。core i7の7700Kが2400くらいなので、2017年の一般向けハイエンドCPUと比較して3倍以上の性能がある事になるが、ZEN2世代CPUが出た現在ではけっして高いスコアではない(12コアの3900Xが7200前後、16コアの3950Xが9300前後。2月発売の3990X(64コア)はリークで25000オーバーという数値が出ている(リンク))

 

*4

GUIツールやofficeファイルの操作、小さな計算、このブログを書くことなどにこのmac miniを使っている。4Kモニタと組み合わせているが、レスポンスは良好。他のマシンにはSSH悪いポートを開けてアクセスしている。

 

*5

http://www.hpc.cmc.osaka-u.ac.jp/wp-content/uploads/2016/07/2016-09-02_ParallelComputingGuide-for-print.pdf

 

biocondaでコンパイルされたバイナリをそのまま使っているのでAVX512が利用できる(link)とは思っていなかったが、実はONになっているのか?と考えたが、AVX512が利用できるminimap2はminimap2のレポジトリのbranchバージョンである。AVX512は使えないと思われる。

  

timeについて

ピークメモリも測定する場合(mac)

/usr/local/Cellar/gnu-time/1.9/bin/gtime -v  

 

ubuntu

apt update && apt install time してから

/usr/bin/time -v

 

さらにランタイムも測定する

time (/usr/bin/time -v echo hello)

 

logも含めて保存する

time (/usr/bin/time -v echo hello) 2>>  time

 

オルガネラゲノムのアノテーションを行うwebサービス AGORA

 

 生物学のビッグデータにとって、次世代シーケンシング(NGS)テクノロジーは注目すべき時代であり、生物学のハイスループットゲノムデータの蓄積につながっている。生物学者がさまざまな生物から高スループットのゲノムデータを取得できたとしても、ゲノムアセンブリと機能的な遺伝子アノテーションは依然として困難な問題である。この研究では、ほぼすべての真核生物におけるオルガネラゲノムの遺伝子アノテーションの新しいアプリケーションの開発に焦点を当てた。オルガネラゲノムの遺伝子アノテーションのための5つの自動化ツールが利用可能である:DOGMA(Wyman et al、2004)、GeSeq(Tillich et al、2017)、CpGAVAS(Liu et al、2012)、Mitofy(Alverson et al、2010 )およびVerdant(McKain et al、2017)。自動化された遺伝子アノテーションのためのほとんどすべてのツールは、主に陸上植物の葉緑体ゲノム用に最適化されている。ミトコンドリアゲノム用に開発されたツールは、植物用のMitofyと動物用のDOGMAの2つだけである。さらに、ほとんどすべてのツールは、リファレンスシーケンスの選択に制限されている。カスタマイズされたデータをリファレンスシーケンスとして使用したり、他の真核生物のリファレンスを変更したりすることは許可されていない。既存のプログラムの結果は、少数の遺伝子と一致していたか、原生生物の色素体およびミトコンドリアゲノムの遺伝子にアノテーションを付けるには不十分であった。

我々(著者ら)はannotator for genes of organelle from the reference sequence analysis (AGORA)というwebのオルガネラ遺伝子アノテーターを紹介する。このプログラムは、ユーザーが選択したリファレンス配列に対して BLAST検索に基づく遺伝子アノテーション用に設計された。 AGORAは、真核生物のほぼすべてのミトコンドリアおよび色素体ゲノムについてテストされた。このプログラムは、選択されたゲノムと生物のカテゴリー(つまり、植物、動物、藻類、または原生生物)に対してユーザーに最適化されている。リファレンス配列のユーザー変更により、プログラムは、真核生物の核ゲノムおよびnucleomorphゲノムまたは細菌ゲノムのすべての遺伝子にアノテーションを付けることができる。このプログラムは、遺伝子内または逆方向反復領域のエキソン-イントロ構造を持つゲノムでも、アノテーション付きの遺伝子に対して高品質の結果を生成する。このプログラムは、機能的タンパク質コーディング遺伝子、tRNAおよびrRNAを含む各遺伝子の開始および終止ヌクレオチド、リファレンスゲノムと比較した各遺伝子情報のBLAST結果、遺伝子内容の視覚化、およびORDRAWによるアレンジメント(Lohse et al 、2013、2007)を提供する。

 

Document

https://bigdata.dongguk.edu/gene_project/AGORA/menu/document.php

 

webサービス

https://bigdata.dongguk.edu/gene_project/AGORA/にアクセスする。

f:id:kazumaxneo:20200111142618p:plain

 

オルガネラゲノムのNCBI accesion IDを指定する(テスト時はエラーになった)。

f:id:kazumaxneo:20200111143126p:plain

chloroplastかmitochondriaか指定する。またGenetic codeがstandardでないなら変更する。

ローカルの配列は下のQueryからアップロードする。

f:id:kazumaxneo:20200111143129p:plain

 

テスト配列の出力

f:id:kazumaxneo:20200111143704p:plain

テスト配列を使用した際は数秒でジョブが終わり、結果が出力された。

 

引用

AGORA: organellar genome annotation from the amino acid and nucleotide references
Jaehee Jung, Jong Im Kim, Young-Sik Jeong, Gangman Yi
Bioinformatics, Volume 34, Issue 15, 01 August 2018, Pages 2661–2663

 

関連


 

 

ロングリードやショートリードのRNA seq情報をもとに転写領域をアセンブリして出力する StringTie2

2020 7/1 インストール方法追記, コマンド追記

2020 7/2  タイトル修正

2020 7/27  merge追記

2022/06/09 論文引用

2022/12/10, 12/28追記

2023/01/21 レポジトリURL修正

 

 RNAシーケンス(RNAシーケンス)データセット内の転写産物の量を測定することは、細胞の働きを理解するための強力な方法である。リードをリファレンスゲノムに合わせるだけで、遺伝子の平均発現の大まかな推定値が得られ、スプライス部位の異なる使用のヒントが得られる[ref.1]。オルタナティブスプライシングは真核生物では非常に一般的であり、推定90%のヒトマルチエクソンタンパク質コーディング遺伝子と30%の非コーディングRNA(ncRNA)遺伝子が複数のアイソフォームを持っている[ref.2、3]。アノテーション付きのヒトタンパク質をコードする遺伝子の数は、過去10年間ほぼ一定のままだが、ncRNA遺伝子およびタンパク質をコードするアイソフォームの数は増加し続けている[ref.4]。

 イルミナなどの第2世代のシーケンサーは、数億のショート(〜100 bp)RNA-seqリードを生成できる。この長さのリードは、非常に小さなエクソンの場合を除いて、通常、2エクソンを超えない。短いリードをアセンブリすることにより、完全長の転写物を再構築し、新規遺伝子と遺伝子アイソフォームを特定できる。トランスクリプトームアセンブリには、de novoとreference-guidedの2つの主なアプローチがある。 Trinity [ref.5]やOases [ref.6]などのde novoトランスクリプトームアセンブラは、リードをゲノムにアラインせずに、リード間のオーバーラップを見つけて、それらを完全なトランスクリプトにチェーンしようとする。このタスクは、相互に大きく重複する多くのアイソフォームを持つパラロガス遺伝子および転写物の存在により複雑になり、その結果、このアプローチは非常に断片化され、エラーを起こしやすいトランスクリプトームを生成する。 Cufflinks [ref.7]、Bayesembler [ref.8]、StringTie [ref.9]、TransComb [ref.10]、Scallop [ref.11]などのリファレンスガイド付きアセンブラは、スプライシングを使用してRNA-seqリードがHISAT [12]やSTAR [13]などのアライナーでアラインされる既存のゲノムを利用する。これらのアセンブラは、アラインメントに基づいてスプライスグラフ(または他のデータ構造)を構築し、それらのグラフを使用して個々のトランスクリプトを構築する。一部のリファレンスガイド付きアセンブラは、オプションのガイドとして既知の転写産物のエクソン-イントロンアノテーションを使用することもでき、可能な場合は既知の遺伝子を優先させることができる。最近の研究[ref.14]では、より多くの転写産物を正確かつ高精度にアセンブルすることにより、StringTieがCufflinksとBayesemblerの両方よりも優れていることがわかった。

StringTieおよびその他のトランスクリプトームアセンブラーは、各トランスクリプトに割り当てられたアライメントされたリード数に基づいて、トランスクリプトの存在量を推定する。最近では、Sailfish [ref.15]、Salmon [ref.16]、Kallisto [ref.17]などの代替方法により、正確なk-merマッチングに基づいて既知の転写産物にリードをアサインすることで豊富さを推定できることが示された。しかし、これらのアラインメントフリーの方法は、新規遺伝子またはアイソフォームを検出することができず、アラインメントベースのパイプラインと比較して、低存在量および低分子量RNA定量化のパフォーマンスが低下する[ref.18]。

 StringTieの最初のリリースでは、全ゲノムアセンブリ用に開発されたスーパーリードの構築を介して、de novoトランスクリプトームアセンブリの限定バージョンを使用する方法を提案した[ref.19]。概念的に、スーパーリードは、k-merルックアップテーブルに基づいた一意の拡張がある限り、ショートリードの各端を拡張することによって構築される。これにより、ショートリードのエラー率が低い合成のロングリードコレクションが作成される。それらはより長いため、ゲノムに一意にアラインする可能性が高く、遺伝子のスプライスグラフが単純化される可能性がある。スーパーリードは、StringTie 1.0(以降StringTie1)の限られた容量で使用され、ペアエンドリード間のギャップを埋めるだけであった。その限られた実装では、リードのペアを置き換えるためにスーパーリードが使用され、単一のペアでないリードのように処理できるようになった。スーパーリードを使用する際の難しさの1つは、ゲノムアセンブリ用にスーパーリードを作成するために使用されるアルゴリズムにエラー修正ステップが含まれていることである。もう1つの問題は、完全なスーパーリードに多数のショートリードが含まれている可能性があるため、定量化ステップ中に単一のリードとしてカウントできないことである。そのため、スーパーリード間でリード範囲を分散するために、期待値最大化(EM)アルゴリズムを開発した。 

 第二世代のシーケンサーは非常に多くのリードを生成するが、それらのリード長は通常非常に短く、ほとんどのRNA-seq実験では75〜125 bpの範囲である。これらの短いリードは、多くの場合、複数の場所にアラインし、このようなリードを「マルチマッピング」と指定する。ショートリードには、2つ以上のエクソンにまたがることがめったにないという制限もある。これらの問題は、Pacific Biosciences(PacBio)やOxford Nanopore Technologies(ONT)などの第3世代のシーケンシングテクノロジーによって軽減できる。 10,000µbpを超えるリード長を生成できるこれらのロングリードテクノロジーは、全ゲノムアセンブリを劇的に改善し[ref.20]、RNAシーケンス実験に使用すると、アイソフォーム同定、そして発見の精度が大幅に向上する可能性がある[ref.21,22,23]。第3世代のシーケンサーによって生成されるリードの一部はRNA転写産物の全長をカバーするが、多くは必然的に部分的な転写産物のみをキャプチャする。これはさまざまな理由で発生する。たとえば、(1)RNAは急速に分解され、シーケンシングのためにキャプチャされるまでに全長より短くなることがある。 (2)長い分子は、ライブラリの準備中に壊れる可能性がある。または(3)cDNAシーケンスでは、逆転写ステップで完全なRNA分子を捕捉できない場合がある。したがって、トランスクリプトを完全にカバーするリードのみを考慮する計算ツールは、多くのリードを破棄せざるを得ず、感度が大幅に低下する可能性がある。しかし、これまでは、ロングリードはトランスクリプトームアセンブリに広く採用されていなかった。これは、部分的にエラー率がはるかに高く(通常8〜10%以上)、アライメントが困難であるためである[ref.24、25]。またシーケンサースループットがはるかに低いため、最高発現の遺伝子以外のすべての正確な定量は不可能である。

(一部略)

ここでは、StringTie2トランスクリプトアセンブラの主要な新しいリリースであるStringTie2を紹介する。これは、ショートリードとロングリードの両方、および完全な長さのスーパーリードをアセンブリできる。 33のIllumina RNA-seqデータセットに関する結果は、StringTie2がより正確であることを示している。

 

 

インストール

condaを使いubuntu18.04とmacos10.14でテストした(anaconda3.7環境)。

依存

  • StringTie is compatible with a wide range of Linux and Apple OS systems

本体 Github

#get latest release
git clone https://github.com/mpertea/stringtie2
cd stringtie2
make release

#bioconda (link)
mamba install -c bioconda -y stringtie

> stringtie

$ stringtie 

StringTie v2.1.2 usage:

 stringtie <input.bam ..> [-G <guide_gff>] [-l <label>] [-o <out_gtf>] [-p <cpus>]

  [-v] [-a <min_anchor_len>] [-m <min_tlen>] [-j <min_anchor_cov>] [-f <min_iso>]

  [-C <coverage_file_name>] [-c <min_bundle_cov>] [-g <bdist>] [-u] [-L] [-e]

  [--ptf <f_tab>] [-x <seqid,..>] [-A <gene_abund.out>] [-h] {-B | -b <dir_path>} 

Assemble RNA-Seq alignments into potential transcripts.

 Options:

 --version : print just the version at stdout and exit

 --conservative : conservative transcriptome assembly, same as -t -c 1.5 -f 0.05

 --rf assume stranded library fr-firststrand

 --fr assume stranded library fr-secondstrand

 -G reference annotation to use for guiding the assembly process (GTF/GFF3)

 --ptf load point-features from a given 4 column feature file <f_tab>

 -o output path/file name for the assembled transcripts GTF (default: stdout)

 -l name prefix for output transcripts (default: STRG)

 -f minimum isoform fraction (default: 0.01)

 -L use long reads settings (default:false)

 -R if long reads are provided, just clean and collapse the reads but do not assemble

 -m minimum assembled transcript length (default: 200)

 -a minimum anchor length for junctions (default: 10)

 -j minimum junction coverage (default: 1)

 -t disable trimming of predicted transcripts based on coverage

    (default: coverage trimming is enabled)

 -c minimum reads per bp coverage to consider for multi-exon transcript

    (default: 1)

 -s minimum reads per bp coverage to consider for single-exon transcript

    (default: 4.75)

 -v verbose (log bundle processing details)

 -g maximum gap allowed between read mappings (default: 50)

 -M fraction of bundle allowed to be covered by multi-hit reads (default:1)

 -p number of threads (CPUs) to use (default: 1)

 -A gene abundance estimation output file

 -B enable output of Ballgown table files which will be created in the

    same directory as the output GTF (requires -G, -o recommended)

 -b enable output of Ballgown table files but these files will be 

    created under the directory path given as <dir_path>

 -e only estimate the abundance of given reference transcripts (requires -G)

 -x do not assemble any transcripts on the given reference sequence(s)

 -u no multi-mapping correction (default: correction enabled)

 -h print this usage message and exit

 

Transcript merge usage mode: 

  stringtie --merge [Options] { gtf_list | strg1.gtf ...}

With this option StringTie will assemble transcripts from multiple

input files generating a unified non-redundant set of isoforms. In this mode

the following options are available:

  -G <guide_gff>   reference annotation to include in the merging (GTF/GFF3)

  -o <out_gtf>     output file name for the merged transcripts GTF

                    (default: stdout)

  -m <min_len>     minimum input transcript length to include in the merge

                    (default: 50)

  -c <min_cov>     minimum input transcript coverage to include in the merge

                    (default: 0)

  -F <min_fpkm>    minimum input transcript FPKM to include in the merge

                    (default: 1.0)

  -T <min_tpm>     minimum input transcript TPM to include in the merge

                    (default: 1.0)

  -f <min_iso>     minimum isoform fraction (default: 0.01)

  -g <gap_len>     gap between transcripts to merge together (default: 250)

  -i               keep merged transcripts with retained introns; by default

                   these are not kept unless there is strong evidence for them

  -l <label>       name prefix for output transcripts (default: MSTRG)

 

Error: no input file provided!

 

 

実行方法

ショートリードかロングリードをリファレンスゲノムにマッピングして得たbamを指定する。複数指定も可能。

stringtie -p 12 -o out.gtf sample1.bam sample2.bam sample3.bam 2> log
  • -o    output path/file name for the assembled transcripts GTF (default: stdout)
  • -p    number of threads (CPUs) to use (default: 1)

  

bamとリファレンスのgtf(ガイド)を指定する。

stringtie -L -G human-chr19_P.gff -p 12 -o long_reads_guided.out.gtf long_reads.bam 2> log
  • -G   reference annotation to use for guiding the assembly process (GTF/GFF3)
  • -L   use long reads settings (default:false)

 

複数のgtfをマージするオプションもあります。

stringtie --merge in1.gtf in2.gtf in3.gtf ... -o output.gtf 2> log

#reference gtf/gff3も指定する(*1)
stringtie --merge in1.gtf in2.gtf in3.gtf ... -G ref.gff3 -o output.gtf 2> log
  • stringtie --merge [Options] { gtf_list | strg1.gtf ...}

 

追記

こちらの論文ではPacBioのLong RNA sequencinng単独をStringTie2で解析したり、illuminaのRNA seqのbamとPacBioのLong RNA sequencinngからのbamをマージしてStringTie2で解析するのワークフローよりも、StringTie2の新しいmixオプション(2022年の論文#2)のオプションを使う事(PacBioとilluminaそれぞれのbamを個別に指定する)を推奨しています。詳しくはこちらExample 3: Using StringTie --mixの部分)。

 

引用

#1

Transcriptome assembly from long-read RNA-seq alignments with StringTie2

Sam Kovaka, Aleksey V. Zimin, Geo M. Pertea, Roham Razaghi, Steven L. Salzberg,  Mihaela Pertea
Genome Biology volume 20, Article number: 278 (2019)

 

#2

2022/06/09

”ハイブリッドリードアセンブリを実行するStringTieの新リリースを紹介する。ロングリードとショートリードの両方の長所を利用することで、StringTieによるハイブリッドリードアセンブリは、ロングリードのみ、またはショートリードのみのアセンブリよりも精度が高く、データセットによっては、正しくアセンブリされた転写産物の数を2倍以上に増やすことができ、ロングリードデータのアセンブリのみよりも大幅に高い精度を得ることができる”

Improved transcriptome assembly using a hybrid of long and short reads with StringTie
Alaina Shumate,Brandon Wong,Geo Pertea,Mihaela Pertea 

https://doi.org/10.1371/journal.pcbi.1009730

Published: June 1, 2022

関連


 

 

 

*1

テストした限り、基本的にリファレンスの注釈は維持され、そこにはない注釈が追加される。リファレンスとオーバーラップしていて微妙に長さが異なる部位などはリファレンスが優先される。

真核生物のゲノムプロジェクトにおいて共同研究者と共にアノテーションを効率的に進めるためのwebサービス GenSAS

2020 1/9 タイトル修正

2020 7/19 追記

2020 7/23 追記

 

 Genome Sequence Annotation Server(GenSAS、https://www.gensas.org)は、構造的および機能的アノテーション、および手動キュレーションのための安全なWebベースのゲノムアノテーションプラットフォームである。 GenSASは、ユーザーによるインストールを必要とせず、一般的なコマンドラインベースのアノテーションツールを単一の使いやすいオンラインインターフェイスに統合する。 GenSASはJBrowseとApolloを統合しているため、ユーザーはアノテーションデータを表示し、遺伝子モデルを手動でキュレートできる。埋め込まれた指示とより詳細なGenSASユーザーガイドにより、はアノテーションプロセスを段階的に、ユーザーにガイドする。ゲノムアセンブリファイルに加えて、ユーザーはアノテーションプロセスで使用するために、生物固有の転写産物、タンパク質、およびRNAシーケンスのエビデンスをアップロードすることもできる。 NCBI RefSeq転写産物およびタンパク質データベースの最新バージョンと、SwissProtおよびTrEMBLタンパク質データベースがすべてのユーザーに提供されている。 GenSASプロジェクトを他のGenSASユーザーと共有して、共同アノテーションを有効にすることができる。アノテーションが完了すると、GenSASはアノテーション付き遺伝子モデルの最終ファイルを一般的なファイル形式で生成し、他のアノテーションツールで使用したり、リポジトリに投稿したり、publicationsで使用したりできる。

 

Available Tools

https://www.gensas.org/tools

 

GenSAS tutorial Jan 2015

 

注意

GenSASはchromosomeレベルのアセンブリを期待しているため、ショートリードから得たcontig配列などでは動作しません(長い配列だけ取り出せば可能)。注意して下さい。

 

アノテーションreadyなgenome配列であるかどうかがもっとも重要になります。

https://f1000research.com/articles/7-148あたりを読んで、十分な品質のデータになっているか確認して下さい。

 

webサービス

https://www.gensas.orgにアクセスする。

f:id:kazumaxneo:20200109015710p:plain

初回はユーザー登録が必要。アカウント申請してから、連絡が来るでしばらくかかる。アカウント申請時はどのようなゲノムプロジェクトなのか、種名などを記載する必要がある。

 

ログインしてプロジェクトページにアクセスする。

f:id:kazumaxneo:20200109015639p:plain

 左から右のタブに順番に進めるようになっている。

  

配列を決定したゲノム配列をアップロードする。Sequecneタブを選択。

f:id:kazumaxneo:20200109120439p:plain

配列のタイプを選択する。

f:id:kazumaxneo:20200109120444p:plain

 

f:id:kazumaxneo:20200109120819p:plain

アセンブリバージョンも指定する。

f:id:kazumaxneo:20200109121252p:plain

 

配列をアップロードし終えたら、processingジョブがスタートする。終わるまでしばらく時間がかかる。

f:id:kazumaxneo:20200109121924p:plain

完了した。

 

projectタブに進む。Begin a new projectをクリックし、プロジェクトの詳細を記載していく。

f:id:kazumaxneo:20200109121837p:plain

 

GFF3 

すでに予測済みの遺伝子情報、RNA seqから取得した転写領域情報があれば、GFF3形式でアップロードする。それ以外の例えばrepetitive regionのGFF3などがあれば、それらもアップロードできる。

f:id:kazumaxneo:20200109233427p:plain

 

Evidence

自身の種のESTや完全長cDNAデータがあったり、closely relatedな種のタンパク質情報があるなら指定する。アラインして使用される。

f:id:kazumaxneo:20200109233113p:plain

 

Rpeats

複雑性の低い配列や反復配列をNでマスクする。

f:id:kazumaxneo:20200109233123p:plain

パラメータを指定して実行する。

f:id:kazumaxneo:20200109233135p:plain

Repeat options

f:id:kazumaxneo:20200109233140p:plain

Masking

Repeat Maskerのジョブが終わったら、このタブでマスクを実行する。

f:id:kazumaxneo:20200109233155p:plain

 

Align

既存のデータベースのタンパク質配列やRNA seqデータアラインメントする。

f:id:kazumaxneo:20200109234748p:plain

BLAST、BLAT、PASAを使ってNCBI nrのcDNAをアラインする。

f:id:kazumaxneo:20200109234755p:plain

f:id:kazumaxneo:20200109234836p:plain

f:id:kazumaxneo:20200109235218p:plain

またはTophatやHISAT2を使ってRNA seqデータをマッピングする。

f:id:kazumaxneo:20200109233226p:plain

HISAT2はSRAのデータを直接使用することもできる。

f:id:kazumaxneo:20200109235318p:plain

 

Structural

ab initioの遺伝子予測。全てのプログラムを走らせる。

f:id:kazumaxneo:20200110001126p:plain

 

Consensus

予測結果について、EvidenceModelerを使ってコンセンサスセットを作成する。

f:id:kazumaxneo:20200110001226p:plain


OGS

f:id:kazumaxneo:20200110002207p:plain

 

ラン中のjobや終わったjobは右端にまとめられる。また登録したメールアドレスにjob完了のメールが届くようになっている。

f:id:kazumaxneo:20200110104024p:plain

クリックすると別のタブで開かれる。
f:id:kazumaxneo:20200110104209p:plain

この右のジョブリンクをクリックすることで、ジョブ結果を見たり、各ファイルをダウンロードできる。

 

また、そのジョブそのものを消すことができる。クリックして開いてDelete this jobsを選択。

f:id:kazumaxneo:20200719194223p:plain

 

 

最後にpublishボタンを押す事でpublishingのジョブが開始される。

f:id:kazumaxneo:20200724161058p:plain

 

引用

Structural and Functional Annotation of Eukaryotic Genomes with GenSAS
Humann JL, Lee T, Ficklin S, Main D

Methods Mol Biol. 2019;1962:29-51

 

関連