シーケンシングデータを解析するためのバイオインフォマティクスソリューションは数多く存在するが、系統樹の作成を最終目的とした全ゲノムシーケンス（WGS）データからの標的配列検索のためのオプションはほとんど存在しない。利用可能なツールは、特に深い系統レベルで苦労しており、アミノ酸空間検索を必要とするため、偽陽性結果の割合が増加する可能性がある。また、多くのツールはインストールが難しく、十分なユーザーリソースがない可能性がある。ここでは、自由に利用できる相同性検索ツールを用いて、アセンブルされたWGSデータから相同性検索を行うプログラムについて説明する。ゲノムアノテーションが存在する2つの分岐した昆虫種（>200 My）と、高度に保存された遺伝子座とより多様な遺伝子座のそれぞれの大規模セットについて、他の一般的なバイオインフォマティクスツールと比較して、その性能を評価した。本ソフトウェアは、すべてのデータセットにおいて、最も多くの遺伝子を最大限のカバレッジで回収し、偽陽性の割合が低いことから評価すると、他のソフトウェアと同等以上に、十分にキュレートされたまたは未アノテーションの、低デプスまたは高デプスのショットガンおよびターゲットキャプチャアセンブリからオルソログを検索することができる。この基準の組み合わせを評価すると、ALiBaSeqは、テストしたすべてのタイプのデータで最も包括的で正確な系統樹アライメントを収集するツールとして、しばしば最高の評価を得ている。ソフトウェア（Pythonで実装）、チュートリアル、マニュアルは、https://github.com/AlexKnyshov/alibaseq で自由に利用できる。
Githubより
ソフトウェアの中核であるalibaseq.pyは、コンティグを含むFASTAファイル（例：NGSリードアセンブリファイル）から相同領域を検索するように設計されている。検索は、BLAST、HMMER、LASTZ検索のタブ区切り出力表やSAM / BAMアライメントファイルを読み込み、その結果をアセンブリファイルから検索することに基づいて行われる。系統推論用に遺伝子領域をコンパイルする（処理中の全分類群を遺伝子座ごとにグループ化し、このデータを与えられた遺伝子座ファイルに付加する）ように設計されているが、これは必須ではなく、別の出力構造を選択することができる。オプションとして、リバースサーチ（レシプロカルなベストヒットチェック）テーブルとリファレンスサーチ（ベイトが由来する分類群の完全アセンブリ/プロテオームに対して検索）テーブルを提供することができる。
tutorial
https://github.com/AlexKnyshov/alibaseq/tree/master/phylogenetic_tutorial
インストール
python2環境でテストした。
依存
- Python 2 or 3, Biopython, Pysam (only if SAM / BAM used as input)
mamba create -n alibaseq python=2.7 -y
conda activate alibaseq
mamba install -c bioconda -y biopython pysam blast
#本体
git clone https://github.com/AlexKnyshov/alibaseq.git
cd alibaseq/
> python alibaseq.py -h
usage: alibaseq.py [-h] -b table -f {S,M,SM} -x {n,s,a,b} [-t assembly]
[-q query] [-o output] [-s logsuffix]
[--om {query,target,combined}] [-e N] [-i N] [-B N] [-c N]
[--fl N] [--lr {none,actual,range}] [--is] [--hit-ovlp N]
[--ctg-ovlp N] [--recip-ovlp N]
[--bt {blast,hmmer22,hmmer18,hmmer15,lastz,sam,bam}]
[--btR {blast,bed}]
[--ac {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}]
[--acr {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}]
[--acR {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}]
[-r file/folder] [-R file]
[-m {e/b-i,e-b-i,b-e-i,i-b-e,i-e-b,b-i-e,e-i-b}]
[--rescale-metric] [--metric-merge-corr N] [--no-hs]
[--ref-hs] [--keep-strand] [--rm-rec-not-found]
[--hmmer-global] [--amalgamate-hits] [--max-gap N]
[--cname] [--both-strands {1,0}] [--srt N]
[--samScore metric] [--dd {all,none,random}] [--log-header]
[--synteny {1,0}] [-d namedelim]
ALiBaSeq (Alignment-Based Sequence extraction)
required arguments:
-b table alignment table file (default: None)
-f {S,M,SM} file / folder mode (default: None)
-x {n,s,a,b} extraction type: n (whole contig), s (only best hit
region), a (extract all hit regions and join them), b
(extract region between two outmost hit regions)
(default: None)
optional arguments:
-t assembly assembly file (default: None)
-q query query file(s) to which extracted results are to be
appended; if not specified, sequences are extracted
into blank files (default: None)
-o output output folder for modified files with extracted
sequences (default: alibaseq_out)
-s logsuffix output log suffix (default: default)
--om {query,target,combined}
output mode: group in files per bait [query], per
sample [target], or combine in a single file
[combined] (default: query)
-e N evalue cutoff (default: 0.01)
-i N identity cutoff (default: 0.0)
-B N bitscore cutoff (default: 0.0)
-c N number of contigs to extract; if set to 0, then
extract all contigs; if set to -1, then extract the
best and all close matches (default: 1)
--fl N flanks on each side in bp (default: 0)
--lr {none,actual,range}
local reciprocator setting (default: range)
--is perform contig stitching (default: False)
--hit-ovlp N allowed hit overlap on query, >= 1 in bp, or relative
0 < N < 1 (default: 0.1)
--ctg-ovlp N allowed contig overlap on query, >= 1 in bp, or
relative 0 < N < 1 (default: 0.2)
--recip-ovlp N contig overlap on query for reciprocator selection, >=
1 in bp, or relative 0 < N < 1 (default: 10)
--bt {blast,hmmer22,hmmer18,hmmer15,lastz,sam,bam}
alignment table type (default: blast)
--btR {blast,bed} reference alignment table type (default: blast)
--ac {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}
alignment coordinate type (default: dna-dna)
--acr {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}
reciprocal alignment coordinate type (default: dna-
dna)
--acR {dna-dna,tdna-aa,aa-tdna,aa-aa,tdna-tdna}
reference alignment coordinate type (default: dna-dna)
-r file/folder reciprocal search output file or folder (default:
None)
-R file bait to reference contig correspondence file (default:
None)
-m {e/b-i,e-b-i,b-e-i,i-b-e,i-e-b,b-i-e,e-i-b}
order of metrics to use to select best matches (e -
evalue, b - bitscore, i - identity) (default: e/b-i)
--rescale-metric divide metric value by length of hit region (default:
False)
--metric-merge-corr N
modify combined metric by this value (default: 1.0)
--no-hs do not run hit stitcher (default: False)
--ref-hs run hit stitcher on reciprocal table (slow) (default:
False)
--keep-strand keep original contig direction (default: False)
--rm-rec-not-found remove hits without matches in reciprocal search
(default: False)
--hmmer-global use HMMER contig score instead of domain score
(default: False)
--amalgamate-hits combine score for different hits of the same contig
(default: False)
--max-gap N max gap between HSP regions in either query or hit,
use 0 for no filtering (default: 0)
--cname append original contig name to output sequence name
(default: False)
--both-strands {1,0} allow both strands of the same contig region to be
considered (default: 1)
--srt N score ratio threshold, greater which the hits
considered be close matches (default: 0.9)
--samScore metric metric to use for scoring matches (default: MAPQ)
--dd {all,none,random}
in case hit matches several query with exactly equal
score, assign such hit to [all queries / none of the
queries / at random to only one] (default: none)
--log-header add a header to the table-like log files (default:
False)
--synteny {1,0} only stitch hits that are in synteny to query
(default: 1)
-d namedelim sample/contigID delimiter (default: |)
テストラン
１、リファレンスゲノムのBLASTnデータベースの作成。
cd alibaseq/phylogenetic_tutorial/
makeblastdb -in reference/reference_genome.fasta -dbtype nucl -parse_seqids
２、リファレンスゲノムに由来するベイトをリファレンスゲノムに対して検索。
mkdir refblast_results
bash ../blast_wrapper.sh baits/ reference/reference_genome.fasta 1e-10 blastn 4 n
mv *.blast refblast_results
３、BLASTするアセンブリのリストを作成（提供されたBLASTおよび相互検索ラッパースクリプトで使用するために必要）。
mkdir refblast_results
bash ../blast_wrapper.sh baits/ reference/reference_genome.fasta 1e-10 blastn 4 n
mv *.blast refblast_results
上のコマンドで指定したベイト配列
baits/
４、サンプルアセンブリの BLAST データベースの作成。
for f in assemblies/*.fasta
do
makeblastdb -in $f -dbtype nucl -parse_seqids
done
５、サンプルアセンブリのベイトを探す。
mkdir blast_results
bash ../blast_wrapper.sh baits/ assemblies/ 1e-10 dc-megablast 4 n assembly_list.txt
mv *.blast blast_results
６、baits-vs-sample 検索で見つかったコンティグをリファレンスゲノムに対して検索。
bash ../reciprocal_search.sh blast_results assemblies \
reference/reference_genome.fasta dc-megablast 4 n ../reciprocal_get_contigs.py \
assembly_list.txt
７、alibaseqを実行する（python3環境ならalibaseqPy3.pyを使用）。
python ../alibaseq.py -x a -f M -b blast_results -t assemblies -q baits -s tutorial_logs -o tutorial_results \
-e 1e-10 --is --amalgamate-hits -r blast_results -R refblast_results/reference_genome.fasta.blast
出力
tutorial_results/
引用
New alignment-based sequence extraction software (ALiBaSeq) and its utility for deep level phylogenetics
Alexander Knyshov, Eric R.L. Gordon, Christiane Weirauch
PeerJ. 2021 Mar 31;9:e11019
