stoat test - Pa-Tou/stoat GitHub Wiki

This subcommand tests the association between genotypes retrieved by stoat graph or stoat vcf and a phenotype.

Input files

Snarl genotypes

A snarl_genotypes.tsv.gz file produced by stoat graph or stoat vcf

Phenotype

A TSV file with the phenotype for all samples. It must have two columns names SAMPLE and PHENO. For binary phenotypes, the PHENO column must have 0 or 1. Quantitative phenotypes can be provided too.

A binary phenotype file might look like:

SAMPLE	PHENO
samp1	0
samp2	1
...

A quantitative phenotype file might look like:

SAMPLE	PHENO
samp1	-0.6
samp2	3.8
...

Gene expression for an eQTL analysis

For an expression QTL study, the phenotype file is a table with

  • a first column gene_name
  • one column per sample with the gene expression in that sample

For example:

gene_name	samp1	samp2
geneA	12.708787454384092	7.945852683899044
geneB	11.130996245885658	6.959378045490106

A TSV with the gene position is also required. It must have the following columns: gene_name, chr, start, end. For example:

gene_name	chr	start	end
geneA	chr1	0	10000
geneB	chr2	100	10100

Covariates (optional)

Covariates can also be provided and will be included in the regression models.

The input TSV file must have a SAMPLE column. Other columns must match the names specified with -C/--covar-name.

For example:

SAMPLE	SEX	PC1	PC2	PC32
samp1	1	1.2	562.25	42.25
samp2	0	5.2	359.25	65.24

Usage

The type of test is specified with the -m/--method parameter. It can be:

  • exact: a naive selection of the snarls where samples perfectly segregate by binary phenotype.
  • chi2: chi-squared and Fisher's exact test, for binary phenotypes without covariates.
  • logreg: logistic regression, for binary phenotypes with or without covariates. (default)
  • linreg: linear regression for quantitative phenotypes with or without covariates.

A typical command would look like:

# chi-squared test on a binary phenotype
stoat test -g snarl_genotypes.tsv.gz -p case_control.tsv -m chi2 -o stoat_chi2

# logistic regression on a binary phenotype
stoat test -g snarl_genotypes.tsv.gz -p case_control.tsv -m logreg -o stoat_logreg

# linear regression on a quantitative phenotype
stoat test -g snarl_genotypes.tsv.gz -p quantitative_peno.tsv -m linreg -o stoat_quant

In those examples, we also specify a different output directory with -o.

Including covariates

For the linreg and logreg regression-based methods, we can include covariates in the model.

stoat test -g snarl_genotypes.tsv.gz -p phenotype.tsv -m linreg -c covariates.tsv -C SEX,PC1

To only include SEX and PC1 taken from the covariates.tsv file.

Expression QTL study

Including the gene position file (-P) will trigger the eQTL mode.

stoat test -g snarl_genotypes.tsv.gz -p gene_expression.tsv -P gene_positions.tsv

The -w arguments controls how far from the gene should variants be tested. The default is 1 Mb (i.e. -w 1000000).

We use the linear regression model to test association between snarls and gene expression. Hence, we can also include covariates, for example:

stoat test -g snarl_genotypes.tsv.gz -p gene_expression.tsv -P gene_positions.tsv -c covariates.tsv -C SEX,PC1,PC2,PC3

Variant filters

Snarls where we expect low power can be excluded from the tests using two filters:

  • -M/--maf the minimum allele frequency. In the case of multiallelic snarls, we use the second most frequent allele. Hence, it means at least two alleles have frequency above the specified threshold. The default is 0.05.
  • -I/--min-individuals the minimum number of individuals with at least one allele traversing that snarl. The default is 0 but it is usually recommended to set to it to ~10 to ensure stable tests (especially for regression models).

stoat test full list of arguments

Stoat v0.0.4
Usage: stoat test [options]
Output: stoat.assoc.pvalues.tsv.gz written to the specified output directory

  -g, --genotype FILE             Path to the genotype file from stoat graph or stoat vcf
  -m, --method STR                Which test method to use: chi2 (Fisher/Chi-Squared), linreg (linear regression), 
                                                            logreg (logistic regression), exact (exact phenotype/genotype match)
  -p, --phenotype FILE            Path to the phenotype file
  -P, --gene-position FILE        Path to the gene position file (activates the eQTL testing mode)
  -w, --max-gene-distance INT     Include snarls up to this distance from the gene when looking for eQTLs [1000000]
  -c, --covariate FILE            Path to the covariate file
  -C, --covar-name NAME           Covariate column name(s) used (default: use all covariates in file)
  -I, --min-individuals INT       Minimum number of individuals with at least one allele in the snarl [2]
  -M, --maf FLOAT                 Minimum allele frequency threshold [0.05]
  -t, --threads INT               Number of threads to use [1]
  -V, --verbose INT               Verbosity level (0=error, 1=warn, 2=info, 3=debug, 4=trace) [2]
  -o, --output FILE               Output directory name [stoat_output]
  -u, --no-bgzip                  Don't compress the output file with bgzip
  -a, --ascii                     Print the STOAT ascii art banner
  -h, --help                      Print this help message

Output

Files are written in the output directory specified with -o. A log will also be saved in this directory. This directory will be created if necessary.

The main output file is a TSV with the pvalues of each snarl in stoat.assoc.pvalues.tsv.gz. The format is described here.

Visualization

Work in progress.

The StoatPlot R package is currently being developed to help explore STOAT's results. It could be used to prepare Manhattan plots, QQ Plots, or to visualize specific snarls.

Experimental sequenceTubeMap visualization

We are also experimenting with outputting paths in the pangenome to highlight associated snarls. This experimental option is not available yet. The idea is to generate a GAF file in the output directory, which can be used with the sequenceTubeMap to visualize the associated path in the pangenome.

Legend: Colors represent the different paths group (red: group 1 & blue: group 0) and opacity represente the number of samples in that paths (number of samples passing trought each paths % 60).