macでインフォマティクス

macでインフォマティクス

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

カスタムアノテーションを使った GO enrichment解析の例

2022/01/08追記, 1/13インストール追記

2022/12/25 説明追記

 

タイトルの通りの内容です。

質問があったのでそれに対応した記事になります。

 

1、アノテーションファイルの準備

TrinotateかeggNOG mapperを使ってGO termをアサインしていることを想定している。

 

A - Trinotateのアノテーション結果

Trinotate(Trinityやtransdecoderを開発したBrian Haas参考)が作ったアノテーションツール)のアノテーション結果のSQLiteファイルを所有しているなら、orf_to_go.rbスクリプトレポジトリ)を使って遺伝子(転写産物)とGO termのリストを抽出できる。

git clone https://bitbucket.org/njcaruana/trinotate-proteomics.git
cd trinotate-proteomics/
./orf_to_go.rb Trinotate.sqlite > orf2go.map

TrinotateはTrinityと連携しており色々連携可能だが、モデル生物から距離がある場合、GO termのアサイン率はeggNOG mapperと比べて劣る可能性がある。

 

B - eggNOG mapperのアノテーション結果

eggNOG mapperのアノテーション結果(.tsv)を所有しているなら、既にGO termやKO termがアサインされた結果を所有していることになる。コメント行を捨て、必要な列だけ保存する。

#先頭と末尾のコメント行(##)を捨てる
grep -v "##" eggnog_result.tsv > filtered.tsv

#GOは10列目(12列目はKO)
cut -f 1,10 filtered.tsv > orf2go.map

*下の補足1のような手順で1行に遺伝子名<tab>GO1つのみ、に変換してから使うことを考えるかもしれないが、topGOはGOの階層関係を読み込む。親子関係の情報を無くしてしまうと結果が変化するので気を付ける。また、GOがアサインされていない遺伝子は除かないとエラーになるので注意。

このようなファイルを用意する。ここではタブ区切り。1行目は#query<tab>GOs

orf2go.map

 

 

2、クエリの準備

トランスクリプトーム解析やプロテオーム解析などで得た関心のある遺伝子リスト(典型的には発現変動遺伝子のリスト)を準備する。1のアノテーション結果と同じ遺伝子名になっていないといけない。ここではファイル名”input”で保存した。

f:id:kazumaxneo:20220107235623p:plain

遺伝子名のみを1行1遺伝子ずつ記載したファイル。これで準備O.K

 

 

3、GOエンリッチメント解析

準備が出来たらGOエンリッチメント解析を行う。カスタムアノテーションに対応したいくつかの実装があり、webでもagriGO(紹介)のようなカスタムアノテーションが利用できるサービスがあるが、ここではカスタムアノテーションの使い方の詳しい説明があるBioconductorのtopGOパッケージを使う。

Bioconductor - topGO

https://bioconductor.riken.jp/packages/3.3/bioc/html/topGO.html

topGO manual

https://bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf

topGOをインストールする。これは自分の主観だが、Bioconductor のパッケージはBiocManagerで導入できるものの、ローカル環境にパッケージが増えてくると時間がかかる上に割と失敗したりと、非効率的な面が残っている。実はBioconductor のパッケージはbiocondaでも提供されているので、今回はmambaで仮想環境を作って素早く導入する(全てのパッケージ対応しているかどうかは不明)。base環境に入れないように注意する。

mamba create -n topGO -y
conda activate topGO

#conda (bioconda)
mamba install -c conda-forge -c bioconda bioconductor-topgo -y
mamba install -c conda-forge -c bioconda bioconductor-rgraphviz -y

 

準備が出来たらtopGOで解析する。流れはtopGOのマニュアルに則っている。

library(topGO)
library(Rgraphviz)

#readMappingsを使ってバックグラウンド(step1で準備)を読み込む。
geneID2GO <- readMappings(file = "orf2go.map", sep = "\t", IDsep = ",")
#check
str(head(geneID2GO))
#gene name
geneNames <- names(geneID2GO)
#DEGsなどのクエリのリスト(step2)を読み込む
list <- unlist(as.list(read.csv("input", header = TRUE)))
geneList <- factor(as.integer(geneNames %in% list ))
names(geneList) <- geneNames

#これで準備ができた。GODataオブジェクトを構築(ここではBPだがMFとCCも試す)
#nodeSize = 10 は、10未満の語彙からGO階層を刈り取るために使用される。
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,annot = annFUN.gene2GO, gene2GO = geneID2GO, nodeSize = 10)

#topGOdataクラスのオブジェクトを得たらエンリッチメント分析を始めることができる。
#ここでは遺伝子数に基づくフィッシャーの正確検定を行う
#runTest関数3つの引数( topGOdata クラスのオブジェクト,GOグラフ構造を扱う方法,
検定統計量)を指定する
resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

#(ほかにも遺伝子スコアに基づいて濃縮度を計算するKolmogorov-Smirnov ライクな検定が利用可能)classicかelimが利用可能。
resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")

#エンリッチメントテスト実施後、結果を分析し解釈するため、GenTableを使う。GenTable は最も有意な GO タームとそれに対応する p 値を分析する。
GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)

 

エンリッチメントテスト実施後、結果を分析し解釈するため、GenTableが用意されている。GenTable は最も有意な GO タームとそれに対応する p 値を分析する。

GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)

f:id:kazumaxneo:20220108140251p:plain

 GO termの順位とp値を、classicとelim間で比較したりもできる。

GenTable(GOdata, classicFisher = resultFisher, classic = resultFis, elimKS = resultKS.elim, orderBy = "weight", ranksOf = "classic", topNodes = 10)

マニュアルでは各検定結果のp 値を比較するための流れも用意されている(マニュアル3.3)。

 

有意なGO termがどのようなものかGOグラフ上で調べることができる。

allRes <- GenTable(GOdata, classic = resultFis,orderBy = "weight", ranksOf = "classic", topNodes = 20)

#showSigOfNodesは誘導されたサブグラフを表示
showSigOfNodes(GOdata, score(resultFis), firstSigNodes = 5, useInfo = 'all')

#PDF保存(printGraphはshowSigOfNodesのワープ機能)
printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)

重要なノードは長方形で表される。プロットされたグラフは、これらの有意なノードによって生成される上位グラフである。色は相対的な重要度を表す。濃い赤(最も重要)から明るい黄色(最も重要でない)までの範囲。各ノードの中の文字は、最初の 2 行がGO 識別子とトリミングされた GO 名、3 行目には生の p-value 、4行目は有意な遺伝子数とその遺伝子にアノテーションされた総遺伝子数を示している。

f:id:kazumaxneo:20220108003029p:plain

 

topGOで利用できるどのような検定、メソッドを使うかについては、マニュアルのRunning the enrichment testを確認して下さい。

 

引用

Alexa A, Rahnenfuhrer J (2021). topGO: Enrichment Analysis for Gene Ontology. R package version 2.46.0.

 

Trinotate Proteomics

https://bitbucket.org/njcaruana/trinotate-proteomics/src/master/

 

参考

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

 

Enrichment analysis with topGO

https://gist.github.com/iracooke/73c58d11c921fae4fa7c

 

 

https://rpubs.com/aemoore62/TopGo_colMap_Func_Troubleshoot

 

補足1

#1遺伝子に複数GOがアサインされているが、これを1行に1遺伝子1語彙になるように仕分ける(変換後は同じ遺伝子が複数回登場する)。

eggnogmapper <- read.csv("filtered.tsv", header = TRUE, sep = "\t")
head(eggnogmapper ) #check

#GOは10列目、KOなら12列目
go <- eggnogmapper[c(1,10)]
head(go,3) #check

#変換
library(data.table)
go <- data.table(go)


go_sort <- go[, list(GOs = unlist(strsplit(GOs , ","))), by = X.query]
# 書き出す
write.table(go_sort, file = "orf2go.map", sep = "\t", quote=F)

 

補足2

上のリンクからダウンロードできるtopGOのソースコードにはカスタムアノテーションのデモデータとしてtopGO/inst/geneid2go.mapが含まれています。マニュアルの4.3でも使用されているファイルと思われます。

 

関連


 

 

ナノポアシグナルデータを効率的かつ並列に解析するための SLOW5フォーマットを扱う slow5tools

 

 現在、最もポピュラーなシグナルレベル解析は、ソフトウェアNanopolish/f5cを用いたDNAメチル化プロファイリングである。この使用例をもとに、ハイパフォーマンス・コンピューティング(HPC)システムでのFAST5データ解析について解析した(論文補足説明2)。FAST5 は、ONT が定義した特定のスキーマを持つ階層型データフォーマット 5 (HDF5) ファイルである。HDF5は、1998年に開発された単一のソフトウェアライブラリを用いてのみ読み書き可能な、大規模データ保存用の汎用ファイル形式である。分析の結果、以下のことがわかった。(1) CPUの並列スレッド数を増やすと、典型的なメチル化コールジョブの実行時間が比較的短くなること (論文Extended Data Fig. 1a) 、 (2) これは非効率なデータ処理ではなく、非効率なデータアクセス(ファイル読み込み)によること (論文Extended Data Fig. 1b) 、(3)根本的なボトルネックは、HDF5ファイルを読み込むためのソフトウェアライブラリの制限で、複数のCPUスレッドからの並列入出力(I/O)要求が直列化されてしまい、並列CPUリソースを効率的に使用できないことに原因があった(論文Extended Data Fig.)。

 並列計算は大規模データセットのスケーラブルな解析を可能にし、現代のゲノミクスの中心となっている。しかし残念ながら、FAST5フォーマットには固有の非効率性があり、高度なHPCシステムを利用しても、ナノポアシグナルデータの解析は法外に遅くなるであろう(論文図1b)。例えば、オーストラリアのNational Computing Infrastructure(世界最大級の学術用スーパーコンピュータ、論文補足表2-HPC-Lustre参照)で利用可能な最大リソース割り当てでは、約30×ヒトゲノムデータセットのゲノムワイドDNAメチル化プロファイリングに14日以上かかった。さらに、全体の実行時間の大部分(90%以上)がFAST5ファイルの読み込みに費やされていることを考えると、さらなるソフトウェアの最適化による性能向上は、ファイルの読み込みにかかる時間に比べてわずかであると思われる。

 そこで、FAST5形式固有の制限を克服するために、ナノポアシグナルデータの効率的かつスケーラブルな解析のために設計されたファイル形式、SLOW5を作成した(論文図1b)。SLOW5は、FAST5に含まれるすべての情報をエンコードしているが、FAST5ファイルの読み込みに必要なHDF5ライブラリには依存していない。SLOW5フォーマットの可読バージョンはタブ区切り値(TSV)ファイルで、1行につき1つのナノポアリードのメタデータと時系列信号データをエンコードし、グローバルメタデータはファイルヘッダに格納されている(論文表1および補足説明3)。並列ファイルへのアクセスは、SLOW5ファイル内の各リードの位置(バイト単位)を指定するバイナリインデックスファイルによって容易になる(論文補足注3)。SLOW5は、人間が読めるASCII形式か、配列アラインメントを保存するための代表的なSAM/BAM形式に類似したコンパクトで効率的なバイナリ形式、BLOW5でエンコードすることができる。バイナリ形式は、zlibおよび'vbz' (Z-standard + StreamVByte) アルゴリズムによる圧縮をオプションでサポートしており、並列アクセスを可能にしながらストレージフットプリントを最小化する。

中略

FAST5データアクセスの非効率性は、遅延と費用を生み出し、研究および臨床ゲノミクスの多くのアプリケーションにおけるONTシーケンスの実現可能性を制限している。また、このような摩擦は、ナノポアシグナルデータに直接アクセスするバイオインフォマティクスソフトウェアの開発を阻む要因にもなっている。これは、2009年に開発された、シンプルで効率的なオープンソースのSAM/BAM配列アライメントフォーマットとは対照的である。

 SLOW5フォーマットは、ナノポアシグナルデータを効率的かつ並列に解析するためのフレームワークを提供し、あらゆる用途に対応する。SLOW5の読み書きは、C言語(slow5lib)およびPython(pyslow5)言語用の効率的なソフトウェアアプリケーションプログラミングインターフェース(API)によって管理されている(論文Methods)。これにより、既存のFAST5 APIを置き換えることで、既存のパッケージを含むサードパーティソフトウェアへのSLOW5の統合が容易になる。注目すべきは、サードパーティソフトウェアSigmapがSLOW5を採用するために必要なコードはわずか〜70行であるのに対し、同じツール内のFAST5アクセスに必要なコードは〜2,600行である点であろう。これは、FAST5を読むために必要なHDF5ライブラリに依存しない、完全なオープンソースであるSLOW5 APIのシンプルさを示している。SLOW5形式のシンプルで直感的な構造とともに、ナノポアデータ解析のための活発でオープンなソフトウェア開発をサポートする。

 

Githubより

Slow5toolsは、SLOW5形式のデータの変換(FAST5 <-> SLOW5)、圧縮、閲覧、索引付け、操作のためのシンプルなツールキットです。
SLOW5は、オックスフォード・ナノポア・テクノロジーズ(ONT)社のデバイスからのシグナルデータを保存するための新しいファイルフォーマットです。SLOW5は、効率的でスケーラブルな解析を妨げ、開発者にとって多くの頭痛の種となっている標準的なFAST5シグナルデータフォーマットの固有の制限を克服するために開発されました。SLOW5は、人が読めるASCIIフォーマット、またはよりコンパクトで効率的なバイナリフォーマット(BLOW5)でエンコードすることができます - これはDNA配列アラインメントを保存するための代表的なSAM/BAMフォーマットに類似しています。BLOW5バイナリ形式は、zlib(DEFLATE)圧縮やその他の圧縮方法をサポートし、効率的な並列アクセスを可能にしながらも、データの保存領域を最小化することができます。詳細なベンチマーク実験により、SLOW5形式はFAST5形式と比較して1桁以上速く、サイズも大幅に小さいことを示しています。

 

 

slow5tools - Full documentation

https://hasindu2008.github.io/slow5tools/getting_started.html

SLOW5 Specification(v0.2)

SLOW5 Specification | slow5specs

 

インストール

ubuntu18でテストした。

Github

#conda (bioconda)
mamba install -c bioconda slow5tools -y

#binary(VERSION=v0.3.0)
wget "https://github.com/hasindu2008/slow5tools/releases/download/$VERSION/slow5tools-$VERSION-x86_64-linux-binaries.tar.gz" && tar xvf slow5tools-$VERSION-x86_64-linux-binaries.tar.gz && cd slow5tools-$VERSION/
./slow5tools

> slow5tools -h

Usage: slow5tools [OPTIONS] [COMMAND] [ARG]

Tools for using slow5 files.

 

OPTIONS:

    -h, --help       Display this message and exit.

    -v, --verbose    Verbosity level.

    -V, --version    Output version information and exit.

 

COMMANDS:

    f2s or fast5toslow5   convert fast5 file(s) to SLOW5/BLOW5

    s2f or slow5tofast5   convert SLOW5/BLOW5 file(s) to fast5

    merge                 merge SLOW5/BLOW5 files

    split                 split SLOW5/BLOW5 files

    index                 create a SLOW5/BLOW5 index file

    get                   display the read entry for each specified read id

    view                  view the contents of a SLOW5/BLOW5 file or convert between different SLOW5/BLOW5 formats and compressions

    stats                 prints statistics of a SLOW5/BLOW5 file to the stdout

    cat                   quickly concatenate SLOW5/BLOW5 files of same type [experimental]

    quickcheck            quickly checks if a SLOW5/BLOW5 file is intact

 

ARGS:    Try 'slow5tools [COMMAND] --help' for more information.

レポジトリでは自分の計算機環境でビルドする手順(必要なライブラリの導入方法や対応linuxなども)も説明されています。

 

 

実行方法  (manualの通り)

fast5ディレクトリをBLOW5ファイルに変換する

slow5tools f2s fast5_dir -d blow5_dir

出力例

input

f:id:kazumaxneo:20220107141529p:plain

output

f:id:kazumaxneo:20220107141546p:plain

非常に高速に動作する。テストデータ(80MB fast5 => BLOW5変換では0.5秒かからなかった)

 

1 つの fast5 ファイルを SLOW5 ASCII に変換する

slow5tools f2s file.fast5 -o file.slow5

 

fast5ファイルのディレクトリをBLOW5ファイルに変換する

slow5tools f2s fast5_dir -d blow5_dir -c zstd -s svb-zd

 

ディレクトリ内の全てのBLOW5ファイルを1つのBLOW5ファイルにマージする

slow5tools merge blow5_dir -o file.blow5

 

BLOW5 ファイルを SLOW5 ASCII に変換して標準出力

slow5tools view file.blow5

 

 BLOW5ファイルをSLOW5 ASCIIに変換

slow5tools view file.blow5 -o file.slow5

 

SLOW5 ファイルを BLOW5 に変換

slow5tools view file.slow5 -o file.blow5

 

 Blow5 ファイルのディレクトリを fast5 に変換

slow5tools s2f blow5_dir -d fast5

出力例(上の画像と同じデータ)

f:id:kazumaxneo:20220107141640p:plain

 

slow5/blow5 ファイルをインデックス化する

slow5tools index file.blow5

 

コメント

論文中では、FAST5から圧縮BLOW5に切り替えることで小さなデータでも読み込みが数十倍高速化された事、GPUと組み合わせると30×ヒトゲノムのメチル化プロファイリングが数十倍高速化された事、Amazon AWSやポータブルコンピューティング用の小型デバイスなどどのようなシステムでも改善が見られたこと、PromethION内蔵コンピューティングシステムに導入してシークエンシング中にリアルタイムにFAST5=>SLOW5変換できたこと(つまり待ち時間ゼロ)などが書かれています。

引用
Fast nanopore sequencing data analysis with SLOW5
Hasindu Gamaarachchi, Hiruna Samarakoon, Sasha P. Jenner, James M. Ferguson, Timothy G. Amos, Jillian M. Hammond, Hassaan Saadat, Martin A. Smith, Sri Parameswaran & Ira W. Deveson

nature  nature biotechnology  brief communications, Published: 03 January 2022

 

本当かどうかわかりませんが、SRAもこのフォーマットになるという話も出ています。

メタゲノムデータに適用可能な機械学習モデル SignalP 6.0

 

 シグナルペプチド(SP)は、すべての生物において、タンパク質の分泌や移動を制御する短いアミノ酸配列である。SPは配列データから予測することができるが、既存のアルゴリズムでは既知のSPの種類を全て検出することはできない。本稿では、5種類のSPをすべて検出し、メタゲノムデータに適用可能な機械学習モデル、SignalP 6.0を紹介する。

 

 SPは、真核生物では分泌経路(Sec)、原核生物では細胞膜(内膜)を通過するタンパク質の標的となるN末端の短いアミノ酸配列である。実験的にSPを網羅的に同定することは困難であるため、計算機によるSPの予測は細胞生物学の研究にとって重要な意味を持つ。SP予測ツールは、一般的な分泌経路やtwin-arginine translocation (Tat) 経路をたどるタンパク質を同定し、シグナルペプチダーゼ (SPase) がSPを切断する配列上の位置を予測することができる。SignalP 5.0は、SPase I (Sec/SPI) またはSPase II (Sec/SPII, 原核生物リポタンパク質) によって切断されるSec基質、およびSPase I (Tat/SPI) によって切断されるTat基質を予測することができる.。しかし、アノテーションデータがないため、SignalP 5.0はSPase IIによって切断されたTat基質やSPase III (prepilin peptidase, 時にはSPase IV2とも呼ばれる) によって処理されたSec基質を検出することができない。このようなSec/SPIII SPは、原核生物において接着、運動、DNAの取り込みに重要な役割を果たすIV型ピリン様タンパク質のトランスロケーションを制御していることが分かっている。さらに、SignalP 5.0は、SPの生物学的機能の根幹をなす部分領域(N末端のn領域、疎水性のh領域、C末端のc領域)を定義できないため、SP構造に関して不可知論的であった。

 SignalP 6.0は、全生物領域にわたる数百万の未注釈タンパク質配列の情報を用いたタンパク質言語モデル(LM)6、7、8、9をベースにしている。LMは、タンパク質の生物学的特性や構造を理解するための意味的な表現を作成する。SignalP 6.0は、これらのタンパク質表現を用いて、従来のバージョンでは検出できなかった新たなタイプのSPを予測し、モデル作成に使用したタンパク質と遠縁のタンパク質や起源不明のメタゲノムデータをより適切に外挿することができる。さらに、SPのサブリージョンを特定することも可能である。

 

Github


webサービス(オンライン)

https://services.healthtech.dtu.dk/service.php?SignalP-6.0

f:id:kazumaxneo:20220105205416p:plain

mirror;https://dtu.biolib.com/SignalP-6

 

(HPの説明より)SignalP 6.0サーバーは、古細菌グラム陽性菌グラム陰性菌、真核生物のタンパク質におけるシグナルペプチドの存在とその切断部位の位置を予測します。BacteriaとArchaeaでは、SignalP 6.0は以下の5種類のシグナルペプチドを識別することができます。

  • Sec/SPI: Sec/SPI:Secトランスロコンによって輸送され、シグナルペプチダーゼI (Lep) によって切断される「標準的」な分泌シグナルペプチド
  • Sec/SPII:Secトランスロコーンで輸送され、シグナルペプチダーゼII(Lsp)で切断されるリポ蛋白質シグナルペプタイド
  • Tat/SPI。Tatトランスロコンにより輸送され、シグナルペプチダーゼI(Lep)により切断されたTatシグナルペプチド
  • Tat/SPII: Tatリポ蛋白のシグナルペプチドがTatトランスロコンにより輸送され、シグナルペプチダーゼII(Lsp)により切断されたもの。
  • Sec/SPIII:ピリンおよびピリン様シグナルペプチドがSecトランスロコンを介して輸送され、シグナルペプチダーゼIII(PilD/PibD)によって切断されたもの。

さらに、SignalP 6.0はシグナルペプチドの領域を予測します。シグナルペプチドの種類に応じて、n-、h-、c-領域の位置や、その他の特徴的な領域を予測することができます。

真核生物のタンパク質 シグナルペプチドの有無がタンパク質の局在の全てではないことをお忘れなく! 真核生物のタンパク質の局在についてもっと知りたい場合は、タンパク質細胞内局在予測ツールDeepLocをお試しください。また、シグナルペプチドを持つタンパク質がGPIアンカーを持ち、細胞膜の外側に付着しているかどうかをNetGPIという予測ツールでチェックすることもできます。

 

 タンパク質配列を貼り付けるかアップロードする(10アミノ酸以上)。タンパク質の最大数は5000だが、ユーザーのエントリが100件を超えるとタイムアウトすることがある。

f:id:kazumaxneo:20220105212714p:plain

パラメータを確認してsubmitする。

 

exaamle outputファイルの出力

f:id:kazumaxneo:20220105213100p:plain

f:id:kazumaxneo:20220105213113p:plain

 

ローカルマシンへのインストール

fastモデルとslowモデル(HPによると、fastより精度が高いが6倍ほど遅い)がダウンロードできる。

https://services.healthtech.dtu.dk/software.php

 

感想

論文図2にはv5と比較した時の結果が提示されています。切断部位の予測精度と感度の増加と、配列同一性が低い時の感度と精度(マシューズ相関係数で評価)の増加が著しいようです。

引用

SignalP 6.0 predicts all five types of signal peptides using protein language models

Felix Teufel, José Juan Almagro Armenteros, Alexander Rosenberg Johansen, Magnús Halldór Gíslason, Silas Irby Pihl, Konstantinos D. Tsirigos, Ole Winther, Søren Brunak, Gunnar von Heijne & Henrik Nielsen

nature biotechnology, Brief Communication, Published: 03 January 2022

 

関連


複数の生物をサポートする機能的エンリッチメント解析ツール GeneSCF

 

 ChIP-sequencing、RNA-sequencing、DNA sequencing、定量的メタボロミクスなどのハイスループット技術により、膨大な量のデータが生成される。研究者は、これらのハイスループット研究から影響を受けた遺伝子の生物学的意義を解釈するために、しばしばfunctional enrichment toolsに依存している。しかし、現在利用可能なfunctional enrichment toolsは、機能データベースリポジトリからの新しいエントリに適応するために頻繁に更新される必要がある。そこで、KEGG、Reactome、Gene Ontology などのソースデータベースから直接更新された情報を使用して、機能エンリッチメント解析を実行できる簡便なツールが必要とされている。

 本研究では、GeneSCF (Gene Set Clustering based on Functional annotations) と呼ばれるコマンドラインツールを設計し、遺伝子セットの機能関連情報をリアルタイムで予測することに焦点を当てた。このツールは、KEGG, Reactome, Gene Ontologyなどの著名な機能データベースから、4000以上の生物の情報を扱うことができるように設計されている。本ツールを2つの公開データセットに適用し、生物学的に関連する機能情報を予測することに成功した。

 GeneSCFは、リアルタイムに参照機能データベースを用いてエンリッチメント解析を行うことができるため、他のエンリッチメントツールと比較して信頼性が高いと言える。また、ハイスループット・データの下流解析に利用できる他のパイプラインとの統合が容易なツールである。さらに重要なことは、GeneSCFは複数の遺伝子リストを異なる生物種で同時に実行することができるため、ユーザーの時間を節約することができる。このツールはすぐに使えるように設計されているため、複雑なコンパイルやインストール手順は必要ない。

 

特徴

  • 多くの生物種に対応。
  • 複数のソースデータベースを用いて、複数の遺伝子リストのエンリッチメントを一度に解析することが可能(GO、KEGGのみ対応)。
  • GOの語彙/パスウェイ/ファンクションと関連する遺伝子をシンプルな表形式でプレーンテキストファイルでダウンロードできる。

 

インストール

ubuntu18でテストした。

依存

GeneSCFはLinuxシステムでのみ動作し、Ubuntu, Mint, Cent OSで正常に動作することが確認されている。他のLinuxディストリビューションでも動作する可能性はある。

For graphical output or plots

  • R (version > 3.0) and 'ggplot2' is required

Github

git clone https://github.com/genescf/GeneSCF.git
cd GeneSCF/
export PATH=$PATH:$PWD

> ./geneSCF -h

GeneSCF USAGE: 

 

./geneSCF -m=[update|normal] -i=[INPUT_PATH/INPUT_FILE] -t=[gid|sym] -o=[OUTPUT_PATH/FOLDER/] -db=[GO_all|GO_BP|GO_MF|GO_CC|KEGG|REACTOME|NCG] -p=[yes|no] -bg=[#TotalGenes] -org=[see,org_codes_help]

 

==========

Options: ALL PARAMETERS MUST BE SPECIFIED

==========

 

[-m= | --mode=]        For normal mode use 'normal' and for update mode use 'update' without quotes.

 

[-i= | --infile=]    Input file contains list of Entrez GeneIDs or OFFICIAL GENE SYMBOLS.

            The genes must be new lines seperated (One gene per line).

 

[-t= | --gtype=]    Type of input in the provided list either Entrez GeneIDs 'gid'

            or OFFICIAL GENE SYMBOLS 'sym' (Without quotes, Example for 

            human 'sym' => HUGO gene symbols).

 

[-db= | --database=]    Database you want to find gene enrichment which is either 

            geneontology 'GO_all' or geneontology-biological_process 

            'GO_BP' or geneontology-molecular_function 'GO_MF' or 

            geneontology-cellular_components 'GO_CC' or kegg 'KEGG' or 

            reactome 'REACTOME' or Network of Cancer Genes 'NCG' (Without quotes).

 

[-o= | --outpath=]    Existing directory to save output file. The output will be with saved in the 

            provided location as {INPUT_FILE_NAME}_{database}_functional_classification.tsv 

            (tab-seperated file).

 

[-bg= | --background=]    Total background genes to consider (Example : ~20,000 for human).

 

[-org= | --organism=]    Please see organism codes(For human in KEGG ->'hsa' in Geneontology -> 'goa_human').

            For REACTOME and NCG use 'Hs' (human).

 

[-p= | --plot=]        For additional graphical output use 'yes' or 'no'.This requires R version > 3.0 and 

            'ggplot2' R package to be pre-installed on the system.

 

[-h | --help]        For displaying this help page.

 

 

prepare_database USAGE: 

 

./prepare_database -db=[GO_all|GO_BP|GO_MF|GO_CC|KEGG|REACTOME|NCG] -org=[see,org_codes_help]

 

==========

Options: ALL PARAMETERS MUST BE SPECIFIED

==========

 

[-db= | --database=]    Options:[GO_all|GO_BP|GO_MF|GO_CC|KEGG|REACTOME] 

 

[-org= | --organism=]    Options:[see,org_codes_help]

> ./prepare_database -h

prepare_database USAGE: 

 

./prepare_database -db=[GO_all|GO_BP|GO_MF|GO_CC|KEGG|REACTOME|NCG] -org=[see,org_codes_help]

 

==========

Options: ALL PARAMETERS MUST BE SPECIFIED

==========

 

[-db= | --database=]    Options:[GO_all|GO_BP|GO_MF|GO_CC|KEGG|REACTOME] 

 

[-org= | --organism=]    Options:[see,org_codes_help]

 

 

 

データベースの準備

デフォルトでヒト用の遺伝子オントロジーKEGG、Reactome、NCGからなるデータベースが用意されているが、ほかのモデル生物だと用意する必要がある。そのためのコマンドが準備されている。

#GO
prepare_database -db=GO_all -org=goa_human

#KEGG
prepare_database -db=KEGG -org=hsa

#REACTOME
prepare_database -db=REACTOME -org=Hs

#NCG (catalogue of known and candidate cancer genes)
prepare_database -db=NCG -org=Hs

コマンドを実行すると、指定したデータベースがダウンロードされる(GeneSCF/class/lib/db/)。

f:id:kazumaxneo:20220104210428p:plain

-p=yesのオプションは機能しなかった。マニュアルの間違い?

 

 

実行方法

遺伝子リストに対してエンリッチメント解析を行う。

必要な遺伝子リストは、公式遺伝子シンボル(例:ヒトのHGNC、--gtype=sym)またはEntrez GeneID(--gtype=gid)形式の遺伝子リストを含むプレーンテキスト。入力遺伝子リストは、発現変動遺伝子だったりエピジェネティック修飾に異常がある遺伝子などの個々の研究において影響を受けると予測される遺伝子リストとなる。このほか、ランには、-mで指定するupdateモードとノーマルモード(-modeパラメータ)のどちらかを指定する必要がある。updateモードでは、対応するデータベースと生物からリアルタイムに機能情報を取得してエンリッチメント解析を行う。通常モードでは、updateモードにより過去に取得・利用した機能データベースに対してエンリッチメント解析を行う。

#テストラン
cd GeneSCF/

mkdir outdir1
geneSCF -m=update -i=test/H0.list -o=outdir1/ -t=sym -db=GO_all -bg=20000 --plot=yes -org=goa_human
  • -m=normal | update    For normal mode use normal and for update mode use *update* without quotes
  • -i= <INPUT-TEXT-FILE>    Input file contains list of Entrez GeneIDs or OFFICIAL GENE SYMBOLS. The genes must be new lines seperated (One gene per line)
  • -t=gid | sym    Type of input in the provided list either Entrez GeneIDs gid or OFFICIAL GENE SYMBOLS sym (default: *gid*)
  • -db=GO_all | GO_BP | GO_CC | GO_MF | KEGG | REACTOME | NCG
  • -o= <OUTPUT-DIRECTORY>   Existing directory to save output file (Don't forget to use trailing slash at end of the directory name). The output will be saved in the provided location as {INPUT_FILE_NAME}_{database}_functional_classification.tsv (tab-seperated file). 
  • -bg= <Total background>    Total number of genes to be considered as background (Example : ~20,000 for human). It is important to choose the background appropriately. Sometimes your samples do not express all the genes. For example, if you are using differentially expressed genes for gene set enrichment analysis, you must choose total number of genes detected in your experiment including control and treatment samples irrespective of their differential status as your background (NOT the total number of genes in the annotation of the corresponding organism you are working with). All the genes from the annotation can be used when working with genes found in genome-wide studies (example, ChIP-seq, WGBS, etc.,).
  • -org=<see organism codes>
  • -p=yes|no    For additional graphical output use yes or no.This requires R version > 3.0 and ggplot2 R package to be pre-installed on the system

GeneSCF は、分割表(論文図1C)から、PERL の統計モジュールを利用し、標準的な多重検定補正法を用いた Fisher's Exact 検定を用いて、エンリッチメントの有意性を計算する。多重仮説補正法には、Benjamini and Hochberg、Hommel一段階処理、Bonferroni一段階処理、Hochbergステップアップ処理、Benjamini and Yekutieli多重仮説補正係数がある。GeneSCF の出力は、ユーザが提供した遺伝子リストからのヒット数に基づいてランク付けされた関数のリストがタブ区切りのテーブルになったもの。また、分割表を用いたフィッシャーの正確検定によるp値と、異なる多重仮説補正法によるFDR)値の列が含まれる(論文より)。

出力例

f:id:kazumaxneo:20220104230720p:plain

H0.list_GO_all_goa_human_functional_classification.tsv

f:id:kazumaxneo:20220104230919p:plain

H0.list_GO_all_goa_human_enrichment_plot.pdf

f:id:kazumaxneo:20220104230737p:plain

KEGG

mkdir outdir2
mkdir -p mapping/DB/
geneSCF -m=normal -i=test/H0.list -o=outdir2/ -t=sym -db=KEGG -bg=20000 --plot=yes -org=hsa

出力例

f:id:kazumaxneo:20220104231031p:plain

 

f:id:kazumaxneo:20220104231056p:plain

 

Githubではバッチモードの使い方についての説明もあります。アクセスして確認して下さい。

引用

GeneSCF: a real-time based functional enrichment tool with support for multiple organisms
Santhilal Subhash & Chandrasekhar Kanduri 
BMC Bioinformatics volume 17, Article number: 365 (2016) 

 

 

動的に生成されるRスクリプトを用いてバルクRNA-seqの自動探索と可視化を行う Searchlight2

2022 1/5 複数比較の例追記、コマンドの誤字修正

 

 バルクRNA-seqデータが処理されると、すなわちアラインメントされ、発現および差分表が作成されると、生物学的性質の探索、視覚化および解釈が行われる重要なプロセスが残る。可視化・解釈パイプラインを使用しないと、このステップは時間と手間がかかり、多くの場合Rを使用して完了する。市販の可視化・解釈パイプラインは包括的だが、自由に使用できるパイプラインは現在限られている。ここでは、自由に利用できるバルクRNA-seq可視化および解釈パイプラインであるSearchlightを紹介する。Searchlightは、グローバル、パスウェイ、単一遺伝子レベルに焦点を当てた包括的な統計および視覚的分析、3つのワークフローによる生物または実験の複雑さに関係なくほとんどの差分実験デザインとの互換性、レポート、およびユーザーフレンドリーなRスクリプトとShinyアプリによるプロットの下流ユーザー修正のサポートなどを提供する。Searchlightは、現在のベストツール(VIPERとBioJupies)よりも優れた自動化を提供することを示す。標準的なバルクRNA-seq処理パイプラインと並行して、Searchlightを使用すると、論文原稿品質の図まで含めてバルクRNA-seqプロジェクトやtimed re-analysis studyを3時間未満で完了できることを実証した。

 

Githubより

 バルクRNA-seqデータの処理、すなわちアライメント、発現および発現変動遺伝子のテーブルの作成が完了すると、生物学の探究、視覚化、解釈を行う重要なプロセスが残ります。最終的には、報告書、論文、原稿の図表を作成します。この下流工程を完了するための典型的な方法は、手動でコード化したコマンドラインとRベース(または類似の)の分析です。しかし、この方法では、数日から数週間を要することもあり、大変な労力を要します。

 Searchlight2は、バルクRNA-seqの探索、可視化、解釈のパイプラインで、下流域の解析を自動化することを目的としています。標準的なアライメントおよび処理パイプライン(Star2、Hisat2、Kallisto、DEseq2、EdgeRなど)と共に使用すれば、研究者はバルクRNA-seqプロジェクトを数時間のうちに、。Rによる手動解析と区別がつかないほどの水準で、最小限の労力で完了させることが可能です。

複雑なパイプラインを使用したり、理解したりする必要はありません。その強みは

  • 強力で広く利用されている幅広い分析・可視化手法。
  • 発現、差分発現、シグネチャー解析をカバーする独立したワークフローを使用できる。これらは、デザイン、複雑さ、生物学、生物に関係なく、実験との互換性を提供する。また、解析を簡素化する。
  • 包括的なレポートを作成。
  • すべてのプロットにRとR Shinyを使用しているため、バイオインフォマティシャンやウェットラボの科学者がすべてのプロットを視覚的に簡単に変更することができる。
  • R-snippetデータベースにより、ユーザーはすべてのプロットのデフォルトの外観を自分の好みに合わせて変更することができる。
  • Searchlight2は100%自動化されている。
  • Searchlight2は、サンプルシート、発現マトリックス、任意の数の発現差テーブルなど、典型的なRNA-seqダウンストリーム解析の入力を受け付ける。
  • バイオインフォマティシャン、RNA-seqサービスプロバイダー、ベンチサイエンティストが最小限の労力でバルクRNA-seq研究プロジェクトを迅速に進め、さらなる詳細解析や別の解析アプローチにリソースを解放できるように設計されている。
  • Searchlight2はアライメントやリードのカウント、発現量や差分発現量の計算を行わないため、処理パイプラインではないことに留意する。これらのステージは、Searchlight2を使用する前に完了している必要がある。
  • 発現量の行列(TPM、RPKM、Rlogなど)と少なくとも1つの差分発現テーブル(DEseq2、EdgeRなど)があれば、どんな処理パイプラインでもSearchlight2を使用できる。

 

インストール

ubuntu18.04でcondaを使ってR4,0の環境を作ってテストした。

Github

#依存
mamba create -n searchlight2 python=2.7 -y
conda activate searchlight2
mamba install Scipy Numpy -y
mamba install -c bioconda r -y

#Rの依存
install.packages("ggplot2")
install.packages("reshape")
install.packages("amap")
install.packages("amap")
install.packages("gridExtra")
install.packages("gtable")
install.packages("ggally")
install.packages("network")
install.packages("sna")

#shinyを使う場合(一番下で紹介)の依存;Rstudio上でインストールする
install.packages("shiny")
install.packages("shinyFiles")
install.packages("fs")
install.packages("shinycssloaders")
install.packages("graphics")
install.packages("dplyr")
install.packages("GGally")

#Searchlight2本体
cd Searchlight2-master/software/

python Searchlight2.py -h

#####################################

#####       Searchlight 2       #####

#####################################

 

Automated bulk RNA-seq data exploration and visualisation using dynamically generated R-scripts

John J. Cole & Carl S. Goodyear, et al.

University of Glasgow (Scotland) 2021

 

v2.0.0

option -h not recognized

 

 

 

実行方法

Searchlight2は1つのコマンドで実行される。まず、入力ファイルを検証し、それらを1つの「マスター遺伝子テーブル」に結合し、そこから下流の解析が行われる。次に、各ワークフローを繰り返し、中間ファイル、統計解析結果ファイル、プロットごと、ワークフローごとのRスクリプト、プロット、HMTLによるレポート、そして最後にShinyアプリが生成される。

Searchlight2のランには最小構成で4つのファイルが必要。入力のフォーマットには厳密で、ファイル間の対応関係が違うとエラーを起こす。例えば、ファイル間でサンプル名が微妙に間違っているケースは許容されない。

 

必要なファイル1- 発現量のマトリックスTPM、RPKM、Rlogなど)

最初の列は遺伝子ID(Ensembl、Refseqなど)とする。1行目は、最初のセルが"ID "、残りはサンプル名とするヘッダー行。サンプル名は数字で始めることはできず、数字、文字、アンダースコア(_)のみ使用可能。

レポジトリの例

f:id:kazumaxneo:20220104021448p:plain

 

必要なファイル2- Differential expression table

DESeq2やEdgeRなどを使った標準的な発現変動遺伝子のテーブル。遺伝子は行単位で、列は遺伝子ID、log2 fold change、p-value、adjusted p-value(この順)のみに切り詰めたもの。1行目にヘッダ行が必要で、ヘッダは正確に”ID”、”Log2Fold”、”p”、”p.adj”でなければならない。大文字・小文字は区別しない、引用符は無視する。IDの型は、発現マトリックスと同じでなければならない。

レポジトリの例

f:id:kazumaxneo:20220104022039p:plain

 

必要なファイル3- サンプルシート (SS)

標準的なタブ区切りのサンプルシートで、1列目に各サンプル名、2列目に所属するサンプルグループを記載する。ヘッダ行があり、ヘッダは正確に"sample", "sample_group"でなければならない。大文字と小文字の区別はなく、引用符も無視される。サンプルグループの階層が複数ある場合(細胞の種類、治療法、年齢など)、3列名以降に、任意のヘッダー(先頭行)から始めてさらにグループを記載できる。サンプル名とサンプルグループ名は数字で始めることはできず、数字、文字、アンダースコア(_)のみを使用することができる。

f:id:kazumaxneo:20220104022553p:plain

 

必要なファイル4- バックグラウンドファイル(BG)

全遺伝子がリストアップされている典型的なバックグラウンドアノテーションファイル。EnsemblsのBiomartから簡単に作成することができる。遺伝子は行単位で、特定のアノテーションは列単位で記述する。ファイルには以下のアノテーションが必要。Gene ID、Gene Symbol、Chromosome、Start position、Stop position、Biotype(遺伝子タイプ)。ヘッダ行には、”ID”、”Symbol”、”Chromosome”、”Start”、”Stop”、”Biotype”ヘッダを正確に記述する。大文字と小文字の区別はなく、引用符は無視する。遺伝子記号が不明な場合は、IDを使用する。Biotypeがわからない場合は、その列のすべてのセルに "gene "と記入する。

f:id:kazumaxneo:20220104095157p:plain

 

任意で必要なファイル5- 遺伝子セットデータベースファイル

GMT形式の有効な遺伝子セットデータベースファイル(GO、string、KEGGなど)を指する(例;Searchlight2/gene_set_databases/)。1つの遺伝子セットにつき1行で、最初のセルが遺伝子セット名(細胞周期など)、2番目のセルが任意のメモ(ない場合は単にNAと入力)、そして遺伝子セット内の各遺伝子で構成されている。

GOデータベースの例

https://github.com/Searchlight2/Searchlight2/blob/master/gene_set_databases/GO_bp_human.tsv

f:id:kazumaxneo:20220104100157p:plain

 

テストラン

必要な4つのオプションを指定する。過剰発現解析も行うなら--oraオプションで必要なファイル5を指定する。--deと--oraは複数のパラメータを指定する。複数パラメータ書くオプションは、パラメータ間にスペースを入れてはならない

python Searchlight2-master/software/Searchlight2.py --out path=results \
--bg file=Searchlight2/backgrounds/mouse_GRCm38.p6.tsv \
--em file=Searchlight2/sample_datasets/EM.tsv \
--ss file=Searchlight2/sample_datasets/SS.tsv \
--de file=Searchlight2/sample_datasets/WT_vs_KO.tsv,numerator=KO,denominator=WT
--ora file=Searchlight2/gene_set_databases/GO_bp_mouse.tsv,type=GO_bp

(注;すべてのパスは「root」からの完全なものでなければならない)

過剰発現解析が必要ないなら--oraオプションを消す。上流調節因子解析も含める場合は、コマンドに --ura パラメータを追加する。

 

出力例

results/

f:id:kazumaxneo:20220103224136p:plain

results/all_genes/

f:id:kazumaxneo:20220103224216p:plain

results/all_genes/de_workflows/KO_vs_WT/

f:id:kazumaxneo:20220103224252p:plain

 

report.htmlを開く。

 

f:id:kazumaxneo:20220103224104p:plain

f:id:kazumaxneo:20220103234019p:plain

f:id:kazumaxneo:20220103234030p:plain

f:id:kazumaxneo:20220103234045p:plain

f:id:kazumaxneo:20220103234059p:plain

f:id:kazumaxneo:20220103234109p:plain

f:id:kazumaxneo:20220103234116p:plain

f:id:kazumaxneo:20220103234136p:plain

f:id:kazumaxneo:20220103234145p:plain

f:id:kazumaxneo:20220103234158p:plain

f:id:kazumaxneo:20220103234214p:plain

f:id:kazumaxneo:20220103234228p:plain

f:id:kazumaxneo:20220103234240p:plain

f:id:kazumaxneo:20220103234306p:plain

f:id:kazumaxneo:20220103234340p:plain

f:id:kazumaxneo:20220103234403p:plain

f:id:kazumaxneo:20220103234422p:plain

 

 

Rに不慣れなユーザーのためにR shinyの環境も用意されている。利用するには出力のshiny/下にあるserver.rをRstudioで開く。

f:id:kazumaxneo:20220103222344p:plain

server.rが読み込まれたら、Rstudio左上のパネルの右上にあるスタートボタンをクリックする。

f:id:kazumaxneo:20220103222516p:plain

 

ブラウザが立ち上がる。

f:id:kazumaxneo:20220103222701p:plain

 

遺伝子、ワークフロータイプを選択する。

f:id:kazumaxneo:20220103222758p:plain

分析方法を選択するとそのまま視覚化される。

f:id:kazumaxneo:20220103222835p:plain

PCA;各主成分の寄与率の推移

f:id:kazumaxneo:20220103222916p:plain

 

作図パラメータは細かい部分まで調整可能。図のサイズを大きくした。

f:id:kazumaxneo:20220103223405p:plain

 

レポートhtml付きサンプル解析結果がsample_datasets/results.zipに含まれています。興味ある方はアクセスしてみて下さい。

Rの依存関係が正しく導入されていない場合、各図の出力には失敗しますが、Searchlight2の処理は最後まで実行されます。その場合、全体のログを見てもどこが失敗したのかわかりにくいわけですが、Searchlightは出力の各ディレクトリにRの作図コードが含まれているので、Rで実行してどんなエラーが出るか確認してみると原因ははっきりします(例えばcd xxx/xxx/PCA/ してRを立ち上げ、source("pca.R"))。

 

 

2022 1/5 追記

2種類以上のサンプルを持つために複数の比較が必要な場合、Searchlight2はこれらのデータセットを完全にサポートしている。これらは、単に「アップ」または「ダウン」制御を超えたシグネチャを含むことができるため、「シグネチャベース」と呼ばれている。

 

独立した2群間比較を複数ペアで行う。

複数の差分発現ファイル(例:WT vs KO、KO vs KO_rescue)がある場合、その分だけ--deパラメータを追加するだけです。これだけで他の比較も同時に実行できます。

サンプルデータセットだとWT vs KO、KO vs KO_rescueで比較したいので、2つ--deを指定しています。以下のコマンドになります。

python Searchlight2/software/Searchlight2.py --out path=results \
--bg file=Searchlight2/backgrounds/mouse_GRCm38.p6.tsv \
--em file=Searchlight2/sample_datasets/EM.tsv \
--ss file=Searchlight2/sample_datasets/SS.tsv \
--de file=Searchlight2/sample_datasets/WT_vs_KO.tsv,numerator=KO,denominator=WT \
--de file=Searchlight2/sample_datasets/KO_rescue_vs_KO.tsv,numerator=KO_rescue,denominator=KO

出力例

f:id:kazumaxneo:20220105103606p:plain

 

Githubでは、このような独立したいくつかの2群間比較に加えて、3つのグループを一緒に調べる手順についても説明されています。具体的にはmultiple differential expression (MDE)ワークフローを使用することで実現できるようになっており、-mdeパラメータを使用します。-mdeパラメータは、2つ以上の差分比較を同時に行うことを可能にするものです。コマンドが複雑になりますので、レポジトリのIncluding a formal signature analysisのセクションを読んでみて下さい。

引用

Searchlight: automated bulk RNA-seq exploration and visualisation using dynamically generated R scripts
John J Cole, Bekir A Faydaci, David McGuinness, Robin Shaw, Rose A Maciewicz, Neil A Robertson, Carl  Goodyear

BMC Bioinformatics. 2021 Aug 19;22(1):411

 

参考


Uniprotのパンプロテオーム

明けましておめでとうございます。

今年もよろしくお願い致します。

 

 パンプロテオームとは、closely related (高度に関連した)生物群(例えば、同じ細菌種の複数の株)によって発現されると考えられるタンパク質の完全な集合のこと(panはギリシャ語で"whole")。分類群内の全配列の代表的なセットを提供し、そのグループのリファレンスプロテオームだけでは見られないユニークな配列を捕捉することができる。

Uniprotが提供しているUniProtKBパンプロテオームは、主に種レベルの全ての非冗長プロテオームを網羅したもので、系統比較、ゲノム進化や遺伝子多様性の研究に利用できる。UniProt FTPの Pan proteomes/からダウンロードできるようになっている。

Pan proteomes

https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/pan_proteomes/

f:id:kazumaxneo:20211231231731p:plain

拡大

f:id:kazumaxneo:20211231231810p:plain

表示されているのはUNiprot Proteome identifier(UPID)という識別子UPIDは、UniProtのプロテオームを構成するタンパク質コレクションに割り当てられる一意の識別子で、UP'という文字と9桁の数字からなる。Uniprotの説明には、リリース間で安定しているため、UniProtプロテオームを引用する際に使用可能とある。

 

UPIDのIDはUniprot proteomesで検索できる。

UP00000012と検索。

f:id:kazumaxneo:20211231234241p:plain

リファレンスプロテオームとしてヒットしている。

 

 

関連

 

 

補足

あるプロテオームがより大きなパンプロテオームの一部であるタンパク質を持つ場合、プロテオームページの 'Pan proteome' の行にその旨が表示される。

 

 

 

 

 

ATAC-Seq、ChIP-Seq、WESなどのcDNA汚染の検出と除去を行う cDNA-detector

 

 意図的または偶然に実験システムに導入された外因性cDNAは、そのシステムから得られた次世代シーケンサーライブラリーにおいて、その遺伝子に対するリードカバレッジの追加として現れることがある。適切に認識・管理されない場合、この外来シグナルによるクロスコンタミネーションは、研究結果の解釈を誤らせる原因となる。しかし、この問題は、現在のシーケンス処理パイプラインでは日常的に対処されていない。ここでは、DNAシーケンス実験における外来cDNAの汚染を特定し、除去するための計算ツールであるcDNA-detectorを紹介する。cDNA-detectorがアライメントファイルからcDNAを迅速かつ正確に同定できることを実証した。cDNA-detectorは、検出されたcDNAからアライメントを汚染しないようにするためのメカニズムを提供する。シミュレーションの結果、cDNA-detectorは高感度かつ特異的であり、既存のツールを凌駕していることが示された。また、cDNA-detectorを複数の公共データベース(TCGA、ENCODE、NCBI SRA)に適用し、汚染遺伝子がシーケンシング実験に出現し、誤ったカバレッジピークコールを導くことを示した。

 cDNA-detector は、NGS ライブラリー中の cDNA 検出と除去を行う、使いやすく正確なツールである。この2段階のデザインは、候補を手動で確認することができるため、真のバリアント除去のリスクを低減する。意図的および偶発的に導入されたcDNAによる汚染は、広く利用されているコンソーシアムのデータセットにおいてさえ、過小評価されている問題であり、偽の結果につながる可能性があることが分かった。本発見は、ダウンストリーム解析の前に、NGSライブラリーから汚染されたcDNAを高感度に検出し、除去することの重要性を強調するものである。

 

 ChIP-seq、ATAC-seq、全エキソーム、全ゲノムなどの超並列DNAシーケンス実験が広く研究手法として用いられている。これらの手法は、DNA関連タンパク質の結合部位、ヒストン修飾状態、クロマチンアクセシビリティ生殖細胞および体細胞ゲノム変異の同定に役立っている。一般的な機能研究では、目的の野生型遺伝子やバリアントをDNAベクターで実験系に導入し、シークエンスライブラリーに登場させることができる。目的の外来遺伝子に由来するシーケンスリードは、アライメントステップにおいて内在性遺伝子座にマッピングされる。さらに、同じ研究室で行われた他の実験や共有の装置で得られた微量のベクターDNAがDNAライブラリーを汚染し、アライメントに外来遺伝子が含まれる可能性がある。これらの汚染遺伝子からのシグナルは、例えば偽のコピー数増加、偽の生殖細胞または体細胞バリアントコール、誤ったカバレッジピークコールなど、下流の結果や解釈に影響を与える可能性がある。特に、ベクター由来のシグナルがエクソームやクロマチンプロファイリング由来の真のシグナルと重なる場合、アライメントにおけるこのような外来遺伝子のシグナルの検出は困難な場合がある。NGSライブラリのベクター汚染を検出する方法は、これまでほとんど開発されていない。SeqClean や VecScreenのように、既知のクローニングベクターからのシーケンスリードを同定するものもあるが、cDNA の挿入を検索するものではない。Vecuumは、汚染遺伝子を特定するために開発されたが、このツールは、既知のベクターバックボーン配列に依存している。

 クローン化されたcDNAにはイントロンや生理的UTR領域がないため、この性質を利用して汚染される可能性のある遺伝子を同定することができる。これらのcDNAからのリードは、エキソン境界で部分的にゲノムにアライメントし、マッピングされていない配列はベクター(5'および3'末端)または隣接するエキソンにマッチする。ゲノムアライメントでは、これらのリードは「切り取られた」ように見えるため、ゲノムに完全にマッピングされる真のシグナルリードと区別することができる(論文Fig. 1A)。

(実際の検出手順は論文のImplementationのセクションで説明されています)

 

インストール

導入手順はレポジトリに記載されている。ここでは、著者らによって公開されているdockerイメージを使ってテストした。

Github

#dockerhub
docker pull rheinbaylab/cdna-detector

cdna-detector.py

usage: cdna-detector.py prepare --annotation <gtf/bed> --genome <genome sequence file>  [options]

 

required arguments:

  --annotation          gene annotation files

                        input format: gtf/bed

  --genome              genome fasta file

                        input format: fa/fasta

 

optional arguments:

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

  --format              input annotation format: gtf or bed

                         If format is "bed", only 3 or 4 columns will be used.

                        - default: gtf

  --output_dir          output directory

                        - default: .

  --chr_exclude         exclude chromosomes, multiple chromosomes are separated by ","

                        - conflict with --chr_include

                        - default: chrM

                        - example: --chr_exclude chrX,chrY,chrM

  --chr_include         only include chromosomes, multiple chromosomes are separated by ","

                        - conflict with --chr_exclude

                        - example: --chr_include chr1,chr2

  --source              the program that generated this feature, it is located in 2nd column of gtf. multiple sources are separated by ","

                        - default: all source

                        - example: --source havana

  --feature_type_cds    feature name for identifying CDS regions. It is located in 3rd column of gtf. Multiple features are separated by ","

                        - default: CDS

  --featureid_gene_id   attribute to show gene id. 

                        - default: gene_id

  --featureid_gene_name 

                        attribute to show gene name. 

                        - default: gene_name

  --featureid_cds_ranking 

                        attribute to show exon number.

                        - default exon_number

  --featureid_transcript_id 

                        attribute to show transcript id. 

                        - default transcript_id

 

 

テストラン

cdna-detector detect - あらかじめ定義された遺伝子モデルファイルを用いて、BAM形式の配列アラインメント中のcDNA汚染の可能性を検索する。

解析対象のbamファイル、遺伝子モデル、出力ディレクトリを指定する。

cdna-detector detect --bam example_data/bam_file/tmp.sample.bam --sample_id tmp_sample --gene_model hg19 --output_dir output_dir --n_threads 12
  • --bam   The input file is in BAM format. Recommend results from software bwa
  • --sample_id    Identifier. Used as output prefix
  • --gene_model    Link to gene model.
    - INPUT: hg19 hg38 | mm10 (gene model generated from subcommand "prepare")
  • --output_dir    output directory. - default: "."
  • --n_threads     number of threads. - integer - default: 1

 

入力にcDNAが検出されなかった場合、cdna-detector detectはSAMPLE_ID.logのログファイルのみ生成する。cDNAが検出された場合、cdna-detector detectは、汚染候補exonと遺伝子に関するオリジナルの統計、高信頼性フィルタリングされたexonと遺伝子、クリーンステップで使用するための汚染の疑いがある遺伝子の位置を含む座標ファイルなど、いくつかのファイルを出力する(マニュアルより)。

出力例

output_dir/

f:id:kazumaxneo:20211231112027p:plain

 

cdna-detector clean -  前ステップcdna-detector detectで検出された汚染配列のクリーニング。

cdna-detector detectで生成されたSAMPLE_ID.merge_region.tsvとソースBAMファイルを指定する。

cdna-detector clean --bam example_data/bam_file/tmp.sample.bam  --sample_id tmp_sample --region output_dir/tmp_sample.merge_region.tsv --output_dir output_clean

候補領域から同定したコンタミリードが削除され、クリーンなBAMファイルが生成される。

出力例

output_clean/

f:id:kazumaxneo:20211231121601p:plain

引用

cDNA-detector: detection and removal of cDNA contamination in DNA sequencing libraries

Meifang Qi, Utthara Nayar, Leif S. Ludwig, Nikhil Wagle & Esther Rheinbay 
BMC Bioinformatics volume 22, Article number: 611 (2021)

 

2021年はこの記事で最後にします。

今年もお世話になりました。

来年もよろしくお願い致します。