macでインフォマティクス

macでインフォマティクス

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

アノテーションパイプライン PASA

2003 Nucleic Acids Research annotation eukaryotic genome annotation database docker

2020 10/4 コマンドの間違い修正

2020 10/5 アップデートのコマンド修正

2023/01/05, 01/12追記 アップデートのコマンド追記

 

　ゲノム配列に対する発現配列データのスプライスアラインメントは、真核生物ゲノムにおける遺伝子の包括的なアノテーションにおいて重要なツールであることが証明されている。これにより、利用可能なすべてのトランスクリプトデータを包括的に取り込み、微妙なスプライシングのバリエーションを捉えることができる。この方法で同定された完全および部分的な遺伝子構造は、シロイヌナズナゲノムアノテーション（TIGR release v.4.0）を改善するために使用された。アライメントアセンブリは、いくつかの新規遺伝子と1000以上のオルタナティブスプライシングバリエーションの自動モデリングを可能にし、約27,000のアノテーションされたタンパク質コード化遺伝子の約半分の更新（UTRアノテーションを含む）を可能にした。PASA(Program to Assemble Spliced Alignments)ツールのアルゴリズムと、シロイヌナズナ遺伝子アノテーションの自動更新の結果を紹介する。

  PASAは、ESTとcDNAのアラインメントに基づき、互換性のある重複するアラインメントのセットをマージしてユニークな最大のアラインメントを作成する。このアルゴリズムは、シロイヌナズナのFL-cDNAの非翻訳領域（UTR）を統合して最大化し、FL-cDNAのアラインメントが利用できなかった領域のESTエビデンスを統合するために採用され、シロイヌナズナゲノムアノテーションの更新と改善のための最大の基質を提供した。

 

Documentation

https://github.com/PASApipeline/PASApipeline/wiki

 

インストール

condaでも導入できるが、依存が多いため、公式のdockerイメージを使ってテストした。

依存

Github

docker

PASA_Docker · PASApipeline/PASApipeline Wiki · GitHub

docker pull pasapipeline/pasapipeline:latest

#こちらも動作することを確認(link)
docker pull hamiltonjp/pasapipeline:2.5.2

#bioconda (link)
conda install -c bioconda pasa

> perl ../Launch_PASA_pipeline.pl

# perl ../Launch_PASA_pipeline.pl

 

############################# Options ###############################

# 

#   * indicates required

#

#

# --config|-c * <filename>  alignment assembly configuration file

#

# // spliced alignment settings

# --ALIGNERS <string>   aligners (available options include: gmap, blat... can run both using 'gmap,blat')

# -N <int>              max number of top scoring alignments (default: 1)

# --MAX_INTRON_LENGTH|-I  <int>         (max intron length parameter passed to GMAP or BLAT)  (default: 100000)

# --IMPORT_CUSTOM_ALIGNMENTS_GFF3 <filename> :only using the alignments supplied in the corresponding GFF3 file.

# --trans_gtf <filename>      :incorporate cufflinks or stringtie--generated transcripts

#

#

# // actions

# --create|-C               flag, create database

# --replace|-r               flag, drop database if -C is also given. This will DELETE all your data and it is irreversible.

# --run|-R               flag, run alignment/assembly pipeline.

# --annot_compare|-A               (see section below; can use with opts -L and --annots)  compare to annotated genes.

# --ALT_SPLICE     flag, run alternative splicing analysis

 

# // input files

# --genome|-g * <filename>  genome sequence FASTA file (should contain annot db asmbl_id as header accession.)

# --transcripts|-t * <filename>  transcript db 

# -f <filename>    file containing a list of fl-cdna accessions.

# --TDN <filename> file containing a list of accessions corresponding to Trinity (full) de novo assemblies (not genome-guided)

#

# // polyAdenylation site identification  ** highly recommended **

# -T               flag,transcript db were trimmed using the TGI seqclean tool.

#    -u <filename>   value, transcript db containing untrimmed sequences (input to seqclean)

#                  <a filename with a .cln extension should also exist, generated by seqclean.>

#

#

#

# Misc:

# --TRANSDECODER   flag, run transdecoder to identify candidate full-length coding transcripts

# --CPU <int>      multithreading (default: 2)

# --PASACONF <string> path to a user-defined pasa.conf file containing mysql connection info

#                      (used in place of the $PASAHOME/pasa_conf/conf.txt file)

#                      (and allows for users to have their own unique mysql connection info)

#                      (instead of the pasa role account)

#

# -d               flag, Debug 

# -h               flag, print this option menu and quit

#

#########

#

# // Transcript alignment clustering options (clusters are fed into the PASA assembler):

#

#       By default, clusters together transcripts based on any overlap (even 1 base!).

#

#    Alternatives:

#

#        --stringent_alignment_overlap <float>  (suggested: 30.0)  overlapping transcripts must have this min % overlap to be clustered.

#

#        --gene_overlap <float>  (suggested: 50.0)  transcripts overlapping existing gene annotations are clustered.  Intergenic alignments are clustered by default mechanism.

#               * if --gene_overlap, must also specify --annots  with annotations in recognizable format (gtf, gff3, or data adapted) (just examines 'gene' rows, though).

#

#

#

# --INVALIDATE_SINGLE_EXON_ESTS    :invalidates single exon ests so that none can be built into pasa assemblies.

#

#

# --transcribed_is_aligned_orient   flag for strand-specific RNA-Seq assemblies, the aligned orientation should correspond to the transcribed orientation.

#

#

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

#

#  // Annotation comparison options (used in conjunction with -A at top).

#   

#  -L   load annotations (use in conjunction with --annots)

#  --annots <filename>  existing gene annotations in recognized format (gtf, gff3, or custom adapted).

#  --GENETIC_CODE (default: universal, options: Euplotes, Tetrahymena, Candida, Acetabularia)

#

###################### Process Args and Options #####################

> scripts/Load_Current_Gene_Annotations.dbi

# scripts/Load_Current_Gene_Annotations.dbi

 

############################# Options ###############################

#

# -c configFile

# -g genome fasta file

# -d debug

#

# ## Annotation database settings:

# -P param string for data adapter: HOOK_GENE_STRUCTURE_UPDATER, as specified in the pasa_conf/conf.txt file.

#

#                Typically, 'param' is just either a gtf or gff3 file containing gene & transcript annotations. (ie. from Gencode)

#

#        (decades ago...)     At TIGR, this would be "SYBTIGR,Sybase,user,password,annotdb[,pred_type]"

#                             Where pred_type is some prediction type (ie. genscan), or current annotations [default]

#

# -h print this option menu and quit

# -v verbose

#

###################### Process Args and Options #####################

 

 

 

テストラン

ここでは/home/kazuに作業ディレクトリtest/を作り実行する。

cd /home/kazu/test/
git clone https://github.com/PASApipeline/PASApipeline.git

#入力ゲノムを解凍、カレントにコピー
gzip -dv PASApipeline/sample_data/genome_sample.fasta.gz
cp PASApipeline/sample_data/genome_sample.fasta .
cp PASApipeline/sample_data/sqlite.confs/alignAssembly.config .
cp PASApipeline/sample_data/all_transcripts.fasta.clean .


#run (PASApipeline/sample_dataに移動してからランしている)
docker run --rm -it -v /tmp:/tmp -v /home/kazu/test:/home/kazu/test \
 pasapipeline/pasapipeline:latest \
 bash -c 'cd /home/kazu/test && \
 /usr/local/src/PASApipeline/Launch_PASA_pipeline.pl \
 -c alignAssembly.config -C -R --ALIGNER gmap \
 -g genome_sample.fasta \
 -t all_transcripts.fasta.clean '
  • --config|-c <filename>   alignment assembly configuration file
  • --run|-R flag    run alignment/assembly pipeline.
  • --create|-C flag    create database
  • --ALIGNERS <string>    aligners (available options include: gmap, blat... can run both using 'gmap,blat')
  • --CPU  <int>  multithreading (default: 2)
  • --genome|-g  <filename>   genome sequence FASTA file (should contain annot db asmbl_id as header accession.)
  • --transcripts|-t <filename>  transcript db

注；　alignAssembly.configのDATABASE=/tmp/sample_mydb_pasa.sqliteが/tmpに作成されるDB名だが、同時に全ての出力アノテーションファイルのprefixにもなる。ラン前に変更しておく（例；/tmp/mygenome）。

 

 出力

f:id:kazumaxneo:20201004153342p:plain

 

 

追記

既存のアノテーションのアップデート（完成度の高い遺伝子モデルができた時に、ESTもしくはRNA seqのassembly（cDNA）を使ってUTRのアノテーションを付加するなど。）

＃まずアップデートするGFF3ファイルにエラーがないか検証（メッセージがでない場合問題なし）
/usr/local/src/PASApipelinemisc_utilities/pasa_gff3_validator.pl original_annotation.gff3

#データベースへのロード
/usr/local/src/PASApipeline/scripts/Load_Current_Gene_Annotations.dbi \
-c alignAssembly.config \
-g genome_sample.fasta \
-P original_annotation.gff3

#PASAの実行。品質の高いfull-length cDNAを使う。上と同じconfigファイルが必要。
/usr/local/src/PASApipeline/Launch_PASA_pipeline.pl --CPU 20 \
-c annotCompare.config -A \
-g genome_sample.fasta \
-t transcripts_clean.fasta

configファイル中のprefix名のGFF3ファイルが出力される。実際にランするにはMySQL DBの名前やパスなどを設定したPASAのconfigファイルと、alignAssembly.config ファイルをいくつか編集してからランする必要がある。詳細はこちらで説明されています。

http://wfleabase.org/release1/PASA_gene_annotation.html

 

 

推奨される解析手順については、Documentのwalk through を確認して下さい。 

 

メモ

  • --TRANSDECODERオプション付きで実行すると、新規遺伝子を追加し、完全長転写物と推定されるものに基づいて、より積極的に注釈を修正する（EVM #55）。
  • misc_utilitiesにはサポートスクリプトが配置されている（link）。
  • GAAPの論文でPASAの使い方も含めた真核生物のアノテーションの手順（structural annotation & functional annotation）について詳しく説明されている（リンク）。

https://github.com/PASApipeline/PASApipeline/tree/master/misc_utilities

引用

Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies
Brian J. Haas, Arthur L. Delcher, Stephen M. Mount,1 Jennifer R. Wortman, Roger K. Smith, Jr, Linda I. Hannick, Rama Maiti, Catherine M. Ronning, Douglas B. Rusch,2 Christopher D. Town, Steven L. Salzberg, Owen White

Nucleic Acids Res. 2003 Oct 1; 31(19): 5654–5666

 

関連