macでインフォマティクス

macでインフォマティクス

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

ゲノムアノテーションとゲノム多様性解析のためのオールインワンワークフロー EASYstrata

2025 NAR Genomics and Bioinformatics fungi eukaryotic genome annotation automatic pipeline recombination suppression dN/dS

 

 生命のツリー全体にわたって新たなリファレンスゲノムとトランスクリプトームがますます利用可能になり、刺激的な疑問に取り組む新たな道が開かれている。しかしながら、ゲノムのアノテーションと進化プロセスの推論には依然として課題があり、方法論の標準化も欠如している。本稿では、これらの課題を克服し、組換え抑制の検出と、その結果として生じるリアレンジメントおよび転移因子の蓄積を容易にする、進化解析用の新しいワークフローを提案する。このワークフローを実現するために、複数のバイオインフォマティクスステップを単一の使いやすいワークフローに統合した。最先端ツールを組み合わせ、転移因子の検出、ゲノムのアノテーション、遺伝子相同性関係の推論、配列間の分岐の計算、進化階層（すなわち、組換え抑制の段階的拡大のフットプリント）とその構造リアレンジメントの推論を行い、結果を視覚化する。EASYstrataと呼ばれるこのワークフローは、ミクロボトリウム属真菌の公開ゲノム42個に再アノテーションを施すために適用された。植物と動物のさらなる事例では、前述のものと同じ層を復元できることを示す。このツールは性染色体または交配型染色体間の分岐を推定することを目的として開発されたが、分岐パターンが関心のある任意のハプロタイプのペアに適用できる。このワークフローは、新たに配列決定されたphased二倍体ゲノムが利用可能になりつつある非モデル種の研究を促進するだろう。

 

Installation instructions

https://github.com/QuentinRougemont/EASYstrata/blob/main/INSTALL/INSTALL.md

 

補足:

組換え抑制:本来なら相同染色体間で起きる組換えが、何らかの理由により､特定の染色体領域で長期的に起きなくなる現象｡組換えが止まると、その領域は1つの塊として世代を超えて受け継がれ､相同染色体間で塩基配列の差が蓄積する｡時間経過によってアミノ酸を変えない同義置換率 dSが高くなったり､機能喪失や擬遺伝子化､トランスポゾンの蓄積が起きやすくなる(論文図1参照)｡

 

インストール

Github

git clone https://github.com/QuentinRougemont/EASYstrata
cd EASYstrata 
#まずcoondaで導入可能なパッケージを導入
mamba env create -f INSTALL/annotation_env.yml   
#for busco:
mamba env create -f INSTALL/busco_env.yml
#for repeatmodeller:
mamba env create -f INSTALL/repeatmodeler.yml

#続いて非conda対応パッケージを導入
mamba create -n superannot -y -c conda-forge -c defaults python=3.11
conda activate superannot
#指示に従って勧める("ALTERNATIVELY: install each non-conda software one by one:"以降の部分)
https://github.com/QuentinRougemont/EASYstrata/blob/main/INSTALL/INSTALL.md

 

実行方法

ランする際は､カレントから見てconfig/configファイルが自動で認識され､この設定をもとに解析フローが進められる｡また､-oで実行するステップを選択できる｡

configを編集する｡

例: https://github.com/QuentinRougemont/EASYstrata/blob/main/example_data/example.config

# config file
#--- COMPULSORY MINIMAL LEVEL OF INFORMATION REQUIRED -----
path=$(pwd) #replace this by your own full path to your data!

genome1="$path"/"example_data/genome1.fa.gz"     #full path to current genome1 assembly (haplotype1 - compressed or not)
genome2="$path"/"example_data/genome2.fa.gz"     #full path to current genome2 assembly (haplotype2)

haplotype1="genome1" #name1 [name of haplotype1 - must match the scaffold/chromosome basename]
haplotype2="genome2" #name2 [name of haplotype2 - must match the scaffold/chromosome basename]

#----- optional --------
ancestral_genome="$path"/"example_data/Mlag129A1.fa.gz" #full path to the ancestral genome to be used whenever possible
ancestral_gff="$path"/"example_data/Mlag129A1.gtf.gz"   #full path to the associated gff to be used whehever possible


#--- annotate or not #
annotate="YES"  #a string (YES/NO)? if annotation = YES then annotation of the genomes will be performed
             #else gtf and fastafiles are expected and only the paml/ds/genespace etc will be performed
RelatedProt="" #/path/to/related_protein.fa" #a full path to a set of external protein data (fasta format) for braker
fungus="YES" #YES or NO set to YES if species is a fungus

## orthoDB ###
orthoDBspecies="Fungi" #one of : "Metazoa" "Vertebrata" "Viridiplantae" "Arthropoda" "Eukaryota" "Fungi" "Alveolata" "Stramenopiles"
#if a species is given then it will be used for annotation (in addition to eventual RNAseq and orthProteins) 
 
#if annotate = NO then gtf should be provided: 
gtf1=""
gtf2=""

#RNASeq data ?
rnaseq="YES"   #YES/NO a string stating wether rnaseq data is available or not
RNAseqlist="$path"/"example_data/rnaseq.list.txt" #"/full/path/to/rnaseq.list.txt" #list of rnaseq data 
#OR BAM file
#bamlist1="/full/path/to/bam.aligned_on_genome1.txt" #list of bam on genome1
#bamlist2="/full/path/to/bam.aligned_on_genome2.txt" #list of bam on genome2


#---- COMPULSORY INFOS ------
#TE INFO:
annotateTE="YES" #YES/NO
TEdatabase="$path"/"example_data/TE.fa.gz"  #a full path to a custom database of species/genus species TE
#NCBI species for de-novo TE:
ncbi_species="fungi"
#if annotateTE="NO" you may want to provide bed formatted TE elements to display on circos plots:
#TEgenome1=""
#TEgenome2=""
#TEancestral=""

#BUSCO SPECIFIC ARGUMENTS:
#busco_lineage name:
busco_lineage="basidiomycota_odb10" #see busco lineage to have the list

#option for Dn/Ds, plot, etc:
scaffold="$path"/"scaffold.txt" #scaffold.txt #a tab separated file providing the genome name in first column, second column the scaffold for the region of interest (X/Y chromosome, supergenes, etc) third colum a string (normal or reversed) expliciting wether the scaffold orientation should be reversed or not


#option for running interproscan at the quality checks steps (can take several days for big dataset)
interpro="NO" #YES #by default set to no

 

このexample.configを使う｡

cp example_data/example.config config/config

 

-o 2だとstep1のゲノムアノテーションと品質チェック､step2のGeneSpaceによる共線性解析だけが実行される｡

./master.sh -o 2 2>&1 |tee log

example.configをテンプレートにすると比較的簡単に実行できたが､細かい修正が多数ある(*1-6)｡

 

コメント

ワークフロー自体はとても良いと思うのですが、実際に正常動作させるには細かい点で相当な修正が必要だという印象を受けました。修正して動作する Docker image を作成し、Docker Hub に公開できればとしばらく試みましたが(そもそもライセンス的に一部ソフトは非対応)、今回は一旦諦めます。熱意があれば、working code としてプルリクエストを出す形の方がずっと建設的かもしれませんね。

引用

EASYstrata: an all-in-one workflow for genome annotation and genomic divergence analysis
Quentin Rougemont, Elise Lucotte, Loreleï Boyer, Alexandra Jalaber de Dinechin, Alodie Snirc, Tatiana Giraud, Ricardo C Rodríguez de la Vega

NAR Genom Bioinform. 2025 Aug 27;7(3):lqaf110

 

関連

 

*1

テストラン用のRNAは､ゲノムがMicrobotryum violaceumが一致もしくは近い系統だったので､M. violaceumのRNA seqであるSRR24741845を使用した｡以下の手順で配置した｡

cd EASYstrata/
mkdir -p example_data/RNA
fastq-dump --outdir example_data/RNA --gzip --split-files SRR24741845
mv example_data/RNA/SRR24741845_1.fastq.gz example_data/RNA/input1.R1.fastq.gz
mv example_data/RNA/SRR24741845_2.fastq.gz example_data/RNA/input1.R2.fastq.gz

 

#2

12行目ancestral_gff=からancestral_gtf=に修正

 

#3 RNAのデータパスを記載したrnaseq.list.txtが必要

> cat /root/EASYstrata/example_data/rnaseq.list.txt

/root/EASYstrata/example_data/RNA/input1.R1.fastq.gz

/root/EASYstrata/example_data/RNA/input1.R2.fastq.gz

で認識した｡

 

#4

インストール手順にないがjavaも必要

apt-get install -y default-jre

 

$5

mkdir -p /root/EASYstrata/haplo1/example_data/
mkdir -p /root/EASYstrata/haplo2/example_data/

cp /root/EASYstrata/example_data/rnaseq.list.txt /root/EASYstrata/haplo1/example_data/
cp /root/EASYstrata/example_data/rnaseq.list.txt /root/EASYstrata/haplo2/example_data/

 

*6

./master.sh -o 2の実行時､gsnapによるRNA seqのショートリードのゲノムへのマッピングでセグメンテーション違反が起き､改善できなかった(自分でビルドし直したSIMD命令を使わないgsnapでも同じエラーが発生)｡ GMAPには本家とfolkがあり､EASYstrataではfolkを使うようになっていたので､Thomas Wu氏の本家ソースコードの最新版をhttp://research-pub.gene.com/gmap/からダウンロードして: Version 2025-07-31 (v2) ､これをビルドしたところ､より進捗したが､繰り返し配列へのマッピングと思われるリードの扱いで再びセグメンテーション違反が発生した｡進捗が見られたので､03_gsnap_PE.shのgsnapのオプション( -M 2 -n 10 -N 1 -w 200000 --pairmax-rna=200000 -E 1 -B 2)を消して､単純にgsnap --gunzip -t 8 -A sam  --dir=03_genomeとすると最後のリードまでsam出力された｡具体的には03_gsnap_PE.shのマッピング部分を以下の通りに変更した｡ただしこれだとイントロンをまたいでアラインメントされないので､ただランできただけとなる｡

gsnap --gunzip -t "$NCPUS_GSNAP" -A sam \
      --dir="$GENOMEFOLDER" -d "$GENOME" \
      "$DATAINPUT"/"${base}"_R1.paired.fastq.gz \
      "$DATAINPUT"/"${base}"_R2.paired.fastq.gz \
      > "$DATAOUTPUT"/"$base".sam