macでインフォマティクス

macでインフォマティクス

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

微生物パンゲノム解析のスコア付けを行う Scoary

 

 ゲノムワイド関連研究(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)

  • ete3
  • six

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の出力を想定している。他のパンゲノムツールの出力を使う場合、1列目にユニークな遺伝子名、15列目以降(excel等ではO列以降)に各ゲノムについて1列ずつ遺伝子の有無が書いてあれば、C-N列の情報は白紙でも受け付けることを確認した。

O列以降の1行目は各ゲノム名。

もう1つは表現型マトリックスのファイル。

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ファイルを正しく扱えるとは限らないため、バグがある場合は報告すること)

 

 

 

論文の表2にサンプルサイズを変えて検出力がnative P値、FDR、経験的p値でどう変化するか示した結果があります。また、Githubには他の色々な使用例が書かれています。確認してください。

引用

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

 

関連