インタラクティブなRNA seq解析webアプリケーション iDEP

2019 5/23 tweet追記、9/26 動画追加、11/30 ツイート追記、12/22 統合TVリンク追加

2020 2/2 ８章補足資料リンク追加、12/15  ツイート追加

2024/04/04 論文追加引用

　RNAシークエンシング（RNA-Seq）[1]は、ゲノムワイドな発現解析のための日常的な技術となった。ますます低コストで、ライブラリー構築およびシーケンシングはしばしば標準的なプロトコルに従って実施することができる。多くの研究者、特にバイオインフォマティクスの経験がない研究者にとって、この技術を十分に活用するためのボトルネックは、発現プロファイルを実用的な洞察に変換する方法である。典型的な分析ワークフローは多くのステップを含み、各ステップは異なるツールを必要とする。これらのツールを正しく学習し、調整し、繋ぐのは時間がかかる。もう1つの障害は、多様な種類の遺伝子IDを持つ散在したアノテーションデータベースである。これらの問題を軽減するため、研究者は、RNA-Seqデータを分析するのに必要な時間と労力を大幅に削減できるアプリケーションの開発を目指している。

　RNA-Seqデータ解析は、クオリティ管理、前処理、マッピング、そしてrawシークエンスリードのsummarizingから始まる。これらのステップは、従来のTuxedo Suite [ref.2、3]またはより高速でアライメントフリーの定量方法[ref.4、5]のいずれかを使用して完了したと想定する。これらのツールは、スタンドアロンでも、GenePattern [ref.6]、Galaxy [ref.7]、CyVerse [ref.8]（紹介）のようなプラットフォームでも使用できる。

　リードマッピングの後、よくあるのは、遺伝子レベルのリードカウントまたは正規化された発現レベルの行列（Fragments Per Kilobase Million、またはFPKM）を取得することである。Rは、そのDNAマイクロアレイデータのような表形式データの強力な可視化と統計分析ツールである。さらに、differentially expressed genes（DEG）および変更されたパスウェイを同定するために、多くの専用のRおよびBioconductor [ref.9]パッケージが開発されてきた。 DESeq2 [ref.10]のようないくつかのパッケージは、特にリードカウントの統計的モデリングのために開発され、広く使われている。しかし、これらのパッケージは時間がかかり、コーディング経験がないと研究者にとって手が届かないことさえある。

  発現データを分析するいくつかのウェブアプリケーションが開発されている（論文 表１）。 STARTアプリ（shinnyトランスクリプトーム解析リソースツール）は、階層的クラスタリング、主成分分析（PCA）、遺伝子レベルのボックスプロット、およびDEGを実行するshinnyアプリである[ref.11]。別の同様のツールであるDegust [ref.12]はEdgeR [ref.13]またはlimma-voom [ref.14]を使って発現解析を実行し、インタラクティブに結果をプロットできる。他のツールには、Sleuth [ref.15]とShinyNGS [ref.16]がある。shinnyでないアプリケーションも開発された。これにはDEIVA [ref.17]とVisRseq [ref.18]が含まれる。いくつかのツールはパスウェイ分析のいくつかの能力を組み込んでいる。定量された発現データについては、ASAP（自動化シングルセル分析パイプライン）[ref.19]により、Gene Ontology（GO）[ref.20]およびKEGG [ref.21]データベースに基づいて正規化、フィルタリング、クラスタリング、およびエンリッチメント分析を行うことができる。 EXPath Tool [ref.22]を使用すると、ユーザーはパスウェイ探索、GOエンリッチメント解析および共発現解析を実行できる。 IRISなどの他のいくつかのShinyベースのツール[ref.23]も開発されている。過去数年間におけるこれらのツールの開発は、RNA-Seqデータの解釈を容易にした。

  本研究では、（1）広範囲にわたる自動遺伝子ID変換、（2）植物と動物の両方に対する包括的な遺伝子アノテーションとパスウェイデータベース、（3）いくつかの方法、詳細なEDAおよびパスウェイ解析、（4）アプリケーションプログラミングインタフェース（API）を介したKEGG [ref.21]、STRING-db [ref.24]などのWebサービスへのアクセス、および（5）スタンドアロン解析用のRスクリプトの生成による再現性の向上、を含むwebアプリケーションを開発する。iDEPを使用して、2つのサンプルデータセットを分析し、論文表1と図1以外のすべての図と表を生成した。(以下略)

iDEP overview 

https://idepsite.wordpress.com/

FAQ

https://idepsite.wordpress.com/faq/

データフォーマット

Data format – iDEP: Gain Insights from RNA-seq

Github

2022/04/26

Cut the crap and let biologists interact with their data!

With iDEP 0.951 update https://t.co/JaBKVawVzZ, you can now easily visualize the expression of all genes on all KEGG pathways.

Just upload your data, define comparisons, and select any pathway at the KEGG tab. pic.twitter.com/MUIVkUxm5T

— Steven Ge (@StevenXGe) April 24, 2022

So happy that I improved the speed of data loading and converting in iDEP(https://t.co/kaZKnaRsF2). Much faster! Initially, I thought it is the database. But it turns out that I only need to examine a small portion of your genes to determine which species your data is from.

— Steven Ge (@StevenXGe) April 26, 2022

2021 10/20

New features for RNA-seq analysis tool iDEP (https://t.co/LlIaeQ71gL)
1. By default enrichment analysis is done using background defined by genes passed the filter.
2. Interactive network enabling users to easily visualize the relatedness of significant pathways EnrichmentMap pic.twitter.com/iFLI6FPPrP

— Steven Ge (@StevenXGe) October 19, 2021

2021 1/19

A new video demo for iDEP focusing on the biological questions.
iDEP web application for RNA-Seq data analysis https://t.co/HBgFkmznDR via @YouTube

— Steven Ge (@StevenXGe) 2021年1月19日

Our iDEP paper is finally out. In addition to the source code in Github, all the supporting databases are available on Zenodo.

iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data https://t.co/9D1J1eX04A #bmcbioinformatics

— Steven Ge (@StevenXGe) December 20, 2018

2019/03/21

Big update to iDEP!
Besides the 315 genomes in Ensembl, we incorporated the annotation for 115 archaeal, 1678 bacterial, and 238 eukaryotic species from STRING-db. iDEP is designed to help smaller labs making sense of transcriptomic data.
Please help spread the word! pic.twitter.com/6rw0X7tFNw

— Steven Ge (@StevenXGe) May 21, 2019

Re-analyze public RNA-seq data can be done easily with interactive tools like iDEP and uniformly processed read counts data from ARCHS4 and DEE2. See how I did it in under 12 minutes.
No more excuses!
We just need good questions and critical thinkers!https://t.co/CZB28e0U0Y

— Steven Ge (@StevenXGe) September 23, 2019 

使い方

local環境でも簡単に実行できるが、ここではidepサイトにアクセスして動作を確認する。

http://bioinformatics.sdstate.edu/idep/ 

v1.1 test mode

http://bioinformatics.sdstate.edu/idep11/

デモデータで読み込むフォーマットを確認する。左上のClick Hereをクリック。

 読み込まれた。

１列目がgene IDになる。Ensembl IDを使うが、一般的なgene IDなら自動認識してEnsembl IDに内部変換してくれる。１行目にはサンプル名を記載する。Control、 TreatmentA、TreatmentBの３条件でそれぞれbiological/technical replicatesが２つずつあるなら、Ctrl_1, Ctrl_2, TrtA_1, TrtA_2, TrtB_1, TrtB_2のように記載する。これでreplicatesとして認識される。

実際に読み込めるデータは、RNA seqのリードカウントデータ（正規化前）、またはFPKMで正規化したデータ、マイクロアレイのデータなどになる。

他の入力例（データフォーマットの解説より）

Gene expression matrix ex.1

Gene expression matrix ex.2

fold-change and P values

Public RNA-Seq and ChIP-Seq dataでは、ARCHS4でマイニングされた7000以上のSRAサンプルの定量データ（CSVファイル）（paper）を解析できる。

DownloadからARCHS4でリードカウントされたデータのCSVをダウンロードして使ってみる。1行目はこのようにした。

左上のリードカウントからCSVをアップロード。

アップロード完了。フォーマットにエラーがあると赤字のメッセージが出る。gene IDは、Ensembl IDに変換され、以後の解析が行われる。IDはEnsembl release 92に基づいており、220の生物種に対応している（link）（2018/12/29現在）。

１、Pre-process（説明）: データの正規化

replicatesとして認識されたサンプルはグラフで同じ色がアサインされている。動作は非常に俊敏で、データのアップロード完了から１０秒くらいで結果は可視化される。

左上のウィンドウからパラメータを変更できる。デフォルトでは1サンプルで0.5 counts per million (CPM) 以上の遺伝子が分析対象となる。 pseudo countとして追加される値はデフォルト4になっている。

２、Heatmap（説明）階層的クラスタリングおよびデンドログラムとヒートマップによる可視化

defaultでは全サンプル間のSD top1000が対象となる。左上のgene SD distributionで可視化。

左上のinteractive heatmapボタンからスケールをリアルタイムに変更できる。50遺伝子に絞った。

グラフはマウスで操作できる。

図と対象遺伝子のCSVは左下からダウンロードできる。

３、k-means（説明）非階層クラスタリングおよびデンドログラムとヒートマップによる可視化

defaultでは全サンプル間のSD top 2000がクラスター分析対象となる。クラスター数のdefaultは4。

t-SNEによる次元圧縮結果

エンリッチされた転写因子（TF）結合モチーフ

エンリッチされたpathway

デンドログラムで可視化（左のVisualize enrichmentボタンより）

４、PCA（説明）主成分分析

最初の次元は、サンプルを最もよく分離し、データの変動の最も大きな割合を説明する倍率変化を表す。それ以降の次元は効果が小さく、その前の次元と直交している。実験計画で複数の因子が含まれる場合、各因子を複数の次元で調査する。ある次元において、ある因子によってサンプルがクラスタリングされる場合、その因子は発現の違いに寄与していることが示唆される（線形モデリングに含める価値がある）。一方、ほとんどあるいは全く効果を示さない因子は、ダウンストリーム解析から除外できる。

５、DEG（説明）発現変動遺伝子の検出

デフォルトではDESeq2が使われる。FDR cutoffは0.1、Min fold changeは2。３つ以上のグループがある場合、全てのペアワイズ比較が実施される。これまでと同様に、結果は左のメニューからダウンロードできる。

ベン図（Venn Diagramボタンより実行できる）

３サンプルの場合

DEGのタブは２つある。DEG2のタブでは、３群以上のデータでも、組み合わせを選び、それぞれ２群比較を実行できる。結果はヒートマップ、MA-plot、Scatter plotなどで可視化される。また、検定されてエンリッチされていると判定されたpathwayが下にp-value付きで表示される。

操作パネル左下にはShinyGO（prepirnt）へのリンクもある。

ShinyGOはユーザーが指定した遺伝子リストを元に、エンリッチされた系を可視化する。機能はiDEPとかなり重複しているが、iDEPにないツールもある。

２群間比較結果の可視化

Scatter plot

MA plot

Volcano plot

６、Pathway（説明）GOエンリッチメント解析

pathway解析はDEGで指定してフィルタリングした条件ではなく、全データからDESeq2 / limmaで出力したfold change値を使って行われることに注意する。有意水準を決めるFDR cutoffはpathway解析パネルにもあるので、こちらで厳密さは調節する。

有名な、いくつかのGOエンリッチメント解析ツールを利用できる。

GAGE（Generally Applicable Gene-set Enrichment ）（pubmed）（遅い）

GSEA (Gene Set Enrichment Analysis) (preranked fgsea)（link）（解説HP）

PGSEA: PAGE (Parametric Analysis of Gene Set Enrichment) の実装 （PDF link）２群間比較

PGSEA w/all sample: 全サンプル間比較

上記の分析は全てbuint-inのデータを使っているが、最後のReactomePA (Reactome Pathway Analysis) は、 Reactomeの遺伝子セットデータベースを使う（Reactome 統合TV解説）。

Reactome (HP) は歴史あるデータベースで、ヒト以外のモデル生物種にも対応している。ペーパーは多数出ている (Google scholar検索結果)。

たくさんのgenesetを利用できる。一覧はマニュアル参照（link）。KEGGに切り替えてエンリッチと判定されたpathwayを可視化する。

Down regulationが35で、NES（Normalised enrichment score for the given gene set）がもっとも低かったHedgehog signaling pathway

可視化するKEGG pathwayは図の上のメニューから選ぶ。

７、Genomes　DGEのGenome上の位置を可視化する

８、Biclustering（説明）サンプル数の多いラージデータセットの解析で（>10）、相関関係にある遺伝子のグループを検出する

biclust R、BCCC、BCXmotifs、BCPlaid、BCSpectral、BCBimax、BCQuest、QUBICなどのRパッケージを利用できる。解説は上の説明リンク参照。

９、Network（説明）共発現解析およびネットワークの可視化。かなり巨大なデータセットで（>15）相関関係がありそうな遺伝子セットを探すために使う（小さなデータセットではクラスタリング解析を行う）。

ネットワークは巨大になるので、全ネットワークの可視化には、localでcytoscape等を使うことが推奨されている。

感想

インタラクティブにパラメータを変更して、結果を見比べながら進められる素晴らしいツールですね。salmonやkallistoなどと組み合わせれば、rawデータが手に入ってから数時間以内にラフな結果は出せるのではないでしょうか（*1）。もちろん、方法について理解していないと使いこなせない訳で、誰でもRNA seq解析できるとは言いません。しかし、RNA seq解析の敷居が下がっているのは確かで、これだけ完成度の高いアプリケーションが出てくると、リードカウントはCyverse（紹介）、DEG検出などはiDEP、と使い分けることで、コンソールが一切使えない環境でも結果を出せてしまえますね。

2019 11/30追記

羊土社さんからRNA-Seqデータ解析の専門書が発売されました。その中の１章で、iDEPの使い方を紹介しています。よろしければ手にとってご確認下さい。どの章もたいへん見応えのある内容になっています。

【新刊】『RNA-Seqデータ解析　WETラボのための鉄板レシピ』発売！ 食材調達（データ入手）から盛付（論文投稿）までシェフ直伝！「内容の更新情報が届く」という驚きの付録つき！ #RNAseqRecipe https://t.co/bPPUTmA2tt

— 羊土社 実験医学 (@Yodosha_EM) 2019年11月29日

実験医学別冊：RNA-Seqデータ解析 WETラボのための鉄板レシピ - 羊土社

８章補足資料

2020 12/14

#RNAseqRecipe の Annual Update が出たようですゾ “実験医学別冊RNA-Seqデータ解析 WETラボのための鉄板レシピ 坊農秀雅／編 2019年11月29日発行 AB判 255ページ 付録：Annual Updateサービスの登録コード ISBN 978-4-7581-2243-6” @Yodosha_sheep / “実験医学別冊：RNA-S…” https://t.co/YJe6vJKIIL

— bonohu / ぼうのう ひでまさ (@bonohu) 2020年12月14日

2019 12/22

ShinyGO

2020 4/21

統合TVでBioJupiesの使い方が紹介されています。

BioJupiesを使ってウェブブラウザ上でRNA-seqデータ解析を行う

2020 5/26

2022/10

少し前にメールで連絡をいただきましたが、完全に新しくなったiDEP 1.0がテストモードで公開されています。

20/27

iDEPにファイルをアップするとすぐに切断されるという現象を経験しましたが、原因は、自分の不注意で、アップしたテキストファイルの先頭に不要な空白のタブが挿入されていたためでした。タブの後ろに隠れスペースがあってもこの現象が発生します。注意してください。

引用

iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data

Steven Xijin Ge, Eun Wo Son, Runan Yao

BMC Bioinformatics 2018 19:534

追記

iDEP Web Application for RNA-Seq Data Analysis

Xijin Ge

Methods Mol Biol. 2021:2284:417-443.

参考

slideshare バイオインフォマティクスによる遺伝子発現解析

https://www.slideshare.net/sesejun/ss-24923282

*1

75-bpのシングルエンドシーケンシングと組み合わせれば、ラン開始から24時間で解析結果まで得られることになる。

pre-processing

mapping

read count

関連ツール

複数フローセル比較にも対応したONTの分析ツール MinIONQC

　Oxford Nanopore Technologies（ONT）の小型で携帯可能な機器MinIONは、DNAシークエンシングに革命をもたらした。それはユーザーがサンプルから数時間でシーケンスまで進めることを可能にし、また非常に長いDNA分子をシーケンスすることができ、そして各フローセルからmany giga basesのデータを提供する。このため、多くの研究グループや企業が社内および現場でのシーケンスのためにこの装置を採用している。

　ここではMinIONQCを紹介する。MinIONのシーケンスデータのクオリティ管理および診断分析を行うための高速、軽量、および診断分析型スクリプトである。 MinIONQCは、主に大量のシークエンシングの迅速で複製可能な比較に焦点が当てられているという点で、関連するツール（De Coster et al、2018; Loman and Quinlan、2014; Stewart and Watson、2017; Watsonら、2015）と異なる。 MinIONQCは、新しいユースケース（新しい設定など）でのMinIONシーケンスの適用や、集計が必要な大規模ゲノムプロジェクトなど、複数のフローセルからのデータ、多くのフローセルからのデータの分析（Austin et al、2017; Jain et al、2017; Jansen et al、2017; Schmidt et al、2017; Tan、2018）を迅速かつ繰り返し比較する必要がある場合に役立つ。

インストール

mac os10.14のR 3.5.1環境でテストした。

依存のインストール 

#Rにて
install.packages(c("data.table", 
 "futile.logger",
 "ggplot2",
 "optparse",
 "plyr",
 "readr",
 "reshape2",
 "scales",
 "viridis",
 "yaml"))

本体 Github

git clone https://github.com/roblanf/minion_qc.git
cd minion_qc/

> Rscript MinIONQC.R -h

$ Rscript MinIONQC.R -h

Usage: MinIONQC.R [options]

Options:

-h, --help

Show this help message and exit

-i INPUT, --input=INPUT

Input file or directory (required). Either a full path to a sequence_summary.txt file, or a full path to a directory containing one or more such files. In the latter case the directory is searched recursively.

-o OUTPUTDIRECTORY, --outputdirectory=OUTPUTDIRECTORY

Output directory (optional, default is the same as the input directory). If a single sequencing_summary.txt file is passed as input, then the output directory will contain just the plots associated with that file. If a directory containing more than one sequencing_summary.txt files is passed as input, then the plots will be put into sub-directories that have the same names as the parent directories of each sequencing_summary.txt file

-q QSCORE_CUTOFF, --qscore_cutoff=QSCORE_CUTOFF

The cutoff value for the mean Q score of a read (default 7). Used to create separate plots for reads above and below this threshold

-p PROCESSORS, --processors=PROCESSORS

Number of processors to use for the anlaysis (default 1). Only helps when you are analysing more than one sequencing_summary.txt file at a time

-s SMALLFIGURES, --smallfigures=SMALLFIGURES

TRUE or FALSE (the default). When true, MinIONQC will output smaller figures, e.g. suitable for publications or presentations. The default is to produce larger figures optimised for display on screen. Some figures just require small text, and cannot be effectively resized.

テストラン

入力はナノポア のAlbacoreでbase callしたときにできるsequencing_summary.txtファイル

cd minion_qc/
Rscript MinIONQC.R -i example_input_minion -o my_example_output_minion -p 2 

 出力

RB7_D3

combinedQC

引用

MinIONQC: fast and simple quality control for MinION sequencing data
R Lanfear M Schalamun D Kainer W Wang B Schwessinger
Bioinformatics. 2018 Jul 23

PCR duplicationにタグをつけたりエラーを取り除く gencore

2018/12/22 タイトル修正

　HIgh depthの次世代シークエンス（NGS）は、癌の精密な診断と治療に広く使用されている。このようなディープシーケンシングデータから、体細胞突然変異を検出して、パーソラナイズされた標的療法または免疫療法のガイドにすることができる。最近、circulating tumor DNA （ctDNA）シーケンシングが癌治療およびモニタリングのための有望なバイオマーカーとして認識されている。腫瘍由来のDNAは、通常、血中の細胞フリーなDNAの一部であるため、ctDNAシークエンシングデータから検出されるバリアントのmutant allele frequency （MAF）は非常に低い（0.1％程度）。このような低スペクトラムのバリアントを検出するために、通常、シーケンシングデプスを増加させる（10,000xを超える可能性がある）。しかし、NGSライブラリーとシークエンシングを作成するプロセスには、エラーがないわけではない。特に、PCR技術を用いたライブラリー増幅は、多くのエラーを引き起こす可能性があり、その結果、NGSデータ解析の結果にいくつかの偽陽性突然変異を引き起こす。
　ライブラリー増幅の結果、NGSデータは重複する可能性がある。シーケンシングデプスが多いほど、データの重複が多くなる。伝統的には、duplicated reads をマークし、ダウンストリーム分析の前にそれらを削除する。デプスの少ないペアエンドのNGSデータの場合、同じ開始および終了マッピング位置のペアリードは、同じ元のDNA断片に由来するduplicated reads として扱うことができる。一緒にクラスタリングされたリードは、次に1つのリードに併合することができる。通常はエラーがランダムに発生するという性質から、クラスター化されたリードグループ内の不一致を除去してコンセンサス読み取りを生成できる。
　しかし、非常に深いシーケンシングのために、異なるDNA断片から同じポジションを有する2つのペアリードが誘導される可能性がある。この可能性は、DNAフラグメントがより短い場合に高くなり得る。例えば、無細胞DNAは、通常、約167bpのピーク長を有し、通常断片化されたゲノムDNAのピーク長よりもはるかに短い。異なるDNA断片由来のシーケンシングリードをよりよく同定するために、unique molecular identifier （UMI）と呼ばれる技術が開発されている。 UMI技術を用いて、各DNAフラグメントは、DNA増幅プロセス前に独特のランダムなバーコードと連結される。 UMIは、リードの正確なクラスタリングに使用できる。
　現在、Picard MarkDuplicatesのような従来の重複除去ツールは、コンセンサスリードを良好に実行することができず、UMI統合データを処理することができない。さらに、Picard MarkDuplicatesの速度が遅すぎ、メモリ使用量が多すぎるため、クラウドベースの展開には適していない。これらの満たされていない要件により、gencoreという新しいツールを開発した。これは、エラーを排除し、コンセンサスリードを生成することによって重複を削除する。

　gencoreには、ソートされたBAMファイルとリファレンスゲノムFASTAファイルの入力が必要になる。 FASTQデータにUMIがある場合は、UMIを読み取りシーケンスからUMIに移動するためにfastp（Chen、Zhou、Chen、＆Gu、2018）を使用して前処理することができる。gencoreは、各クラスタに対して、マッピング位置とUMI（該当する場合）によってペアリードからコンセンサスリードを生成する。

Githubより

インストール

mac os10.14でテストした。

本体 Github

#linux向けバイナリ
wget http://opengene.org/gencore/gencore
chmod a+x ./gencore

#ビルドする
git clone https://github.com/OpenGene/gencore.git
cd gencore
make

 >  ./gencore -h

$ gencore -h

undefined short option: -h

Version: 0.10.0, usage: gencore --ref=string [options] ... 

options:

  -i, --in                   input sorted bam/sam file. STDIN will be read from if it's not specified (string [=-])

  -o, --out                  output bam/sam file. STDOUT will be written to if it's not specified (string [=-])

  -r, --ref                  reference fasta file name (should be an uncompressed .fa/.fasta file) (string)

  -u, --umi_prefix           the prefix for UMI, if it has. None by default. Check the README for the defails of UMI formats. (string [=])

  -s, --supporting_reads     only output consensus reads/pairs that merged by >= <supporting_reads> reads/pairs. The valud should be 1~10, and the default value is 2. (int [=2])

  -a, --ratio_threshold      if the ratio of the major base in a cluster is less than <ratio_threshold>, it will be further compared to the reference. The valud should be 0.5~1.0, and the default value is 0.8 (double [=0.8])

  -c, --score_threshold      if the score of the major base in a cluster is less than <score_threshold>, it will be further compared to the reference. The valud should be 1~20, and the default value is 6 (int [=6])

      --high_qual            the threshold for a quality score to be considered as high quality. Default 30 means Q30. (int [=30])

      --moderate_qual        the threshold for a quality score to be considered as moderate quality. Default 20 means Q20. (int [=20])

      --low_qual             the threshold for a quality score to be considered as low quality. Default 15 means Q15. (int [=15])

  -j, --json                 the json format report file name (string [=gencore.json])

      --debug                output some debug information to STDERR.

      --quit_after_contig    stop when <quit_after_contig> contigs are processed. Only used for fast debugging. Default 0 means no limitation. (int [=0])

  -?, --help                 print this message

実行方法 

ウルトラディープシーケンシングでduplicated readsも多数予想される時、クリーニングされたbamを出力する。

gencore -i in.bam -o out.bam -r hg19.fa -s 2

• -s, --supporting_reads     only output consensus reads/pairs that merged by >= <supporting_reads> reads/pairs. The valud should be 1~10, and the default value is 2. (int [=2])

PCR duplicationにタグをつける。Picard（ピカード）のmarkduplicateの代わりに使える（s=1でリードは除去しない）。

gencore -i in.bam -o out.bam -r hg19.fa -s 1

クオリティの低いリードを除きクリーニングする。

gencore -i in.bam -o out.bam -r hg19.fa -s 1 --score_threshold=8

 UMIt tag付きデータの扱いについてはGithubで確認してください。

引用

gencore: an efficient tool to generate consensus reads for error suppressing and duplicate removing of NGS data
Shifu Chen, Yanqing Zhou, Yaru Chen, Tanxiao Huang, Wenting Liao, Yun Xu, Zhihua Liu, Jia Gu

bioRxiv preprint first posted online Dec. 19, 2018

BED、VCF、GTFをユーザー定義の方法でソートする gsort

gsortはゲノムファイルをソートするためのツール。たとえば、何らかの理由でVCFを並べ替えて、X、Y、2,1,3、などの順序で並べ替えることができる。他のソートツールでは不可能だったGATK order（1 ... X、Y、MT）に一致するようなソートもできる。ソートは、ユーザーが指定したTSVファイルに従って行われる。ソートできるのは、BED、VCF、GTFになる。

インストール

mac os10.14のminiconda2-4.0.5環境でテストした。

 本体 Github

#anaconda環境ならcondaで導入できる
conda install -c bioconda -y gsort

>gsort -h

$ gsort -h

Usage: gsort [--memory MEMORY] [--parent] PATH GENOME

Positional arguments:

  PATH                   a tab-delimited file to sort

  GENOME                 a genome file of chromosome sizes and order

Options:

  --memory MEMORY, -m MEMORY

                         megabytes of memory to use before writing to temp files. [default: 1300]

  --parent, -p           for gff only. given rows with same chrom and start put those with a 'Parent' attribute first

  --help, -h             display this help and exit

実行方法

vcfを指定した ファイル順にソートする

gsort --memory 1500 input.vcf.gz input.genome | bgzip -c > my.crazy-order.vcf.gz

 引用

https://github.com/brentp/gsort

seqkitに新しく追加されたコマンドを確認する

2019 8/7 誤字修正

2023/01/20 translate help更新

seqkitを以前ブログで紹介した時は0..6.0でしたが、１年半近く経ち、2018年12月20日現在ではバージョンが0.9.4まで上がっています。ありがたいことに、bug fixだけでなく、新しいコマンドが複数追加されています。v0.6.1以降に追加されたコマンドについてまとめます。

インストール 

ダウンロード

http://bioinf.shenwei.me/seqkit/download/

解凍するとbinaryファイルできる。

実行権をつけ、それを/usr/local/bin/などに移動して使う。

chmod u+x seqtk #実行権
ls -l seqtk #確認
#パスが通ったディレクトリに移動
mv seqtk /usr/local/bin

#anacondaを使っているならcondaで入る。例えば20181219時点で最新のv0.9.3を入れる
conda install -c bioconda -y seqkit==0.9.3

0.6.0の時点であったコマンドを復習しておく。

アルファベット順

• common　 IDや配列が共通する配列抽出
• dup   　　   配列を指定回、複製
• fq2fa　　　fastqをfastaに変換
• fx2tab　　  fasta/fastqをtab形式に変換し、さらにGCなどの情報を表示
• grep　  　   特定の配列を検索
• head　　　指定数だけ先頭からリードを抽出
• locate　　  配列をリプレース
• rename　   duplicateしたIDをユニークな名前に修復
• rmdup　    duplication除去
• restart　    環状ゲノムのスタート位置変更
• sample　  リードをランダムサンプリング（10%だけランダムに取り出す）
• seq　　     配列を変形（逆相補鎖、DNA<=>RNA、指定サイズ以上/以下）
• shuffle　   ランダムな順番にする
• sliding　    配列を指定の単位ずつずらしながら取り出す
• sort　　    長さでソート
• split　　   分割（fastqを均等にn個に分けて出力）
• stats　      fasta/fastqの簡単なstatistics
• subseq      特定のリードや配列だけ抽出 

新しいコマンド（アルファベット順）

concat （リンク）     同じIDの配列を指定長ずつ合体

$ seqkit concat -h

concatenate sequences with same ID from multiple files

Example: concatenating leading 2 bases and last 2 bases

    $ cat t.fa

    >test

    ACCTGATGT

    >test2

    TGATAGCTACTAGGGTGTCTATCG

    $ seqkit concate <(seqkit subseq -r 1:2 t.fa) <(seqkit subseq -r -2:-1 t.fa)

    >test

    ACGT

    >test2

    TGCG

Usage:

  seqkit concat [flags]

Aliases:

  concat, concate

Flags:

  -h, --help   help for concat

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

上の例のように、１つ目の配列の先頭1-2文字と、２つ目の配列の後ろから1-2文字が結合（同じIDのもの同士）。

seqkit concat <(seqkit subseq -r 1:2 t.fa) <(seqkit subseq -r -2:-1 t.fa)

convert （リンク）       fastqのクオリティエンコードをSolexa形式（1.3）、サンガーイルミナ（1.8）形式で相互変換

$ seqkit convert -h

convert FASTQ quality encoding between Sanger, Solexa and Illumina

Usage:

  seqkit convert [flags]

Flags:

  -d, --dry-run                         dry run

  -f, --force                           for Illumina-1.8+ -> Sanger, truncate scores > 40 to 40

      --from string                     source quality encoding. if not given, we'll guess it

  -h, --help                            help for convert

  -n, --nrecords int                    number of records for guessing quality encoding (default 1000)

  -N, --thresh-B-in-n-most-common int   threshold of 'B' in top N most common quality for guessing Illumina 1.5. (default 4)

  -F, --thresh-illumina1.5-frac float   threshold of faction of Illumina 1.5 in the leading N records (default 0.1)

      --to string                       target quality encoding (default "Sanger")

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

クオリティフォーマット確認

seqkit convert input.fq.gz |seqkit head -n 1

このデータは[Sanger Illumina 1.8

$ seqkit convert SRR5035895_1.fastq |seqkit head -n 1

[INFO] possible quality encodings: [Sanger Illumina-1.8+]

[INFO] guessed quality encoding: Sanger

[INFO] converting Sanger -> Sanger

[WARN] source and target quality encoding match.

@SRR5035895.1 1 length=228

AAAAAGTATCTACGACACGCGATTTAGTGTCGAGAAAGGCGACTAAGACGGTAGCTTCGTCGGTTGCTTCGCTTCCTGCA+

CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGデフォルトの設定でランすると、イルミナ1.8に変換するようになっている。イルミナ1.5に変換。

seqkit convert Illimina1.8.fq.gz --to Illumina-1.5+ > out.fq

duplicate （リンク）      配列を複製

$ seqkit duplicate -h

duplicate sequences N times

You may need "seqkit rename" to make the the sequence IDs unique.

Usage:

  seqkit duplicate [flags]

Aliases:

  duplicate, dup

Flags:

  -h, --help        help for duplicate

  -n, --times int   duplication number (default 1)

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

全配列ではなく、fastaの先頭の配列だけ１回複製する。seqkit renameでユニークなヘッダーにして出力。

seqkit head -n 1 input.fa | seqkit duplicate -n 2 | seqkit rename > out.fa

range （リンク）     指定間のリードを抽出 (headでもtailでも単独では対応できない時に)

$ seqkit range -h

print FASTA/Q records in a range (start:end)

Usage:

  seqkit range [flags]

Flags:

  -h, --help           help for range

  -r, --range string   range. e.g., 1:12 for first 12 records (head -n 12), -12:-1 for last 12 records (tail -n 12)

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

user-no-MacBook-Pro-2:Illumina_MiSeq_metagenome_of_Kelp_associated_biofilm user$ $

先頭から101番目 ~ 150番目の50リード出力

seqkit range -r 101:150 input.fq > out.fq

 後ろから-1番目 ~ -100番目の50リード出力

seqkit range -r -100:-1 input.fq > out.fq

faidx （リンク）    FASTAのindexを作成したり、特定の配列を抽出（samttools faidx互換だが、より高機能）

$ seqkit faidx -h

create FASTA index file and extract subsequence

This command is similar with "samtools faidx" but has some extra features:

  1. output full header line with flag -f

  2. support regular expression as sequence ID with flag -r

Usage:

  seqkit faidx [flags] <fasta-file> [regions...]

Flags:

  -f, --full-head     print full header line instead of just ID. New fasta index file ending with .seqkit.fai will be created

  -h, --help          help for faidx

  -i, --ignore-case   ignore case

  -r, --use-regexp    IDs are regular expression. But subseq region is not suppored here.

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

seqkit faidx input.fa

 input.fa.faiができる。

restart （リンク）　環状ゲノムのスタート位置を変える

$ seqkit restart -h

reset start position for circular genome

Examples

    $ echo -e ">seq\nacgtnACGTN"

    >seq

    acgtnACGTN

    $ echo -e ">seq\nacgtnACGTN" | seqkit restart -i 2

    >seq

    cgtnACGTNa

    $ echo -e ">seq\nacgtnACGTN" | seqkit restart -i -2

    >seq

    TNacgtnACG

Usage:

  seqkit restart [flags]

Flags:

  -h, --help            help for restart

  -i, --new-start int   new start position (1-base, supporting negative value counting from the end) (default 1)

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

 入力環状

ゲノムのポジション10000からリスタート

seqkit restart -i 10000 input.fa > output.fa

split2 （リンク）　入力配列を分割する。

$ seqkit split2 -h

split sequences into files by part size or number of parts

This command supports FASTA and paired- or single-end FASTQ with low memory

occupation and fast speed.

The file extensions of output are automatically detected and created

according to the input files.

Usage:

  seqkit split2 [flags]

Flags:

  -p, --by-part int      split sequences into N parts

  -s, --by-size int      split sequences into multi parts with N sequences

  -f, --force            overwrite output directory

  -h, --help             help for split2

  -O, --out-dir string   output directory (default value is $infile.split)

  -1, --read1 string     read1 file

  -2, --read2 string     read2 file

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

FASTAを10000リードずつに分ける。

seqkit split2 hairpin.fa.gz -s 10000 -f

FASTAを4パートに分ける。

seqkit split hairpin.fa.gz -p 4 -f

 シングルエンドFASTQを2パートに分ける。

seqkit split2 reads_1.fq.gz -p 2 -O out -f

 ペアエンドFASTQを2パートに分ける。

seqkit split2 -1 reads_1.fq.gz -2 reads_2.fq.gz -p 2 -O out -f

tab2fx （リンク）　FASTA/FASTQをタブ区切り（TSV）フォーマットに変換

$ seqkit tab2fx -h

convert tabular format (first two/three columns) to FASTA/Q format

Usage:

  seqkit tab2fx [flags]

Flags:

  -p, --comment-line-prefix strings   comment line prefix (default [#,//])

  -h, --help                          help for tab2fx

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^([^\\s]+)\\s?")

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

TSV形式で10リードだけ出力

seqkit fx2tab input.fq |head -n 10 > output.tsv

配列長、GC率だけ表示（*csvtkはcondaで導入可）

seqkit fx2tab contig.fa.gz -l -g -n -i -H | head -n 10 | csvtk -t -C '&' pretty

• -g    print GC content
• -H   print header line
• -l     print sequence length
• -n    only print names (no sequences and qualities)
• -i     print ID instead of full head

#name                                seq   qual   length   GC

NODE_1_length_254790_cov_21.168655                254790   47.85

NODE_2_length_246122_cov_21.354689                246122   47.88

NODE_3_length_214379_cov_20.356971                214379   48.91

NODE_4_length_200639_cov_20.843756                200639   47.45

NODE_5_length_186136_cov_20.060989                186136   48.74

NODE_6_length_183861_cov_21.459733                183861   47.74

NODE_7_length_180715_cov_21.655153                180715   47.09

NODE_8_length_141474_cov_21.429290                141474   48.11

NODE_9_length_130080_cov_20.966290                130080   47.76

translate （リンク）  DNA/RNAをアミノ酸配列に翻訳

$  seqkit translate -h

translate DNA/RNA to protein sequence (supporting ambiguous bases)

Note:

  1. This command supports codons containing any ambiguous base.

     Please switch on flag -L INT for details. e.g., for standard table:

        ACN -> T

        CCN -> P

        CGN -> R

        CTN -> L

        GCN -> A

        GGN -> G

        GTN -> V

        TCN -> S

        MGR -> R

        YTR -> L

Translate Tables/Genetic Codes:

    # https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes

     1: The Standard Code

     2: The Vertebrate Mitochondrial Code

     3: The Yeast Mitochondrial Code

     4: The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code

     5: The Invertebrate Mitochondrial Code

     6: The Ciliate, Dasycladacean and Hexamita Nuclear Code

     9: The Echinoderm and Flatworm Mitochondrial Code

    10: The Euplotid Nuclear Code

    11: The Bacterial, Archaeal and Plant Plastid Code

    12: The Alternative Yeast Nuclear Code

    13: The Ascidian Mitochondrial Code

    14: The Alternative Flatworm Mitochondrial Code

    16: Chlorophycean Mitochondrial Code

    21: Trematode Mitochondrial Code

    22: Scenedesmus obliquus Mitochondrial Code

    23: Thraustochytrium Mitochondrial Code

    24: Pterobranchia Mitochondrial Code

    25: Candidate Division SR1 and Gracilibacteria Code

    26: Pachysolen tannophilus Nuclear Code

    27: Karyorelict Nuclear

    28: Condylostoma Nuclear

    29: Mesodinium Nuclear

    30: Peritrich Nuclear

    31: Blastocrithidia Nuclear

Usage:

  seqkit translate [flags]

Flags:

  -x, --allow-unknown-codon                     translate unknown code to 'X'. And you may not use flag --trim which removes 'X'

  -F, --append-frame                            append frame information to sequence ID

      --clean                                   change all STOP codon positions from the '*' character to 'X' (an unknown residue)

  -f, --frame strings                           frame(s) to translate, available value: 1, 2, 3, -1, -2, -3, and 6 for all six frames (default [1])

  -h, --help                                    help for translate

  -M, --init-codon-as-M                         translate initial codon at beginning to 'M'

  -l, --list-transl-table int                   show details of translate table N, 0 for all (default -1)

  -L, --list-transl-table-with-amb-codons int   show details of translate table N (including ambigugous codons), 0 for all.  (default -1)

  -T, --transl-table int                        translate table/genetic code, type 'seqkit translate --help' for more details (default 1)

      --trim                                    remove all 'X' and '*' characters from the right end of the translation

Global Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^(\\S+)\\s?")

      --infile-list string              file of input files list (one file per line), if given, they are appended to files from cli arguments

  -w, --line-width int                  line width when outputting FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. can also set with environment variable SEQKIT_THREADS) (default 4)

 １フレームで翻訳

seqkit translate --frame 1 --transl-table 1 input.fa > output.fa

#The Bacterial, Archaeal and Plant Plastid Code
seqkit translate --frame 1 --transl-table 11 input.fa > output.fa

• -T, --transl-table int      translate table/genetic code, type 'seqkit translate --help' for more details (default 1)

 --trimをつけると*を除ける。

codon tableリンク

他にもbashに対応したシェル自動補完コマンドがあります。

引用

SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation

Wei Shen, Shuai Le, Yan Li , Fuquan Hu 

PLOS ONE. 2016. doi:10.1371/journal.pone.0163962.

UMIタグつきraw シーケンシングリードをクラスタリングする calib

　次世代シーケンシングにより、シーケンシングエラーの処理を含む多くの課題が発生する大規模なゲノムデータセットが利用可能になった。これは特にガンゲノミクスに関連する。循環腫瘍DNAからの低い対立遺伝子頻度変動を検出するために使用される。ユニークな分子識別子（UMI）を有するDNA分子のバーコードタグは、シークエンシングエラーを緩和しようと試みる。 UMIタグ分子のPCRコピーは独立にシーケンシングされる。しかしながら、PCRおよびシークエンシングプロセスではエラーを生成し得る。 UMIタグ付きシーケンシングデータを分析するには、同じUMIタグ分子のPCR duplicatesからシーケンシングされたリードを単一のクラスターにグループ化することを目的とした初期クラスタリングステップが必要であり、現在のデータセットのサイズにはこのクラスタリングプロセスがリソース効率が必要である。
結果
　イルミナのような置換エラーがドミナントなシーケンシングプラットフォームによって生成されたUMIタグ付きシークエンシング実験からのペアエンドリードをクラスタリングする計算ツールであるCalibを紹介する。 Calibでのクラスタは、バーコードの類似性とリードシーケンス類似性の両方で定義されているエッジにより連結したグラフとして定義される。グラフは、 locality sensitive hashingとMinHashingテクニックを使用して効率的に構築される。 Calibのデフォルトのクラスタリングパラメータは、Calibとパッケージ化されたシミュレーションモジュールを使用して、UMIとリード長が異なるため、経験的に最適化されている。他のツールと比較して、Calibは、妥当なランタイムとメモリフットプリントを維持しながら、シミュレートされたデータに対して最高の精度を実現する。実際のデータセットでは、Calibはアライメントベースの方法よりもはるかに少ないリソースで実行され、そのクラスターは下流のコールで偽陽性の数を減らす。

calibに関するツイート

インストール

mac os10.14のminiconda2-4.0.5環境でテストした。

依存

• pyfaidx
• numpy
• scipy
• scikit-learn
• biopython
• pandas
• ART Illumina (version 2.5.8)

• plotly

本体　GIthub

git clone https://github.com/vpc-ccg/calib.git calib
cd calib

#To install Calib clustering module:
make

#To install Calib error correction module:
make -C consensus/
cd ..

 > calib

$ calib 

Barcode length must be a positive integer!

Calib: Clustering without alignment using LSH and MinHashing of barcoded reads

Usage: calib [--PARAMETER VALUE]

Example: calib -f R1.fastq -r R2.fastq -o my_out. -e 1 -l 8 -m 5 -t 2 -k 4 --silent

Calib's paramters arguments:

-f --input-forward                 (type: string;   REQUIRED paramter)

-r --input-reverse                 (type: string;   REQUIRED paramter)

-o --output-prefix                 (type: string;   REQUIRED paramter)

-s --silent                        (type: no value; default: unset)

-q --no-sort                       (type: no value; default:  unset)

-l --barcode-length                (type: int;      REQUIRED paramter)

-p --ignored-sequence-prefix-length (type: int;      REQUIRED paramter)

-m --minimizer-count               (type: int;      REQUIRED paramter)

-k --kmer-size                     (type: int;      REQUIRED paramter)

-e --error-tolerance               (type: int;      REQUIRED paramter)

-t --minimizer-threshold           (type: int;      REQUIRED paramter)

-c --threads                       (type: int;      default: 1)

-h --help

実行方法

ペアエンドのfastqとUMIタグの長さを指定する。

calib -f pair_1.fq -r pair_2.fq -l <barcode_tag_length> -o output

タグの長さは、ペアエンドの1つのメイトのバーコードタグ長を記載する。

他にもPCRサイクルによるエラー導入をシミュレートするツールなどが付属しています。

引用

Alignment-free clustering of UMI tagged DNA molecules

Baraa Orabi Emre Erhan Brian McConeghy Stanislav V Volik Stephane Le Bihan Robert Bell Colin C Collins Cedric Chauve Faraz Hach

Bioinformatics, bty888, Published: 23 October 2018

関連ツール

SPAdesとUnicyclerでlarge k-merを使う part2 (SPAdesのテスト)

2019 12/8 誤字修正

127以上のk-merを使うために、SPAdesとUnicyclerをビルドし直した（リンク）。今回は、実際に127以上のk-mer値でアセンブリを行い、アセンブリ性能がどのように変化するか簡単にテストした結果を書く。 Real dataの傾向が知りたいので、GAGE-BのMiseqでシーケンシングされたVibrio cholerae CO 0132データセット(250bp x 2、x100)を使う（GAGE-B paperの表１参照）。

HP:  https://ccb.jhu.edu/gage_b/datasets/index.html

 テストOS

• ubuntu 18.04、Anaconda3.5.0

結果

Vibrio cholerae CO 0132の250bp x 2のペアエンドfastqをHPからダウンロードした。

１、merged.fqを作成

bbmerge.sh in1=reads1.fq in2=reads2.fq out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq ihist=ihist.txt

#リード長を分析
seqkit stats *fq.gz

file                 format  type  num_seqs      sum_len  min_len  avg_len  max_len

GAGEB_reads_1.fq.gz  FASTQ   DNA    796,812  199,999,812      251      251      251

GAGEB_reads_2.fq.gz  FASTQ   DNA    796,812  199,999,812      251      251      251

merged.fq.gz         FASTQ   DNA    755,374  149,383,904       39    197.8      493

mergedの平均サイズが短い。あまり長いライブラリは作られていない。

２、fastpによるクオリティチェック（出力指定なしだとfastpはQT前後の図だけ出力する）

fastp -i GAGEB_reads_1.fq.gz -I GAGEB_reads_2.fq.gz 

GAGEB_reads_1 (before)

GAGEB_reads_2 (before)

いつものように、R1は3'末端1bpだけクオリティが低い。R2はより悪く、200bpあたりでQ30を下回っている。エラー率がアセンブリへ与える影響を調べるためには、トリミングの有り・無しでアセンブリ結果を比較すべきだが、今回は時間の関係でQTあり/無しの比較は行わない。トリミングなしで進め、k-merの選択による違いだけ評価する。

３、spadesパッケージのBayesHammer(link)を使ってerror correction

spades.py -1 GAGEB_reads_1.fq.gz -2 GAGEB_reads_2.fq.gz --merged.fq.gz --only-error-correction -o error_correction


#出力をカレントにコピー
mv error_correction/corrected/*gz .

#rename
mv GAGEB_reads_1.00.0_0.cor.fastq.gz GAGEB_1_cor.fq.gz
mv GAGEB_reads_2.00.0_0.cor.fastq.gz GAGEB_2_cor.fq.gz
mv merged.00.0_1.cor.fastq.gz merged_cor.fq.gz

４、de novo assembly

spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21 -o spades_k21 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31 -o spades_k21-31 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41 -o spades_k21-41 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51 -o spades_k21-51 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61 -o spades_k21-61 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71 -o spades_k21-71 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81 -o spades_k21-81 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91 -o spades_k21-91 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,101 -o spades_k21-101 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111 -o spades_k21-111 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121 -o spades_k21-121 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131 -o spades_k21-131 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141 -o spades_k21-141 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151 -o spades_k21-151 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161 -o spades_k21-161 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171 -o spades_k21-171 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181 -o spades_k21-181 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191 -o spades_k21-191 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201 -o spades_k21-201 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211 -o spades_k21-211 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221 -o spades_k21-221 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231 -o spades_k21-231 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231,241, -o spades_k21-241 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231,241,249 -o spades_k21-249

 （xeon E5 v3 x2サーバーでトータル10時間ほどかかった）

５、評価

quast k21.scaffolds.fasta k21-41_scaffolds.fasta ...（省略） \
-t 24 -o quast_result -1 GAGEB_reads_1.fq.gz -2 GAGEB_reads_2.fq.gz \
-R GCF_000006745.1_ASM674v1_genomic.fna

k21は極端に悪い。その後、k-merサイズが増えるにつれ、アセンブリの連続性は徐々によくなる。k100以上になるとほぼ変化がない。

 Quast report1

リファレンスがあるので、それをgold standardとして、アセンブリ精度に関する分析も行う。

 Quast report2

icarus.html

 Quast report2の重要な項目のみ、表にまとめ直した(*2)。

1.  scaffolds数はk-mer値が増えるにつれ減少
2.  10kb以上のscaffolds数もk-mer値が増えるにつれ減少
3.   NG50はk-mer値180付近で最大になり、それ以上k-merを増やすと減少
4.   genome fraction（カバー率）はk-mer値180-200付近で最大になり、それ以上k-merを増やすと0.2%ほど減少 (10000-bp前後)
5. ミスアセンブリはk-mer値が小さい方が少ない。k-mer値が190を越えるとミスアセンブリ数は10を越える。その大半がrelocations。
6. ミスマッチ、indel数は大差なし

k-mer値が増えるにつれscaffolds数が減少しているが、NG50もk-merが増えると減少している。アセンブリの連続性が改善されてscaffolds数が減少しているなら、直感的にはNG50も増えるように思える。最初は理由がわからなかったが、アセンブリのGFAファイルをbandageで可視化してみると、何が起きているのか理解できてきた。6に続く。
 
 

６、Bandageを使いGFAファイルを可視化する。端末から呼び出し、SVG出力。

K21

Bandage image spades_k21/assembly_graph.fastg k21.svg

 K21-71

Bandage image spades_k21-71/assembly_graph.fastg k21-71.svg

 K21-101

Bandage image spades_k21-101/assembly_graph.fastg k21-101.svg

K21-131

Bandage image spades_k21-131/assembly_graph.fastg k21-131.svg

K21-171

Bandage image spades_k21-171/assembly_graph.fastg k21-171.svg

K21-191

Bandage image spades_k21-191/assembly_graph.fastg k21-191.svg

K21-231

Bandage image spades_k21-231/assembly_graph.fastg k21-231.svg

　Vibrio choleraeのゲノムはchromosome２つだけからなる。k-merが大きくなるについてgraphはシンプルになるが、k-mer190くらいからはorphanなcontigが増えてくる。k-merが小さいうちはこのような断片化は起きなかったことも加えて考えると、断片化は、k-merピークが得られるか得られないのボーダーギリギリまでk-merの値を増やすことで、k-merカバレッジが足りない領域ができて発生したと推測される。

　この仮説は、先ほどの表で、scaffolds数の減少とNG50減少が同時に起こったことも説明できる。すなわち、k-mer値を増やすと、より長いunitigを作り出せるが、同時にk-merカバレッジが不足する領域が確率的に発生し、そこでアセンブリが切れてしまう。k-merを減らしすぎると、断片化は防げ、計算も早くなるが、今度は短いリピートの解消もできなくなる。例えば、上のデータでは、k-merが100から200くらいがバランスがいいように感じられるが、これは、リード長はいくつか、どれだけのカバレッジが得られたか、どのくらいのエラー率か、どのくらい複雑なゲノムか（*3）、またどれだけ計算リソースを使えるか、などによって変化しうる（GAGE-Bの結論とだいたい同じ）。de Bruijn Graphアセンブリでより長いcontigを作る上でのベストなk-merサイズ選択の難しさを示している（*4）。

結論

一度考えましたが考察が甘かったので破棄しました。もう少し考えて追記します。

 以上、１例報告であまり参考になるようなものでなかったかもしれませんが、簡単なテストでした。

今回使用したツール

引用

SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing

Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev,Pavel A. Pevzner, 

J Comput Biol. 2012 May; 19(5): 455–477.

参考

*1

バーコードを付けた相乗りシーケンスなら、１サンプルあたりのコストは大きく抑えられるが、そうでないなら、ロングリードシーケンスのコストは、バクテリア向けには高い。

*2

上では20merずつ変化させたが、ほんとは10merずつ増やしている。見えないかもしれませんが、全データを載せておきます。

 

 こちらはGAGE-B 当時のアセンブリstatistics。

*3

極端な話、全くリピートがないスモールゲノムなら、短いk-merでも完全なアセンブリが可能になる。

*4

https://github.com/rrwick/Unicycler

UnicyclerのGithub HPより転載