macでインフォマティクス

macでインフォマティクス

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

ロングリードのマルチプルシーケンスアラインメントを行う abPOA (POAのSIMD実装拡張)

 

 マルチプルシーケンスアラインメント(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やコンセンサスコールツールとして機能するほか、任意のロングリードシーケンシングエラー訂正やアセンブリワークフローに簡単に統合することができる。

 

 

インストール

ubuntu18.04LTSとmacos10.14でテストした。

ビルド依存

  • gcc (>=6.4.0) and zlib

Github

#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
  • -r    output result mode [0]
     0: consensus (FASTA format)
     1: MSA (PIR format)
     2: both 0 & 1

出力

 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