2024/04/09 追記
ゲノムワイド関連研究（GWAS）は、ヒトの医学やゲノミクスにおいて不可欠なものとなっているが、細菌を対象とした研究はほとんど行われていない。本発表では、パンゲノムの構成要素について、観察された表現形質との関連を、集団の階層性を考慮しながら、進化過程の仮定を最小限にしてスコア化する、超高速で使いやすく、広く応用できるソフトウェアツール、Scoaryを紹介する。著者らは、この手法を従来の一塩基多型（SNP）ベースのGWASと区別するために、 pan-GWASと呼んでいる。ScoaryはPythonで実装され、オープンソースGPLv3ライセンスの下、https://github.com/AdmiralenOla/Scoary で入手可能である。
インストール
依存
- Python (Tested with versions 2.7, 3.4, 3.5, 3.6 and 3.6-dev)
- SciPy (Tested with versions 0.16, 0.17, 0.18)
If you supply custom trees (Optional)
Github
#PyPI(link)
pip install scoary==1.6.16
#option
pip install ete3 six
#conda (link)
mamba install -c bioconda scoary
> scoary -h
usage: scoary [-h] [-t TRAITS] [-g GENES] [-n NEWICKTREE] [-s START_COL]
[--delimiter DELIMITER] [-r RESTRICT_TO] [-o OUTDIR] [-u]
[-p P_VALUE_CUTOFF [P_VALUE_CUTOFF ...]]
[-c [{I,B,BH,PW,EPW,P} [{I,B,BH,PW,EPW,P} ...]]] [-m MAX_HITS]
[--include_input_columns GRABCOLS] [-w] [--no-time] [-e PERMUTE]
[--no_pairwise] [--collapse] [--threads THREADS] [--test]
[--citation] [--version]
Scoary version 1.6.16 - Screen pan-genome for trait-associated variants
optional arguments:
-h, --help show this help message and exit
Input options:
-t TRAITS, --traits TRAITS
Input trait table (comma-separated-values). Trait
presence is indicated by 1, trait absence by 0.
Assumes strain names in the first column and trait
names in the first row
-g GENES, --genes GENES
Input gene presence/absence table (comma-separated-
values) from ROARY. Strain names must be equal to
those in the trait table
-n NEWICKTREE, --newicktree NEWICKTREE
Supply a custom tree (Newick format) for phylogenetic
analyses instead instead of calculating it internally.
-s START_COL, --start_col START_COL
On which column in the gene presence/absence file do
individual strain info start. Default=15. (1-based
indexing)
--delimiter DELIMITER
The delimiter between cells in the gene
presence/absence and trait files, as well as the
output file.
-r RESTRICT_TO, --restrict_to RESTRICT_TO
Use if you only want to analyze a subset of your
strains. Scoary will read the provided comma-separated
table of strains and restrict analyzes to these.
Output options:
-o OUTDIR, --outdir OUTDIR
Directory to place output files. Default = .
-u, --upgma_tree This flag will cause Scoary to write the calculated
UPGMA tree to a newick file
-p P_VALUE_CUTOFF [P_VALUE_CUTOFF ...], --p_value_cutoff
P_VALUE_CUTOFF [P_VALUE_CUTOFF ...]
P-value cut-off / alpha level. For Fishers,
Bonferronis, and Benjamini-Hochbergs tests, SCOARY
will not report genes with higher p-values than this.
For empirical p-values, this is treated as an alpha
level instead. I.e. 0.02 will filter all genes except
the lower and upper percentile from this test. Run
with "-p 1.0" to report all genes. Accepts standard
form (e.g. 1E-8). Provide a single value (applied to
all) or exactly as many values as correction criteria
and in corresponding order. (See example under
correction). Default = 0.05
-c [{I,B,BH,PW,EPW,P} [{I,B,BH,PW,EPW,P} ...]], --correction
[{I,B,BH,PW,EPW,P} [{I,B,BH,PW,EPW,P} ...]]
Apply the indicated filtration measure. Allowed values
are I, B, BH, PW, EPW, P. I=Individual (naive)
p-value. B=Bonferroni adjusted p-value. BH=Benjamini-
Hochberg adjusted p. PW=Best (lowest) pairwise
comparison. EPW=Entire range of pairwise comparison
p-values. P=Empirical p-value from permutations. You
can enter as many correction criteria as you would
like. These will be associated with the
p_value_cutoffs you enter. For example "-c I EPW -p
0.1 0.05" will apply the following cutoffs: Naive
p-value must be lower than 0.1 AND the entire range of
pairwise comparison values are below 0.05 for this
gene. Note that the empirical p-values should be
interpreted at both tails. Therefore, running "-c P -p
0.05" will apply an alpha of 0.05 to the empirical
(permuted) p-values, i.e. it will filter everything
except the upper and lower 2.5 percent of the
distribution. Default = Individual p-value. (I)
-m MAX_HITS, --max_hits MAX_HITS
Maximum number of hits to report. SCOARY will only
report the top max_hits results per trait
--include_input_columns GRABCOLS
Grab columns from the input Roary file. and puts them
in the output. Handles comma and ranges, e.g.
--include_input_columns 4,6,8,16-23. The special
keyword ALL will include all relevant input columns in
the output
-w, --write_reduced Use with -r if you want Scoary to create a new gene
presence absence file from your reduced set of
isolates. Note: Columns 1-14 (No. sequences, Avg group
size nuc etc) in this file do not reflect the reduced
dataset. These are taken from the full dataset.
--no-time Output file in the form TRAIT.results.csv, instead of
TRAIT_TIMESTAMP.csv. When used with the -w argument
will output a reduced gene matrix in the form
gene_presence_absence_reduced.csv rather than
gene_presence_absence_reduced_TIMESTAMP.csv
Analysis options:
-e PERMUTE, --permute PERMUTE
Perform N number of permutations of the significant
results post-analysis. Each permutation will do a
label switching of the phenotype and a new p-value is
calculated according to this new dataset. After all N
permutations are completed, the results are ordered in
ascending order, and the percentile of the original
result in the permuted p-value distribution is
reported.
--no_pairwise Do not perform pairwise comparisons. Inthis mode,
Scoary will perform population structure-naive
calculations only. (Fishers test, ORs etc). Useful for
summary operations and exploring sets. (Genes unique
in groups, intersections etc) but not causal analyses.
--collapse Add this to collapse correlated genes (genes that have
identical distribution patterns in the sample) into
merged units.
Misc options:
--threads THREADS Number of threads to use. Default = 1
--test Run Scoary on the test set in exampledata, overriding
all other parameters.
--citation Show citation information, and exit.
--version Display Scoary version, and exit.
by Ola Brynildsrud (olbb@fhi.no)
実行方法
GUI上でも実行できるが（scoary_GUI）、ここではコマンドライン上での基本的な実行方法を説明する。
gene presence absence tableのCSVファイルと形質を1/0で示したCSVファイルを指定する。gene presence absence tableのフォーマットはroaryの出力もしくはLS-BSRの出力を想定している。他のパンゲノムツールの出力を使う場合、１列目にユニークな遺伝子名、15列目以降（excel等ではO列以降）に各ゲノムについて１列ずつ遺伝子の有無が書いてあれば、C-N列の情報は白紙でも受け付けることを確認した。
O列以降の１行目は各ゲノム名。
もう１つは表現型マトリックスのファイル。
gene presence absence tableのゲノム名と同じ名前で形質の有無を記載したCSVファイルを作成する。形質列は複数列あってもよい。
v1.6.9以降、欠損データを扱えるようになった。欠損値は "NA", ".", "-" のいずれかで指定する（Scoary はこれらの欠損値に対して実際には何らかの不確実性モデルを指定せず、単にさらなる分析から除外することに注意（マニュアルより））。
オプションでサンプルの系譜を定義する系統樹も使える。提供されない場合は、入力遺伝子型ファイルの分離ハミング距離を通して内部で計算される。
準備ができたらscoaryを実行する。macなどで作成した場合、改行形式に注意する（linux推奨）。
scoary -g gene_presence_absence.csv -t traits.csv
出力例。
traits ファイル内の trait 毎に 1 つの csv ファイルが出力される。有意性の順にソートされ、デフォルトではnaive p-value < 0.05 のすべての遺伝子が報告される（オプションで変更可）。
サンプルサイズ、オッズ比、native P値（Fisher’s exact test）、Bonferroni法での補正後の調整後P値(FWER)、Benjamini-Hochbergで補正後の調整後P値（FDR）、すべての検定推定量をランク付けした後の経験的p値（Empirical_p）などが表示される。
バージョン 1.6.12 から、scoary は VCF ファイルを Roary/Scoary 形式に変換する機能を備えている。これにより、バリアント(SNPs, indel, structural variants など)を入力に使用できる。
vcf2scoary myvariants.vcf
(すべてのVCFファイルを正しく扱えるとは限らないため、バグがある場合は報告すること)
論文の表２にサンプルサイズを変えて検出力がnative P値、FDR、経験的p値でどう変化するか示した結果があります。また、Githubには他の色々な使用例が書かれています。確認してください。
追記
- macでエラーが出る時は入力CSVの改行コードがLFになっているか確認する。
引用
Rapid scoring of genes in microbial pan-genome-wide association studies with Scoary
Ola Brynildsrud, Jon Bohlin, Lonneke Scheffer, Vegard Eldholm
Genome Biol. 2016 Nov 25;17(1):238
関連
