macでインフォマティクス

macでインフォマティクス

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

de novo transcriptome 解析のツールまとめ(主にde novo assemblerについて)

随時更新

2020 2/4 help更新

 

 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

 

 

追記

Pre-processing

De novo assembly

Hybrid de novo assembly

Genome guide assembly (よく似たゲノムが利用できる場合) 

Non-coding RNA

  • 真核と原核のSSU rRNA再構成にも使えるphyloFlash
  • 入力配列lからrRNAを検出するBarrnap
  • rRNAを除くSortMeRNA

Long read 

Clustering

Pan-transcriptome

Coding region

Merge

Annotation

GO enrichment analysis

Functional enrichment analysis

Pathway analysis

Evaluation

More comprehensive transcriptome

Visualization

Integrated Analysis Platform (web service)

Other

Database

 

 

 

ここでは代表的なRNA seqのアセンブルプログラムの導入と使い方についてまとめる。

1、Trinity 

Releases · trinityrnaseq/trinityrnaseq · GitHub

#bioconda (link
mamba create -n trinity python=3.9
conda activate trinity
#バージョン指定しないと古いバージョンが入ることがあるので、バージョンを指定するのが無難。2022年5月現在ならレポジトリは2.14が最新、conda配布版は2.13.2が最新(上のリンクから確認)。
mamba install -c bioconda -y trinity=2.13.2

#bowtie1も必要
conda install -c bioconda -y bowtie

#homebrew (試してません)
brew install Trinity

#dockerhub (link)
docker pull trinityrnaseq/trinityrnaseq:latest

Trinity

$ Trinity

 

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

#

# Required:

#

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

#

#  If paired reads:

#

#      --left  <string>    :left reads

#      --right <string>    :right reads

# 

#  Or, if unpaired reads:

#

#      --single <string>   :single reads

#

#

#  --output <string>     :name of directory for output (will be created if it doesn't already exist) 

#                                 default( "trinity_out_dir" )

#

#  if strand-specific data, set:

#  

#      --SS_lib_type <string>  :if paired: RF or FR,  if single: F or R

#

#  

#

#  Butterfly-related options:

# 

#    --bfly_opts <string>            :parameters to pass through to butterfly (see butterfly documentation).

#

#    --bflyHeapSpace <string>        :java heap space setting for butterfly (default: 1000M) => yields command java -Xmx1000M -jar Butterfly.jar ... $bfly_opts

#

#    --no_run_butterfly              :stops after the Chrysalis stage. You'll need to run the Butterfly computes separately, such as on a computing grid.

#                                       Then, concatenate all the Butterfly assemblies by running:

#                                             find trinity_out_dir/ -name "*allProbPaths.fasta" -exec cat {} + > trinity_out_dir/Trinity.fasta

#

#   

#  Inchworm-related options:

#

#    --no_meryl                      :do not use meryl for computing the k-mer catalog (default: uses meryl, providing improved runtime performance)

#    --min_kmer_cov <int>            :min count for K-mers to be assembled by Inchworm (default: 1)

#

#

#   Chyrsalis-related options:

#

#     --max_reads_per_graph <int>     :maximum number of reads to anchor within a single graph (default: 20000000)

#                                     

# Misc:

#

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

#

#  --min_contig_length <int> :minimum assembled contig length to report (def=200)

#

#  --paired_fragment_length <int>  :maximum length expected between fragment pairs (aim for 90% percentile)  (def=300)

#

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

#

#

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

 

 

ラン

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

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

A. Genome guided 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

B. 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 40 --output outdir

注意:Trinityはアップデートが続けられており、古いバージョンより最新のバージョンの方が連続性、精度共に向上する可能性が高いです。

Trinity --version

コマンドを打って、古いバージョンを使っていないか定期的に確認してください。

追記


 

 

 

2、Trans-ABySS

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

依存

#bioconda (link
conda install -c bioconda -y transabyss


#condaを使わないなら依存から順番に入れていく。python2環境でうごくので、pyenvやconda createでpython2.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 -h

$ transabyss -h

usage: transabyss [-h] [--version] [--se PATH [PATH ...]]

                  [--pe PATH [PATH ...]] [--SS] [--outdir PATH] [--name STR]

                  [--stage {dbg,unitigs,contigs,references,junctions,final}]

                  [--length INT] [--cleanup {0,1,2,3}] [--threads INT]

                  [--mpi INT] [-k INT] [-c INT] [-e INT] [-E INT] [-q INT]

                  [-Q INT] [-n INT] [-s INT] [--gsim INT] [--indel INT]

                  [--island INT] [--useblat] [--pid FLOAT] [--walk FLOAT]

                  [--noref]

 

Assemble RNAseq with Trans-ABySS.

 

optional arguments:

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

  --version             show program's version number and exit

 

Input:

  --se PATH [PATH ...]  single-end read files

  --pe PATH [PATH ...]  paired-end read files

  --SS                  input reads are strand-specific

 

Basic Options:

  --outdir PATH         output directory

                        [/Users/kazu/1/transabyss_2.0.1_assembly]

  --name STR            assembly name [transabyss] (ie. output assembly:

                        'transabyss-final.fa')

  --stage {dbg,unitigs,contigs,references,junctions,final}

                        run up to the specified stage [final]

  --length INT          minimum output sequence length [100]

  --cleanup {0,1,2,3}   level of clean-up of intermediate files [1]

 

ABySS Parameters:

  --threads INT         number of threads ('j' in abyss-pe) [1]

  --mpi INT             number of MPI processes ('np' in abyss-pe) [0]

  -k INT, --kmer INT    k-mer size [32]

  -c INT, --cov INT     minimum mean k-mer coverage of a unitig [2]

  -e INT, --eros INT    minimum erosion k-mer coverage [c]

  -E INT, --seros INT   minimum erosion k-mer coverage per strand [0]

  -q INT, --qends INT   minimum base quality on 5' and 3' ends of a read [3]

  -Q INT, --qall INT    minimum base quality throughout a read

  -n INT, --pairs INT   minimum number of pairs for building contigs [2]

  -s INT, --seed INT    minimum unitig size for building contigs [k]

 

Advanced Options:

  --gsim INT            maximum iterations of graph simplification [2]

  --indel INT           indel size tolerance [1]

  --island INT          minimum length of island unitigs [0]

  --useblat             use BLAT alignments to remove redundant sequences.

  --pid FLOAT           minimum percent sequence identity of redundant

                        sequences [0.95]

  --walk FLOAT          percentage of mean k-mer coverage of seed for path-

                        walking [0.05]

  --noref               do not include reference paths in final assembly

 

Written by Ka Ming Nip.

Copyright 2018 Canada's Michael Smith Genome Sciences Centre

Report bugs to <trans-abyss@googlegroups.com>

 

 ランチュートリアル

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

 

 

3、SOAPdenovo-Trans

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

ソースコード

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

BInary (linux)

SOAP :: Short Oligonucleotide Analysis Package

#bioconda (link
conda install -c bioconda -y soapdenovo-trans


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のランタイムはかなり短いが、メモリ要求量はかなり多め。


 

 4、Oases

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

https://github.com/dzerbino/oases

#bioconda (link)
conda install -c bioconda -y oases

#またはソースから。velvetのインストールについては、こちらのリンクのvelvetの項目を参照。
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 = 4

MAXKMERLENGTH = 191

OPENMP

LONGSEQUENCES

 

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

 

ラン

アセンブル。 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'

 

 

5、IDBA-tran

$ idba_tran

not enough parameters

IDBA-Tran - Iterative de Bruijn Graph Assembler for next-generation transcriptome sequencing data.

Usage: idba_tran -r read.fa -o output_dir

Allowed Options: 

  -o, --out arg (=out)                   output directory

  -r, --read arg                         fasta read file (<=600)

  -l, --long_read arg                    fasta long read file (>600)

      --mink arg (=20)                   minimum k value (<=312)

      --maxk arg (=60)                   maximum k value (<=312)

      --step arg (=10)                   increment of k-mer of each iteration

      --inner_mink arg (=10)             inner minimum k value

      --inner_step arg (=5)              inner increment of k-mer

      --prefix arg (=3)                  prefix length used to build sub k-mer table

      --min_count arg (=2)               minimum multiplicity for filtering k-mer when building the graph

      --min_support arg (=1)             minimum supoort in each iteration

      --num_threads arg (=0)             number of threads

      --seed_kmer arg (=30)              seed kmer size for alignment

      --min_contig arg (=200)            minimum size of contig

      --min_transcript arg (=300)        minimum size of transcript

      --similar arg (=0.95)              similarity for alignment

      --max_mismatch arg (=3)            max mismatch of error correction

      --no_local                         do not use local assembly

      --no_coverage                      do not iterate on coverage

      --no_correct                       do not do correction

      --pre_correction                   perform pre-correction before assembly

      --max_isoforms arg (=3)            maximum number of isoforms

      --max_component_size arg (=30)     maximum size of components

 

1、 マージしたfastaの作成

fastqはfastaとして与える必要がある。ペアエンドデータの場合、マージする必要もあるため、IDBAのfq2faコマンドを使ってマージしながらfastaに変換するマージする。fq2faコマンドはgzip圧縮fastqは受け付けないので、解凍してから指定する。

#ペアエンド
fq2fa --merge --filter pair_1.fq pair_2.fq read.fa
#シングルエンド
fq2fa single.fq read.fa

2、IDBA-tranのラン

24スレッド指定、k値は25から125まで5ずつ増やす。 相同性98%、サイズ300bp以上とする。precorrection実行。

idba_tran -r read.fa --num_threads 24 --pre_correction -o out_dir --mink 25 --maxk 125 --step 10 --min_transcript 300 --similar 0.98

 

 

 引用

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

 

参考資料

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

 

追記

評価のペーパーが出ています。コマンドや手順も載っているのでぜひ読んで見て下さい。


 こちらも読んでみて下さい。


A consensus-based ensemble approach to improve de novo transcriptome assembly

https://www.biorxiv.org/content/10.1101/2020.06.08.139964v1

 

f:id:kazumaxneo:20200915192446p:plain