最近のロングリードシークエンシング技術の著しい進歩にもかかわらず、二倍体ゲノムのアセンブルは依然として困難な課題である。主な障害は、高度にヘテロ接合性のある領域を表すオルタナティブなコンティグを区別することである。プライマリなコンティグとセカンダリなコンティグが適切に識別されない場合、一次アセンブリはゲノムの大きさと複雑さの両方を過剰に表し、scaffolds形成などの下流の解析を複雑にしてしまう。
HapSoloは、セカンダリコンティグを識別し、複数のペアワイズコンティグアライメント指標に基づいて一次アセンブリを定義するものである。HapSoloはBUSCOスコアを用いて候補のプライマリアセンブリを評価し、それからコスト関数を用いて候補のアセンブリを区別する。コスト関数はユーザーが定義することができるが、デフォルトではアセンブリ内の欠失、重複、単一BUSCO遺伝子の数を考慮する。HapSoloは何千もの候補アセンブリのコストを最小化するために最急降下法を行う。シャルドネブドウ（Vitis vinifera）、蚊（Anopheles funestus）、韓国のマッドスキッパー（Periophthalmus magnuspinnatus）の3種のゲノムデータを用いて、HapSoloの性能を示した。
HapSoloは候補アセンブリを迅速に同定し、ゲノムサイズの縮小やN50スコアの改善などアセンブリの指標を劇的に改善する。N50スコアは、縮小前のプライマリアセンブリと比較して、シャルドネ、蚊、マッドスキッパーでそれぞれ26％、8.5％、21％向上した。HapSoloの利点はダウンストリーム解析によって増幅され、Hi-Cデータを用いたスキャフォールドによって説明された。例えば、HapSolo を適用する前は、染色体の数に相当する最大の 19 のスキャフォールドに、シャルドネのゲノムの 39%しか捕捉されていなかったことがわかった。HapSoloの適用後、この値は77％に増加した。
インストール
condaで仮想環境を作ってテストした。
依存
conda create -n hapsolo python=2.7 -y
conda activate hapsolo
pip install pandas
conda install -c bioconda ucsc-fatotwobit
conda install -c bioconda busco
本体 Github
git clone https://github.com/esolares/HapSolo.git
cd HapSolo/
#パスを通す
export PATH=$PWD:$PATH
export PATH=$PWD/scripts:$PATH
> preprocessfasta.py -h
$ preprocessfasta.py -h
usage: preprocessfasta.py [-h] -i INPUT [-m MAXCONTIG]
Preprocess FASTA file and outputs a clean FASTA and seperates contigs based on
unique headers. Removes special chars
optional arguments:
-h, --help show this help message and exit
-i INPUT, --input INPUT
Input FASTA file
-m MAXCONTIG, --maxcontig MAXCONTIG
Max size of contig in Mb to output in contigs folder.
> hapsolo.py -h
$ hapsolo.py -h
usage: hapsolo.py [-h] -i INPUT -p PSL -b BUSCOS [-m MAXZEROS] [-t THREADS]
[-n NITERATIONS] [-B BESTN] [-S THETAS] [-D THETAD]
[-F THETAF] [-M THETAM]
Process alignments and BUSCO"s for selecting reduced assembly candidates
optional arguments:
-h, --help show this help message and exit
-i INPUT, --input INPUT
Input Fasta file
-p PSL, --psl PSL BLAT PSL alignnment file
-b BUSCOS, --buscos BUSCOS
Location BUSCO output directories. i.e. buscoN/
-m MAXZEROS, --maxzeros MAXZEROS
Max # of times cost function delta can consecutively
be 0. Default = 10
-t THREADS, --threads THREADS
# of threads. Multiplies iterations by threads.
Default = 1
-n NITERATIONS, --niterations NITERATIONS
# of total iterations to run per gradient descent.
Default = 1000
-B BESTN, --Bestn BESTN
# of best candidate assemblies to return using
gradient descent. Default = 1
-S THETAS, --thetaS THETAS
Weight for single BUSCOs in linear fxn. Default = 1.0
-D THETAD, --thetaD THETAD
Weight for single BUSCOs in linear fxn. Default = 1.0
-F THETAF, --thetaF THETAF
Weight for single BUSCOs in linear fxn. Default = 0.0
-M THETAM, --thetaM THETAM
Weight for single BUSCOs in linear fxn. Default = 1.0
実行方法
1、preprocessing
HapSoloのランにはコンティグごとにBlatアライメントファイルとbusco出力を格納したbuscoディレクトリが必要となる。この前処理を行うスクリプトpreprocessfasta.pyを実行する。
preprocessfasta.py -i contig.fasta
- -i Input FASTA file
- -m Max size of contig to output in contigs folder. By default the maxcontig size is set to 10000000 or 10Mb and is a not a required parameter.
このスクリプトは現在のパスに contigs ディレクトリを作成し（）、contigs/に各 contig の個別の FASTA ファイルを保存する（contig配列ごとにファイルに分けて保存する）。また、このスクリプトは、BUSCOまたはMUMmerの結果のチェック時にエラーの原因になるFASTAヘッダの不正な文字を削除する。またBlatアライメントのために配列を2bitに変換する。ジョブが終わると、contigs/以外にカレントにcontig_new.fastaが出力される。
２、blat 2bit (randomly accessible) format ファイルの作成
faToTwoBit contig_new.fasta contig_new.2bit
３、buscoのラン
４、hapsoloのラン
hapsolo.py -i contig_new.fasta -p contig_new.2bit -b busco_out_dir
連続性が低い、つまりコンティグ数が多いアセンブリを使う場合、膨大なファイルを作成するためか、I/Oにかなり負担がかかるようです。連続性が低いアセンブリ配列を使用する事はあまりないかもしれませんが、注意して下さい。
引用
HapSolo: An optimization approach for removing secondary haplotigs during diploid genome assembly and scaffolding
Edwin A. Solares, Yuan Tao, Anthony D. Long, Brandon S. Gaut
bioRxiv, Posted June 30, 2020
関連
