A method for allele specific analyses across cell states and conditions. stratAS implements population-scale tests and predictive models for allele-specific activity (typically RNA or ATAC data).

• A .vcf file containing phased individual genotypes and an "AS" field listing the allelic counts. This gets turned into a flat file to be loaded in the --input flag (see below).
• A samples file with the columns ID and CONDITION that is in the same order as the vcf input. CONDITION is currently restricted to 0/1. This file goes into the --samples flag.
• A .bed file listing the physical positions of the features to be tested for allelic imbalance. This file goes into the --peaks flag.
• Global/local parameter files generated by params.R (see below). At minimum a global overdispersion parameter file with the columns ID and PHI with entries in any order goes into the --global_params flag.
• R, and the VGAM, optparse libraries installed.

• Phase and impute your genotype data with EAGLE on the Haplotype Reference Consortium imputation server.
• Process all sequence data with the WASP mapping pipeline.
• Use the GATK ASEReadCounter to count reads then convert into vcf (see pipeline/ scripts).

params.R infers the read count distribution (beta binomial) parameters and needs to be run on a whole genome vcf one individual at a time. The input is an --inp_counts file which must contain the headers CHR POS HAP REF.READS ALT.READS where HAP is the 0|1 or 1|0 vcf haplotype code.

A vcf file containing one individual is converted to counts as follows:

zcat $VCF \
| grep -v '#' \
| cut -f 1,2,10 \
| tr ':' '\t' \
| awk '$3 != "1|1" && $3 != "0|0"' \
| cut -f 1-3,7 \
| tr ',' '\t' \
| sed 's/chr//' \
| awk 'BEGIN { print "CHR POS HAP REF.READS ALT.READS" } $4 + $5 > 0' \
 > $OUT.counts

For tumor data, local CNV-specific parameters are also estimated, and the --inp_cnv file must be provided with headers CHR P0 P1 listing the boundaries of CNV regions. This file can additionally contain a CNV column listing the CNV estimate (centered to zero as in TCGA calls) for inclusion as a covariate in the final analysis.

inference is then performed by running:

Rscript params.R --min_cov 5 --inp_counts $OUT.counts --inp_cnv $CNV --out $OUT

A $OUT.global.params file is generated containing the parameters (with header PHI MU N) and (optionally) an $OUT.local.params file is generated containing the positions and parameters for each CNV (with header CHR P0 P1 PHI MU N).

stratas.R computes the AS statistics from a VCF of all individuals and the prior parameters estimated above.

A vcf file containing all individuals is converted to counts and split into batches as follows:

zcat $VCF\
| grep -v '#' \
| cut -f 1-5,9- \
| tr ':' '\t' \
| awk '{ for(i=1;i<=5;i++) { printf "%s ",$i; } for(i=6;i<=10;i++) { if($i == "GT") gtnum=i-6; if($i=="AS") asnum=i-6; } for(i=11;i<=NF;i++) { if( (i-11) % 5 == asnum || (i-11) % 5 == gtnum ) printf " %s",$i; } print ""; }' \
| sed 's/[|,]/ /g' | split -d -l 15000 - $OUT.MAT.split.

Each batch is then processed as follows:

Rscript stratas.R \
--input $OUT.MAT.split.00 \
--samples SAMPLES.ID \
--peaks gencode.protein_coding.transcripts.bed \
--global_param GLOBAL.params \
--local_param LOCAL.params \

stratas.R can compute multi-variate predictive models for downstream integration with GWAS data, using both allelic and total activity.

Inputs are similar to the tests for allele specificity described above and additionally require a total activity matrix, covariates, and the predict flags:

Rscript stratas.R \
--input $OUT.MAT.split.00 \
--samples SAMPLES.ID \
--peaks gencode.protein_coding.transcripts.bed \
--global_param GLOBAL.params \
--local_param LOCAL.params \
--total_matrix KIRC.gexp \
--covar KIRC.gexp.cov \
--predict_snps HM3.extract \
--predict --predict_only > $OUT.profile

This generates predictive model files for every feature in the TWAS/FUSION format, which can then be integrated with GWAS data using the FUSION software. The generated model types inside each model file are "lasso", "lasso.as", "lasso.combined", "top1.as", "top1", "top1.combined"; where lasso versus top corresponds to full locus penalized versus top SNP models; and *, *.as, *.combined indicates models using total activity, allele-specific activity, or both jointly. An $OUT.profile file is generated which lists model performance characteristics.

Important: This analysis differs from the individual allele-specific tests in two ways: (1) the CONDITION value is not used and one model is trained on all samples; (2) models are still trained even if no allele-specific informatino is available. To restrict to models with both sources of signal (which are likely to be more accurate as a group), run the allele-specific tests described above and retain only peaks that have non-NA values.

Example total activity input files (KIRC.gexp, KIRC.gexp.cov, and HM3.extract) based on TCGA data are provided in the example/ directory.

By default, the outputs are printed to screen with each line containing the following entries:

Enabling the --total_matrix flag additionally runs a standard eQTL test (and interaction term), producing the following columns:

Enabling the --combine flag additionally estimates the (inverse-variance weighted) combined statistics from the allele-specific and total models, producing the following columns:

Enabling the --binom flag additionally runs a standard binomial test (and Fisher's test for interaction), producing the following columns:

Enabling the --bbreg flag additionally runs a beta binomial regression with CNV as covariate, and produces the following columns:

Enabling the --indiv flag additionally produces the following columns:

An example locus with significant AS associations can be run by calling:

Rscript stratas.R \
--input example/ENSG00000075240.12.mat \
--samples example/KIRC.ALL.AS.PHE \
--peaks example/ENSG00000075240.12.bed \
--global_param example/KIRC.ALL.AS.CNV \
--local_param=example/KIRC.ALL.AS.CNVLOCAL

All parameters can be displayed by running --help.