macでインフォマティクス

macでインフォマティクス

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

真菌ゲノムのアノテーションパイプライン FunGAP

2017 Bioinformatics fungi annotation docker automated pipeline eukaryotic genome annotation

 

　ゲノム解析が成功するかどうかは遺伝子予測の質にかかっている。fungalゲノムの解読とアセンブルは容易になったが、そのアノテーション手順はまだ標準化されていない。
FunGAP は、真菌ゲノムアセンブリ中のタンパク質をコードする遺伝子を予測するプログラムである。高品質な遺伝子モデルを得るために、複数の遺伝子予測プログラムを実行し、予測されたすべての遺伝子を評価し、既知の配列との相同性から高い支持を得られる遺伝子モデルをアセンブルする。そのために、既知のタンパク質やドメインの相同性に基づいて、各遺伝子モデルの適合性を推定するスコアリング関数を構築した。
FunGAPはPythonスクリプトで書かれており、GitHubで公開されている(https://github.com/CompSynBioLab-KoreaUniv/FunGAP)。このソフトウェアは、非営利目的のユーザーに限り、自由に利用することができる。

 

Githubより

ステップ１

ユーザーから提供されたmRNAリードは、Trinityプログラムによってアセンブルされます。Hisat2リードアライナーとSamtoolsフォーマットコンバータ（SAMファイルからソートされたBAMファイル）により、ゲノムガイド付きアセンブリ用のBAMフォーマットファイルが生成されます。菌類の転写産物アセンブリでは、遺伝子密度が高いとアセンブリー中にUTRオーバーラップが生じるため、Trinityのオプションパラメータ--jaccard_clipを使用しています。このオプションは、隣接する転写産物の融合を避けるのに役立ちます。Hisat2の--max-intronlenオプションで最大イントロン長を2000bpに設定しています。

ステップ2：遺伝子予測
FunGAPでは3つの遺伝子予測ツール；Augustus, Braker, Makerを使用します。予測結果はGFF3ファイルとFASTAファイルに保存され、次のエビデンススコアの計算に使用されます。

Step3：遺伝子モデルの評価とフィルタリング
前のステップでは、3つの遺伝子予測ツールにより、予測された遺伝子のセット（以下、遺伝子モデルとする）が生成されました。FunGAPでは、すべての遺伝子モデルを評価し、最も良いスコアのモデルのみを残すことで、「重ならない」コーディング配列を生成します。評価は BLASTp、BUSCO、InterProScanの3つのツールで行います。アラインメントのビットスコアは、長い遺伝子モデルほど高いアライメントスコアが得られる可能性があるため、長さのカバレッジを乗じています。3つのビットスコアの合計が、各遺伝子モデルのエビデンススコアとなります。最後に、フィルターをかけて、最終的な遺伝子モデルのセットを作成します。

アセンブルされたトランスクリプトームとの塩基配列の類似性は、予測された遺伝子の信頼性を直接証明することができる。FunGAPは各予測遺伝子に対して、Trinityでアセンブルされた転写産物とのBLASTnを実行する。長さの範囲も考慮されます。

上記4つの情報源（ BLASTp、BUSCO、InterProScan、BLASTn）から得られた3つのビットスコアを合計し、各遺伝子モデルのエビデンススコアを算出する（スコアリング関数はGithub参照）。

最後のフィルター処理

FunGAPは少なくとも1つの塩基対でオーバーラップする遺伝子モデルの集合として定義される「遺伝子ブロック」を見つけます。FunGAPは1つの遺伝子ブロック内の全ての遺伝子モデルの組み合わせを取得し、evidence scoreの合計を算出します。エビデンススコアが最も高いブロック内の遺伝子モデルが、その領域の最終遺伝子として選択されます。コーディング配列の短いオーバーラップ（コーディング配列長の10％以下）は許容されます。

 

 

インストール

依存

Required softwares (and tested versions)

  • Hisat2 v2.2.1
  • Trinity v2.12.0
  • RepeatModeler v2.0.1
  • Maker v3.01.03
  • GeneMark-ES/ET v4.65_lic
  • Augustus v3.4.0
  • Braker v2.1.5
  • BUSCO v5.1.2
  • Pfam_scan v1.6
  • BLAST v2.11.0
  • Samtools v1.10
  • Bamtools v2.5.1

Required database

  • Pfam release 34.0

Github

Docker

FunGAP docker image · GitHub

#ここではtanigutiさんのdockerイメージを使う
docker run --rm -it -w /fungap_workspace --rm -v $(pwd):/fungap_workspace docker pull taniguti/fungap-base:v1.1.1 bash

#GeneMarkはライセンスの関係でHPより直接ダウンロードする必要がある (link)。プログラム本体とキーをダウンロードして、プログラムを解凍。中に入る。
cd GeneMark_dir
#Install the key
zcat gm_key_64.gz > ~/.gm_key
#test run 
bash check_install.bash
#パスを通す
export GENEMARK_PATH=/your/path/to/GeneMark-ET/

> /workspace/FunGAP/fungap.py -h

# /workspace/FunGAP/fungap.py -h

usage: fungap.py -g <genome_assembly> -12UA <trans_read_files> -o <output_dir> -a <augustus_species> -b <busco_dataset> -s <sister_proteome>

 

optional arguments:

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

-o [OUTPUT_DIR], --output_dir [OUTPUT_DIR]

Output directory (default: fungap_out)

-1 [TRANS_READ_1], --trans_read_1 [TRANS_READ_1]

Paired-end read1 "<prefix>_1.fastq"

-2 [TRANS_READ_2], --trans_read_2 [TRANS_READ_2]

Paired-end read2 "<prefix>_2.fastq"

-U [TRANS_READ_SINGLE], --trans_read_single [TRANS_READ_SINGLE]

Single read "<prefix>_s.fastq"

-A [TRANS_BAM], --trans_bam [TRANS_BAM]

BAM file (RNA-seq reads alignment to a genome assembly

-g GENOME_ASSEMBLY, --genome_assembly GENOME_ASSEMBLY

Genome assembly file in FASTA format

-a AUGUSTUS_SPECIES, --augustus_species AUGUSTUS_SPECIES

AUGUSTUS species

-b BUSCO_DATASET, --busco_dataset BUSCO_DATASET

BUSCO lineage dataset (run "busco --list-datasets" for the list). Commonly used are fungi_odb10, ascomycota_odb10, and basidiomycota_odb10

-s SISTER_PROTEOME, --sister_proteome SISTER_PROTEOME

Sister proteome sequences in .faa

-c [NUM_CORES], --num_cores [NUM_CORES]

Number of cores to be used (default: 1)

-v, --version show program's version number and exit

--no_braker_fungus No --fungus flag in BRAKER for non-fungus genomes

--no_jaccard_clip No --jaccard_clip flag in Trinity for non-fungus genomes

--no_genemark_fungus No --fungus flag in GeneMark for non-fungus genomes

-M [MAX_INTRON], --max_intron [MAX_INTRON]

Max intron length (Default: 2000 bp)

-t [TRANSLATION_TABLE], --translation_table [TRANSLATION_TABLE]

Translation table (default: 1)

 

 

実行方法

1、タンパク質データベースのダウンロード

python /workspace/FunGAP/download_sister_orgs.py \
 --taxon "Saccharomyces" \
 --email_address byoungnammin@lbl.gov
zcat sister_orgs/*faa.gz > prot_db.faa

 

２、Augustus speciesの確認

python /workspace/FunGAP/get_augustus_species.py \
 --genus_name "Saccharomyces" \
 --email_address your.mail.address

f:id:kazumaxneo:20210811202521p:plain

今回はsaccharomyces_cerevisiae_S288Cを使う（ステップ３）。

 

３、fungap.pyの実行。fastqファイルは解凍して".fastq"指定する。

python /workspace/FunGAP/fungap.py \
 -1 SRR1198667_sampled_1.fastq -2 SRR1198667_sampled_2.fastq \
 -g GCF_000146045.2_R64_genomic.fna \
 -a saccharomyces_cerevisiae_S288C \
 -s prot_db.faa \
 -o fungap_out \
 -b fungi_odb10 \
 -c 12
  • -b   BUSCO lineage dataset (run "busco --list-datasets" for the list). Commonly used are fungi_odb10, ascomycota_odb10, and basidiomycota_odb10 
  • -o　Output directory (default: fungap_out)
  • -1　Paired-end read1 "<prefix>_1.fastq"
  • -2    Paired-end read2 "<prefix>_2.fastq"
  • -A    BAM file (RNA-seq reads alignment to a genome assembly
  • -g    Genome assembly file in FASTA format
  • -a    AUGUSTUS species
  • -s    Sister proteome sequences in .faa
  • -c    Number of cores to be used (default: 1)

fungap.confがないと言われてランできなかった。解決できたら追記します。

 

 

  • Makerは4回実行され、反復的なSNAP遺伝子モデルトレーニングが行なわれる。真菌ゲノムは遺伝子密度が高いため、mRNA において隣接する転写産物の融合を修正するために correct_est_fusion オプションを使用する。split_hit オプションで最大イントロン長を 5000 bp に設定する。50アミノ酸以上のシングルエクソン遺伝子は、single_exonおよびsingle_lengthオプションを設定することで予測される。

  • FunGAP はユーザーが指定した augustus_species パラメーターで Augustus を実行する。リピートマスキングはソフトマスキングされたアセンブリを生成するため、オプション --softmasking がオンになっている。CDS予測の重複を可能にするため、FunGAPは-singlestrandオプションをオンにする。出力はGFF3で、翻訳されたタンパク質配列は簡単なパーススクリプトによってFASTAで生成される。

  • Brakerは、GeneMark-ETとAugustusを用いて、教師なしのRNAシーケンスベースのゲノムアノテーションを行う。リピートマスキングによりソフトマスキングされたアセンブリが生成されるため、オプション --softmasking をオンにする。

  • BLASTPによる系統的に近いゲノムの遺伝子との配列の類似性は、予測された遺伝子が実際の遺伝子であることの証拠となる。系統的に近い生物のプロテオームを提供するには、-sister_proteome 引数を使用する。FunGAPには、NCBIから指定された分類群のタンパク質配列をダウンロードするdownload_sister_orgs.pyというスクリプトが用意されている（上記の１）。

  • BUSCOは真菌類の全ゲノムに保存されているシングルコピーのオルソログの隠れマルコフモデルを提供する。BUSCOのEvidence scoreは、hmmer出力の「full sequence score」とlength coverage [min (query length, target length)/max (query length, target length)]を掛け合わせて算出される。

  • Pfamは、手動で精査されたタンパク質ファミリーのデータベースを提供している。Pfamドメインでアノテーションされた遺伝子モデルは、実際の遺伝子である可能性が高いと想定される。PfamのEvidenceスコアは、InterProScanのXML出力（-f XMLオプション）のhmmer3-matchスコアで直接提供される。1つの遺伝子モデルに複数のドメインが存在する場合は、スコアの合計が使用される。

  • アセンブルされたトランスクリプトームとの配列類似性は、予測された遺伝子の信頼性を直接証明する。FunGAPは各予測遺伝子に対して、Trinityでアセンブルされた転写産物とのBLASTnを実行する。長さの範囲も考慮される。

     

引用

FunGAP: Fungal Genome Annotation Pipeline using evidence-based gene model evaluation
Byoungnam Min, Igor V Grigoriev, In-Geol Choi
Bioinformatics, Volume 33, Issue 18, 15 September 2017, Pages 2936–2937