Skip to content

stoat test

Jump to bottom
Xian Chang edited this page May 12, 2026 · 3 revisions

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

Usage: stoat test [options]

  -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
  -I, --min-individuals INT       Minimum number of individuals per snarl [0]
  -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
  -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.

It contains the following columns are common to all tests:

  • CHR: chromosome name
  • START_POS/END_POS: start/end position of the snarl. If it's a nested snarl not on the reference, this will be the coordinates of the first parent snarl on the reference.
  • SNARL: identifier for the snarl
  • PATH_LENGTHS: comma separated list of the lengths of each allele. For complex variants (which contain nested variants), there could be multiple possible lengths so they are represented as a range (min/max). A SNP will appear as 1,1[,1...]. A simple INDEL will have one zero-length allele and the other non-zero.
  • DEPTH: how deeply nested is the variant? 1 for top-level variant, 2 for nested within a top-level variant, etc.

When -m chi2 is used (binary phenotype, no covariates), the output will include:

  • P_FISHER/P_CHI2: P-values calculated using Fisher's exact test (if two alleles) or the Chi-squared test.
  • ALLELE_COUNT_PER_PHENO encodes the allele distribution across the two phenotype groups for each path in a snarl. Each entry is formatted as X:Y,X':Y',... where X:Y represents one path (allele), with X being the count of samples in group 0 and Y in group 1. Commas separate multiple paths (alleles) within the same snarl.

Otherwise, the output will include:

  • P: P-values calculated using the regression model.
  • ALLELE_COUNT: for each allele, how many samples take this allele. Comma separated list.

For an eQTL study, an additional GENE column will be present.

Example of an output file

Example of the output for a binary phenotype (provided with -b):

#CHR	START_POS	END_POS	SNARL	PATH_LENGTHS	P_FISHER	P_CHI2	GROUP_PATHS	DEPTH
ref	247411	247411	4271_4260	0,616	0.368	0.3172	107:97,93:103	1
ref	247411	247411	4270_4267	1,1	0.1532	0.1447	14:24,79:79	2

Example of the output for a quantitative phenotype (-q option):

#CHR	START_POS	END_POS	SNARL	PATH_LENGTHS	P	RSQUARE	ALLELE_PATHS	DEPTH
ref	244656	244796	4238_4236	0,140	0.5142	2.1519e-03	22,378	1
ref	244553	244553	4236_4234	0,182	0.7975	3.3338e-04	348,52	1

Example of the output for an eQTL analysis (-e option):

#CHR	START_POS	END_POS	SNARL	PATH_LENGTHS	GENE	P	RSQUARE	ALLELE_PATHS	DEPTH
1	100000	100001	2_5	1,1	gene_80	0.2039	8.1409e-03	201,199	1

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).

Clone this wiki locally