マルチプルシーケンスアラインメント(MSA)問題を解決するために、Leeら(2002)によって最初に導入されたのがPartial order alignment(POA)である。POAでは、MSAをdirected acyclic graph(DAG)(有向無閉路グラフ)として表現し、動的計画法(wiki)(DP)を用いて配列をDAGに反復的にアラインメントする。その後、アラインメントグラフに最も重いバンドル化アルゴリズムを適用することで、複数のコンセンサス配列が生成される(Lee, 2003)。最近、Pacific Biosciences(PacBio)およびOxford Nanopore Technologies(ONT)プラットフォーム上でのロングリードシーケンシング技術の出現と人気の高まりに伴い、POAの再評価と関心が高まっており、このアルゴリズムは現在、エラーが発生しやすいロングリードのエラー訂正とアセンブリに広く使用されている(Loman et al、2015; Vaser et al、2017; Volden et al、2018; Gao et al、2019; Ruan and Li、2019)。
POAは古典的なMSAアルゴリズム(Lassmann and Sonnham-mer, 2002)よりもはるかに高速であるが、最先端のロングリードシーケンシング・プラットフォームによって生成される大規模なデータセットが大きな課題となっている。この課題に対処するため、オリジナルのPOAアルゴリズムを高速化するSIMD(Single Instruction Multiple Data)実装が使用された(Vaser et al.) POAのこのSIMDバージョンであるSPOAは、複数の要素を並列に処理する現代のプロセッサにおけるより広いSIMDレジスタを利用する。SIMDベクタは、DP行列の各行の複数の連続したセルのスコアを格納するために使用され、SIMD命令を使用して処理され、各ベクタに格納されたすべてのスコアのために並列更新される。SIMD並列化に加えて、もう一つのアクセラレーション戦略である「banded DP」もまた、配列対配列のアライメントツールにおいて広く用いられている(Chao et al. 具体的には、DP行列の各行または列において、特定の「バンド」内のセルのみを塗りつぶす必要がある。最近、この戦略のグラフバージョンがGraphAligner(Rautiainen and Marschall, 2019)で探求されたが、この場合、DPマトリックスの各行のバンドは、その行の最小スコアセルに基づいて動的に定義される。他のセルは、スコアがその最小スコアから一定の距離内にある場合にのみ、バンド内にあるとみなされる。GraphAlignerは他のグラフアライメントツールよりも高速だが、その戦略は、エラー率の高いロングリードのMSAやコンセンサスコールには適していないかもしれない。これは、GraphAlignerが一般的なスコアリング関数ではなく編集距離をアライメント指標として使用しているためであり、PacBioやONTのシークエンシングデータによく見られる挿入・欠失エラーの多い領域では、誤ったアライメントが発生する可能性がある。
本研究では、SIMDを実装した適応型バンドDPを実行するPOAの拡張版であるabPOAを開発した。abPOAはスタンドアロンのMSAやコンセンサスコールツールとして機能するほか、任意のロングリードシーケンシングエラー訂正やアセンブリワークフローに簡単に統合することができる。
Check out abPOA, an extended version of Partial Order Alignment (POA) based on adaptive banded dynamic programming with SIMD implementation. Generate high-quality consensus seq from error-prone long reads, with significant speed gain over existing tools https://t.co/rbTANkZlAr
— Yi Xing (@YiXing77) 2020年5月6日
インストール
ubuntu18.04LTSとmacos10.14でテストした。
ビルド依存
- gcc (>=6.4.0) and zlib
#bioconda
conda install -c bioconda abpoa -y
#from source
git clone https://github.com/yangao07/abPOA.git
cd abPOA; make
#dot出力も行うならgraphvizも必要。
sudo apt-get install graphviz
> abPOA
$ abPOA
abpoa: adaptive banded Partial Order Alignment
Usage: abpoa [options] <in.fa/fq> > cons.fa/msa.out
Options:
Alignment:
-m --aln-mode INT alignment mode [0]
0: global, 1: local, 2: extension
-M --match INT match score [2]
-X --mismatch INT mismatch penalty [4]
-O --gap-open INT(,INT) gap opening penalty (O1,O2) [4,24]
-E --gap-ext INT(,INT) gap extension penalty (E1,E2) [2,1]
abPOA provides three gap penalty modes, cost of a g-long gap:
- convex (default): min{O1+g*E1, O2+g*E2}
- affine (set O2 as 0): O1+g*E1
- linear (set O1 as 0): g*E1
Adaptive banded DP:
-b --extra-b INT first adaptive banding parameter [10]
set b as < 0 to disable adaptive banded DP
-f --extra-f FLOAT second adaptive banding parameter [0.01]
the number of extra bases added on both sites of the band is
b+f*L, where L is the length of the aligned sequence
Input/Output:
-l --in-list input file is a list of sequence file names [False]
each line is one sequence file containing a set of sequences
which will be aligned by abPOA to generate a consensus sequence
-o --output FILE ouput to FILE [stdout]
-r --result INT output result mode [0]
- 0: consensus (FASTA format)
- 1: MSA (PIR format)
- 2: both 0 & 1
-g --out-pog FILE dump final alignment graph to FILE (.pdf/.png) [Null]
-h --help print this help usage information
-v --version show version number
Parameters under development:
-a --cons-alg INT algorithm to use for consensus calling [0]
- 0: heaviest bundling
- 1: heaviest in column
-d --diploid input data is diploid [False]
-a/--cons-alg will be set as 1 when input is diploid
and at most two consensus sequences will be generated
-q --min-freq FLOAT min frequency of each consensus for diploid input [0.30]
テストラン
git clone https://github.com/yangao07/abPOA.git
cd abPOA/test_data/
abpoa seq.fa > cons.fa
入力の配列
$ cat seq.fa
>1
CGTCAATCTATCGAAGCATACGCGGGCAGAGCCGAAGACCTCGGCAATCCA
>2
CCACGTCAATCTATCGAAGCATACGCGGCAGCCGAACTCGACCTCGGCAATCAC
>3
CGTCAATCTATCGAAGCATACGCGGCAGAGCCCGGAAGACCTCGGCAATCAC
>4
CGTCAATGCTAGTCGAAGCAGCTGCGGCAGAGCCGAAGACCTCGGCAATCAC
>5
CGTCAATCTATCGAAGCATTCTACGCGGCAGAGCCGACCTCGGCAATCAC
>6
CGTCAATCTAGAAGCATACGCGGCAAGAGCCGAAGACCTCGGCCAATCAC
>7
CGTCAATCTATCGGTAAAGCATACGCTCTGTAGCCGAAGACCTCGGCAATCAC
>8
CGTCAATCTATCTTCAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATC
>9
CGTCAATGGATCGAGTACGCGGCAGAGCCGAAGACCTCGGCAATCAC
>10
CGTCAATCTAATCGAAGCATACGCGGCAGAGCCGTCTACCTCGGCAATCACGT
上記のコマンドだとコンセンサス配列が出力される。
$ cat cons.fa
>Consensus_sequence
CCACGTCAATCTATCGAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATCAC
PIR fomatのマルチプルシークエンスアラインメント(多重整列)を出力。
abpoa seq.fa -r1 > cons.out
出力
cat cons.out
>Multiple_sequence_alignment
---CGTCAAT-CTA-TC--G---AAGCA---TACG--CGGGC-AGAGCC--GAA---GACCTCGG-CAATCCA--
CCACGTCAAT-CTA-TC--G---AAGCA---TACG--C-GGC---AGCC--GAACTCGACCTCGG-CAATCAC--
---CGTCAAT-CTA-TC--G---AAGCA---TACG--C-GGC-AGAGCCCGGAA---GACCTCGG-CAATCAC--
---CGTCAATGCTAGTC--G---AAGCA---GCTG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATCAC--
---CGTCAAT-CTA-TC--G---AAGCATTCTACG--C-GGC-AGAGCC--------GACCTCGG-CAATCAC--
---CGTCAAT-CTA-----G---AAGCA---TACG--C-GGCAAGAGCC--GAA---GACCTCGGCCAATCAC--
---CGTCAAT-CTA-TC--GGTAAAGCA---TACGCTC-TGT---AGCC--GAA---GACCTCGG-CAATCAC--
---CGTCAAT-CTA-TCTTC---AAGCA---TACG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATC----
---CGTCAAT-GGA-TC--G----AG-----TACG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATCAC--
---CGTCAAT-CTAATC--G---AAGCA---TACG--C-GGC-AGAGCC--G---TCTACCTCGG-CAATCACGT
abPOAはCのAPIsとしても機能する。詳細はGIthubとペーパーを確認して下さい。
引用
abPOA: an SIMD-based C library for fast partial order alignment using adaptive band
Yan Gao, Yongzhuang Liu, Yanmei Ma, Bo Liu, Yadong Wang, Yi Xing
bioRxiv, posted May 10, 2020