macでインフォマティクス

macでインフォマティクス

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

polyploidのラージゲノムのアセンブラ Meraculous2

 

 ヒトや他のギガベース規模のゲノムの正確なディープショットガンシーケンスは、今や控えめなコストで容易に利用可能になっている。これらのシーケンシングスループットの増加により、大規模かつ複雑なゲノム用のショットガンシーケンスを構築するための新しい世代の計算アルゴリズムが開発された(Preprintより ref.1-3で概説されている)。これらのアプローチには、通常、Euler [ref.4]とVelvet [ref.5]によるショートリードアセンブリで先駆けされたde Bruijn graphアプローチが組み込まれている。いくつかのグループが、これらの方法でヒトおよび他の哺乳動物ゲノムを組み立てた[6-11]。このようなアセンブリは、ゲノムのサイズと反復性のためだけでなく、一塩基変異(SNV)、スモールおよびラージスケールの挿入および欠失(INDEL)、およびより大きな構造変化を包含する異種交配種の本質的なヘテロ接合性のためにも困難である。 大規模で反復的な多型ゲノムを組み立てるには、通常、大規模な共有メモリシステムが必要であり、実行には1週間以上かかり、1つのアセンブリを作成するため、かなりのコンピューティングリソースが必要になる。いくつかの大きなゲノムは純粋に短い(<150bp)ペアエンドシーケンスからアセンブリされているが、どちらのゲノムが良質のアセンブリを可能にする構造を有するか、および/またはどのデータの組み合わせがこれを促進するかは不明である [11、13、14]「Assemblathon」コンペは、アルゴリズムの直接比較を容易にするための共通データセットを提供する新しいアプローチのテストおよび実装に役に立つ[ref.11,13]。 Genome 10Kのような10,000個の脊椎動物ゲノムをアセンブルするようなラージゲノムぷZロジェクトでは[ref.15]、非常に効率的かつ精度の高いアセンブリツールの発達が不可欠である。

 以前著者らは、ハイブリッドk-mer / read-basedアセンブラとして"Meraculous"アルゴリズムを発表した。簡単に言うと、Meraculousは、簡略化されたde Bruijn k-mer graphを効率的に構築してトラバースすることによって、最初にゲノムの固有の領域を予備的な未解決の「UU」コンティグに組み立てる。次に、これらのコンティグをペアエンドデータとアライメントさせることによって連結し、残ったscaffoldsのギャップは、関連するローカルアセンブリを用いて埋められる。 Meraculousは現在のIlluminaシーケンスの高精度シーケンスを活用することで明示的なエラー訂正ステップを避ける。エラー訂正ステップは、アセンブリプロセスにとって冗長なプロセスであると著者らは考えている。最初のバージョンのMeraculousはpolymorphicな二倍体ゲノムに対応しておらず、並列化されているにもかかわらず、シミュレーションデータセットでは約1500万bpの一倍体カビゲノムしかレポートできず、ギガクラスラージゲノムには適用できなかった。

 ここでは、これらの以前の制限を克服して、Meraculous Assemblathon IIエントリーに組み込まれた変更を拡張する、Meraculousの改善について説明する。新機能には、(1)de Bruijn graph内の「バブル」構造の線形連鎖を用いた対立遺伝子多型の明示的な処理、(2)ケーススタディに基づくギャップ閉鎖の改善、および(3)より完全なアセンブリを生成する改良されたscaffoldingアルゴリズム、となる。新しい並列実装によって処理速度と帯域幅の効率が大幅に改善され、JGI Genepoolクラスタ上でリアルタイムに24時間以内にヒトゲノムのアセンブルが可能になった(この計算を実行するために使用されたリソースの詳細はPreprint 表1を参照) )。 Meraculous2のこれらの機能を調べるために、ここでは、ゲノムが決定され、1000ゲノムプロジェクトによって広範囲に分析されたヨーロッパ系祖先の女性である NA12878ヒトデータセットアセンブリへの適用について説明する。Trioフェージングによって(すなわち、NA12878と彼女の両親の配列を組み合わせること)、フェージングされた母系および父系のハプロタイプ配列が推定された(一部略)。 NA12878のゲノムは、以前にALLPATHS-LG [ref.9]とSGA [ref.10]によってアセンブリされているので、Meraculous2の性能をこれらの2つの最先端のアセンブラと比較した。我々(著者ら)は、このヒトセンブリにおいて以前には記載されていないいくつかの方法を提示し、長距離リンケージの正確さを測定するためのいくつかの測定基準について議論する。

  

PDFマニュアル

http://1ofdmq2n8tc36m6i46scovo2e.wpengine.netdna-cdn.com/wp-content/uploads/2014/12/Manual.pdf

以下のような特徴を持つ(マニュアルより)。

  1. Currently Meraculous works with Illumina data only. It relies on Illumina naming conventions and Phred-like sequence quality scores. Long-read/low-depth sequencing platforms are not supported at this time.
  2. An overall mean depth of read coverage of at least 30x is strongly recommended. Low-coverage datasets will likely result in a highly fragmented assembly or an aborted process altogether.
  3. Meraculous deals with genomic diploidy by creating a pseudo-haploid assembly where haplotypes are "squashed", i.e., a contig is formed with a single majority allele. However, the higher the polymorphic rate the less effective this process is. As a result, genomes with polymorphism rates of over 0.05 are better to assemble as haploid, letting Meraculous keep both haplotypes as distinct contigs, in essence imitating a meta- genome.
  4. Although it is capable of assembling small bacterial genomes, it may not be the most resource-efficient choice for these scenarios.
  5. Meraculous relies heavily on distributed and threaded computing and will perform best on a multiple-core server or in a cluster environment. For more on this, see sections 'Operating System requirements' and 'Hardware considerations'

メモリ使用量(マニュアルより)。

f:id:kazumaxneo:20180617193811j:plain

 

インストール

cent os6でテストした。

依存

  • Meraculous can run on any 64-bit Linux system
  • cmake >= 2.8
  • GCC g++ >= 4.4.7
  • GNU make 3.81
  • Boost C++ library >= 1.50.0
  • Perl (>= 5.10)
  • Log4perl.pm (>= 1.31 )
  • gnuplot (>= 3.7)
  • qqacct (optional but highly recommended for Grid Engine cluster environments)

Log4perl.pmは手っ取り早く"cpanm Log4perl.pm"で導入した。

本体 SourceFroge

https://sourceforge.net/projects/meraculous20/

tar -xvf Meraculous-v2.2.5.1-1-ga103cd6.tar
cd Meraculous-v2.2.5.1-1-ga103cd6/
install.sh <installation directory>

./run_meraculous.sh

$ ./run_meraculous.sh

Smartmatch is experimental at /home/uesaka/test2/Meraculous-v2.2.5.1-1-ga103cd6/bin/meraculous.pl line 929.

Smartmatch is experimental at /home/uesaka/test2/Meraculous-v2.2.5.1-1-ga103cd6/bin/meraculous.pl line 2789.

 

Command line arguments for meraculous.pl (Version 2.2.5.1):

 

meraculous.pl 

 

  Required:

  -c|config <config file>         : user config file 

 

  Optional:

  -label   <label>                : provide a label name for new runs ( Default: 'run' )

  -dir     <directory>            : provide a run directory name  ( Default: latest run )

  -restart                        : restart a previously failed assembly

  -resume                         : restart but preserve any partial results    

  -step                           : execute one stage and stop

  -start   <stage>                : re-run starting with this stage

  -stop    <stage>                : stop after this stage

  -archive                        : save any old stage directories (valid only with -restart)

  -cleanup_level [0|1|2]          : decide how agressively the pipeline should clean up intermediate data ( Default: 1)

                                      0 - do not delete any intermediate outputs (disk space footprint may be huge)

                                      1 - delete files that are not used in any of the subsequent stages and that are generally not informative to the user

                                      2 - delete files as soon as possible.  WARNING!!! You will not be able to rerun the

                                          stages individually once they have completed!

 

  -h|help                         : you guessed it: this usage page

  -v|version                      : about this program

 

 

The default configuration file is 'meraculous.params', which must be present

 

The default label name is  <genus>_<species>_[strain] if these are defined in 

the configuration file, and 'run' otherwise;

 

-resume/-restart : If no directory is given, the most recently run dir. is used.

 

Invalid command line combinations:

  -restart with     -resume

  -label   with    -restart or -resume

  -start   without -restart or -resume

  -archive without -restart

 

 

Please contact Eugene Goltsman at egoltsman@lbl.gov if you encounter any problems.

  

 

 

ラン

テストランを実行する。

cd /etc/meraculous/test/pipeline bash
bash /bin/run_meraculous.sh -c meraculous.config

configファイルを指定している。lib_seqのところでfastqとそのパラメータを指定している。詳細はPDFマニュアルのBasic assembly parameters:の項を参照。

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

#

#  Meraculous params file

#

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

 

 

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

#

# Basic parameters

#

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

 

 

 

# Describe the libraries ( one line per library )

# lib_seq [ wildcard ][ prefix ][ insAvg ][ insSdev ][ avgReadLen ][ hasInnieArtifact ][ isRevComped ][ useForContigging ][ onoSetId ][ useForGapClosing ][ 5pWiggleRoom ][3pWiggleRoom]

#

 

lib_seq frags.fastq.25K  FRA  180  10   101 0 0   1 1 1  0 0

lib_seq jumps.fastq.25K  JMP  3000 500  101 1 1   1 2 1  0 0

 

 

genome_size   0.0002

 

diploid_mode 0

 

mer_size 21

 

min_depth_cutoff 7

 

num_prefix_blocks 4

 

 

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

#

# Advanced parameters 

#

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

 

no_read_validation 1

 

use_cluster 0

 

local_num_procs       4

 

local_max_retries   0

 

 

 

10分程度でテストランはおわる。"ran to completion successfully" が表示されれば、正常にランできている。

> ls -alt

f:id:kazumaxneo:20180617192118j:plain

ステージごとにサブフォルダが作られる。

 

実際にアセンブリを行うには、テストと同様にconfigファイルを準備すればよい。

run_meraculous.sh -c <config file>

 

 

 

 

引用

Meraculous2: fast accurate short-read assembly of large polymorphic genomes

Jarrod A. Chapman, Isaac Y. Ho, Eugene Goltsman, Daniel S. Rokhsar

Submitted on 2 Aug 2016 (v1), last revised 7 Nov 2017 (this version, v2)

https://arxiv.org/abs/1608.01031

 

Meraculous-2D: Haplotype-sensitive Assembly of Highly Heterozygous genomes

Eugene Goltsman, Isaac Ho, Daniel Rokhsar

(Submitted on 29 Mar 2017)

https://arxiv.org/abs/1703.09852

 

Meraculous: de novo genome assembly with short paired-end reads.

Chapman JA, Ho , Sunkara S, Luo S, Schroth GP, Rokhsar DS.

PLoS One. 2011;6(8):e23501.

 

DACCOR

 

 シーケンシングリードからのゲノムの再構築は、デノボアセンブリによって達成でき、重複したリードが同定され、コンティグと呼ばれるより長い連続配列に拡張される。あるいは、highly closedなリファレンスゲノムが利用可能であれば、このゲノムに対してリードをマッピングでき、再構成はリードのコンセンサスに基づく。

 マッピングベースのアプローチでは、BWA(Li&Durbin、2009)やBowtie2(Langmead&Salzberg、2012)などのプログラムを使用して、次世代シーケンシング(NGS)技術によって生成されたリードを既知リファレンスゲノムにマッピングする(Veeramah&Hammer、2014)。次いで、サンプルが十分なカバレッジ深度でシーケンスされたと仮定して、アライメントされたリードのコンセンサス配列を使用して、新たに配列決定されたサンプルのゲノム配列を生成する。これにより、短い挿入、欠失、一塩基変異(SNV)または一塩基多型(SNP)の迅速な同定が可能になる。マッピングプログラムは、通常、アラインメントクオリティに相当する各リードのアライメントスコアを計算する(Li、Ruan&Durbin、2008)。スコアは、リードが正しいゲノムポジションに置かれている確率を定量化する。低いマッピングクオリティのリードはコンタミや低いクオリティでシーケンスされたリードに由来しており、フィラタリングで除外できる(Smith、Xuan&Zhang、2008)。低クオリティのリードの他に、リピート領域へのマッピングも、ユニークなポジションにマッピングできない場合、低クオリティのスコアを生成する可能性がある。マッピングクオリティが低いリードのフィルタリングには、これらのリードも含まれる。このフィルタリングはancient DNA(aDNA)(Bos et al、2016)で行われることが多く、このようなサンプルの場合、それぞれの再構成ゲノムのリピート領域が一般に影響を受け、多くの未知(N)が生じる。また、デノボアセンブリアプローチは、リピート領域を再構成する際に問題があり、少なくともリード長よりもはるかに長いものがある(Simpson&Durbin、2012)。

 しかしながら、リピート領域は多くのゲノムにおいて必須の役割を果たす(Shapiro&von Sternberg、2005)。数百から数千ものそのような領域が原核および真核染色体に存在する(Treangen et al、2009)。ヒトゲノムは、例えば、約50%がリピート領域からなる(International Human Genome Sequencing Consortium、2001)。バクテリアのタンデムリピート領域は、外膜タンパク質と関連しているようであり、病原体が宿主に適応するのを助けることを示唆している(Denoeud&Vergnaud、2004)(一部略)。

 Illumina SLRプラットフォーム、PacBio、Oxford Nanoporeなどの新しいシーケンシング技術は、ほとんどの反復領域をカバーできるロングリードを作成することができる(Huddleston et al、2014)。しかし、これらの技術をDNAサンプルに適用することは必ずしも可能ではない。例えば、aDNAプロジェクトでは、DNAの平均抽出断片長は約44〜72塩基対である(Sawyer et al、2012)。 Illuminaのショートリード技術でも断片全体をシーケンスできるため、ロングリード技術でこれらの断片をシーケンシングしても情報は得られない。短いDNA断片はaDNAに限定されない。これらはまた、培養困難なT.pallidumのような配列にも現れる(Arora et al、2016 pubmed)。

 ショートリードを使用してリピート領域をよりよく解決するために、研究者は最初に、マッピングの前に重複領域と低複雑領域をマスクする(Frith、Hamada&Horton、2010)。 1つの例は、マスクされたバージョンが既に利用可能なヒトリファレンスゲノムである(University of California Santa Cruz(UCSC)、2014)。マスクされたゲノムがない場合、RepeatMasker(Smitt、Hubley&Green、1996)のようなプログラムでリピート領域を特定し、マスクされたリファレンスを作成することができる。 RepeatMaskerは、既知のリピート領域のライブラリを使用し、入力シーケンスと比較する。これはゲノムのマスキングを可能にするが、リピート領域の新規同定は不可能である。

 リピート領域を新しく識別できるプログラムの1つはVMatch(Kurtz、2003)である。それはリピート領域を特定するためにprefix配列(Weiner、1973)を使用する。 VMatchは、マスキングタスク(Assuncao et al。、2010)と同様に、繰り返し領域のアノテーション(例えば、Lindow&Krogh(2005))のために複数のゲノムプロジェクトに適用されている。もちろん、GalaxyのウェブツールであるRepeatExplorer(Novak et al、2013)やシーケンスリード基づくRepARK(Koch、Platzer&Downie、2014)のようなデノボの繰り返し検索プログラムもある。現在までユーザーがNGSデータからリピート同定とゲノムの再構成を組み合わせることができるツールは存在しない。

 著者らのアプローチの一般的な考え方は、最初に、与えられたリファレンスゲノム中のすべてのリピート領域の新規同定から始まる。 NGSデータ中のリードは、同定されたすべてのリピート領域と同様に、リファレンスゲノムに対して個別にマッピングされ、再構成する。各試料の再構成ゲノムの塩基対分解能を高めるために、個々に再構成されたリピート領域を用いて、ゲノムの対応する再構築されたリピート領域を置換する。これらゲノム中のリピート領域の検出、特徴付けおよび再構成すべてのステップは、DACCORという全自動ツールで結びつけられる。

 

f:id:kazumaxneo:20180615213456j:plain

ワークフロー。論文より転載。

  

インストール

cent os6とmac os10.12でテストした。

 本体 Github

https://github.com/Integrative-Transcriptomics/Daccor

git clone https://github.com/Integrative-Transcriptomics/Daccor.git
cd Daccor/Releases/

java -jar DACCOR-1.0.0.jar 

$ java -jar DACCOR-1.0.0.jar 

Program: DACCOR

Version: 1.0.0

 

usage: Daccor <tool> <options>

--HELP--

Tool Description

 

identify Analyze a reference file and identify repetitive regions

reconstruct Reconstruct each region separately with EAGER

combine Combine the reconstructed regions with the reconstructed references

pipeline Pipeline from the identification of repetitive regions

to the combined reconstruction of the genome

--End of HELP--

java -jar DACCOR-1.0.0.jar identify

$ java -jar DACCOR-1.0.0.jar identify

usage:   [-g <arg>] [-h] -i <arg> [-k <arg>] [-m <arg>] [-M <arg>] -o <arg> [-p

       <arg>] [-rl <arg>] [-S] [-t <arg>] [-v <arg>] [-ws]

 

--HELP--

     -g,--gff <arg>           Path of GFF file

     -h,--help                Prints options

     -i,--input <arg>         input filename, REQUIRED

     -k,--kmersize <arg>      size of initial kmer [readlength/2 OR 17]

     -m,--mismatches <arg>    number of mismatches allowed [0]

     -M,--margin <arg>        Margin around repeat to extract [0]

     -o,--output <arg>        output filename, REQUIRED

     -p,--processes <arg>     Number of threads per genome [1]

     -rl,--readlength <arg>   readlength [17]

     -S,--separately          Analyze each sequence separately [false]

     -t,--threshold <arg>     Min length of displayed results [readlength OR 51]

     -v,--vmatch <arg>        don't analyze but parse vmatch output

     -ws,--writeSeparately    write entries into separate files [false]

--End of HELP--

> java -jar DACCOR-1.0.0.jar reconstruct

$ java -jar DACCOR-1.0.0.jar reconstruct

usage:   [-e <arg>] -f <arg> [-h] -i <arg> -o <arg> -r <arg>

 

--HELP--

     -e,--eager <arg>     config file for reconstructed samples already

                          reconstructed with EAGER

     -f,--fastqs <arg>    config file for fastq files, REQUIRED

     -h,--help            Prints options

     -i,--input <arg>     input filename (reference), REQUIRED

     -o,--output <arg>    output folder, REQUIRED

     -r,--repeats <arg>   config file for repeat input files, REQUIRED

--End of HELP--

>java -jar DACCOR-1.0.0.jar combine

$ java -jar DACCOR-1.0.0.jar combine

usage:   -g <arg> [-h] -o <arg> -r <arg>

 

--HELP--

     -g,--genomes <arg>   config file for reconstructed genomes, REQUIRED

     -h,--help            Prints options

     -o,--output <arg>    output folder, REQUIRED

     -r,--repeats <arg>   config file for reconstructed repeats, REQUIRED

--End of HELP--

 

 

ラン

1、Identify repeats de novo

リファレンスゲノムGCF_000008605.1_ASM860v1_genomic.fnaのリピートをde novo検出する(examples/data/identify/GCF_000008605.1_ASM860v1_genomic.fna)。

java -jar DACCOR-1.0.0.jar identify -i GCF_000008605.1_ASM860v1_genomic.fna -k 31 -m 5 -o Nichols_repeats -p 10R -t 101

head Nichols_repeats_k31_Mis5_Thres101.bed

 $ head Nichols_repeats_k31_Mis5_Thres101.bed 

NC_000919.1 134890 136735 repeat_0 0 +

NC_000919.1 135149 135365 repeat_1 0 +

NC_000919.1 151090 152935 repeat_0 0 +

NC_000919.1 151349 151565 repeat_1 0 +

NC_000919.1 229895 230295 repeat_2 0 +

NC_000919.1 230294 231763 repeat_3 0 +

NC_000919.1 231827 235110 repeat_4 0 +

NC_000919.1 278328 278728 repeat_2 0 +

NC_000919.1 278728 280197 repeat_3 0 +

NC_000919.1 280271 283554 repeat_4 0 +

> head -n 15 Nichols_repeats_k31_Mis5_Thres101_Summary.csv |column -t

 $ head -n 15 Nichols_repeats_k31_Mis5_Thres101_Summary.csv |column -t

>NC_000919.1                                                                                               Treponema   pallidum    subsp.      pallidum    str.       Nichols     complete  genome

category                                                                                                   value

Total                                                                                                      number      of          repetitive  regions     23

Total                                                                                                      number      of          different   repetitive  regions    11

Minimum                                                                                                    length      115

Maximum                                                                                                    length      3283

Average                                                                                                    length      965

Number                                                                                                     of          repetitive  regions     containing  mutations  2

Number                                                                                                     of          repetitive  regions     also        contained  in          other     sequences  0

fastaName                                                                                                  repeatName  startindex  endIndex    length      number     mismatches  indices   of         mismatches

NC_000919.1Treponemapallidumsubsp.pallidumstr.Nicholscompletegenome_pos_[134890,151090]_length_1845        repeat_0    134890      136735      1845        0

NC_000919.1Treponemapallidumsubsp.pallidumstr.Nicholscompletegenome_pos_[135149,151349,671306]_length_216  repeat_1    135149      135365      216         0

NC_000919.1Treponemapallidumsubsp.pallidumstr.Nicholscompletegenome_pos_[134890,151090]_length_1845        repeat_0    151090      152935      1845        0

NC_000919.1Treponemapallidumsubsp.pallidumstr.Nicholscompletegenome_pos_[135149,151349,671306]_length_216  repeat_1    151349      151565      216         0

NC_000919.1Treponemapallidumsubsp.pallidumstr.Nicholscompletegenome_pos_[229895,278328]_length_400         repeat_2    229895      230295      400         0

——

 

2、Reconstruct identified repetitive regions individually

スクリプト中のfastq-dumpのパスを修正してからランする。

./runReconstructExample.sh

以下のような実行内容になっている。

cat runReconstructExample.sh 

$ cat runReconstructExample.sh 

#!/bin/sh

 

# location of the fastq-dump program from the SRA toolkit(https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/)

fastqdump=/usr/local/bin/fastq-dump

 

# download the syphilis sample SRR3584843

echo "downloading fastq files"

 

if ! -f data/SRR3584843_R1.fastq.gz || ! -f data/SRR3584843_R2.fastq.gz ; then

$fastqdump --gzip -O data --split-3 SRR3584843

mv data/SRR3584843_1.fastq.gz data/SRR3584843_R1.fastq.gz

mv data/SRR3584843_2.fastq.gz data/SRR3584843_R2.fastq.gz

fi

 

# setting the variables

genome=../examples/data/identify/GCF_000008605.1_ASM860v1_genomic.fna

fqConfig=../examples/data/exampleFQConfig.txt

outPrefix=SRR3584843

repeatConfig=../examples/data/exampleRepeatConfig.txt

 

echo "starting reconstruction"

java -jar Daccor-0.0.1.jar reconstruct -f $fqConfig -i $genome -o $outPrefix -r $repeatConfig

#if ! -f  のところは存在すれば真(参考)。

前半はfastqの調達作業。中盤は変数やパス定義。最終的に実行しているのは以下のコマンドになる。

java -jar Daccor-0.0.1.jar reconstruct -f exampleFQConfig.txt -i GCF_000008605.1_ASM860v1_genomic.fna-o SRR3584843 -r exampleRepeatConfig.txt

mapping base assembly toolのeagerの読み込みでエラーになる。修正できたら追記します。 

 

3、Combine reconstructed sequences with reconstructed genome

./runCombineExample.sh

 

 

引用

DACCOR-Detection, characterization, and reconstruction of repetitive regions in bacterial genomes

Seitz A, Hanssen F, Nieselt K

PeerJ. 2018 May 29;6:e4742.

 

 

 

クロロプラストゲノムの自動アセンブリパイプライン Fast-Plast

 

 Fast-Plastは、既存および新規のプログラムを活用して、葉緑体ゲノム全体を迅速にアセンブリし、検証するパイプライン。 十分なデータを持つほとんどのデータセットについて、Fast-Plastは自動で完全長の葉緑体ゲノムアセンブリを生成できる。 Fast-Plastは、葉緑体配列生成に加え、存在する葉緑体遺伝子を同定する。Fast-PlastはIlluminaデータを使う。

 Fast-Plastは、de Bruijn graphベースのSPAdesとafinで実装された低いカバレッジでコンティグのギャップを埋めるiterative なシードベースのアセンブリを組み合わたアプローチを使っている。 パイプラインは、quadripartiteな葉緑体ゲノム構造から領域を同定し、慣習に従い順番をアサインする。 最終的なアセンブリ品質を評価するためにカバレッジ分析が実行される。

ゲノムスキミング(従来のゲノムプロジェクトとの比較(pubmed) )やlow pass シーケンス解析(Biostars)を想定して構築されたプログラムになっている。キャプチャしてエンリッチしてシーケンスしたデータを入力とする際は、bbnormやkhmerでノーマイライズしてから使うことが推奨されている。

 

インストール

cent os6でテストした。

依存

 本体 Github

https://github.com/mrmckain/Fast-Plast

git clone https://github.com/mrmckain/Fast-Plast.git
cd Fast-Plast/
perl INSTALL.pl #依存を導入するか聞いてくる。質問に"All"と答えると全導入される。

perl fast-plast.pl 

$ perl fast-plast.pl 

ERROR: Missing reads file(s).

Usage:

        fast-plast.pl [-1 <paired_end_file1> -2 <paired_end_file2> || -single <singe_end_file>] -name <sample_name> [options] 

    or

        fast-plast.pl -help

            -1 <filenames>          File with forward paired-end reads. Multiple files can be designated with a comma-delimited list. 

                                    Read files should be in matching order with other paired end files.

            -2 <filenames>          File with reverse paired-end reads. Multiple files can be designated with a comma-delimited list. 

                                    Read files should be in matching order with other paired end files.

            -s <filenames>          File with unpaired reads. Multiple files can be designated with a comma-delimited list.

 

            PAIRED END AND SINGLE END FILES CAN BE PROVIDED SIMULTAENOUSLY.

 

            -n <sample_name>        Name for current assembly. We suggest a species name/accession combination as Fast-Plast will use 

                                    this name as the FASTA ID in the final assembly.

 

    Advanced options:

 

            --threads               Number of threads used by Fast-Plast.  [Default = 4]

            --adapters              Files of adapters used in making sequencing library. Users can select "Nextera" for Nextera adapters, "TruSeq" for TruSeq adapters, leave the default (NEB), or provide their own. [Default = NEB-PE]

            --bowtie_index          Order for sample to draw references for mapping. If order exists, then all available samples for that order will be used. 

                                    If order does not exist in default set or the terms "all" or "GenBank" are given, one exemplar from each available order is used 

                                    to build the Bowtie2 indicies. [default="All"]

            --user_bowtie           User supplied bowtie2 indices. If this option is used, bowtie_index is ignored.

            --posgenes              User defined genes for identification of single copy/IR regions and orientation. Useful when major rearrangments are present in user plastomes.

            --coverage_analysis     Flag to run the coverage analysis o

 

ラン

参照とする生物種の目(order)を指定する。イネ目ならPoales。以下のリストが現在指定可能になっている。

f:id:kazumaxneo:20180616212416j:plain

ペアエンドfastqを指定してランする。--bowtie_indexでゲノムを指定する。Genome skimmngを想定して素早く自動で結果を出すことを念頭に置いており、fastqのアダプターやクオリティトリミングは未実行を想定してプロセスは進められる。

perl fast-plast.pl -1 pair1.fastq -2 pair2.fastq --name output --bowtie_index Poales --coverage_analysis --clean light --threads 12

> ls -alth output/

$ ls -alth output/

total 36K

drwxr-xr-x  2 uesaka user 4.0K Jun 16 21:56 Final_Assembly

drwxr-xr-x  8 uesaka user 4.0K Jun 16 21:56 .

drwxr-xr-x  2 uesaka user 4.0K Jun 16 21:56 4_Afin_Assembly

drwxr-xr-x  2 uesaka user 4.0K Jun 16 21:56 5_Plastome_Finishing

drwxr-xr-x  3 uesaka user 4.0K Jun 16 21:52 3_Spades_Assembly

-rw-r--r--  1 uesaka user  211 Jun 16 21:52 output3_Plastome_Summary.txt

drwxr-xr-x  2 uesaka user 4.0K Jun 16 21:52 2_Bowtie_Mapping

drwxr-xr-x  2 uesaka user 4.0K Jun 16 21:51 1_Trimmed_Reads

drwxr-xr-x 12 uesaka user 4.0K Jun 16 21:51 ..

 

解析ステップごとにサブフォルダができる。最終出力のFASTAはFinal_assembly/に保存される。1つのcontigまでアセンブリされると、5_Plastome_Finishing/ができる。出力の詳細はGithubで確認してください。

https://github.com/mrmckain/Fast-Plast

 

bin/でビルドされたspades.pyがエラーを起こしたので、手っ取り早くspadesのversion3.12のシンボリックリンクspades.pyと置き換えてランしました。応急処置でしたが、一応うまくランできます。

 

引用

https://github.com/mrmckain/Fast-Plast

 

de novo transcriptomeのcontigクラスタリングツール Corset

 

 RNA-seqは、トランスクリプトームの様々な側面を研究するための強力な技術である。それは、遺伝子発見、選択的スプライシングイベントの検出、継時的発現分析、融合の検出、SNPおよび転写後エディティングなどの変異の同定を含む広範囲の用途を有する[ref.1]、[ref.2]。マイクロアレイのような従来技術に比べてRNA-seqの利点の1つは、データの生成と分析にリファレンスゲノムとアノテーションが必要ないため、非モデル生物のトランスクリプトーム全体の分析が可能になることである。リファレンスゲノムが利用できない場合、トランスクリプトームはRNA-seqのリードから直接de novo assemblyされる[ ref.3]。新たなトランスクリプトームアセンブリのためのいくつかのプログラムが存在する:Velvet [ref.6]およびAbyss [ref.7]ゲノムアセンブラを拡張するOases [ref.4]およびTrans-abyss [ref.5]ならびにTrinity [ref.8] ]。これらのプログラムは、コンティグ(contig)と呼ばれる転写配列へ数百万のショートリードをアセンブリすることができる。

 RNA-seqの一般的かつ生物学的に重要な用途の1つは、2つ以上の条件の間で発現が異なる遺伝子を同定することである[ref.9]。しかし、デノボアセンブリされたトランスクリプトームに対するdifferential expression解析を行うことは、遺伝子あたり複数のコンティグが報告されているため困難である。トランスクリプトームアセンブラーは、同じ遺伝子のアイソフォームを区別し、それぞれを別々に出力するため、共通の配列を持つ複数のコンティグが生じる。さらに、しばしば、異なるアイソフォームを真に代表してはいないコンティグを作るが、これは、二倍体またはプールされた集団内のシークエンシングエラー、リピートのカバレッジ変動や変異などのアーティファクトから生じている。その結果、トランスクリプトームのアセンブラは、転写物の断片化されたバージョン、またはSNPかindelだけ異なるコンティグを繰り返し報告することがよくある。驚くべきことに、シークエンシングエラー、SNP、または選択的スプライシングを伴わないデータのアセンブリでさえ、遺伝子ごとに複数のコンティグを生成することがシミュレーションによって示されている[ref.10]。したがって、デノボアセンブリによって生成されるコンティグの数は、典型的には大きい。例えば、8,000万リードのアセンブリから数十万のコンティグが生じる[ref.11]。

 de novoトランスクリプトームアセンブリによって生成される必然的に長いコンティグのリストは、differential expression解析のためにいくつかの問題を引き起こす:i)リードは重複配列に対して明確にアラインできないため、リードのアライメント起源の決定には誤りが多くなる。 ii)より多くのコンティグ間でリードが割り当てられなければならないので、コンティグあたりの平均カウントが減少しdifferential expressionの統計試験の能力が低下する、; iii)複数テストの調整がより厳しくなる。そしてiv)differential expressionのコンティグが同定されても、多数の遺伝子がリストに複数回出現するので、解釈が困難となる。Differential expression解析をコンティグよりも遺伝子を使い行えば、これらの困難は克服されるであろう。しかしながら、de novo アセンブリされたコンティグセットから遺伝子レベルの発現を評価するための手順は、これまでの文献では十分に検討されていない。

 De novoアセンブリされたトランスクリプトームからdifferential expression genes (DEG) を同定するには、いくつかの段階がある(論文 図1(下の図))。最初にアセンブリされ、リードはコンティグへマッピングされる。コンティグはクラスター化され遺伝子にし、その後、各遺伝子の発現レベルを解析する(図1)。DEGを検出するために統計試験を行う。

 いくつかの研究では、この分析パイプラインの個々のステップを比較している。例えば、様々なアセンブラの検討や、アセンブリ前の品質管理のメリットが検討されている[ref.12]〜[ref.15]。同様に、リードカウントベースの統計的な検定を行う方法も比較評価されている[ref.16]、[ref.17]。しかし、トランスクリプトームアセンブリから遺伝子レベルのカウントを得るための系を比較したり、提案したりする研究はほとんどない[ref.18]、[ref.19]。(一部略)

 この論文では、任意のde novoトランスクリプトームアセンブリから遺伝子レベルのカウントを得るための方法とソフトウェアであるCorsetを紹介する。Corsetは、新しくアセンブリされたトランスクリプトーム配列にマルチマッピングされたリードのセットを受け取り、共通するリードの比率および発現パターンに基づいてコンティグを階層的にクラスタリングする。発現パターンが異なる場合、リードの配列がシェアされているパラログなどの遺伝子間を区別可能である。マップされたリードを使用して、Corsetは遺伝子レベルのカウントを出力する。次いで、遺伝子レベルのカウントは、edgeRおよびDESeqなどのカウントベースのフレームワークを使用して、differential expressionの試験を容易に実行できる。著者らは、Corsetが、複数の評価基準に従って、ほかのクラスタリング手法と比較して一貫して優れた性能を発揮することを示している。Corsetはアセンブラに依存しないメソッドであるため、さまざまなソースからコンティグや転写産物を組み合わせることができる。また、Corsetはクラスタリングとカウントのステップが1回実行で行えるため、使うのが簡単である。

f:id:kazumaxneo:20180614221350j:plain

非モデル生物でリードカウントベースのDEG検出を行う流れ。論文より転載。

 

簡単に言うと、同じ転写産物には、contigが複数に分かれていてもリードが冗長にマッピングされうることを利用してクラスタリングする手法です。複数mappingが全て報告されているbamを使います。

 

インストール

cent os6とmac os10.12でテストした。

 本体 Github

https://github.com/Oshlack/Corset

git clone https://github.com/Oshlack/Corset.git
cd Corset/
./configure
make
make install


#condaでも導入できる(リンク)。
conda update corset

> ./corset

$ ./corset

 

Running Corset Version 1.06

No input files specified

 

Corset clusters contigs and counts reads from de novo assembled transcriptomes.

 

Usage: corset [options] <input bam files>

 

Input bam files:

The input files should be multi-mapped bam files. They can be single, paired-end or mixed

and do not need to be indexed. A space separated list should be given.

e.g. corset sample1.bam sample2.bam sample3.bam

or just: corset sample*.bam

 

If you want to combine the results from different transcriptomes. i.e. the same reads have 

been mapped twice or more, you can used a comma separated list like below:

corset sample1_Trinity.bam,sample1_Oases.bam sample2_Trinity.bam,sample2_Oases.bam ...

 

Options are:

 

-d <double list> A comma separated list of distance thresholds. The range must be

                  between 0 and 1. e.g -d 0.4,0.5. If more than one distance threshold

                  is supplied, the output filenames will be of the form:

                  counts-<threshold>.txt and clusters-<threshold>.txt 

                  Default: 0.3

 

-D <double>      The value used for thresholding the log likelihood ratio. The default 

                  value will depend on the number of degrees of freedom (which is the 

                  number of groups -1). By default D = 17.5 + 2.5 * ndf, which corresponds 

                  approximately to a p-value threshold of 10^-5, when there are fewer than

                  10 groups.

 

-m <int>         Filter out any transcripts with fewer than this many reads aligning.

                  Default: 10

 

-g <list>        Specifies the grouping. i.e. which samples belong to which experimental

                  groups. The parameter must be a comma separated list (no spaces), with the 

                  groupings given in the same order as the bam filename. For example:

                  -g Group1,Group1,Group2,Group2 etc. If this option is not used, each sample

                  is treated as an independent experimental group.

 

-p <string>      Prefix for the output filenames. The output files will be of the form

                  <prefix>-counts.txt and <prefix>-clusters.txt. Default filenames are:

                  counts.txt and clusters.txt

 

-f <true/false>  Specifies whether the output files should be overwritten if they already exist.

                  Default: false

 

-n <string list> Specifies the sample names to be used in the header of the output count file.

                  This should be a comma separated list without spaces.

                  e.g. -n Group1-ReplicateA,Group1-ReplicateB,Group2-ReplicateA etc.

                  Default: the input filenames will be used.

 

-r <true/true-stop/false> 

                  Output a file summarising the read alignments. This may be used if you

                  would like to read the bam files and run the clustering in seperate runs

                  of corset. e.g. to read input bam files in parallel. The output will be the

                  bam filename appended with .corset-reads.

                  Default: false

 

-i <bam/corset>  The input file type. Use -i corset, if you previously ran

                  corset with the -r option and would like to restart using those

                  read summary files. Running with -i corset will switch off the -r option.

                  Default: bam

 

-l <int>         If running with -i corset, this will filter out a link between contigs

                  if the link is supported by less than this many reads. Default: 1 (no filtering)

 

-x <int>         If running with -i corset, this option will filter out reads that align to more

                  than x contigs. Default: no filtering

 

Citation: Nadia M. Davidson and Alicia Oshlack, Corset: enabling differential gene expression 

          analysis for de novo assembled transcriptomes, Genome Biology 2014, 15:410

 

Please see https://github.com/Oshlack/Corset/wiki for more information

 

 

 

ラン

最低bamだけ指定すればランできる。

corset input.bam

cluster.txtが出力される。

 

 ここでは公式の例に従いテストランを行う。データのダウンロードからアセンブリマッピングクラスタリング、DEG検出まで行う本格的なものになっている。yeastのRNA seqデータを使う。

 

1、ダウンロード

ダウンロードしてFastqに変換。

cat >sra_run_ids.txt <<EOF 
SRR453566
SRR453567
SRR453568
SRR453569
SRR453570
SRR453571
EOF

#ダウンロード
prefetch --option-file sra_run_ids.txt
mv ~/ncbi/public/sra/SRR4535* sra_files/

 forループでfastqに変換。

#fastqディレクトリを作成
mkdir fastq

for srr in $(cat sra_run_ids.txt)
do
fastq-dump --split-3 --skip-technical --readids
--defline-seq '@$sn/$ri' --defline-qual '+' sra_files/${srr}.sra
-O fastq
done
  • --split-3    Legacy 3-file splitting for mate-pairs:  First biological reads satisfying dumping conditions are placed in files *_1.fastq and *_2.fastq If only one biological read is  present it is placed in *.fastq Biological reads and above are ignored. 

> ls -alh fastq/

$ ls -alh fastq/

total 16G

drwxr-xr-x 2 uesaka user 4.0K Jun 15 12:16 .

drwxr-xr-x 5 uesaka user 4.0K Jun 15 12:10 ..

-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:13 SRR453566_1.fastq

-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:13 SRR453566_2.fastq

-rw-r--r-- 1 uesaka user 1.7G Jun 15 12:14 SRR453567_1.fastq

-rw-r--r-- 1 uesaka user 1.7G Jun 15 12:14 SRR453567_2.fastq

-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:15 SRR453568_1.fastq

-rw-r--r-- 1 uesaka user 1.3G Jun 15 12:15 SRR453568_2.fastq

-rw-r--r-- 1 uesaka user 916M Jun 15 12:15 SRR453569_1.fastq

-rw-r--r-- 1 uesaka user 916M Jun 15 12:15 SRR453569_2.fastq

-rw-r--r-- 1 uesaka user 1.5G Jun 15 12:16 SRR453570_1.fastq

-rw-r--r-- 1 uesaka user 1.5G Jun 15 12:16 SRR453570_2.fastq

-rw-r--r-- 1 uesaka user 1.4G Jun 15 12:17 SRR453571_1.fastq

-rw-r--r-- 1 uesaka user 1.4G Jun 15 12:17 SRR453571_2.fastq

目に見える問題は起きてない。

 

2、クオリティトリミング

trimmomatic(brewで導入可)を使っているが、ペアエンドを扱えるなら他のツールでもOK。

for srr in $(cat sra_run_ids.txt)
do
trimmomatic PE -phred33 fastq/${srr}_1.fastq fastq/${srr}_2.fastq fastq/${srr}_P1.fastq fastq/${srr}_U1.fastq
fastq/${srr}_P2.fastq fastq/${srr}_U2.fastq
LEADING:20 TRAILING:20 MINLEN:50
done

Pはペアエンド、Uはアンペア。アンペアになったリードは使わない。

 

3、de novoアセンブリ

全データをマージし、Trinityを使いde novoアセンブリを実行する。

#マージ。公式では、このstepでsedを使い、Trinityのランに必要な/1と/2を付加している。今回はダウンロード時に付加済み
cat fastq/*_P1.fastq > fastq/all_1.fastq
cat fastq/*_P2.fastq > fastq/all_2.fastq

#アセンブリを実行
Trinity --seqType fq --left fastq/all_1.fastq --right fastq/all_2.fastq --max_memory 100G --CPU 20 --full_cleanup
  • --full_cleanup   only retain the Trinity fasta file, rename as ${output_dir}.Trinity.fasta
  • --CPU   number of CPUs to use, default: 2
  • --max_memory   suggested max memory to use by Trinity where limiting can be enabled. (jellyfish, sorting, etc)

メモリは100GB指定している(CPU指定が増えると比例してメモリ使用量も増加)。カレントにtrinity_out_dir.Trinity.fastaが出力される。

参考: Trinityがエラーを起こし解決が難しければ、Docker環境で動かすのが楽かもしれません。Dockerならオーバーヘッドが比較的少なくランできます。例えばcentosのコンテナ入れて、そこにpyenv入れてcondaを導入し、conda installでTrinityをインストールする(なるべく手を抜く)。

 

4、マッピング

salmon(紹介)やbowtie2を使う。全てのアライメントをレポートするフラグを立てる必要がある。

#index
bowtie-build trinity_out_dir.Trinity.fasta Trinity

#mapping
for srr in $(cat sra_run_ids.txt)
do
bowtie -p 20 --all -S Trinity -1 fastq/${srr}_P1.fastq -2 fastq/${srr}_P2.fastq |samtools sort -@ 20 -O BAM -o fastq/${srr}.bam -
done
  •  --allをつけ、全てのアライメントを出力するようにしている。fastq/にbamが出力される。salmonを使う例は公式を参考にしてください。実際に動かして確認しましたが、流れは同じです。

 

5、クラスタリング

Corsetにbamの付属情報を与える。今回の6サンプルは、SRR453566-68の3つがbatch growth、SRR453569-71の3つがchmostat 条件になっている。条件が同じか違うかどうかは-gフラグで与える。次のコマンドでは、12は違うグループと見なされる。

cd fastq/

#サンプルごとにクラスタリング。6サンプルを並列実行。子プロセスが終わるのを待つwaitをつけている。
for FILE in `ls *.bam` ; do
corset -r true-stop $FILE &
done
wait

#統合
corset -g 1,1,1,2,2,2 -n A1,A2,A3,B1,B2,B3 -i corset *.corset-reads

clusters.txt、count.txtが出力される。 

head clusters.txt

$ head clusters.txt

TRINITY_DN1053_c0_g1_i1 Cluster-0.0

TRINITY_DN1053_c0_g1_i2 Cluster-0.0

TRINITY_DN2782_c0_g1_i1 Cluster-1.0

TRINITY_DN1713_c0_g9_i1 Cluster-2.0

TRINITY_DN1710_c0_g11_i2 Cluster-3.0

TRINITY_DN1836_c0_g1_i1 Cluster-4.0

TRINITY_DN180_c0_g1_i1 Cluster-5.0

TRINITY_DN947_c0_g2_i1 Cluster-6.0

TRINITY_DN865_c0_g1_i1 Cluster-7.0

1列目が転写産物のID、2列目がアサインされたcluster ID。

 >head count.txt

$ head counts.txt

A1 A2 A3 B1 B2 B3

Cluster-0.0 361 541 397 361 458 562

Cluster-1.0 0 1 2 2 3 8

Cluster-2.0 0 3 3 3 4 4

Cluster-3.0 0 5 1 5 0 0

Cluster-4.0 10 5 4 9 15 25

Cluster-5.0 35 25 19 867 1663 1949

Cluster-6.0 7 6 2 3 3 10

Cluster-7.0 1650 2430 1856 608 702 829

Cluster-8.0 1041 1575 1214 734 849 1132

1列目がクラスターIDで、それ以降の列がそのクラスターにアサインされたリードのカウント。

 

6、edgeRを使ったDEG判定

> R #Rを立ち上げる

library(edgeR)
count_data<-read.delim("counts.txt",row.names=1)
group<-factor(c(1,1,1,2,2,2))
dge<-DGEList(counts=count_data,group=group)

dge<-DGEList(counts=count_data,group=group)
keep <- rowSums(cpm(dge)>1) >= 3
dge <- dge[keep,,keep.lib.sizes=FALSE]

dge<-calcNormFactors(dge)
dge<-estimateCommonDisp(dge)
dge<-estimateTagwiseDisp(dge)

results<-exactTest(dge)
topTags(results)

Comparison of groups:  2-1 

                  logFC    logCPM        PValue           FDR

Cluster-3549.0 7.781056  8.703321  0.000000e+00  0.000000e+00

Cluster-1221.0 6.381804  6.825066  0.000000e+00  0.000000e+00

Cluster-3976.0 4.793833  7.872592  0.000000e+00  0.000000e+00

Cluster-5221.0 4.487102  7.319362 7.386281e-321 9.764664e-318

Cluster-34.2   7.148952 10.719010 5.103542e-317 5.397506e-314

Cluster-2899.0 4.225118  7.564345 2.930408e-310 2.582666e-307

Cluster-1599.0 5.453314  6.612733 2.553171e-309 1.928738e-306

Cluster-4134.0 8.640848  8.879755 4.277488e-309 2.827419e-306

Cluster-502.0  8.604537  9.306601 2.350392e-296 1.380986e-293

Cluster-2771.0 5.183932  7.621110 2.961453e-291 1.566016e-288

続き。

allTags<-topTags(results,n=dim(results$table)[1])1
#ここではFDR0.05で切る。
deClusters=rownames(allTags)[allTags$FDR<0.05]
#保存
write.table(deClusters,"clusters_of_interest.txt",quote=F,row.names=F,col.names=F)

head clusters_of_interest.txt 

$ head clusters_of_interest.txt 

Cluster-3549.0

Cluster-1221.0

Cluster-3976.0

Cluster-5221.0

Cluster-34.2

Cluster-2899.0

Cluster-1599.0

Cluster-4134.0

Cluster-502.0

Cluster-2771.0

 

 

7、FASTAへの変換

用意されているスクリプト(Corset-tools)を使えば、上のcluster IDをfastaに変換できる。

#Corset-toolsをクローンする。
git clone https://github.com/Adamtaranto/Corset-tools
python Corset-tools/fetchClusterSeqs.py -i trinity_out_dir.Trinity.fasta -t clusters_of_interest.txt -o contigs_of_interest.fasta -c clusters.txt

(python2.7環境で実行)

cluster IDに該当するMulti-FASTAが出力される。

f:id:kazumaxneo:20180616093103j:plain

 

8、アノテーション

Multi-FASTAが得られたら、blast2goやTrinnotate(紹介)を使ってアノテーションをかける。

 

 

 

以前一度紹介しましたが、中途半端だったのでまとめ直しました。

引用

Corset: enabling differential gene expression analysis for de novo assembled transcriptomes

Davidson NM, Oshlack A

Genome Biol. 2014 Jul 26;15(7):410.

 

corset-project - Example.wiki

https://code.google.com/archive/p/corset-project/wikis/Example.wiki

 

fastqをクラスタリングする QCluster

 次世代シーケンシング(NGS)技術によって生成されるデータ量は、現在のコンピュータシステムのストレージおよびデータ処理能力に挑戦しているペースで増加している[ref.1]。現在の技術は1回の実行で5千億本以上のDNAを生産し(論文執筆時点)、今後のシーケンサーはこのスループットを向上させることを約束する。シークエンシング技術の急速な進歩により、ゲノムリシークエンシング、RNA-Seq、ChIP-Seqなど多くの異なるシークエンシングベースのアプリケーションが可能になった[ref.2]。このような大きなファイルの取り扱いと処理は、ほとんどのゲノム研究プロジェクトにとって大きな課題の1つになっている。

 シーケンス間の類似性を確立するために、アライメントに基づく方法がかなりの間使用されている[ref.3]。しかし、アライメント方法を適用することができないか、またはそれらが適していない場合がある。例えば全ゲノムの比較は、トラディショナルな従来のアラインメント技術では実施することはできない[ref.4-6]。高速なアライメントヒューリスティックが存在するが、アライメントには通常時間がかかり、次世代シークエンシング技術(NGS)[ref.7,8]によって生成される大規模シーケンスデータには適していないという欠点がある。これらの理由から、多くのアライメントフリー技術が長年提案されている[ref.9]。
 配列を比較するためのアラインメントフリー方法の使用は、異なる用途において有用であることが判明している。いくつかのアライメントフリーの方法は、異なる生物間の進化的関係を研究するためにパターン分布を使用する[ref.4,10,11]。 ChIP-Seqデータ[ref.12-14]およびエントロピープロファイル[ref.15,16]におけるエンハンサー検出のために、いくつかのアライメントフリー方法が考案されている。もう1つのアプリケーションは、遠タンパク質の遠縁関係の分類であり、これは洗練された計数手順で対処することができる[17,18]。 NGSに基づくゲノムのアセンブリフリー比較は、最近になって研究されている[ref.7,8]。包括的なレビューについては、[ref.9 pubmed]を参照してください。
 この研究では、アラインメントフリーリードクラスタリング能力を探求する。クラスタリング技術は、エラー訂正[ref.19]からマイクロRNAの発見[ref.20]まで、多くの異なるアプリケーションで広く使用されている。(一部略)。
Solovyov et  al (ref.21) はNGSをアライメントフリーでクラスタリングする最初の比較方法の1つを提示した。彼らは、k-mer数に基づいて異なる遺伝子および異なる種に由来するリードをクラスタリングすることに焦点を合わせた。Dタイプの測定(第2節参照)、特にD2 *が、同じ遺伝子または種からのリードを効率的に検出し、クラスタリングすることができることを示した(誤差に焦点を当てているref.20とは対照的である)。この論文では、これらの尺度にクオリティ情報を組み込むことによって、この研究を拡張する。

 NGSプラットフォームによって生成されるクオリティスコアは、NGSデータのさまざまな分析のための基本的なもので多くのことに関係している:リードをリファレンスゲノムにマッピングする、エラー訂正[ref.19]。挿入および欠失の検出[ref.23]および他の多くが挙げられる。さらに、将来の世代のシークエンシング技術では、誤った塩基が多数存在する、ロングリードが生成される[ref.24]。 1回のリードあたりの平均エラー数は15%まで増加するため、アライメント・フリーのフレームワークとデノボ・アセンブリで、高品質の情報を活用することが基本となる。

(以下略)

 

インストール

cent os6でテストした。

公式HP

http://www.dei.unipd.it/~ciompin/main/qcluster.html

ZIPをダウンロードして解凍し、makeする。

cd Cluster/
make

./qCluster

$ ./qCluster

Missing/extra input file name

Centroid based (k-means-like) clustering of sequences in n-mer frequency space

 

usage: /Users/kazumaxneo/Downloads/qCluster/qCluster [-h] [-c num_clusters] [-d dist_type] [-e num_method] [-k nmer_length] [-m max_iterations] [-N num_method] [-n] [-P num_method] [-p pseudocount] [-R] [-r] [-S seed] [-t num_trials] [-v] [-w] fastq_file

 

-c num_clusters (5 clusters by default)

-d dist_type:

a: d2* distance

c: chi square statistic

d: d2 distance

e: regular euclidean (L2) distance; default

k: Kullback-Leibler divergence

s: symmetrized Kullback-Leibler divergence

-e num_method: method for computation of quality expected 

  value of words:

1: average quality of the word over the full dataset (default)

2: average quality of each base in the same word over the full 

  dataset and calculate the expected value with Markovian model

3: average quality of each baseover the full dataset and 

  calculate the expected value with Markovian model

-h print this message and exit

-k nmer_length: length of word (2-mers by default)

-m max_iterations: number of maximum iterations of the algorithm 

  without an improvement (min. 2) 

-N num_method: Divide each quality vector by the taxicab norm so 

  that the sum of the elements is (about) 1. Possible values:

0: no normalization

1: divide by the norm of the vector itself

2: divide by the norm of the frequency vector (default)

3: as method 2 and project the vector on the hyperplane 1*X-1=0

-n normalize frequency matrix to make each column univariant; 

  implies L2 distance

-P num_method: method for computation of expected frequancy of words:

1: uses a markovian model to compute the expected frequency of

  words basing only on the words of the cluster

  (P2Local methond: default)

2: Average frequency of every word over the entire dataset

  (P1Global method)

3: Markovian model based on average frequency of single bases

  on the full dataset (P2Global method)

-p pseudocount: (default: 1); can be fractional. 

  Must be nonzero with KL and simmetrized KL distance.

-R do not redistribute missing quality among other bases. 

-r reverse complement and stack together

-S seed: initial seed for random number generator

-t num_trials: repeat clustering num_trials, choosing the best 

  partitioning. Default: 1

-v output progress messages (repeat for increased verbosity)

-w write sequences from each cluster to a file

 

 

ラン

テストラン。3つのクラスターに分類する。

QCluster -d a -c 3 -k 3 -t 5 -S 0 -w Example/sequences.fastq 
  • -c   num_clusters (5 clusters by default)
  • -d   dist_type:-d dist_type:
  1. a : d2* distance
  2. c : chi square statistic
  3. d : d2 distance
  4. e : regular euclidean (L2) distance; default
  5. k : Kullback-Leibler divergence
  6. s : symmetrized Kullback-Leibler divergence 
  • -t   num_trials: repeat clustering num_trials, choosing the best -t num_trials: repeat clustering num_trials, choosing the best    partitioning. Default: 1

3つのfastqが出力される。

 

 

引用

QCluster: Extending Alignment-Free Measures with Quality Values for Reads Clustering

Matteo CominAndrea LeoniMichele Schimd

WABI 2014: Algorithms in Bioinformatics pp 1-13

 

SPAdesアセンブラ

 ref.1

 人体から海洋までほとんどの環境のバクテリアは研究所でクローン化できないため、既存のNGS(Next Generation Sequencing)技術を使用してシーケンスを決定することはできない。これは、Human Microbiome Project(HMP)(論文より Gill et al、2006)から抗生物質発見(Li and Vederas、2009)までのさまざまなプロジェクトの主要なボトルネックを表している。例えば、HMPの重要な問題は、バクテリアが互いにどのように相互作用するかである。これらの相互作用は、しばしば、他のバクテリアとのコミュニケーションまたはそれらを殺すために産生される様々なペプチドによって行われる。しかしながら、質量分析(このような研究のための重要な技術)は、かなり完全なプロテオームの知識を必要とするため、ヒトマイクロバイオームのpeptidomics研究は今や限られている。一方、新しいペプチド抗生物質の研究は、Non-Ribosomal Peptide Syntetases(NRPS)をコードする遺伝子のDNA配列決定(Sieber and Marahiel、2005)により大きく恩恵を受けるが、既存のメタゲノミクスアプローチは、通常このような例外的に長い遺伝子(60,000ヌクレオチド)をアセンブリすることができない。

 HMPと新しい抗生物質の発見は、Single-Cell Sequencing(SCS)によって革命を起こす多くのプロジェクトの2つの例にすぎない。 SCSの実験的(Ishoey et al、2008; Navin et al、2011; Islam et al、2011)および計算(Chitsaz et al、2011)の両局面における最近の改良により、1細胞からバクテリアゲノムを決定する可能性が開かれた。特に、Chitsaz et al(2011)は、SCSで生物の代謝を推定するのに十分な多数の遺伝子を捕捉できることを示した。多くのアプリケーション(プロテオミクスや抗生物質の発見を含む)では、大多数の遺伝子を捕捉することは、完全なゲノムと同じくらい有用である。

 現在、MDA(Multiple Displacement Amplification)はシングルセル増幅の主要技術である(Dean et al、2001)。しかし、MDAは極端な増幅バイアス(異なる領域間でカバレッジの差が桁違いに大きい)をもたらし、キメラリードとキメラペアリードが生じ、これが次のアセンブリを複雑にする。既存のアセンブラはこれらの問題に対処できないことから、Rodrigue et al. (2009) は、SCSが直面している課題は実験的ではなく計算的なものになっていると指摘した。最近の論文(Marcy et al、2007; Woyke et al、2010; Youssef et al、2011; Blainey et al、2011; Grindberg et al、2011)は、SCSにおけるフラグメントアセンブリの課題を示している。

 Chitsaz et al (2011)はE+V-SC assemblerを導入し、modified EULER-SR とmodified velvet を組み合わせ、SCSデータのフラグメントアセンブリの大幅な改善を達成した。しかし、Chitsaz et alの共著者である我々(この論文の著者ら)(G.T.およびP.A.P.)は、 (2011)では、SCSの可能性を最大限に活用するために、(既存のツールを変更するだけでなく)アルゴリズム設計を変更する必要があることに気付いた。

 著者らは、SPAdesアセンブラを発表し、SCSと標準(マルチセル)バクテリアの両方のデータセットのために、数多くの新しいアルゴリズムソリューションを導入し、最先端のアセンブラを改良した。フラグメントアセンブリは、そのk-mersのセットから文字列を再構成する問題として、しばしば抽象化される。この抽象化は、多くのフラグメントアセンブリアルゴリズムの基礎となるアセンブリへのde Bruijnアプローチに自然につながる。しかしながら、NGSデータのより高度な抽象化は、文字列中の距離≒dのk-mersのペア(以下、k-bimersと呼ぶ)から文字列を再構築する問題を考慮する。残念ながら、前者の抽象化には単純なアルゴリズムがある一方で、後者の抽象化の分析は、主にde Bruijn graphの後処理ヒューリスティックになっている。ペアリード(以下、bireads2と呼ぶ)の解析のための多くのヒューリスティックが提案されているが(Pevzner et al、2001; Pevzner and Tang、2001; Zerbino and Birney、2008; Butler et al、2008; Simpson et al、2009 ; Chaisson et al、2009; Li et al、2010)、bireads2の利用はアセンブリのステージの中でおそらく最も不十分な段階であると考えられる。 Medvedev et al(2011a)は、後者の抽象化に適した新しい手法であり、標準的なDe Bruijnグラフより重要な利点を持つPaired de Bruijnグラフ(PDBG)を導入した。しかし、PDBGは理論的には実用的ではなくむしろ理論的に導入されており、主にリード間を固定距離とした非現実的なケースを対象としていた。

 IduryとWaterman(1995)がフラグメントアセンブリのためにde Bruijn graphを導入したとき、多くの人々はサンガーのリードのエラー率が高いため、このアプローチを非実用的と見なしていた(ref.3)。Pevzner et al (2001)は、このボトルネックを解消するためにエラー訂正手順を導入した。同様に、PDBGsは、NGSに特徴的なbiread(and thus k-bimer)距離の変動のために実用的でないように見えるかもしれない。 SPAdesは、k-bimersの正確な距離を明らかにするk-bimer調整を導入し、PDBGに触発されたペアのアセンブリグラフを導入することによって、このボトルネックに対処する。特に、SPAdesはペアリードを利用することができる。 E + V-SCは、リードを使用したが、キメラリードペアのレベルの上昇によって生じるミスアセンブリを回避するために、ペアエンドのpairingを無視した。

 読者はPevzner et al(2004 pubmed)で紹介されたA-Bruijn graphの概念に精通していると仮定する 。この論文のDe Bruijn graph、PDBG、およびその他のgraphは、A-Bruijn graphの特別なケースである。 SPAdesは、最初のde Bruijnグラフを作成するためだけにk-mersを使用し、後でそれらについて「忘れる」という意味で普遍的なA-Bruijnアセンブラである。それ以降のステージでは、k-mersによってラベル付けする必要がないgraphに対してのみgraph理論操作を実行する。操作はgraphトポロジーカバレッジ、およびシーケンス長に基づいているが、シーケンスそのものには基づいていない。最後の段階で、コンセンサスDNA配列が復元される。著者らは、同じフレームワーク内でA-Bruijn graphのいくつかのバリエーション(例えば、Bruijn graphの対と多重化)を実装し、これらのgraphが有用であることが証明されている他のアプリケーションに適用する汎用アセンブラを設計した(Bandeira et al、 2007、2008; Pham and Pevzner、2010)。

 セクション2では、SPAdesのステージの概要を説明する。セクション3では、de Bruijn graph、Multisized de Bruijn graph、およびPDBGを定義する。セクション4-6は、SPAdesのさまざまな段階をカバーしている。セクション7では、シングルセルおよび培養された大腸菌(E.coli)データセットについて、SPAdesおよび他のアセンブラベンチマークする。セクション8では、アセンブリgraphの構築とデータ構造に関する追加の詳細を示す。第9章では、ペアアセンブリgraphを構築する詳細な例を示す。セクション10では、ユニバーサルアセンブラの概念についてさらに説明する。

セクション2

以下では、SCSで特に面倒な問題を扱うSPAdesの4つの段階について説明する。

(1)ステージ1(アセンブリgraph構築)はすべてのNGSアセンブラによって扱われ、しばしばde Bruijn graph簡略化(EULER / Velvetにおけるバルジ/バブル除去)と呼ばれる。この論文では、複数のk-merを使ったde Bruijn graphを使用し、新しいバルジ/チップ除去アルゴリズムを実装し、キメラを検出して除去し、biread情報を距離ヒストグラムに集約し、実行されたgraph操作をバックトラックすることを可能にするアセンブリグラフ構築の新しいアプローチを提案する。
(2)ステージ2(k-bimer調整)は、アセンブリgraphの距離ヒストグラムとパスのジョイント解析を使用して、ゲノム内のk-mer間の正確な距離推定値(アセンブリgraph内のエッジ)を導き出す。
(3)ステージ3は、PDBGアプローチからインスパイアされたペアのアセンブリgraphを構築する。
(4)ステージ4(コンティグ構築)は、サンガーシーケンス時代に十分に研究された(Ewing et al、1998)。 NGSプロジェクトは通常、高いカバレッジを特徴とするため、NGSアセンブラはかなり正確なコンティグを生成する(ただしSCSの精度は低下する)。 SPAdesは、コンティグのDNA配列を構築し、グラフの簡略化(論文 8.6節参照)を目的としてコンティグへリードをマッピングする。
これまでの研究では、さまざまなアセンブラをエラー訂正ツールと組み合わせることで性能が向上することが実証されている(Pevzner et al、2010; Kelley et al、2010; Ilie et al、2010; Gnerre et al、2011)。しかし、 (Quake [Kelley et al、2010]などの)ツールは、暗黙的にほぼ均一なカバレッジを仮定しているため、シングルセルのデータでは機能しない。 Chitsaz et al。(2011)はVelvet-SCとEULERのエラー訂正(Chaisson and Pevzner、2008)を結びつけ、E + V-SCツールとなった。この論文では、SPAdesは、アセンブリ前にエラー修正と品質トリミングのためにHammer(Medvedev et al、2011b)(SCSを対象としています)によるエラー修正を行う。

 

オリジナルのSPAdesは以下のようなフローでアセンブリを実行する。

f:id:kazumaxneo:20180614103714j:plain

プレゼンテーションPDFより転載。 de Brujin graph構築にはショートリード情報が用いられ、ロングリード情報はエラーの多寡に関わらず使われない。

 

公式ページ

http://cab.spbu.ru/software/spades/

3.12マニュアル

http://cab.spbu.ru/files/release3.12.0/manual.html 

 

New Frontiers of Genome Assembly with SPAdes 3.1

質疑応答がちょっと面白いです。

ppt1

f:id:kazumaxneo:20180613212251j:plain

ppt2

f:id:kazumaxneo:20180613212329j:plain

 

インストール

mac os 10.13でテストした。

 バイナリをダウンロードする。Linuxは上のマニュアル参照(バイナリが用意されている)。

curl http://cab.spbu.ru/files/release3.12.0/SPAdes-3.12.0-Darwin.tar.gz -o SPAdes-3.12.0-Darwin.tar.gz
tar -zxf SPAdes-3.12.0-Darwin.tar.gz
cd SPAdes-3.12.0-Darwin/bin/

./spades.py 

$ ./spades.py 

SPAdes genome assembler v3.12.0

 

Usage: ./spades.py [options] -o <output_dir>

 

Basic options:

-o <output_dir> directory to store all the resulting files (required)

--sc this flag is required for MDA (single-cell) data

--meta this flag is required for metagenomic sample data

--rna this flag is required for RNA-Seq data 

--plasmid runs plasmidSPAdes pipeline for plasmid detection 

--iontorrent this flag is required for IonTorrent data

--test runs SPAdes on toy dataset

-h/--help prints this usage message

-v/--version prints version

 

Input data:

--12 <filename> file with interlaced forward and reverse paired-end reads

-1 <filename> file with forward paired-end reads

-2 <filename> file with reverse paired-end reads

-s <filename> file with unpaired reads

--merged <filename> file with merged forward and reverse paired-end reads

--pe<#>-12 <filename> file with interlaced reads for paired-end library number <#> (<#> = 1,2,...,9)

--pe<#>-1 <filename> file with forward reads for paired-end library number <#> (<#> = 1,2,...,9)

--pe<#>-2 <filename> file with reverse reads for paired-end library number <#> (<#> = 1,2,...,9)

--pe<#>-s <filename> file with unpaired reads for paired-end library number <#> (<#> = 1,2,...,9)

--pe<#>-m <filename> file with merged reads for paired-end library number <#> (<#> = 1,2,...,9)

--pe<#>-<or> orientation of reads for paired-end library number <#> (<#> = 1,2,...,9; <or> = fr, rf, ff)

--s<#> <filename> file with unpaired reads for single reads library number <#> (<#> = 1,2,...,9)

--mp<#>-12 <filename> file with interlaced reads for mate-pair library number <#> (<#> = 1,2,..,9)

--mp<#>-1 <filename> file with forward reads for mate-pair library number <#> (<#> = 1,2,..,9)

--mp<#>-2 <filename> file with reverse reads for mate-pair library number <#> (<#> = 1,2,..,9)

--mp<#>-s <filename> file with unpaired reads for mate-pair library number <#> (<#> = 1,2,..,9)

--mp<#>-<or> orientation of reads for mate-pair library number <#> (<#> = 1,2,..,9; <or> = fr, rf, ff)

--hqmp<#>-12 <filename> file with interlaced reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)

--hqmp<#>-1 <filename> file with forward reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)

--hqmp<#>-2 <filename> file with reverse reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)

--hqmp<#>-s <filename> file with unpaired reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9)

--hqmp<#>-<or> orientation of reads for high-quality mate-pair library number <#> (<#> = 1,2,..,9; <or> = fr, rf, ff)

--nxmate<#>-1 <filename> file with forward reads for Lucigen NxMate library number <#> (<#> = 1,2,..,9)

--nxmate<#>-2 <filename> file with reverse reads for Lucigen NxMate library number <#> (<#> = 1,2,..,9)

--sanger <filename> file with Sanger reads

--pacbio <filename> file with PacBio reads

--nanopore <filename> file with Nanopore reads

--tslr <filename> file with TSLR-contigs

--trusted-contigs <filename> file with trusted contigs

--untrusted-contigs <filename> file with untrusted contigs

 

Pipeline options:

--only-error-correction runs only read error correction (without assembling)

--only-assembler runs only assembling (without read error correction)

--careful tries to reduce number of mismatches and short indels

--continue continue run from the last available check-point

--restart-from <cp> restart run with updated options and from the specified check-point ('ec', 'as', 'k<int>', 'mc', 'last')

--disable-gzip-output forces error correction not to compress the corrected reads

--disable-rr disables repeat resolution stage of assembling

 

Advanced options:

--dataset <filename> file with dataset description in YAML format

-t/--threads <int> number of threads

[default: 16]

-m/--memory <int> RAM limit for SPAdes in Gb (terminates if exceeded)

[default: 250]

--tmp-dir <dirname> directory for temporary files

[default: <output_dir>/tmp]

-k <int,int,...> comma-separated list of k-mer sizes (must be odd and

less than 128) [default: 'auto']

--cov-cutoff <float> coverage cutoff value (a positive float number, or 'auto', or 'off') [default: 'off']

--phred-offset <33 or 64> PHRED quality offset in the input reads (33 or 64)

[default: auto-detect]

ほかにmetaspades.py、rnaspades.py、dipspades.py、plasmidspades.py、/truspades.pyがある。以下のような使い分けになる。

  • metaSPAdes – a pipeline for metagenomic data sets (see metaSPAdes options).
  • plasmidSPAdes – a pipeline for extracting and assembling plasmids from WGS data sets (see plasmidSPAdes options).
  • rnaSPAdes – a de novo transcriptome assembler from RNA-Seq data (see rnaSPAdes manual).
  • truSPAdes – a module for TruSeq barcode assembly (see truSPAdes manual).
  • dipSPAdes – a module for assembling highly polymorphic diploid genomes (deprecated, see dipSPAdes manual).

spadesはbrewでも導入できる。

brew install spades

SPAdesは現在でも開発が続けられており、バージョンを変えると末端部位やカバレッジの低い領域の結果が揺れ動くことがあります。アップグレードする時は注意してください。

 

 

ラン

ペアエンドリードfastqのアセンブリ。非圧縮fastqにも対応。

spades.py -1 pair1.fq.gz -2 pair2.fq.gz -k auto --careful -t 12 -o output
  • -1     file with forward paired-end reads-1 <filename> file with forward paired-end reads
  • -2     file with reverse paired-end reads
  • -t     number of threads-t/--threads <int> number of threads [default: 16]
  • -k     comma-separated list of k-mer sizes (must be odd and-k <int,int,...> comma-separated list of k-mer sizes (must be odd and less than 128) [default: 'auto']
  • --careful     tries to reduce number of mismatches and short indels--careful tries to reduce number of mismatches and short indels
  • -o   directory to store all the resulting files (required)-o <output_dir> directory to store all the resulting files (required)

インターレースfastqは--12、シングルエンドfastqは-sで指定する。シングルセルシーケンスなどで複数fastqがある場合は、数値をつけて指定する。例えばインターレースfastqを複数指定するなら"--pe1-12 sample1.fq --pe2-12 sample2.fq --pe3-12 sample3.fq ..."。詳細はマニュアルの Specifying multiple librariesの節を参照。

 

ペアエンドリードfastqとONTのロングリードのハイブリッドアセンブリ

spades.py -1 pair1.fq.gz -2 pair2.fq.gz --nanopore long.fa -k auto --careful -t 12 -o output

 

ランが終わるとk-merごとにアセンブリデータフォルダができる。scaffolds.fastaが最終出力となる。correctedにはエラー訂正されたfq.gzも保存されている。

f:id:kazumaxneo:20180613230424j:plain

 

ケースバイケースですが、うまくパラメータを使えばパフォーマンスが上がる可能性があります。気づいたことを書いておきます。

  1. --carefulをつけると、付属のbwaを使い、アセンブリ後にミスマッチを訂正するステップが挿入される。フラグを立てると、マッピングの時間分だけランタイムが増えるが、スモールゲノムではフラグをつけることが推奨されている(defaultではオフ)。ラージゲノムでは非推奨。純粋な培養株のシーケンスを想定して修復するため、rnaSPADdesとmetaSPAdesにはこのオプションはない。
  2. たとえハイクオリティなPacbioのロングリードでも、ロングリードだけでSPAdesを動かすことは非推奨になってます。
  3. --only-error-correctionをつけると、BayesHammerによるエラー訂正だけが実行されます。エラー訂正されたfastqが欲しいときに使えます。
  4. --only-assemblerをつけると、エラー訂正なしでアセンブリが実行されます。非推奨ですが、ほかのエラー訂正ツールで修復したfastqを使う時にも利用できます。
  5. カビくらい大きなゲノムになるとかなりメモリを食います。論文を読んで、ランできそうか確認してください。defaultでは-m --memoryは最大250GBで、越えるとランを停止します。
  6. 軽度のコンタミネーションは、--cov-cutoff <NUM>でカバレッジcut-offをつけることで抑制できます。アセンブリ後にカバレッジGC、taxonomy profileを分析して、異種生物ゲノムのcontigのコンタミカバレッジで区別できそうなら使用して下さい。

 

 

http://quast.bioinf.spbau.ru

 

  • Unicycler

 

引用

ref.1

SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing

Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev,Pavel A. Pevzner, 

J Comput Biol. 2012 May; 19(5): 455–477.

 

ref.2

ExSPAnder: a universal repeat resolver for DNA fragment assembly

Prjibelski AD, Vasilinetc I, Bankevich A, Gurevich A, Krivosheeva T, Nurk S, Pham S, Korobeynikov A, Lapidus A, Pevzner PA.

Bioinformatics. 2014 Jun 15;30(12):i293-301

 

ref.3

Assembling short reads from jumping libraries with large insert sizes

Vasilinetc I, Prjibelski AD, Gurevich A, Korobeynikov A, Pevzner PA.

Bioinformatics. 05 Oct 5;(0):6-8.

 

ref.4

hybridSPAdes: an algorithm for hybrid assembly of short and long reads

Antipov D, Korobeynikov A, McLean JS, Pevzner PA

Bioinformatics. 2016 Apr 1;32(7):1009-15.

 

ref.5

TruSPAdes: barcode assembly of TruSeq synthetic long reads

Bankevich A, Pevzner PA

Nat Methods. 2016 Mar;13(3):248-50.

 

ref.6

metaSPAdes: a new versatile metagenomic assembler

Nurk S, Meleshko D, Korobeynikov A, Pevzner PA

Genome Res. 2017 May;27(5):824-834.

 

ref.7

BayesHammer: Bayesian clustering for error correction in single-cell sequencing

Nikolenko SI, Korobeynikov AI, Alekseyev MA

BMC Genomics. 2013;14 Suppl 1:S7.

 

 

参考

バクテリアのゲノムアセンブリ

Tychus: a whole genome sequencing pipeline for as- sembly, annotation and phylogenetics of bacterial genomes

https://www.biorxiv.org/content/biorxiv/early/2018/03/16/283101.full.pdf

 

アセンブリの進め方。

http://compgenomics2017.biology.gatech.edu/images/a/ad/Genome_Assembly_Final_Pipeline_and_Results.pdf

 

バクテリアのRNA seq定量ツール EDGE-pro

 

 バクテリアゲノム中の遺伝子の発現を測定することは、感染の治療法の開発から合成ゲノムの作成まで、非常に幅広い用途を有する。バクテリアにおける遺伝子発現研究は、代謝経路を研究し、変異株の特性を同定し、他の点ではそれらのゲノムにおける分子過程をよりよく理解するために用いられてきた(論文より ref1.2)。

 10年以上にわたって、マイクロアレイは遺伝子発現を研究するための主な手段であった。しかしながら、マイクロアレイ技術は、プローブが設計された転写産物のみを捕捉するので、したがって、よく研究された種の既知の株における既知の遺伝子へ適用を制限する。あるいは、定量PCR(qPCR)の使用は、ゲノム中のすべての遺伝子ではなく、特定の遺伝子を定量することを可能にするが、この技術はゲノム全体のスケールではるかに高価である。第2世代シーケンシングの品質、効率、およびコストの近年の改良によって、マイクロアレイ解析を置き換え、RNAを直接捕捉してシーケンスするいわゆるRNA-seq実験の爆発がもたらされた。マイクロアレイ実験では、参照株と新規の株との違いによって、マイクロアレイ上のいくつかのプローブとのハイブリダイゼーションが妨げられる可能性がある。対照的に、RNA-seq実験では、転写された遺伝子が配列決定され、ゲノムにアライメントされる。アラインメントアルゴリズムは、多くのミスマッチを許容することができ、それによって標的ゲノムが参照から分岐した場合でも遺伝子発現の高感度測定が可能になる。さらに、全転写物が配列決定されるので、RNA-seqデータはゲノムのオペロン構造も明らかにする。

 RNA-seq技術の導入以来、サンプル中の遺伝子発現を定量化し、複数の試料間の遺伝子発現の差異を見出すためにソフトウェア方法が開発されている(ref.5-9)。しかし、遺伝子発現を推定するためのツールは、真核生物遺伝子の構造を同定する目的で開発された。これらのツールは、遺伝子内のイントロン領域を見いだし、高等真核生物で一般的な選択的スプライスバリアントを見出す努力の多くに焦点を当てている。逆に、バクテリア遺伝子はイントロンを有さず、選択的にスプライシングされない。したがって、それらの転写産物を分析する際にスプライスバリアントを探す必要はない。

 バクテリアゲノムは非常に密集しており、その多くは互いに重なっている。これまでのRNA-seqソフトウェア法は、一般に、ヒトや他の哺乳動物(ほとんどのRNA-seq実験の主な標的)においては重複が極めてまれであるため、重複する遺伝子を扱う手段を提供していない。これとは対照的に、典型的なバクテリアゲノムの約90%はタンパク質をコードしている(ref.10)。広い進化範囲にまたがる220の原核生物ゲノムの研究(ref.10)から、これらの種の全遺伝子の29%は5 'または3'で互いにオーバーラップがあることが示されている 。これらのオーバーラップは、わずか数塩基対(bp)から100bpをはるかに超える範囲にわたる。重複する遺伝子は、同じ鎖上または反対の鎖上に存在し得る; 従って、ストランド特異的RNA-seqは、せいぜい部分的な解を提供するだけである。オーバーラップ領域内にマップされたRNAシーケンスリードの場合、2つの遺伝子のどちらがリードをもたらしたかを決定することが不可能であり、したがって、各原核生物遺伝子の遺伝子発現を決定する際に課題が生じる。

 さらに、分析のための要求を複雑にするバクテリアRNA-seq技術は、ポリアデニル化不在のために、真核生物RNA seqプロトコルと少なくとも1つの大きな相違点を有することである。真核生物の転写産物上の長いポリAテイルは捕捉するプローブとして使用できるが、バクテリアRNA-seqではランダムプライミング法が必要となる。もう一つの課題は、捕捉された細菌転写産物の80% がrRNAになることである。 rRNAを除去するための方法が開発されているが、RNA-seq実験ではまだ多くのrRNAが依然として出現する。

 真核生物ゲノムとバクテリアゲノムとの間の相違、およびRNA-seqプロトコル間の相違のために、発現解析のための既存のプログラムは、バクテリアゲノムからRNA-seqデータに適用されると、ほとんど機能しないか完全にbreak down する。従って、バクテリアRNA-seqデータにおける遺伝子発現のレベルを推定するためには、新しいバイオインフォマティクス法が必要である。現在のところ(論文執筆時点)、この目的のための独立型ツールは存在しない。複数のバクテリアRNA-seqプロジェクトが公開されているが、これらは遺伝子発現を定量化するためにad hoc(その場限り)な方法を用いており、Bowtie、MAQ、SOAP、BWA、ELAND、Novoalign、などの次世代配列アライナーを使用して、最初に入力ゲノムにリードをアライメントする。これらのアライメントから、各遺伝子にマッピングされたリード数をカウントし、通常は各遺伝子の長さでカウントを正規化する。これらのad hocなアプローチのいくつかのレビューは、Guell et al (ref.27) やVan Vliet.(ref.28)に見つかる。

 標準的なアラインメント手法が直面する課題の1つは、すべてのデータセットにおいて、いくつかのリードがゲノムの複数の場所にアライメントされうることである。これらの複数の場所のどれが真のシーケンス源であるかを決定することは困難であり、時には不可能である場合がある。特に、リピートがある場合はそうである。この問題を回避するために、いくつかの以前の方法では、マルチアラインされたリードを単に破棄する。この戦略は、リピートを含む遺伝子の見かけの発現レベルを有意に(そして間違って)減少させる可能性がある。より深刻な問題は、ファミリー内のどのリードもファミリー遺伝子のすべてのコピーに等しくマッピングされるような遺伝子ファミリーで生じる。マルチアライメントリードをカウントする他の方法では、リードがマップされている各場所に分数を割り当てたり、リードをマルチマップされた場所のいずれかにランダムにアサインしたりする。これらのどれも完全な解決法を提供していないが、我々(この論文の著者ら)が示すように、分数のリードカウントの使用は実際にはうまく機能する。

 バクテリアの遺伝子発現レベルを推定するために現在使用されている臨時の方法は、合理的にうまくいくが、しばしば使用するのは容易ではない。一部のユーザーは、複数のソフトウェアツールを連続して実行する必要があり、あるプログラムの出力は、次のツールの入力が正しくない場合があり、データを再フォーマットするために付随するプログラムが必要になる。本論文では、原核生物ゲノムの遺伝子発現レベルを推定するために特別に設計された最初のスタンドアローン方式であるEDGE-pro(Estimated Degree of Gene Expression in PROkaryots)を紹介する。 EDGE-proは、上記の課題に対する解決策を提供する効率的なソフトウェアシステムである。

 EDGE-proパイプラインは、リファレンスゲノム、そのゲノム内のタンパク質コード遺伝子の座標を含むタンパク質翻訳テーブル(ptt)、tRNAおよびrRNA遺伝子の座標を含む別のテーブル(rnt)、およびRNA-seq 自分自身を読み込む。 pttとrntのテーブルファイルが利用できない場合は、Glimmerを実行してタンパク質コード遺伝子を検索し、tRNAscan-SEとRNAmmerをRNA遺伝子に変換することで、ゲノム配列から別々に生成することもできる。EDGE-proパイプラインは、以下の4つの手順で構成されている。

  1. メインマッピングステップ。すべてのリードがリファレンスゲノムにアライメントされる。
  2. マルチアライメントリードのフィルタリングステップ。
  3. リファレンスゲノムにおける各塩基のカバレッジデプスの計算
  4. カバレッジデプスを、各遺伝子のRPKM(マッピングされた百万単位当たりの遺伝子キロベース当たりのリード量)に変換する。

各ステップの詳細は論文で確認してください(リンク)。このツールはedgeRと直接関係があるわけではありません。このツールは、内部でbowtieを動かしマッピングして、フィルタリングした後にリードカウントを出すプロセスを自動で行うツールです。

 

公式ページ

https://ccb.jhu.edu/software/EDGE-pro/index.shtml

マニュアル

https://ccb.jhu.edu/software/EDGE-pro/MANUAL

 

インストール

cent os6でテストした。

依存

  • Bowtie2

本体中に、linux 64bit向けのコンパイル済みのBowtie2があるので、Bowtie2にパスが通ってなければそれを使ってください。

 

EDGE-pro本体のダウンロード (linux向けバイナリ リンク)。

wget https://ccb.jhu.edu/software/EDGE-pro/EDGE_pro_v1.3.1.tar.gz
tar -vxf EDGE_pro_v1.3.1.tar.gz

 

perl EDGE_pro_v1.3.1/edge.pl -h

Usage

    edge.pl <-g genome> <-p ptt> <-r rnt> <-u reads>

# [OMP_NUM_THREADS=n] PATH/edge.pl <-g genome> <-p ptt> <-r rnt> <-u reads> [-v reads2, if paired-end] [-m minInsert] [-M maxInsert] [-t bowtie threads] [-s path to Source code] [-o prefix of Output files] [-w window] [-i initial transcription site window] [-x similarity] [-l read length] [-c minimum coverage] [-n count type: 0-partial count (default), 1-random count] [-h]

(anaconda3-5.0.0) [uesaka@cyano EDGE_pro_v1.3.1]$ perl edge.pl -h

USAGE:

------

[OMP_NUM_THREADS=n] PATH/edge.perl <-g genome> <-p ptt> <-r rnt> <-u reads> [options]

 

MANDATORY FILES:

----------------

-g genome: fasta file containing bacterial genome. If multiple chromosomes/plasmids exist, they must be combined into one file before running EDGE (see MANUAL)

-p ptt: ptt file with coordinates of coding genes, in Genbank format. If multiple chromosomes/plasmids exist, they must be combined into one file (see MANUAL).

-r rnt: rnt file with coordinates of rRNAs and tRNAs, in Genbank format. If multiple chromosomes/plasmids exist, they must be combined into one file (see MANUAL).

-u reads: fastq file of reads. If multiple fastq files exists, see MANUAL to combine them, if possible.

 

OPTIONAL FILES/PARAMETERS:

--------------------------

-v reads2: fastq file of mates in paired-end data. If this file is not entered, single-end reads are assumed.

-m min: min is an integer value. It is minimum insert size in paired-end library. Default: 0.

-M max: max is an integer value. It is maximum insert size in paired-end library. Default: 500.

-t threads: threads is an integer value. It is the number of threads to be used by Bowtie2. Default: 1.

OMP_NUM_THREADS is an integer environmental optional parameter that specifies the number of threads to be used to count per base coverage. Note that it is entered before the command ./edge.perl. Default: 16.

-s source_dir: It is a string specifying the absolute of relative path to the directory that contains all scripts. Default: working directory.

-o prefix: It is a string specifying the prefix of all output files. Default: out.

-w window: It is an integer specifying the window size close to overlapping region used to approximate the coverage of a gene close to the overlapping region in order to distrbute the coverage of the overlapping region between two overlapping genes. Default: 100.

-i untranslated region: It is an integer specifying the window size of the untranslated region bewteen the initial transcription site and the start codon. Default: 40.

-x similarity: It is a decimal number spcifying the percentage used to determine when two coverage values are considered similar. For example, if the similarity is x, and coverage of a region is C, then another region is considered similarly expressed if its coverage is in the interval [(1-x)*C,(1+x)*C]. Default: 0.15.

-l read length: It is an integer specifying the read length. If read length is not specofied, the first 1000 reads are used to approximate the read length.

-c minimum coverage: It is an integer specifying the minimum average coverage of gene for gene to be considered expressed. Coverage lower than specified is assumed to be noise and gene is considered to not be expressed. Default: 3.

-n count type: It is 0 or 1 specifying how to count reads that map to multiple places. 0 denotes giving a partial count to each place where the read maps. 1 denotes picking randomly one of the places where the read maps and assigning full count to that one place. Default: 0.

-h: display this help

 

 

 

ラン

 Role of CRP K100 positive charge on Escherichia coli global transcriptomeのWTのtripicates(SRR5416993 - SRR5416995)と K100 variant (SRR5416999 - SRR5417001)を使い、2群間比較を行う。SRAにDepositされたデータは3グループあるので、3グループの多群間比較も行える。目的により変わるが、ここではクオリティトリミング後もシーケンスデータは十分量あるものと仮定し、解析は盲目的に最後まで実行する(検証ステップはOptional stepに記載)。

 

 1、シーケンスデータのダウンロード。

以下のURLのSRAデータを使用する。E.coliのWTと変異株のRNA seq解析データとなる。

https://www.ncbi.nlm.nih.gov/bioproject/PRJNA381697

#ヒアドキュメントでダウンロードリスト作成
cat >sra_run_ids.txt <<EOF
SRR5416993
SRR5416994
SRR5416995
SRR5416999
SRR5417000
SRR5417001
EOF

#ダウンロード(srapyを使う(Github))
get-run.py -s -f sra_run_ids.txt -d sra_files/

#またはsra-toolkitを使う
prefetch --option-file sra_run_ids.txt
mv ~/ncbi/public/sra/SRR541* sra_files/

 

2、ペアエンドのfastqに変換。

forループで処理。pfastq-dump(紹介)を使えばより早く終わります。

mkdir fastq

for srr in $(cat sra_run_ids.txt)
do
fastq-dump --split-files --skip-technical --readids \
--defline-seq '@$sn/$ri' --defline-qual '+' sra_files/${srr}.sra \
-O fastq
done

 

3、クオリティフィルタリング。

高速なfastp(紹介)を使い、アダプターとクオリティトリミングしてfastq.gz出力。forループで処理。

for srr in $(cat sra_run_ids.txt)
do
fastp -i fastq/${srr}_1.fastq -I fastq/${srr}_2.fastq \
-o fastq/${srr}_clean1.fastq.gz -O fastq/${srr}_clean2.fastq.gz \
-h fastq/${srr}_report.html -j fastq/${srr}_report.json
done

htmlとjsonのレポートもそれぞれ出している。 統合した1つのレポート出したいなら、afterQC(紹介)などMultiQC(紹介)が対応したツールに切り替えて、MulitiQCに渡すのが楽ですfastp実行時に"-w 8"フラグ立てるとエラーが起きたのでdefaultの"-w 3"で実行した。

 

4、コンカテネート。

EDGE-proはペアエンド情報を使わないので、EDGE-pro入力前にペアエンドをコンカテネートする(間違いが発生しないよう丁寧に書いたが、今回のケースではfastq-dumpのコマンドを--split-filesを付けずに実行してOK)。EDGE-proはfq.gzを扱えないのでgz解凍も行う。

for srr in $(cat sra_run_ids.txt)
do
cat fastq/${srr}_clean2.fastq.gz fastq/${srr}_clean2.fastq.gz > fastq/${srr}_merge.fq.gz
sleep 1s
gzip -dv fastq/${srr}_merge.fq.gz
done

論文より  However, because bacterial RNA-seq analysis does not need to link together exons across splice junctions, paired-end sequencing does not provide as much of an advantage.

 

  

Optional1:  BBsketch(紹介)でシーケンスデータを1つ分析。

sendsketch.sh in=SRR5416995.fastq.gz reads=100k

sendsketch.sh in=SRR5416995.fastq.gz reads=100k

Adding /home/uesaka/.linuxbrew/Cellar/bbtools/37.77/resources/blacklist_refseq_species_300.sketch to blacklist.

Loaded 1 sketch in 0.738 seconds.

 

Query: 1/1 DB: RefSeq SketchLen: 10032 Seqs: 200000 Bases: 27673848 gSize: 3391752 Quality: 0.9365 File: SRR5416995.fastq.gz

WKID KID ANI Complt Contam Matches Unique noHit TaxID gSize gSeqs taxName

47.53% 36.64% 97.51% 72.88% 5.45% 4768 0 4717 1403831 4414395 1 Escherichia coli str. K-12 substr. MC4100

43.88% 30.91% 97.24% 64.02% 9.10% 4402 0 4717 457401 4800053 8 Escherichia sp. 4_1_40B

43.98% 30.00% 97.25% 62.08% 9.00% 4412 0 4717 2082622 4975840 193 Escherichia sp. R7

43.99% 29.79% 97.25% 61.64% 8.99% 4413 0 4717 2082619 5014011 219 Escherichia sp. R4

42.40% 30.81% 97.13% 64.98% 10.58% 4254 0 4717 457400 4659828 47 Escherichia sp. 1_1_43

44.06% 29.39% 97.26% 60.75% 8.92% 4420 0 4717 1806490 5074112 170 Achromobacter sp. ATCC35328

K12のゲノム配列をリファレンスとして使って問題なさそう。

 

Optional2:  統合レポート作成(multiqc紹介)。

fastqc fastq/*merge.fq
multiqc fastq/

html出力。極端におかしなサンプルはない。

f:id:kazumaxneo:20180613131848j:plain

 

step1−4でfastqのpreprocessingは完了。続いてmappingとリードカウントを行う。

 

 

5、リファレンス配列やコード領域情報のダウンロード。

Genome2D(HP)のbacteira(リンク)のMG1655(リンク)から、リファレンスゲノム(.fna)、プロテインテーブル(.ptt)、tRNA&rRNAテーブル(.rnt)をダウンロード(リンク)。

wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.fna 
wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.ptt
wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.rnt

> head NC_000913.rnt 

head NC_000913.rnt 

Escherichia coli str. K-12 substr. MG1655, complete genome. - 1..4641652

178 RNAs

Location Strand Length PID Gene Synonym Code COG Product

16952..17006 + 55 556503834 sokC b4413 - - -

77367..77593 + 227 556503834 sgrS b4577 - - -

189712..189847 + 136 556503834 tff b4414 - - -

223771..225312 + 1542 556503834 rrsH b0201 - - 16S ribosomal RNA of rrnH operon

225381..225457 + 77 556503834 ileV b0202 - - Ile tRNA

225500..225575 + 76 556503834 alaV b0203 - - Ala tRNA

225759..228662 + 2904 556503834 rrlH b0204 - - 23S ribosomal RNA of rrnH operon

head NC_000913.ptt 

$ head NC_000913.ptt 

Escherichia coli str. K-12 substr. MG1655, complete genome. - 1..4641652

4140 proteins

Location Strand Length PID Gene Synonym Code COG Product

190..255 + 21 16127995 thrL b0001 - - thr operon leader peptide

337..2799 + 820 16127996 thrA b0002 - COG0527E,COG0527E Bifunctional aspartokinase/homoserine dehydrogenase 1

2801..3733 + 310 16127997 thrB b0003 - COG0083E,COG0083E homoserine kinase

3734..5020 + 428 16127998 thrC b0004 - COG0498E,COG0498E L-threonine synthase

5234..5530 + 98 16127999 yaaX b0005 - - DUF2502 family putative periplasmic protein

5683..6459 - 258 16128000 yaaA b0006 - COG3022S,COG3022S peroxide resistance protein, lowers intracellular iron

6529..7959 - 476 16128001 yaaJ b0007 - COG1115E,COG1115E putative transporter

 

 

6、EDGE-proによるマッピング定量

EDGE-proでマッピングを実行。forループで処理。不安なら1つずつ。

#最初にダウンロードしたEDGE-proのプログラムディレクトリ(EDGE_pro_v1.3.1)がカレントにあるものとする。
#出力ディレクト
mkdir read_count

#edge.plを実行。
for srr in $(cat sra_run_ids.txt)
do

perl EDGE_pro_v1.3.1/edge.pl -g NC_000913.fna -p NC_000913.ptt -r NC_000913.rnt \
-u fastq/${srr}_merge.fq -o read_count/${srr}.out \
-s EDGE_pro_v1.3.1/ -t 20 -c 3
done

確認できないようなハイスピードで多量のログが標準出力されるのでトラブルではないかと心配になりますが、エラーコードが出てこなければ大丈夫です。

1サンプルの出力を確認。 

$ head SRR5416993.out.rpkm_0

gene_ID                  start_coord       end_coord     average_cov          #reads            RPKM

 

b0001                            190             255            40.5              19              84

b0002                            337            2799             3.5              61               7

b0003                           2801            3733            16.3             109              34

b0004                           3734            5020            13.0             120              27

b0005                           5234            5530            11.3              24              23

b0006                           5683            6459            35.6             198              74

b0007                           6529            7959            19.8             203              41

b0008                           8238            9191          1053.5            7179            2185

空行のほか、大量にスペースができてしまう。次のステップで消す。

 

7、形式を変更。

DESeq2のカウントテーブル形式に変換するスクリプトを走らせる。

cd read_count/

perl ../EDGE_pro_v1.3.1/additionalScripts/edgeToDeseq.perl \
SRR5416993.out.rpkm_0 SRR5416994.out.rpkm_0 SRR5416995.out.rpkm_0 \
SRR5416999.out.rpkm_0 SRR5417000.out.rpkm_0 SRR5417001.out.rpkm_0

 perlスクリプト処理でエラーが起きたので、手作業で修正する。DESeqやedgeRで利用するリードカウントを抽出する。rpkmはいらない。

 

大量のスペース、余分な列(2列目)、そして空行があるのでsedgrepで消す。そしてリードカウントのカラムだけ抽出し、各サンプルのカラムを結合する。rpkmを使うなら、5カラムでなく6カラム目を抽出する。ループにし忘れましたが、できればループ処理してください。

sed -e 's/s+/	/g' SRR5416993.out.rpkm_0 |cut -f 1 | grep -v '^$' > name
sed -e 's/s+/ /g' SRR5416993.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416993_read_count
sed -e 's/s+/ /g' SRR5416994.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416994_read_count
sed -e 's/s+/ /g' SRR5416995.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416995_read_count
sed -e 's/s+/ /g' SRR5416999.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416999_read_count
sed -e 's/s+/ /g' SRR5417000.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5417000_read_count
sed -e 's/s+/ /g' SRR5417001.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5417001_read_count

#出力ファイルのカラムを結合
paste name SRR5416993_read_count SRR5416994_read_count SRR5416995_read_count SRR5416999_read_count SRR5417000_read_count SRR5417001_read_count > count.txt

#vi count.txt #1列目を修正。以下のようにした。
> head -n count.txt 

original count.txt 

$ head -n 1 count.txt 

gene_ID #reads #reads #reads #reads #reads #reads

modified count.txt 

$ head count.txt 

gene_ID WT1 WT2 WT3 K100R1 K100R2 K100R3

b0001 19 21 23 9 109 31

b0002 61 69 80 47 89 127

b0003 109 120 152 132 172 181

b0004 120 129 116 157 125 207

b0005 24 9 11 13 22 26

b0006 198 195 234 241 157 303

b0007 203 197 202 243 194 203

b0008 7179 7028 8557 7635 7808 8414

b0009 345 276 340 337 403 351

 

 

 

8、edgeRやDEseq2で正規化、DEGを抽出。Rで可視化。

本ツールではrpkmを出しますが、現在はrpkmは推奨されていません。最近の論文(pubmed)のように、リードカウントの方を使いedgeRやDEseq2で正規化してDEGを出してください。edgeRの流れについては、手前味噌ですが、2017年に勉強会用の記事をまとめましたので、参考にしてください(シロイヌナズナのRNA seq解析)。DEseq2を使った流れについては、北大の中川先生がデータのダウンロードから分かりやすく説明されています。

https://ncrna.jp/blog/item/388-deseq2-ggplot2

また次のリンクは英語になりますが、こちらはEDGE-proのチュートリアルです。

https://bioinformatics.uconn.edu/edge-pro-tutorial/

データの取り込みからDESeq2の検定まで通して説明されています。ただしDEseq2のアップデートで一部使えないコマンドが出てきている?ようです(未確認)。

 

 

 

 

バクテリアRNA seqを行なっている人は、情報が少なく苦労されてると思います。今更かもしれませんが、参考にしてください。

 

このツールにこだわらなければ、Kallistoでもいいかもしれません。

https://groups.google.com/forum/#!topic/kallisto-sleuth-users/U8KdYSCjPyE

https://groups.google.com/forum/#!topic/kallisto-sleuth-users/Yu_T4jkgSuE

https://www.biostars.org/p/238650/

Nature Biotechnology volume 34, pages 525–527 (2016)

https://www.nature.com/articles/nbt.3519

 

引用

EDGE-pro: Estimated Degree of Gene Expression in Prokaryotic Genomes

Magoc T1, Wood D, Salzberg SL.

Evol Bioinform Online. 2013;9:127-36.

 

参考文献

A survey of best practices for RNA-seq data analysis

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4728800/

The impact of quality filter for RNA-Seq

https://www.sciencedirect.com/science/article/pii/S0378111915003145

Preprint

FADU: A Feature Counting Tool for Prokaryotic RNA-Seq Analysis

https://www.biorxiv.org/content/early/2018/06/03/337600

 

参考HP

RNA seq Blog

https://www.rna-seqblog.com/edge-pro-estimated-degree-of-gene-expression-in-prokaryots/

How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use?

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4878611/

Question: Prokaryotic Differential expression analysis- RNA seq data and software

https://www.biostars.org/p/217160/

 

GUIツール

https://cs.wellesley.edu/~btjaden/Rockhopper/