近年のハイスループットシーケンシング(HTS)技術は、低コストで数百万の短いDNA配列(リードと呼ばれる)を生成するため、デノボアセンブリプロジェクトにとって魅力的である。しかしながら、これらのリードは数百bpの長さしかないため、アセンブラ(例えば、[論文より ref.1,2])がゲノムを再構築することを困難にする。その結果、アセンブリ出力はコンティグ、すなわち長い断片にアセンブリされたゲノムの断片からなることが多い。
しかし、HTS技術は、アセンブリの連続性を高めるために使用できるペアエンドリードを作成するためのプロトコルを提供する。ペアエンドリードは、リード間の距離をインサートサイズと呼ぶ既知の距離と方向でシーケンシングされている。ペアエンドリード内の2つのリードが異なるコンティグCaおよびCbに属する場合、CaとCbの間にリンクが作成される(論文 図1aを参照 link)。このリンクから、CaとCbの相対的な順序、方向、距離を推論することができる。
コンティグのリンクと順序付けのプロセスはscaffoldingと呼ばれている。ペアエンドリードに加え、関連生物のリファレンス配列[ref.3]、制限酵素地図[ref.4]およびRNA-seqデータ[ref.5]などの情報をコンティグリンクに用いることができる。しかし、リファレンスベースのアセンブリは、ほとんどのde novoシーケンシングプロジェクトには適用できない。制限酵素地図はしばしば利用できず、RNA-seqデータは遺伝子をカバーするだけで、リード間の距離に関する情報は含まれていない。これにより、(しばしば唯一適用可能な)ペアエンド情報がスキャフォールディングの情報源として最も一般的に使用されている。
残念なことにペアエンドのscaffoldingは問題に直面する。リードは、シーケンスエラー、ヘテロ接合性、およびゲノムのリピートのため擬似リンクを作成する可能性があり、これらの擬似リンクはコンティグ間の順序付けおよび方向付けをあいまいにする。したがってscaffoldingの問題とは、コンティグの一貫した順序付けおよび配向を見つけるために、正しいリンクを検出および利用すること、とまとめられる。scaffoldingの既存の公式化はNP完全であることが証明されているが、これらの公式化が目的に関して最適な解を見出すときでさえ、実際の(すなわち生物学的な)問題を解決するかどうかはまだ不明である。これらのアプローチは、個々のリンクの正確性を評価することにほとんど重点を置かずに、コンティグリンクによって誘導されるグラフの構造的特性に焦点を当ててきた。著者らのアプローチでは、間違ったリンクを削除し、洗練された統計を使用して、リンクしているリードが基礎となるライブラリの配布からのものであるのか、 ミスアライメントによるものか評価する。
論文の次のセクションでは、scaffoldingとその関連作業の形式化と、作業の概要とこの仕事のmotivationについて説明する。我々(著者ら)のアルゴリズムは、BESST(Bias Estimating Stepwise Scaffolding Tool)と呼ばれるインプリメンテーションで実現され、メソッドのセクションで詳細に説明されている。このアルゴリズムは、Picea abies(20 Gbp)のゲノムプロジェクト[ref.6]の使用で証明されているように、非常に大きく複雑なゲノム上にもスケールする。さらに、より幅広いインサートサイズ分布でも優れたscaffolding能を持つ。
(一段落省略)
我々(著者ら)の結果は、単一のscaffolderがすべてのデータセットで他のスキャフォルダーより優れているわけではないことを示しているが、他のスタンドアロンscaffolderとの比較でBESSTは最も好ましい結果を示している。さらに、ライブラリのインサートサイズ分布の標準偏差が大きい場合、他のスタンドアローンのscaffolderよりも優れている(以下省略)。
BESSTに関するツイート
インストール
cent os7でテストした。
本体 Github
conda install -c bioconda besst
> runBESST
runBESST
usage: BESST [-h] -c CONTIGFILE -f BAMFILES [BAMFILES ...] -orientation
ORIENTATION [ORIENTATION ...] [-r READLEN [READLEN ...]]
[-m MEAN [MEAN ...]] [-s STDDEV [STDDEV ...]] [-z COVCUTOFF]
[-z_min LOWER_COVCUTOFF] [-T THRESHOLD [THRESHOLD ...]]
[-e EDGESUPPORT [EDGESUPPORT ...]] [-k MINSIZE [MINSIZE ...]]
[-filter_contigs CONTIG_FILTER_LENGTH] [--min_mapq MIN_MAPQ]
[--iter PATH_THRESHOLD] [--score_cutoff SCORE_CUTOFF]
[--max_extensions MAX_EXTENSIONS] [-a HAPLRATIO]
[-b HAPLTHRESHOLD] [-K KMER] [-M MMER] [-g] [-o OUTPUT] [-d] [-y]
[-q] [--no_score] [-devel] [-plots] [--separate_repeats]
[--NO_ILP] [--FASTER_ILP] [--print_scores] [--dfs_traversal]
[--bfs_traversal] [-max_contig_overlap MAX_CONTIG_OVERLAP]
[--version]
BESST - scaffolder for genomic assemblies.
optional arguments:
-h, --help show this help message and exit
-y Deactivate pathfinder module for including smaller
contigs.
--no_score Statistical scoring is not performed. BESST instead
searches for paths between contigs.
Required arguments:
-c CONTIGFILE Fasta file containing contigs.
-f BAMFILES [BAMFILES ...]
Path(s) to bamfile(s).
-orientation ORIENTATION [ORIENTATION ...]
Mapped orientation for each library given with -f
option. Valid input is either fr (forward reverse
orientation) or rf (reverse forward orientation).
Library information:
Library parameters that can be set, e.g., read length, insert size mean/std_dev, coverage etc.
-r READLEN [READLEN ...]
Mean read length of libraries.
-m MEAN [MEAN ...] Mean insert size of libraries.
-s STDDEV [STDDEV ...]
Estimated standard deviation of libraries.
-z COVCUTOFF User specified coverage cutoff. (Manually filter out
contigs with coverage over this value)
-z_min LOWER_COVCUTOFF
User specified coverage cutoff. (Manually filter out
contigs with coverage under this value)
Ploidy:
Options involving detection of split allele contigs in diploid assemblies.
-a HAPLRATIO Maximum length ratio for merging of haplotypic
regions.
-b HAPLTHRESHOLD Number of standard deviations over mean/2 of coverage
to allow for clasification of haplotype. Example:
contigs with under mean/2 + 3sigma are indicated as
potential haplotypes (tested later) if -b 3
-g Haplotype detection function, default = off
Parameters/variables/threshold involved in BESST.:
-T THRESHOLD [THRESHOLD ...]
Threshold value filter out reads that are mapped too
far apart given insert size.
-e EDGESUPPORT [EDGESUPPORT ...]
Threshold value for the least nr of links that is
needed to create an edge. Default for all libs:
Inferred by BESST (see value in output in
Statistics.txt).
-k MINSIZE [MINSIZE ...]
Contig size threshold for the library (all contigs
below this size is discarded from the 'large contigs'
scaffolding, but included in pathfinding). Default:
Set to same as -T parameter
-filter_contigs CONTIG_FILTER_LENGTH
Contigs under this size is discarded from all
scaffolding (including pathfinding). they are also not
included in the resulting scaffold output fasta file
and any of the other files as well as all statistics.
Default: All contigs are included
--min_mapq MIN_MAPQ Lowest mapping quality allowed in order to use the
read pair (both reads needs to have equal or bigger
mapq than this value). This value is compared to the
mapping quality column in the BAM file.
--iter PATH_THRESHOLD
The number of iterations performed in breadth first
search for placing smaller contigs.
--score_cutoff SCORE_CUTOFF
Only store paths with score larger than score_cutoff.
--max_extensions MAX_EXTENSIONS
Maximum number of (large) scaffolds to search for
paths extensions with.
--NO_ILP Warning: Do not use this option! This option is
included only for benchmarking purposes of old BESST
algorithm. This option will give poor results as it
mimics earlier versions of BESST.
--FASTER_ILP Faster but worse performing heuristic solution to
solving ILPs. May be used if BESST is too slow.
However, lowering --iter is usually more effective to
reduce scaffolding time.
--print_scores Print BESST scores on edges in the Scaffolding graph.
Various other parameters:
-K KMER k-mer size used in de brujin graph for obtaining
contigs in assembly, default 50
-M MMER m-mer usted for creating connection graph. Should be
set lower than k-mer size
-o OUTPUT Path to output directory. BESST will create a folder
named 'BESST_output' in the directory given by the
path.
-d Deactivate sequencing duplicates detection
-q Parallellize work load of path finder module in case
of multiple processors available.
-devel Run in development mode (bug checking and memory usage
etc.)
-plots Plot coverage, histograms of scores e.t.c.
--separate_repeats Separates contigs classified as repeats by BESST into
a file 'repeats.fa'. They are not included in the main
scaffolding output file with this flag specified.
--dfs_traversal Depth first search with DP in the contig graph
(default).
--bfs_traversal Choose a breadth first search in the contig graph.
Default is the new depth first search with a DP
approch that seems to outperform previous traversals.
This option is kept for evaluation but is not
reccomended to specify.
-max_contig_overlap MAX_CONTIG_OVERLAP
BESST checks for overlapping ends in contigs that are
adjacent in a scaffold. This parameter sets the
maximum identical overlap to search for, default is
200.
--version show program's version number and exit
Source code and manual: http://github.com/ksahlin/BESST
Please cite:
Sahlin K, Vezzi F, Nystedt B, Lundeberg J, Arvestad L (2014) BESST--efficient scaffolding of large fragmented assemblies. BMC Bioinformatics 15, 281
ラン
アセンブリしたcontigと、そのcontigにマッピングしたbamを指定する。複数ライブラリのシーケンシングデータがある時は各々のbamを順番に指定する。
runBESST -c contig.fa -f file1.bam file2.bam ... -o output -orientation {fr/rf}
ランの詳細はGithubのREADMEを確認して下さい。
引用
BESST--efficient scaffolding of large fragmented assemblies
Sahlin K, Vezzi F, Nystedt B, Lundeberg J, Arvestad L
BMC Bioinformatics. 2014 Aug 15;15:281
Assembly scaffolding with PE-contaminated mate-pair libraries.
Sahlin K, Chikhi R, Arvestad L
Bioinformatics. 2016 Jul 1;32(13):1925-32