macでインフォマティクス

macでインフォマティクス

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

ハイブリッドアセンブリのためのアライメントフリー scaffolding graph構築ツール Fast-SG

2018 10/26 タイトル修正 

 

 ゲノム全体のデノボアセンブリの主要な課題は、リピートを解決することである[論文より 1,2]。リピートは、ゲノムの複数の位置で生じるほぼ同一のゲノム配列に対応する。この課題に対処するために、主に2つのタイプのアプローチが提案されている。一つはのショート・リードのペアの情報[ref.3]を使ったアプローチと、それ以外のロング・リード[ref.4]を使ったアプローチである。後者は、ロングリード中にリピートを完全にキャプチャすることが目的である。そのようなロングリードの非反復サフィックスプレフィックスシーケンスを使用してユニークなオーバーラップを計算すると、元のリードをコンティグと呼ばれる大きな配列に曖昧さなく拡張でき、時には(必ずしもそうではないが)完全なゲノム配列を直接推測できる。前者のアプローチは、genome scaffoldingと呼ばれる操作に関連している必要がある。ショートリードは、オーバーラップの計算[rerf.5]またはde Bruijn graph[ref.6]のいずれかによって、上記のように最初にコンティグに組み立てられる。しかし、この場合得られたコンティグは全ゲノムにまたがっていない。ほとんどの場合、ずっと短くなる。それらは、第2のステップで結合される(すなわち、一緒にリンクされる)必要がある。リンキング情報は、一般に、ペアエンドまたはメイトペアシーケンスによって提供される。一般的に、両端がシーケンスされている1kbより大きいゲノムフラグメントはメイトペアライブラリと呼ばれ、そうでなければペアエンドライブラリと呼ばれる。ペアリードを使用するgenome scaffoldingは、コンティグ間にギャップ(すなわち、未知の配列)を導入し、それによってゲノム配列全体には至らず、いわゆるscaffold配列のセットにつながる。したがって、scaffoldは、順序と方向が決まったコンティグのセットを表す。

 Genome scaffoldingの問題は、最初にHusonらによって定式化された[ref.7]。著者らが提案した方法は、ノードがコンティグを表現し、エッジがメイトペア(重み)の数、方向、および2つの異なるコンティグ間の距離を符号化する、scaffold graphと呼ばれるものを構築することから始まった。greedyアルゴリズムを使用して、ヒューリスティックにscaffold配列に対応する最適な経路を得る。

 Husonらの定式化から開発されたスキャフォールディング法の大部分は、同種のグラフを使用しており、超高速のショートリードアライナー[ref.8-10]でscaffoldの基礎を構築している[ref.3]。この分野におけるアルゴリズムの革新は、主に最適経路(通常は最大の重み)を選択する方法に重点を置いており、したがって、精度の高いlarge scaffoldsを得ることができる。ダイナミックプログラミング[ref.11]、breadth-first search [ref.12]、maximum weight matching[ref.13]、branch and bound[14]などの様々なアプローチが提案されている。

 新しいロングリードシークエンシング技術(Pacific Biosciences、Oxford Nanopore)は、高レベルのエラー(現在の平均15%)を含む非常に長いリード(> 10kb)を生成することによって、急激にゲノムアセンブリを変更した。これらの新技術は解決可能なリピート配列の景観を拡大した[ref.15]。現在、このようなロングリード[ref.4,16]を使用するde novoアセンブラは、バクテリアゲノムをFinishさせ、ヒトゲノムの高度に連続した再構成を生成することができる[ref.4,17]。しかし、オーバーラップを計算すること[ref.5]に基づく大規模なゲノムのデノボアセンブリは計算的に激しく[ref.4]、エラーを含むロングリードシーケンスのエラー訂正のためにかなりのカバレッジ(50X)を必要とし、これらの方法の大きなゲノムのデノボアセンブリへの幅広い適用を妨げている[17]。

 しかしながら、新規ライブラリ作成技術の補完的な長距離情報と関連させると、ロングライブラリを用いたデノボアセンブリは染色体レベルにスケーラビリティがあることが証明されている[18,19] [20,21]。このような新規の実験ライブラリは、イルミナ(Illumina)のマシンでシーケンスされ、従来のペアエンドリードに至る。したがってDOVETAILゲノミクス[20](参考ページ)は1-200kbの範囲で有用なリンク情報を生成し、一方10Xゲノミクス[ref.22]は、バーコードを賢く活用することで、最大100kbのリンクリードを生成する。どちらの技術も、アセンブリパイプライン内で長距離情報を使用してscaffolding graphを構築し、scaffolding配列を得るための独自のアルゴリズム解を適用する。両方の技術は、高価で時間のかかる長距離のメイトペアライブラリを作製するために必要な実験プロトコルを、ショートリードシーケンスで置き換えることを目指して考案された[23,24]。

 原則的に、ロングレンジの情報は、後者の、ロングリードの実際のサイズの制限範囲内で直接抽出することができる。このような情報はハイブリッドアセンブリ法を考案するために使用することができ、ショートリードアセンブリから高品質のコンティグがscaffolding graphのノードとして使用され、エッジがロングリードからのリンク情報を使用して作成され、scaffoldはショートリードscaffolderで形成される。しかし、現在、ロングリードからscaffolding graphを構築するためのアルゴリズムが欠けている。そのようなアルゴリズムは、効率的な既存のショートリードアルゴリズムを新しいハイブリッドアセンブリパイプラインを構成するのに再利用することを可能にする。

 既存の効率的なショートリードスキャナーと互換性がある構造標準を維持しながら、適度な計算資源を用いて、ショートリードまたはロングリードから超短時間でこのようなgraphを構築できることは、ここで取り上げる主な課題である。我々(著者ら)が提案しているFAST-SGは、さまざまな配列ソース(イルミナ、パシフィックバイオサイエンス、オックスフォードナノポア)からの情報だけでなく、アライメントフリーのアルゴリズム[ref.25]ストラテジーを使用し、スケーラビリティ、スピード、およびモジュール性を最大化するように考えられた。後者の特徴は、特に、大きなゲノムの効率的なアセンブリを可能にする新規のハイブリッドアセンブリパイプラインを定義することを可能にする。

 FAST-SGは包括的な一連の標準データセット[ref.3、26]とベンチマークを使用して広範にテストされた。我々(著者ら)は、FAST-SGが大きなゲノムのハイブリッドアセンブリを可能にし、特にカバレッジの浅いロングリードカバレッジデータ(5X〜10X)で有効であることを示す。このハイブリッド戦略は、バクテリア人工染色体(BACs、180kb)までのインサートサイズを有することができるいくつかの合成メイトペアライブラリーの構築から成り、長距離のscaffoldを生成するためにショートリードscaffolderと組み合わせることができる。このような戦略は、適度な計算資源でヒトゲノムサイズまでスケーリングできる。さらに、FAST-SGは従来のショート・リード・アライナーよりも高速(7X-15X)で、ショートリードのメイトペアデータを使用したスキャフォールディングの強力な代替手段であることを示す。

 我々(著者ら)は、FAST-SGの効果的なハイブリッドアセンブリのための手順を提供することによって結論を出し、提案する戦略がどのようにロングリードを使用してギャップを埋めるように拡張され、そしてscaffoldsのエラーを訂正するか議論する。

 

Fast-SGに関すつツイート

 

インストール

mac os10.12でテストした。

依存

  •  KMC

本体 Github

git clone https://github.com/adigenova/fastsg
cd fastsg/
make all

#KMCのダウンロード(pre-compile) ./KMC/bin/にないと認識しないので解凍後、KMC/bin/にKMCバイナリを放り込む
#mac
wget https://github.com/refresh-bio/KMC/releases/download/v3.0.0/KMC3.mac.tar.gz
tar -zxvf KMC3.mac.tar.gz
mv KMC3.mac KMC
cd KMC
mkdir bin && mv kmc* bin/ && cd ../

make test

 test結果

ont.I1000.FastSG_K15.sam

ont.I2000.FastSG_K15.sam

ont.I3000.FastSG_K15.sam

End of mapping in thread

Query of 56850535 kmer-pairs  in 34.01s rate 1671387.23 per second

Query of 286428 pairs  in 34.01s

Number of total mapped kmer-pairs 4751139

Number of total mapped reads-pairs 144987

Number of total long reads 1298

Number of total mapped pairs: 144988

touch test.fast-sg.K15.done

make[1]: Leaving directory '/Volumes/4TB-gene_new2/Fast-SG/fastsg'

 

 

ENDING Fri Oct 26 20:06:55 2018 FAST-SG with K=15

 

all done

正常に動いている。手動で行うなら

cd fastsg/ #fastsgのrootに移動
perl FAST-SG.pl -k 15 -l example/ecoli-reads.txt -r example/ecoli-illumina.fa.gz -p test -t 20

ヘルプ

perl FAST-SG.pl   #"FastSG++"をperlのラッパーで呼び出している。

$ perl FAST-SG.pl 

 

 

ERROR: Mandatory options need to be specified -k -l -r -p

 

 

 

Usage: FAST-SG.pl  [<arguments>]

 

 

Basic options:

 

-k *K-mer size [12-256] could be a single value or a range [15-40:5]

-l *Read configuration file

-r *Contig fasta file

-p *Output prefix

-t Number of thread to be used [1]

 

 

Advanced options:

 

Illumina reads options:

-x Vector hit size [10]

-y Mimimum score to report reads [6]

-z Mimimum score for pair rescue [4]

 

The value of -x must be larger than -y and -y larger than -z [-x > -y > -z]. 

To modify these values take into account the k-mer size, the length and error rate of the reads.

Specified values must be larger than defaults values.

 

Long reads options:

-u Vector hit size [20]

-v Mimimum score to report a pair of reads [15]

-w Mimimum score for pair rescue [4]

-s Length of syntethic read [200]

 

The value of -u must be larger than -v and -v larger than -w [-u > -v > -w]. 

To modify these values take into account the k-mer size, the length and error rate of the reads.

Specified values must be larger than defaults values.

 

Other options:

-c Full path to KMC3 directory [FAST-SG-DIR/KMC]

-N Show pipeline commands [-N 1]

-R Delete results to restart FAST-SG [-R 1]

--help Show help message

--version Display current version of FAST-SG

 

Note: Options with description starting with (*) are mandatory.

 

 

テストラン

NA12878のハイブリッドアセンブリを例にフローが説明されている。この流れを確認する。

https://github.com/adigenova/fast-sg/wiki/Hybrid-scaffolding-of-NA12878

 

1、illuminaのショートリードのde novoアセンブリ結果をダウンロード。

cd fastsg/ #fastsgのrootに移動

#DiscoverDENOVOによるilluminaショートリードのアセンブリ結果をダウンロードする
wget ftp://ftp.broadinstitute.org/pub/crd/Discovar/assemblies/51400.newchem/a.lines.fasta -O NA12878.asm.fa

#fa.faiを利用して2kb以上のcontigだけ取り出す。
samtools faidx NA12878.asm.fa
sort -rn -k2,2 NA12878.asm.fa.fai | awk '{if($2 >=2000){print $1}}' > NA12878.asm.2kb.lst
samtools faidx NA12878.asm.fa $(cat NA12878.asm.2kb.lst | xargs) > NA12878.asm.2kb.fa

 

2、ONTのロングリードをダウンロードし、コンカテネートしてFASTAに変換。このリンク先のデータ(link)が使われている(preprint)。

cat > downlod.txt <<EOF
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-16056159-FAF15665.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-17958431-FAF13748.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-2901545329-FAF10039.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3439856925-FAF09968.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3709819546-FAF09277.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-3976726082-FAF14035.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4109802543-FAF15694.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4111860526-FAF09713.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4178920553-FAF18554.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4244782843-FAF15630.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4245291640-FAF09640.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-4249180049-FAF09701.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-82266371-FAF15586.fastq.gz
http://s3.amazonaws.com/nanopore-human-wgs/rel4-nanopore-wgs-87644245-FAF05869.fastq.gz
EOF

#wgetにテキストで与えダウンロード。
mkdir fastq && cd fastq
wget -i downlod.txt

#concatenate
zcat fastq/*.fastq.gz | awk 'BEGIN{h=1;s=2;i=1}{if(NR==h){print ">ONT_"i; h+=4;i++;} if(NR == s){print $0; s+=4;}}' | fold | gzip > ULTRA-LONG-RENAME-FOLD.fa.gz

 

3、configファイルを準備する。Fast-SGを使い、ONTのロングリードを使ってilluminaショートリードアセンブリをscaffoldingして、2kbから180kbまでのsynthetic librariesを作成する場合、以下の書式になる。

cat > ultra-long-conf.txt <<EOF
long ont ULTRA-LONG-RENAME-FOLD.fa.gz 2000,4000,6000,8000,10000,12000,14000,16000,18000,20000,30000,40000,50000,60000,70000,80000,100000,120000,150000,180000 1
EOF

以下のような書式になっている。

 

#shortリードの場合
#type libID Path(fwd) Path(rev) SAM(1:single 2:paired)
short lib1 example/ecoli.ill-sim.fwd.fq.gz example/ecoli.ill-sim.rev.fq.gz 1

#longリードの場合
#type libID Path(long read) Insert-sizes SAM(1:single 2:paired)
long pac example/ECOLI-PACSEQ.subset.fasta.gz 1000,2000,3000,5000 1
long ont example/ECOLI-ONT-1D.subset.fasta.gz 1000,2000,3000,5000 1

 example  https://github.com/adigenova/fast-sg/blob/master/example/ecoli-reads.txt

 

4、実行

perl FAST-SG.pl -k 22 -l ultra-long-conf.txt –r NA12878.asm.2kb.fa.gz -p NA12878-ultra –t 20 

 -k、-l、-r、-pはx必須。複数k-merを指定する際は"-k 15-25:5"のように指定する。エラー訂正されていないロングリードを使う場合は、15-22程度の短いk-merを使うことが推奨されている。

ジョブが終わると、insert librariesで指定したサイズのsamファイルができる。

*私の環境では4日経ってもジョブが終わらなかったので、一旦テストを終了します。  

 

5、samをforwardとreverseに分割する。

awk -f fastSG2scaff.awk -v name=$(basename ultra_ont_raw.I2000.fast-sg_K22.sam) -v k=22 ultra_ont_raw.I2000.fast-sg_K22.sam

 

6、ScaffMatch(リンク)を使ってScaffoldingを実行する。

scaffmatch -m -w UONT-K22-7X -c NA12878.asm.2kb.fa \
-s 201,416,632,847,1063,1279,1495,1710,1925,2141,3216,4291,5364,6434,7495,8555,10675,12781,15946,19133 \
-i 2010,4166,6323,8478,10636,12793,14950,17106,19259,21414,32169,42912,53641,64340,74955,85557,106753,127813,159463,191339 \
-p fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr,fr \
-1 ultra_ont_raw.I2000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I4000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I6000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I8000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I10000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I12000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I14000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I16000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I18000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I20000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I30000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I40000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I50000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I60000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I70000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I80000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I100000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I120000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I150000.Fast-SG_K22.fwd.sam,ultra_ont_raw.I180000.Fast-SG_K22.fwd.sam \
-2 ultra_ont_raw.I2000.Fast-SG_K22.rev.sam,ultra_ont_raw.I4000.Fast-SG_K22.rev.sam,ultra_ont_raw.I6000.Fast-SG_K22.rev.sam,ultra_ont_raw.I8000.Fast-SG_K22.rev.sam,ultra_ont_raw.I10000.Fast-SG_K22.rev.sam,ultra_ont_raw.I12000.Fast-SG_K22.rev.sam,ultra_ont_raw.I14000.Fast-SG_K22.rev.sam,ultra_ont_raw.I16000.Fast-SG_K22.rev.sam,ultra_ont_raw.I18000.Fast-SG_K22.rev.sam,ultra_ont_raw.I20000.Fast-SG_K22.rev.sam,ultra_ont_raw.I30000.Fast-SG_K22.rev.sam,ultra_ont_raw.I40000.Fast-SG_K22.rev.sam,ultra_ont_raw.I50000.Fast-SG_K22.rev.sam,ultra_ont_raw.I60000.Fast-SG_K22.rev.sam,ultra_ont_raw.I70000.Fast-SG_K22.rev.sam,ultra_ont_raw.I80000.Fast-SG_K22.rev.sam,ultra_ont_raw.I100000.Fast-SG_K22.rev.sam,ultra_ont_raw.I120000.Fast-SG_K22.rev.sam,ultra_ont_raw.I150000.Fast-SG_K22.rev.sam,ultra_ont_raw.I180000.Fast-SG_K22.rev.sam

 UONT-K22-7X/Scaffolds.faができる。

テスト中 

 

引用

Fast-SG: an alignment-free algorithm for hybrid assembly

Di Genova A, Ruz GA, Sagot MF, Maass 

Gigascience. 2018 May 1;7(5).