🌟
Helixerによるゲノムアノテーション

2024/06/23に公開
eukaryote
geneprediction
genomeannotation
tech
 Helixer–de novo Prediction of Primary Eukaryotic Gene Models Combining Deep Learning and a Hidden Markov ModelHelixer – ディープラーニングと隠れマルコフモデルを組み合わせた真核生物遺伝子モデルのde novo予測
https://www.biorxiv.org/content/10.1101/2023.02.06.527280v2.full

 はじめに遺伝子構造予測はゲノム配列から生物学的な知見を得るための重要なステップであるが、一方、ゲノムシーケンスプロジェクトにおいては計算量や複雑な処理が課題となっている。
現在、主流のde novo隠れマカロフモデルは生物学的な複雑性をモデル化する能力に限界があり、採用しているパイプラインはリソース集約型であるものが多い。
加えて、結果の品質はRNAseqのデータやタンパク質の配列情報といった利用可能な外部データの有無と充実度合いに左右される。


Helixerのアーキテクチャ
ここでは、遺伝子予測にDeep learningを適用した以前の論文をもとに、DNA配列のみからプライマリな遺伝子モデルを高速に予測するHelixerを紹介する。
予測精度は他のde novoツールのよりもリファレンスに近いスコアであり、予測結果はそのまま後続の解析に使用したり既存のパイプラインに統合して品質向上に寄与できるレベルである。

 インストールSingularityの使用が推奨されているのでSingularityコンテナを使用する。
https://github.com/weberlab-hhu/Helixer
また、HelixerはGPU利用モードも存在する。使用する場合には以下サイトを参考にその要件を確認しておく。
https://github.com/gglyptodon/helixer-docker
Singularityコンテナはdockerhubから取得する。
singularity pull docker://gglyptodon/helixer-docker:helixer_v0.3.4_cuda_12.2.2-cudnn8
次にHelixerは事前学習モデルをベースに遺伝子領域の推定を行う。著者提供のトレーニングファイルが存在するのでfetch_helixer_models.pyで取得する。
singularity exec --nv helixer-docker_helixer_v0.3.4_cuda_12.2.2-cudnn8.sif fetch_helixer_models.py --all
.local/share/Helixer/modelsに必要なファイルがダウンロードされる。また、トレーニングファイルは以下のサイトからダウンロードされている。
https://zenodo.org/records/10836346
デフォルトで利用可能な系統ごとの事前学習済みのモデルは以下の4つ。
land_plant: 陸生植物
vertebrate: 脊椎動物
invertebrate: 無脊椎動物
fungi: 菌類
特定の系統のモデルを取得したい場合には--lineageで指定する。すべてダウンロードする場合には--allをオプションに与える。
もし上記以外の生物種を対象にする場合にはトレーニングファイルを別途生成するといった方法もある。
https://github.com/weberlab-hhu/Helixer/blob/main/docs/training.md

 Helixerで解析するゲノム配列Helixerを実行するにあたり、アノテーション対象のゲノム配列を取得する。小さなゲノムを持つ真核生物から中程度のゲノムサイズを持つものをいくつか選択した。
Helixerのtestデータ
パン酵母 (12.1 Mbp)
アユのゲノム (449.1 Mbp)
マウスのリファレンスゲノム (2.7 Gbp)
まずはwgetで配列情報をダウンロードする。
mkdir -p dat && cd dat

# Testデータ
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-47/fasta/arabidopsis_lyrata/dna/Arabidopsis_lyrata.v.1.0.dna.chromosome.8.fa.gz
# パン酵母
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
# アユ
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/036/571/765/GCA_036571765.1_ASM3657176v1/GCA_036571765.1_ASM3657176v1_genomic.fna.gz
# Mus musculs GRCm39
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.fna.gz

cd ..
アユのゲノム配列は後々も使用するので、扱いやすいようにrenameする。
rename用の対応ファイルはこちらのものを使用。
rename_GCA_036571765.1.tsv
cd dat
seqkit replace -K -j 8 -p '(.+)$' -r '{kv}' -k rename_GCA_036571765.1.tsv GCA_036571765.1_ASM3657176v1_genomic.fna.gz -o GCA_036571765.1_ASM3657176v1_genomic_SM.fna
# 一応半角スペースをinplaceで置換
sed -i "s/ /_/g" GCA_036571765.1_ASM3657176v1_genomic_SM.fna
これで、chrかscaffoldの形式に置換できた。

 Helixerの実行USAGE Helixer.pyusage: Helixer.py [-h] [--config-path CONFIG_PATH] [--compression {gzip,lzf}] [--no-multiprocess] [--version] --fasta-path FASTA_PATH --gff-output-path GFF_OUTPUT_PATH [--species SPECIES] [--temporary-dir TEMPORARY_DIR]
                  [--subsequence-length SUBSEQUENCE_LENGTH] [--lineage {vertebrate,land_plant,fungi,invertebrate}] [--model-filepath MODEL_FILEPATH] [--batch-size BATCH_SIZE] [--no-overlap] [--overlap-offset OVERLAP_OFFSET]
                  [--overlap-core-length OVERLAP_CORE_LENGTH] [--debug] [--window-size WINDOW_SIZE] [--edge-threshold EDGE_THRESHOLD] [--peak-threshold PEAK_THRESHOLD] [--min-coding-length MIN_CODING_LENGTH]

optional arguments:
  -h, --help            show this help message and exit
  --version             show program's version number and exit

Data input and output:
  --config-path CONFIG_PATH
                        Config in form of a YAML file with lower priority than parameters given on the command line.
  --fasta-path FASTA_PATH
                        FASTA input file.
  --gff-output-path GFF_OUTPUT_PATH
                        Output GFF file path.
  --species SPECIES     Species name.
  --temporary-dir TEMPORARY_DIR
                        use supplied (instead of system default) for temporary directory

Data generation parameters:
  --compression {gzip,lzf}
                        Compression algorithm used for the intermediate .h5 output files with a fixed compression level of 4. (Default is "gzip", which is much slower than "lzf".)
  --no-multiprocess     Whether to not parallize the numerification of large sequences. Uses half the memory but can be much slower when many CPU cores can be utilized.
  --subsequence-length SUBSEQUENCE_LENGTH
                        How to slice the genomic sequence. Set moderately longer than length of typical genic loci. Tested up to 213840. Must be evenly divisible by the timestep width of the used model, which is typically 9. (Default is lineage dependent
                        from 21384 to 213840).
  --lineage {vertebrate,land_plant,fungi,invertebrate}
                        What model to use for the annotation.
  --model-filepath MODEL_FILEPATH
                        set this to override the default model for any given lineage and instead take a specific model
  --no-overlap          Switches off the overlapping after predictions are made. Predictions without overlapping will be faster, but will have lower quality towards the start and end of each subsequence. With this parameter --overlap-offset and --overlap-
                        core-length will have no effect.

Prediction parameters:
  --batch-size BATCH_SIZE
                        The batch size for the raw predictions in TensorFlow. Should be as large as possible on your GPU to save prediction time. (Default is 8.)
  --overlap-offset OVERLAP_OFFSET
                        Offset of the overlap processing. Smaller values may lead to better predictions but will take longer. The subsequence_length should be evenly divisible by this value. (Default is subsequence_length / 2).
  --overlap-core-length OVERLAP_CORE_LENGTH
                        Predicted sequences will be cut to this length to increase prediction quality if overlapping is enabled. Smaller values may lead to better predictions but will take longer. Has to be smaller than subsequence_length (Default is
                        subsequence_length * 3 / 4)
  --debug               add this to quickly the code runs throughwithout loading/predicting on the full file

Post processing parameters:
  --window-size WINDOW_SIZE
  --edge-threshold EDGE_THRESHOLD
  --peak-threshold PEAK_THRESHOLD
  --min-coding-length MIN_CODING_LENGTH
GPUを使った遺伝子予測はsingularity execに--nvオプションを付ける。testデータとパン酵母についてはGPUを使って実行した。
また、1稿目にはHelixer_postの存在を忘れており、Helixer.pyの出力を使った追加解析が必要と思っていた。
ただ、Helixer.pyのオプションに、window-sizeやedge-threshold、peak-threshold、min-coding-lengthといったpostステップで利用するオプションが指定できることに気づいたので、一括実行されることを期待してパラメータに加えて実行した。
mkdir -p helixer_out/test helixer_out/pan helixer_out/ayu helixer_out/mouse

# test data
singularity exec --nv helixer-docker_helixer_v0.3.4_cuda_12.2.2-cudnn8.sif \
Helixer.py \
--fasta-path dat/Arabidopsis_lyrata.v.1.0.dna.chromosome.8.fa.gz \
--lineage land_plant \
--gff-output-path helixer_out/test/Arabidopsis_lyrata_chromosome8_helixer.gff3 \
--window-size 100 \
--edge-threshold 0.1 \
--peak-threshold 0.8 \
--min-coding-length 60 \
2>&1 | tee Arabidopsis_lyrata_helixer.log

# パン酵母
singularity exec --nv helixer-docker_helixer_v0.3.4_cuda_12.2.2-cudnn8.sif \
Helixer.py \
--fasta-path dat/GCF_000146045.2_R64_genomic.fna.gz \
--lineage fungi \
--gff-output-path helixer_out/pan/GCF_000146045.2_R64_genomic_helixer.gff3 \
--window-size 100 \
--edge-threshold 0.1 \
--peak-threshold 0.8 \
--min-coding-length 60 \
2>&1 | tee R64_helixer.log
以下の2種については、GPUに搭載されているメモリの容量が足りなかったため--nvオプションをつけずにCPUベースで実行。
RAMを多く使用するので、最低でもアユは80GB程度、Mouseについては126GB程度は必要となるので注意してください。
# ayu 
singularity exec helixer-docker_helixer_v0.3.4_cuda_12.2.2-cudnn8.sif \
Helixer.py \
--fasta-path dat/GCA_036571765.1_ASM3657176v1_genomic_SM.fna \
--lineage vertebrate \
--gff-output-path helixer_out/ayu/GCA_036571765.1_ASM3657176v1_genomic_SM_helixer.gff3 \
--window-size 100 \
--edge-threshold 0.1 \
--peak-threshold 0.8 \
--min-coding-length 60 \
2>&1 | tee GCA_036571765.1_ASM3657176v1_genomic_SM_helixer.log

# Mouse
singularity exec helixer-docker_helixer_v0.3.4_cuda_12.2.2-cudnn8.sif \
Helixer.py \
--fasta-path dat/GCF_000001635.27_GRCm39_genomic.fna.gz \
--lineage vertebrate \
--gff-output-path helixer_out/mouse/GCF_000001635.27_GRCm39_genomic_helixer_inpost.gff3 \
--window-size 100 \
--edge-threshold 0.1 \
--peak-threshold 0.8 \
--min-coding-length 60 \
2>&1 | tee GRCm39_helixer_inpost.log
HelixerPostは、スライディングウィンドウによる評価を使用して遺伝子が含まれる可能性のあるゲノム領域を決定する。
スライディングウィンドウは、一定の長さ(ウィンドウサイズ)の領域を少しずつ移動させて解析する手法のこと。
Start/Stopコドン、RNAの転写に関連する既知の情報を利用して、隠れマルコフモデルによる遺伝子予測モデルを作成します。
遺伝子を含むウィンドウを決定するために、設定された幅（例えば100bp）のスライディングウィンドウが遺伝子間（UTR/Coding/Intron）の含有量について評価。
遺伝子が含まれる候補領域は、ウィンドウ内の平均遺伝子スコアが閾値を超えると開始され、平均遺伝子スコアがそのウィンドウを下回るまで継続。
候補領域は、必要な閾値以上の遺伝子スコアを持つウィンドウを少なくとも1つ含んでいれば受け入れ。といった感じの流れのようです。
テストデータとパン酵母は2分かからないくらい、アユは6時間、Mouseは19.7時間程度を要しました。ゲノムサイズによって推定時間と要求されるPCスペックが異なる感じです。
出力されたgff3を見るとこんな感じ。exon領域や5′UTR, 3′UTR位置の推定結果が確認できます。
また、igvに非圧縮のゲノムfastaとgff3をアップロードするとstrand情報も含めて視覚的に状況を確認することが出来ます。
ゲノムのコードデング領域の推定後は一般的にはQCを行います。QCについては以降の記事でまとめます。

 主要な更新履歴2024/11/28: Singularity containerのバージョン変更、モデルのURL修正

 参考Helixer: cross-species gene annotation of large eukaryotic genomes using deep learning

https://doi.org/10.1093/bioinformatics/btaa1044

 備考自作PC1
OS:  Ubuntu 22.04
CPU: Ryzen 5950X 32 threads
GPU: RTX3070Ti メモリ 8GB
Helixer–de novo Prediction of Primary Eukaryotic Gene Models Combining Deep Learning and a Hidden Markov Model

はじめに

インストール

Helixerで解析するゲノム配列

Helixerの実行

主要な更新履歴

参考

備考

Discussion
