macでインフォマティクス

macでインフォマティクス

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

Trinityのalign_and_estimate_abundance.plスクリプト

2011 Nature Biotechnology de novo transcriptome

2021 12/25追記

 

　Trinityに付属するスクリプトalign_and_estimate_abundance.plは、最初にBowtieアライナーを使用してRNA-SeqリードをTrinity転写物にアラインし、その後、RSEMを使ってアバンダンス推定を実行する。

 RSEMは、isoformなどの一意にリードをマッピングして定量することが難しい配列に対して、isoform 内でユニークにマッピングされる領域からのマッピングを利用して、尤もらしい定量値を得る（正しくはExpectation-Maximization (EM)アルゴリズムを使用してMLのアバンダンス推定値を計算（ペーパーのMethod参照））。Long noncoding RNA のデータセットを使ったベンチマークペーパー(pubmed)で、アラインメントベースのカウント方法の中ではhtseq-count やfeatureCountsよりも定量精度が高いことが示されている。このRSEMを組み込んだalign_and_estimate_abundance.plは、精度を担保することが難しいDe novo transcriptome解析でも数少ない有用なスクリプトである。使い方だけ簡単に紹介しておく。

 

インストール

ubuntu18.04LTSでcondaの仮想環境を作ってテストした。

Github

conda create -n trinity python=3.8
conda activate trinity
conda install -c bioconda -y trinity

#rsemも使う
conda install -c bioconda rsem

> align_and_estimate_abundance.pl -h

$ align_and_estimate_abundance.pl -h

 

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

#

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

#  Essential parameters:

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

#

#  --transcripts <string>           transcript fasta file

#

#  --seqType <string>               fq|fa

# 

#  If Paired-end:

#

#     --left <string>

#     --right <string>

#  

#   or Single-end:

#

#      --single <string>

#

#   or (preferred):

#

#      --samples_file <string>    tab-delimited text file indicating biological replicate relationships.

#                                   ex.

#                                        cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq

#                                        cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq

#                                        cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq

#                                        cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq

#

#                      # if single-end instead of paired-end, then leave the 4th column above empty.

#

#

#

#  --est_method <string>           abundance estimation method.

#                                        alignment_based:  RSEM

#                                        alignment_free: kallisto|salmon

#  

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

#  Potentially optional parameters:

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

#

# --output_dir <string>            write all files to output directory 

#                                  (note, if using --samples_file, output_dir will be set automatically according to replicate name))

#  

#

#  if alignment_based est_method:

#       --aln_method <string>            bowtie|bowtie2 alignment method.  (note: RSEM requires either bowtie or bowtie2)

#                                       

###########

# Optional:

# #########

#

# --SS_lib_type <string>           strand-specific library type:  paired('RF' or 'FR'), single('F' or 'R').

#                                  

# --samples_idx <int>               restricte processing to sample entry (index starts at one)

#

#

# --thread_count                   number of threads to use (default = 4)

#

# --debug                          retain intermediate files

#

#  --gene_trans_map <string>        file containing 'gene(tab)transcript' identifiers per line.

#     or  

#  --trinity_mode                   Setting --trinity_mode will automatically generate the gene_trans_map and use it.

#

#

#  --prep_reference                 prep reference (builds target index)

#

#

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

#

#  Parameters for single-end reads:

#

#  --fragment_length <int>         specify RNA-Seq fragment length (default: 200) 

#  --fragment_std <int>            fragment length standard deviation (defalt: 80)

#

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

#  

#   bowtie-related parameters: (note, tool-specific settings are further below)

#

#  --max_ins_size <int>             maximum insert size (bowtie -X parameter, default: 800)

#  --coordsort_bam                  provide coord-sorted bam in addition to the default (unsorted) bam.

#

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

#  RSEM opts:

#

#  --bowtie_RSEM <string>          if using 'bowtie', default: "--all --best --strata -m 300 --chunkmbs 512"

#  --bowtie2_RSEM <string>         if using 'bowtie2', default: "--no-mixed --no-discordant --gbar 1000 --end-to-end -k 200 "

#                                ** if you change the defaults, specify the full set of parameters to use! **

#

#  --include_rsem_bam              provide the RSEM enhanced bam file including posterior probabilities of read assignments.

#  --rsem_add_opts <string>        additional parameters to pass on to rsem-calculate-expression

#

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

#  kallisto opts:

#

#  --kallisto_add_opts <string>  default:   

#

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

#

#  salmon opts:

#

#  --salmon_add_opts <string>    default: 

#

#

#  Example usage

#

#   ## Just prepare the reference for alignment and abundance estimation

#

#    /home/kazu/miniconda3/envs/trinity/bin/align_and_estimate_abundance.pl --transcripts Trinity.fasta --est_method RSEM --aln_method bowtie --trinity_mode --prep_reference

#

#   ## Run the alignment and abundance estimation (assumes reference has already been prepped, errors-out if prepped reference not located.)

#

#    /home/kazu/miniconda3/envs/trinity/bin/align_and_estimate_abundance.pl --transcripts Trinity.fasta --seqType fq --left reads_1.fq --right reads_2.fq --est_method RSEM --aln_method bowtie --trinity_mode --output_dir rsem_outdir

#

##  ## prep the reference and run the alignment/estimation

#

#    /home/kazu/miniconda3/envs/trinity/bin/align_and_estimate_abundance.pl --transcripts Trinity.fasta --seqType fq --left reads_1.fq --right reads_2.fq --est_method RSEM --aln_method bowtie --trinity_mode --prep_reference --output_dir rsem_outdir

#

#   ## Use a samples.txt file:

#

#    /home/kazu/miniconda3/envs/trinity/bin/align_and_estimate_abundance.pl --transcripts Trinity.fasta --est_method RSEM --aln_method bowtie2 --prep_reference --trinity_mode --samples_file samples.txt --seqType fq  

#

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

 

 

 

 

実行方法

Trinityでアセンブルして得られた転写物配列を指定する。strand specificなら--SS_lib_type をつける（FRとRF）。--coordsort_bamをつけると、IGVで可視化するのに利用できるbamも生成される。

align_and_estimate_abundance.pl --transcripts Trinity.fasta \
--seqType fq --left all-1.fq.gz --right all-2.fq.gz \
 --est_method RSEM \
 --aln_method bowtie2 \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20
 --output_dir RSEM_outdir
  •  --seqType <string>    fq|fa
  • --est_method <string>     alignment_based:  RSEM.  alignment_free: kallisto|salmon
  • --trinity_mode    Setting this mode will automatically generate the gene_trans_map and use it.
  • --prep_reference    prep reference (builds target index)
  • --coordsort_bam   provide coord-sorted bam in addition to the default (unsorted) bam.
  • --thread_count    number of threads to use (default = 4)
  • --SS_lib_type     strand-specific library type:  paired('RF' or 'FR'), single('F' or 'R').

if alignment_based est_method:

  •  --aln_method <string>      bowtie|bowtie2 (note: RSEM requires either bowtie or bowtie2)

 

2つのファイルRSEM.isoforms.resultsとRSEM.genes.resultsが結果となる。

f:id:kazumaxneo:20201204200025p:plain

 

De novo transcriptome解析ではRSEM.genes.resultsが出力。TPMやexpected countが示されている（詳細はRNA seqデータ解析本の３章P.79を参照）。

RSEM.genes.results

f:id:kazumaxneo:20201204200036p:plain

 

 

追記

複数サンプル、レプリケートがあるならfastqとサンプルとの関係を示したリストファイルを作成し、それを"--samples_file <list>"で指定する。リストファイルは、helpに書いてあるように、サンプルグループ名<tab>反復名<tab>ペアエンドR1<tab>ペアエンドR2と書く（シングルエンドなら４列目は不要）。

cond_A    cond_A_rep1    A_rep1_left.fq    A_rep1_right.fq
cond_A    cond_A_rep2    A_rep2_left.fq    A_rep2_right.fq
cond_B    cond_B_rep1    B_rep1_left.fq    B_rep1_right.fq
cond_B    cond_B_rep2    B_rep2_left.fq    B_rep2_right.fq

 

このリストファイルを指定する。strand ppecificなら--SS_lib_typeでペアエンドの向きを指定する（FRとRF）。

#RSEM
align_and_estimate_abundance.pl --transcripts Trinity.fasta \
--samples_file list \
 --est_method RSEM \
 --aln_method bowtie2 \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20 \
--SS_lib_type FR \
 --output_dir RSEM_outdir

#kallisto
align_and_estimate_abundance.pl --transcripts Trinity.fasta \
--samples_file list \
 --est_method kallisto \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20 \
--SS_lib_type FR \
 --output_dir RSEM_outdir

#salmon
align_and_estimate_abundance.pl --transcripts Trinity.fasta \
--samples_file list \
--est_method salmon \
 --trinity_mode \
 --prep_reference \
 --coordsort_bam \
 --thread_count 20 \
--SS_lib_type FR \
 --output_dir RSEM_outdir
  • --samples_file <string>    tab-delimited text file indicating biological replicate relationships.
  • --SS_lib_type     strand-specific library type:  paired('RF' or 'FR'), single('F' or 'R').

 

RSEMについて

別のアラインメントプログラムを手動で実行してSAMフォーマットのアラインメントをrsem-calculate-expressionに提供することもできるが、アライナーは、単一の ベスト アラインメントだけではなく、有効なアラインメント全てを報告する設定になっている必要がある。

引用

Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data

Manfred G. Grabherr, Brian J. Haas, Moran Yassour, Joshua Z. Levin, Dawn A. Thompson, Ido Amit, Xian Adiconis, Lin Fan, Raktima Raychowdhury, Qiandong Zeng, Zehua Chen, Evan Mauceli, Nir Hacohen, Andreas Gnirke, Nicholas Rhind, Federica di Palma, Bruce W. Birren, Chad Nusbaum, Kerstin Lindblad-Toh, Nir Friedman, Aviv Regev

Nat Biotechnol. 2011 Jul; 29(7): 644–652 

 

参考

ftp://ftp.broadinstitute.org/pub/users/bhaas/rnaseq_workshop/rnaseq_workshop_2014/Trinity_workshop_activities.pdf

 


http://cbsuss05.tc.cornell.edu/doc/Trinity_Lecture2.pdf

 

関連


参考


 

補足

こちらの論文では、Trinityでアセンブリした配列から疑わしいtranscriptsをフィルタリングする目的でalign_and_estimate_abundance.plスクリプトが使用されています。