macでインフォマティクス

macでインフォマティクス

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

de novo transcriptome 解析を行うためのRNAのアセンブルツール

 

 RNAの発現は分子によって大きな差があるため、どのようなデータセットでも万能なRNAアセンブルツールというのは存在しない。そのため単一のアセンブルツールでアセンブルするより、複数のアセンブルツールでかつ複数のk-merサイズでアセンブルして、そのプール全体からよりよいmRNAセットを選抜してくるほうが、より完全なRNAのコレクションになりやすい。一般に、より大きいk値は、高い遺伝子発現レベルまたはより長いコンティグで良好に機能し、小さいk値は、低い遺伝子発現レベルまたはより短いコンティグで良好に機能する。したがって、Trans-ABySS、Oases-M(Oaseのmultiple-kバージョン)、IDBA-Tranなどのいくつかのアセンブラでは、複数のアセンブリが異なるk値が合併されてより高い感度が得られるが、同時により多くの偽陽性転写産物もできる傾向にある。

以下のポスターの右側に、その複数のアセンブルツールを使う覚書があるので確認してみてください。

Do and Don't

https://f1000research.com/posters/5-1695

 

したがって、必然的にde novo transcriptome解析では複数のアセンブルツールを使ってアセンブルを行うことになる。ここでは、de novo transcriptomeの代表的なアセンブルプログラムのインストールのリンクや流れををまとめておく。

 

修正 1/31

TransAbySS

 

ツールのインストール 

 

Trinity 

Releases · trinityrnaseq/trinityrnaseq · GitHub

#brewで導入できる。
brew install Trinity

> trinity

$ trinity

 

 

 

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

#

#     ______  ____   ____  ____   ____  ______  __ __

#    |      ||    \ |    ||    \ |    ||      ||  |  |

#    |      ||  D  ) |  | |  _  | |  | |      ||  |  |

#    |_|  |_||    /  |  | |  |  | |  | |_|  |_||  ~  |

#      |  |  |    \  |  | |  |  | |  |   |  |  |___, |

#      |  |  |  .  \ |  | |  |  | |  |   |  |  |     |

#      |__|  |__|\_||____||__|__||____|  |__|  |____/

#

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

#

# Required:

#

#  --seqType <string>      :type of reads: ('fa' or 'fq')

#

#  --max_memory <string>      :suggested max memory to use by Trinity where limiting can be enabled. (jellyfish, sorting, etc)

#                            provided in Gb of RAM, ie.  '--max_memory 10G'

#

#  If paired reads:

#      --left  <string>    :left reads, one or more file names (separated by commas, no spaces)

#      --right <string>    :right reads, one or more file names (separated by commas, no spaces)

#

#  Or, if unpaired reads:

#      --single <string>   :single reads, one or more file names, comma-delimited (note, if single file contains pairs, can use flag: --run_as_paired )

#

#  Or,

#      --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.

#

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

##  Misc:  #########################

#

#  --SS_lib_type <string>          :Strand-specific RNA-Seq read orientation.

#                                   if paired: RF or FR,

#                                   if single: F or R.   (dUTP method = RF)

#                                   See web documentation.

#

#  --CPU <int>                     :number of CPUs to use, default: 2

#  --min_contig_length <int>       :minimum assembled contig length to report

#                                   (def=200)

#

#  --long_reads <string>           :fasta file containing error-corrected or circular consensus (CCS) pac bio reads

#                                   (** note: experimental parameter **, this functionality continues to be under development)

#

#  --genome_guided_bam <string>    :genome guided mode, provide path to coordinate-sorted bam file.

#                                   (see genome-guided param section under --show_full_usage_info)

#

#  --jaccard_clip                  :option, set if you have paired reads and

#                                   you expect high gene density with UTR

#                                   overlap (use FASTQ input file format

#                                   for reads).

#                                   (note: jaccard_clip is an expensive

#                                   operation, so avoid using it unless

#                                   necessary due to finding excessive fusion

#                                   transcripts w/o it.)

#

#  --trimmomatic                   :run Trimmomatic to quality trim reads

#                                        see '--quality_trimming_params' under full usage info for tailored settings.

#                                  

#

#  --no_normalize_reads            :Do *not* run in silico normalization of reads. Defaults to max. read coverage of 50.

#                                       see '--normalize_max_read_cov' under full usage info for tailored settings.

#                                       (note, as of Sept 21, 2016, normalization is on by default)

#     

#  --no_distributed_trinity_exec   :do not run Trinity phase 2 (assembly of partitioned reads), and stop after generating command list.

#

#

#  --output <string>               :name of directory for output (will be

#                                   created if it doesn't already exist)

#                                   default( your current working directory: "/Users/kazumaxneo/Downloads/v0.99/trinity_out_dir" 

#                                    note: must include 'trinity' in the name as a safety precaution! )

#  

#  --full_cleanup                  :only retain the Trinity fasta file, rename as ${output_dir}.Trinity.fasta

#

#  --cite                          :show the Trinity literature citation

#

#  --verbose                       :provide additional job status info during the run.

#

#  --version                       :reports Trinity version (Trinity-v2.4.0) and exits.

#

#  --show_full_usage_info          :show the many many more options available for running Trinity (expert usage).

#

#

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

Trinityにはgenome guide assemblyとde novo asemblyの2つの方式がある。

https://github.com/trinityrnaseq/trinityrnaseq/wiki/Genome-Guided-Trinity-Transcriptome-Assembly

ラン

1、genome guide assembly

#まずindexをつけ、ゲノムにmappingする。
#index
mkdir genome/ #空フォルダを作っておく
STAR --runMode genomeGenerate --genomeDir genome/ --genomeFastaFiles chromosome.fa --runThreadN 12 --sjdbGTFfile reference.gtf #gtfがないゲノムなら、gtfのオプションは指定なしでindex。

#mapping 予めmergeしておいたfastqを使いマッピング(例えば cat *pair1.fq > merge_pair1.fq)
STAR --genomeDir genome/ --readFilesIn merged_pair1.fq merged_pair2.fq --runThreadN 30 --outSAMtype BAM SortedByCoordinate --outFileNamePrefix sample

#bamのアライメント情報を使いtrinityでgenome guide assembly (k-merはdefaultの25)
Trinity --genome_guided_bam sample_mergeAligned.sortedByCoord.out.bam --genome_guided_max_intron 50000 --max_memory 100G --CPU 20

2、De novo assembly

#全サンプルをmergeしたfastqを使い、de novoでアセンブル (mergeはcat *pair1.fq > merge_pair1.fq)
Trinity --seqType fq --left pair1.fastq --right pair2.fastq --max_memory 100G --CPU 24

 

 

  • Trans-ABySS

GitHub - bcgsc/transabyss: de novo assembly of RNA-Seq data using ABySS.

依存

#python2環境でうごくので、pyenv等で2.7に切り替える。

#BLATbrewで導入する。ただしABySSはbrewで導入するとversion2以降が入ってしまう。上のリンクからyosemite用のバイナリをダウンロードしてパスを通す。
brew install blat

#依存のigraphのインストール。
pip install python-igraph

#本体
git clone https://github.com/bcgsc/transabyss.git
cd transabyss/
./transabyss

$ ./transabyss

Trans-ABySS 1.5.5

CMD: ./transabyss

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Found Trans-ABySS directory at: /Users/kazumaxneo/Documents/transabyss

Found Trans-ABySS `bin` directory at: /Users/kazumaxneo/Documents/transabyss/bin

Found script at: /Users/kazumaxneo/Documents/transabyss/bin/skip_psl_self.awk

Found script at: /Users/kazumaxneo/Documents/transabyss/bin/skip_psl_self_ss.awk

Found `blat' at /usr/local/bin/blat

Found `MergeContigs' at /usr/local/bin/MergeContigs

Found `abyss-filtergraph' at /usr/local/bin/abyss-filtergraph

Found `abyss-map' at /usr/local/bin/abyss-map

Found `abyss-pe' at /usr/local/bin/abyss-pe

Found `abyss-junction' at /usr/local/bin/abyss-junction

# CPU(s) available: 8

# thread(s) requested: 1

# thread(s) to use: 1

ERROR: No input reads specified! Use option '--pe' to specify paired-end reads and/or option '--se' to specify single-end reads.

リンクを張るかパスを通しておく。

 ランチュートリアル

https://github.com/bcgsc/transabyss/blob/master/TUTORIAL.txt

 ラン

アセンブル。 

#transabyssでアセンブル
#K=21
transabyss --pe merged_pair1.fq merged_pair2.fq -k 21 --threads 10 --length 100 --outdir k21
#K=31
transabyss --pe merged_pair1.fq merged_pair2.fq -k 31 --threads 10 --length 100 --outdir k31
#K=41
transabyss --pe merged_pair1.fq merged_pair2.fq -k 41 --threads 10 --length 100 --outdir k41
#K=51
transabyss --pe merged_pair1.fq merged_pair2.fq -k 51 --threads 10 --length 100 --outdir k51
#K=61
transabyss --pe merged_pair1.fq merged_pair2.fq -k 61 --threads 10 --length 100 --outdir k61

# transabyss-mergeで複数のk-merアセンブルをmerge
transabyss-merge k21.fa k31.fa k41.fa k51.fa k61.fa --mink 21 --maxk 61 --prefix k21. k31. k41. k51. k61

 

 

 

  • SOAPdenovo-Trans

macではgnuのmakeに失敗したので、深追いせずcent OSに導入した。

ソースコード

SOAPdenovo-Trans - Browse /SOAPdenovo-Trans at SourceForge.net

BInary (linux)

SOAP :: Short Oligonucleotide Analysis Package

mkdir SOAPdenovo && cd SOAPdenovo/
wget https://downloads.sourceforge.net/project/soapdenovotrans/SOAPdenovo-Trans/bin/v1.03/SOAPdenovo-Trans-bin-v1.03.tar.gz
tar zxvf SOAPdenovo-Trans-bin-v1.03.tar.gz

#カレントにパスを通す。
echo export PATH=\$PATH:`pwd`\ >> ~/.bash_profile && source ~/.bash_profile

> SOAPdenovo-Trans-31mer 

$ SOAPdenovo-Trans-31mer 

 

The version 1.03: released on July 19th, 2013

 

 

Usage: SOAPdenovo-Trans <command> [option]

    pregraph     construction kmer-graph

    contig       eliminate errors and output contigs

    map          map reads to contigs

    scaff        scaffolding

    all          doing all the above in turn

SOAPdenovo-Trans-127merもある。

ラン

はじめにconfig_fileを作る。

cat > SOAP.config <<TEXT 
max_rd_len=150
[LIB] avg_ins=192
reverse_seq=0
asm_flags=3
q1=R1_pairedout.fastq
q2=R2_pairedout.fastq
TEXT

config_fileを指定してアセンブル。 

SOAPdenovo-Trans-31mer all -K 21 -p 20 -s config_file -o SOAPdenovo-Trans_output_k21

 SOAPdenovo-Trans_output_k21.scafStatistics にアセンブルのサマリーが出る。アセンブル結果はSOAPdenovo-Trans_output_k21.scafSeq に出る。

SOAPdenovo-Transのランタイムはかなり短いが、メモリ要求量はかなり多め。


 

  • Oases

このOases自体、velvetの複数k-merアセンブルをマージするツールとなる。 

https://github.com/dzerbino/oases

git clone --recursive https://github.com/dzerbino/oases 
cd oases
make

> ./oases

$ ./oases 

oases - De novo transcriptome assembler for the Velvet package

Version 0.2.09

 

Copyright 2009- Daniel Zerbino (dzerbino@soe.ucsc.edu)

This is free software; see the source for copying conditions.  There is NO

warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 

Compilation settings:

CATEGORIES = 2

MAXKMERLENGTH = 64

 

Usage:

./oases directory [options]

 

directory : working directory name

 

Standard options:

--help : this help message

--version : print version and exit

--citation : print citation to Oases manuscript and exit

-ins_length2 <integer> : expected distance between two paired-end reads in the second short-read dataset (default: no read pairing)

-ins_length_long <integer> : expected distance between two long paired-end reads (default: no read pairing)

-ins_length*_sd <integer> : est. standard deviation of respective dataset (default: 10% of corresponding length)

[replace '*' by nothing, '2' or '_long' as necessary]

-unused_reads <yes|no> : export unused reads in UnusedReads.fa file (default: no)

-amos_file <yes|no> : export assembly to AMOS file (default: no export)

-alignments <yes|no> : export a summary of contig alignment to the reference sequences (default: no)

 

Advanced options:

-cov_cutoff <floating-point> : removal of low coverage nodes AFTER tour bus or allow the system to infer it (default: 3)

-min_pair_count <integer> : minimum number of paired end connections to justify the scaffolding of two long contigs (default: 4)

-min_trans_lgth <integer> : Minimum length of output transcripts (default: hash-length)

-paired_cutoff <floating-point> : minimum ratio allowed between the numbers of observed and estimated connecting read pairs

Must be part of the open interval ]0,1[ (default: 0.1)

-merge <yes|no> :Preserve contigs mapping onto long sequences to be preserved from coverage cutoff (default: no)

-edgeFractionCutoff <floating-point> : Remove edges which represent less than that fraction of a nodes outgoing flow

Must be part of the open interval ]0,1[ (default: 0.01)

-scaffolding <yes|no> : Allow gaps in transcripts (default: no)

-degree_cutoff <integer> : Maximum allowed degree on either end of a contigg to consider it 'unique' (default: 3)

 

Output:

directory/transcripts.fa

directory/contig-ordering.txt

パスの通ったディレクトリに移動しておく。 

velvetのインストールについては、こちらのリンクのvelvetの項目を参照。

 

ラン

アセンブル。 velveth => velvetg => oasesと進む。

velveth velvet_k21 21 -fastq merged.fq 
velvetg velvet_k21 -ins_length1 300 -ins_length2 3000 -read_trkg yes
oases velvet_k21/ -min_trans_lgth 100 -ins_length 215

oasesはスクリプトでもランできる。

python scripts/oases_pipeline.py -m 21 -M 23 -d 'data/test_reads.fa'

 

  

以前の記事でも触れていましたが、内容が膨れ上がってわかりにくかったので、視点をアセンブラに絞ってまとめなおしました。


 

 

参考資料

Sildeshareの関連プレゼン資料

 

論文3つ

Combining Transcriptome Assemblies from Multiple De Novo Assemblers in the Allo-Tetraploid Plant Nicotiana benthamiana

https://dl.sciencesocieties.org/publications/tpg/pdfs/8/1/plantgenome2014.10.0064

https://www.ncbi.nlm.nih.gov/pubmed/27054874

論文4  mergeしてパフォーマンスを上げる

https://www.ncbi.nlm.nih.gov/pubmed/28224052

 

 引用

SOAPdenovo-Trans: de novo transcriptome assembly with short RNA-Seq reads.

Xie Y, Wu G, Tang J, Luo R, Patterson J, Liu S, Huang W, He G, Gu S, Li S, Zhou X, Lam TW, Li Y, Xu X, Wong GK, Wang J.

Bioinformatics. 2014 Jun 15;30(12):1660-6.

 

 

Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels.

Schulz MH1, Zerbino DR, Vingron M, Birney E.

Bioinformatics. 2012 Apr 15;28(8):1086-92.

 

Full-length transcriptome assembly from RNA-Seq data without a reference genome.

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

Nature Biotechnology 29, 644–652 (2011)

 

De novo assembly and analysis of RNA-seq data.

Gordon Robertson, Jacqueline Schein[…]Inanc Birol

Nature Methods 7, 909–912 (2010).

 

Velvet: algorithms for de novo short read assembly using de Bruijn graphs.

D.R. Zerbino and E. Birney.

Genome Res. 2008 May;18(5):821-9.

 

De novo transcriptome assembly workflow

https://www.protocols.io/view/de-novo-transcriptome-assembly-workflow-ghebt3e?step=9