共分散に基づく共選択的圧力下での多型の発見やエピスタシスは、近年、集団ゲノム科学において非常に注目されている。染色体上の対立遺伝子の集団レベルでの共分散を統計的にモデル化し、多型のペア間の依存性をモデルフリーで検定することで、細菌集団における選択のパターンを見出すことに成功したことが示されている。ここでは、計算効率が高く、多くの細菌のパンゲノムのスケールで解析が可能なモデルフリー手法SpydrPickを紹介する。SpydrPickは集団構造に対する効率的な補正を組み込んでおり、明示的な系統樹を必要とせずにデータ中の系統的なシグナルを調整することができる。また、ゲノムワイド関連研究で用いられるマンハッタンプロットに似た、新しいタイプの結果の可視化を導入し、同定された共進化のシグナルを迅速に探索することができるようになった。シミュレーションにより、本手法の有用性を示すとともに、この種の解析が最も成功しやすい時期についての知見を得ることができる。この方法を2つの主要なヒト病原体、Streptococcus pneumoniaeとNeisseria meningitidisの大規模な集団ゲノムデータセットに適用したところ、病原性と抗生物質耐性に関する共選択の既知ターゲットと新規ターゲットの両方が発見された。
（中略）
集団の配列データから共進化シグナルを検出するための比較手法は、ここ数十年の間に多くの注目を集めている。その顕著な例として、大規模なタンパク質アラインメントにおける非隣接部位間の共分散の統計的解析が、タンパク質の3次元構造における部位間の接触を予測するのに有効であることが証明されている。タンパク質構造中の接触部位は共通の構造的制約の下で共進化するため、タンパク質アラインメント中に検出可能な相関の痕跡を与える。同様に、共通の選択圧のもとで進化した部位は、適切な表現型データがない場合でも、配列アラインメントから検出可能な共選択パターンを生じさせる可能性がある。そのため、近年、細菌集団のゲノムワイド塩基配列の探索的共変化解析が注目されている。この解析の目的は、共通の選択圧のもとで共進化し、エピスタティック相互作用に関与している可能性のある部位を発見することである。
インストール
#conda (link)
mamba create -n spydrpick -y
conda activate spydrpick
mamba install -c bioconda spydrpick -y
> SpydrPick -h
$ SpydrPick -h
SpydrPick: MI-ARACNE genome-wide co-evolution analysis.
Usage: SpydrPick [options] <alignmentfile> [-o <outputfile>]
Option '--help' will print a list of available options.
-h [ --help ] Print this help message.
--version Print version information.
-v [ --verbose ] Be verbose.
--mi-threshold arg (=-1) The MI threshold value. Experience
suggests that a value of 0.11 is often
reasonable. Zero indicates no threshold
and negative values will trigger
auto-define heuristics.
--mi-values arg (=0) Approximate number of MI values to
calculate from data (default=#samples*1
00).
--mi-pseudocount arg (=0.5) The MI pseudocount value.
--mi-threshold-iterations arg (=10) Number of iterations for estimating
saving threshold.
--mi-threshold-pairs arg (=0) Number of sampled pairs for estimating
saving threshold (0=determine
automatically).
--ld-threshold arg (=0) Threshold distance for linkage
disequilibrium (LD).
--no-aracne Skip ARACNE, only calculate MI.
-t [ --threads ] arg (=-1) Number of threads per MPI/shared memory
node (-1=use all hardware threads that
the OS/environment exposes).
--alignmentfile arg The input alignment filename(s). When
two filenames are specified, only
inter-alignment links will be probed
for.
--include-list arg Name of file containing list of loci to
include in analysis.
--exclude-list arg Name of file containing list of loci to
exclude from analysis.
--mapping-list arg Name of file containing loci mappings.
--sample-list arg The sample filter list input filename.
--sample-weights arg Name of file containing sample weights.
--input-indexing-base arg (=1) Base index for input.
--no-filter-alignment Do not filter alignment columns by
applying MAF and GAP thresholds.
--maf-threshold arg (=0.01) Minor state frequency threshold. Loci
with less than 2 states above threshold
are removed from alignment.
--gap-threshold arg (=0.14999999999999999)
Gap frequency threshold. Positions with
a gap frequency above the threshold are
excluded from the pair-analysis.
--no-sample-reweighting Turn sample reweighting off, i.e. do
not try to correct for population
structure.
--sample-reweighting-threshold arg (=0.90000000000000002)
Fraction of identical positions
required for two samples to be
considered identical.
--rescale-sample-weights Rescale sample weights so that
n(effective) == n.
--output-state-frequencies Write column-wise state frequencies to
file.
--output-sample-weights Output sample weights.
--output-sample-distance-matrix Output triangular sample-sample Hamming
distance matrix.
--output-indexing-base arg (=1) Base index for output.
--output-alignment Write alignment to file.
--output-filtered-alignment Write filtered alignment to file.
--genome-size arg Genome size, if different from input.
Default = 0: detect size from input.
--linear-genome Assume linear genome in distance
calculations.
--begin arg (=1) The first locus index to work on. Used
to define a range.
--end arg (=-1) The last locus index to work on (-1=end
of input). Used to define a range.
-o [ --aracne-outputfile ] arg (=aracne.out)
The ARACNE output filename. This is a
binary file for "plot.r".
--aracne-edge-threshold arg (=2.2204460492503131e-16)
Equality tolerance threshold. Edges
differing by less than this value are
considered equal in strength.
--aracne-block-size arg (=16384) Block size for graph processing.
--aracne-node-grouping-size arg (=16) Grouping size for node processing.
実行方法
FASTA形式の多重整列ファイルを指定する。A、C、G、Tは異なるカテゴリとしてマップさが、その他の文字は1つのギャップカテゴリにマップされる。大文字と小文字は区別されない。
SpydrPick -v core_gene_alignment.aln
- -v Be verbose.
- --ld-threshold <arg> (=0) Threshold distance for linkage disequilibrium (LD).
- --linear-genome Assume linear genome in distance calculations.
出力例
Githubより
- SpydrPickのメイン出力ファイル *.spydrpick_couplings.?-based.*edges は、入力アラインメントの列（左から右）に従って番号付けされたMI値（Mutual information value; 論文で定義された2つの確率変数の相互依存性を表す情報理論的な尺度）と位置インデックスのペア（デフォルトでは1ベースのインデックス使用、--output-indexing-baseでコントロール）の降順リストで空白区切りのリストが含まれている。MI値が外れ値閾値以上で --ld-閾値より離れている位置ペアは、いくつかの追加データと共に *.outliers というファイルでも見つけることができる。
- デフォルトの設定は全体的に合理的な結果を得られるように設計されている。
- 注目すべき例外は、連鎖不平衡（LD）の閾値距離オプション --ld-threshold=
Rスクリプトgwes_plot.rを使うと、結果のedgesファイルから論文のFIg.3AやFIg.4のようなGWES (genome-wide epistasis and co-selection study) Manhattan plotをプロットできる。使用するにはRのコード；gwes_plot.rの10行目input_full_filepathで出力されたedgeファイルをフルパスで指定する。
git clone https://github.com/santeripuranen/SpydrPick.git
cd SpydrPick/
Rscript gwes_plot.r
引用
Genome-wide epistasis and co-selection study using mutual information
Johan Pensar, Santeri Puranen, Brian Arnold, Neil MacAlasdair, Juri Kuronen, Gerry Tonkin-Hill, Maiju Pesonen, Yingying Xu, Aleksi Sipola, Leonor Sánchez-Busó, John A Lees, Claire Chewapreecha, Stephen D Bentley, Simon R Harris, Julian Parkhill, Nicholas J Croucher, Jukka Corander
Nucleic Acids Research, Volume 47, Issue 18, 10 October 2019, Page e112
