Ubuntu Manpage: QTLtools cis - cis QTL analysis

Provided by: qtltools_1.3.1+dfsg-2build2_amd64

NAME

       QTLtools cis - cis QTL analysis

SYNOPSIS

       QTLtools   cis   --vcf  [in.vcf|in.vcf.gz|in.bcf|in.bed.gz]  --bed  quantifications.bed.gz
       [--nominal float | --permute integer | --mapping in.txt] --out output.txt [OPTIONS]

DESCRIPTION

       This mode maps cis (proximal) quantitative trait loci (QTLs) that  affect  the  phenotype,
       using       linear       regression.        The       method      is      detailed      in
       <https://www.nature.com/articles/ncomms15452>.   We  first  regress   out   the   provided
       covariates  from the phenotype data, followed by running the linear regression between the
       phenotype residuals and the genotype.  If --normal and --cov  are  provided  at  the  same
       time,  then  the residuals after the covariate correction are rank normal transformed.  It
       incorporates an efficient permutation scheme to control for differential multiple  testing
       burden  of  each  phenotype.  You can run a nominal pass (--nominal threshold) listing all
       genotype-phenotype associations below a certain threshold, a permutation  pass  (--permute
       no_of_permutations)  to empirically characterize the null distribution of associations for
       each phenotype separately, thus adjusting the nominal p-value of the best association  for
       a  phenotype,  or  a  conditional  analysis pass (--mapping filename) to discover multiple
       proximal QTLs with independent effects on a phenotype.

       As multiple molecular phenotypes can belong to  higher  order  biological  entities,  e.g.
       exons  of genes, QTLtools cis allows grouping of phenotypes to maximize the discoveries in
       such particular cases.  Specifically, QTLtools can either aggregate multiple phenotypes in
       a given group into a single phenotype via PCA (--grp-pca1) or by taking their mean (--grp-
       mean), or directly use all individual phenotypes in an extended  permutation  scheme  that
       accounts  for  their  number  and  correlation structure (--grp-best).  In our experience,
       --grp-best outperforms the other options for expression QTLs (eQTLs).

       The conditional analysis  pass  first  uses  permutations  to  derive  a  nominal  p-value
       threshold  per  phenotype  that  varies  and  reflects the number of independent tests per
       cis-window.  Then, it uses a forward-backward stepwise regression to learn the  number  of
       independent  signals  per  phenotype,  determine the best candidate variant per signal and
       assign all significant hits to the independent signal they relate to.

       Since linear regressions assumes normally distributed data, we highly recommend using  the
       --normal  option  to rank normal transform the phenotype quantifications in order to avoid
       false positive associations due to outliers.

OPTIONS

--vcf [in.vcf|in.bcf|in.vcf.gz|in.bed.gz]
Genotypes in VCF/BCF format, or another molecular phenotype in BED format. If
there is a DS field in the genotype FORMAT of a variant (dosage of the genotype
calculated from genotype probabilities, e.g. after imputation), then this is used
as the genotype. If there is only the GT field in the genotype FORMAT then this is
used and it is converted to a dosage. REQUIRED.

--bed quantifications.bed.gz
Molecular phenotype quantifications in BED format. REQUIRED.

--out output.txt
Output file. REQUIRED.

--cov covariates.txt
Covariates to correct the phenotype data with.

--normal
Rank normal transform the phenotype data so that each phenotype is normally
distributed. RECOMMENDED.

--std-err
Calculate and output the standard error of the beta (slope).

--window integer
Size of the cis window flanking each phenotype's start position. DEFAULT=1000000.
RECOMMENDED=1000000.

--nominal float
Calculate the nominal p-value for the genotype-phenotype associations and print out
the ones that pass the provided threshold. Give 1.0 to print everything. Mutually
exclusive with --permute and --mapping.

--permute integer
Adjust the best nominal p-value for this phenotype accounting for the number of
variants and the linkage disequilibrium in its cis-window. We recommend at least
1000 permutation for the final analysis, and in most cases you will see diminishing
returns when going over 5000. However, if you are doing exploratory analyses like
which/how many covariates to include, you can go as low as 100. Mutually exclusive
with --nominal and --mapping. RECOMMENDED=1000.

--mapping thresholds_filename
The conditional analysis. First you need to run a permutation analysis, then
generate a thresholds file using the runFDR_cis.R script in the script directory.
Mutually exclusive with --nominal and --permute.

--grp-best
Correct for multiple phenotypes within a phenotype group. Mutually exclusive with
--grp-pca1 and --grp-mean.

--grp-pca1
Run PCA on phenotypes within a phenotype group and use PC1 for association testing.
Mutually exclusive with --grp-best and --grp-mean.

--grp-mean
Average phenotypes within a group and use the results for association testing.
Mutually exclusive with --grp-best and --grp-pca1.

--chunk integer1 integer2
For parallelization. Divide the data into integer2 number of chunks and process
chunk number integer1. Chunk 0 will print a header. Mutually exclusive with
--region and --region-pair. Minimum number of chunks has to be at least the same
number of chromosomes in the --bed file.

--region chr:start-end
Genomic region to be processed. E.g. chr4:12334456-16334456, or chr5 Mutually
exclusive with --chunk and --region-pair.

--region-pair chr:start-end chr:start-end
Genomic region for genotypes followed by region for phenotypes to be processed.
Mutually exclusive with --chunk and --region.

OUTPUT FILES

--nominal output file
Space separated output file with the following columns (certain columns are only printed
based on options). We recommend including chunk 0 to print out a header in order to
avoid confusion.

5.2 n_phe_in_grp Only printed if --group-pca1 | --group-mean.
The number of phenotypes in the phenotype
group.
6 n_var_in_cis The number variants in the cis window for this
phenotype.
7 dist_phe_var The distance between the variant and the
phenotype start positions.
8 var_id The variant ID.
9 var_chr The variant chromosome.
10 var_from The start position of the variant.
11 var_to The end position of the variant.
12 nom_pval The nominal p-value of the association between
the variant and the phenotype.
13 r_squared The r squared of the linear regression.
14 slope The beta (slope) of the linear regression.
14.1 slope_se The standard error of the beta. Only printed
if --std-err is provided.
15 best_hit Whether this varint was the best hit for this
phenotype.

--permute output file
Space separated output file with the following columns (certain columns are only printed
based on options). We recommend including chunk 0 to print out a header in order to
avoid confusion.

1 phe_id | grp_id The phenotype ID or if one of the grouping
options is provided, then phenotype group ID
2 phe_chr The phenotype chromosome
3 phe_from Start position of the phenotype
4 phe_to End position of the phenotype
5 phe_strd The phenotype strand
5.1 phe_id | ve_by_pc1 | n_phe_in_grp Only printed if --group-best | --group-pca1 |
--group-mean. The phenotype ID, variance
explained by PC1, or number of phenotypes in
the phenotype group for --group-best, --group-
pca1, and --group-mean, respectively.
5.2 n_phe_in_grp Only printed if --group-pca1 | --group-mean.
The number of phenotypes in the phenotype
group.
6 n_var_in_cis The number variants in the cis window for this
phenotype.
7 dist_phe_var The distance between the variant and the
phenotype start positions.
8 var_id The most significant variant ID.
9 var_chr The most significant variant's chromosome.
10 var_from The start position of the most significant
variant.
11 var_to The end position of the most significant
variant.
12 dof1 The number of degrees of freedom used to
compute the p-values.
13 dof2 Estimated number of degrees of freedom used in
beta approximation p-value calculations.
14 bml1 The first shape parameter of the fitted beta
distribution (alpha parameter). These should
be close to 1.
15 bml2 The second shape parameter of the fitted beta
distribution (beta parameter). This
corresponds to the effective number of
independent tests in the region.
16 nom_pval The nominal p-value of the association between
the most significant variant and the
phenotype.
17 r_squared The r squared of the linear regression.
18 slope The beta (slope) of the linear regression.

18.1 slope_se The standard error of the beta. Only printed
if --std-err is provided.
19 adj_emp_pval Adjusted empirical p-value from permutations.
This is the adjusted p-value not using the
beta approximation. Simply calculated as:
(number of p-values observed during
permutations that were smaller than or equal
to the nominal p-value + 1) / (number of
permutations + 1). The most significant p-
value achievable would be 1 / (number of
permutations + 1).
20 adj_beta_pval Adjusted empirical p-value given by the fitted
beta distribution. We strongly recommend
using this adjusted p-value in any downstream
analysis.

--mapping output file
Space separated output file with the following columns (certain columns are only printed
based on options). We recommend including chunk 0 to print out a header in order to
avoid confusion.

1 phe_id | grp_id The phenotype ID or if one of the grouping
options is provided, then phenotype group ID
2 phe_chr The phenotype chromosome
3 phe_from Start position of the phenotype
4 phe_to End position of the phenotype
5 phe_strd The phenotype strand
5.1 phe_id | ve_by_pc1 | n_phe_in_grp Only printed if --group-best | --group-pca1 |
--group-mean. The phenotype ID, variance
explained by PC1, or number of phenotypes in
the phenotype group for --group-best, --group-
pca1, and --group-mean, respectively.
5.2 n_phe_in_grp Only printed if --group-pca1 | --group-mean.
The number of phenotypes in the phenotype
group.
6 n_var_in_cis The number variants in the cis window for this
phenotype.
7 dist_phe_var The distance between the variant and the
phenotype start positions.
8 var_id The most significant variant ID.
9 var_chr The most significant variant's chromosome.
10 var_from The start position of the most significant
variant.
11 var_to The end position of the most significant
variant.
12 rank The rank of the association. This tells you
if the variant has been mapped as belonging to
the best signal (rank=0), the second best
(rank=1), etc ... As a consequence, the
maximum rank value for a given phenotype tells
you how many independent signals there are
(e.g. rank=2 means 3 independent signals).
13 fwd_pval The nominal forward p-value of the association
between the most significant variant and the
phenotype.
14 fwd_r_squared The r squared of the forward linear
regression.
15 fwd_slope The beta (slope) of the forward linear
regression.
15.1 fwd_slope_se The standard error of the forward beta. Only
printed if --std-err is provided.
16 fwd_best_hit Whether or not this variant was the forward
most significant variant.
17 fwd_sig Whether this variant was significant.
Currently all variants are significant so this
is redundant.

18 bwd_pval The nominal backward p-value of the
association between the most significant
variant and the phenotype.
19 bwd_r_squared The r squared of the backward linear
regression.
20 bwd_slope The beta (slope) of the backward linear
regression.
20.1 bwd_slope_se The standard error of the backward beta. Only
printed if --std-err is provided.
21 bwd_best_hit Whether or not this variant was the backward
most significant variant.
22 bwd_sig Whether this variant was significant.
Currently all variants are significant so this
is redundant.

NOMINAL ANALYSIS EXAMPLES

       o Nominal  pass,  correcting  for  technical  covariates,  rank  normal  transforming  the
         phenotype, and printing out associations with a p-value <= 0.01 on chromosome 22 between
         17000000 and 18000000 bp, and excluding some samples (see QTLtools (1)):

         QTLtools  cis  --vcf  genotypes.chr22.vcf.gz  --bed  genes.50percent.chr22.bed.gz  --cov
         genes.covariates.pc50.txt.gz  --nominal  0.01  --region chr22:17000000-18000000 --normal
         --out nominals.txt --exclude-samples sample_names_to_exclude.txt

       o Nominal pass with parallelization  correcting  for  technical  covariates,  rank  normal
         transforming  the  phenotype,  and printing out associations with a p-value <= 0.01.  To
         facilitate parallelization on compute  cluster,  we  developed  an  option  to  run  the
         analysis into chunks of molecular phenotypes.  For instance, to run analysis on chunk 12
         when splitting the example data set into 20 chunks, run:

         QTLtools  cis  --vcf  genotypes.chr22.vcf.gz  --bed  genes.50percent.chr22.bed.gz  --cov
         genes.covariates.pc50.txt.gz    --nominal    0.01   --chunk   12   20   --normal   --out
         nominals_12_20.txt

       o If you want to submit the whole analysis with 20 jobs on a  compute  cluster,  just  run
         (qsub needs to be changed to the job submission system used [bsub, psub, etc...]):

         for j in $(seq 0 20); do
             echo  "QTLtools  cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz
             --cov genes.covariates.pc50.txt.gz --nominal  0.01  --chunk  $j  20  --normal  --out
             nominals_$j\_20.txt" | qsub
         done

PERMUTATION ANALYSIS EXAMPLES

       o Permutation  pass,  correcting  for  technical  covariates, rank normal transforming the
         phenotype, and running 1000 permutations with a specific random seed  on  chromosome  22
         between 17000000 and 18000000 bp:

         QTLtools  cis  --vcf  genotypes.chr22.vcf.gz  --bed  genes.50percent.chr22.bed.gz  --cov
         genes.covariates.pc50.txt.gz --permute 1000  --region  chr22:17000000-18000000  --normal
         --seed 1354145 --out permutation.txt

       o Permutation  pass  with parallelization correcting for technical covariates, rank normal
         transforming  the   phenotype,   and   running   5000   permutations.    To   facilitate
         parallelization  on  compute  cluster,  we  developed an option to run the analysis into
         chunks of molecular phenotypes.   For  instance,  to  run  analysis  on  chunk  12  when
         splitting the example data set into 20 chunks, run:

         QTLtools  cis  --vcf  genotypes.chr22.vcf.gz  --bed  genes.50percent.chr22.bed.gz  --cov
         genes.covariates.pc50.txt.gz   --permute   5000   --chunk   12   20    --normal    --out
         permutations_12_20.txt

       o If  you  want  to  submit the whole analysis with 20 jobs on a compute cluster, just run
         (qsub needs to be changed to the job submission system used [bsub, psub, etc...]):

         for j in $(seq 0 20); do
             echo "QTLtools cis --vcf genotypes.chr22.vcf.gz  --bed  genes.50percent.chr22.bed.gz
             --cov  genes.covariates.pc50.txt.gz  --permute  5000  --chunk  $j  20 --normal --out
             permutations_$j\_20.txt" | qsub
         done

       o When your phenotypes in the BED file are grouped, you can perform a permutation pass  at
         the phenotype group level in order to discover group-level QTLs:

         QTLtools  cis  --vcf  genotypes.chr22.vcf.gz  --bed  exons.50percent.chr22.bed.gz  --cov
         genes.covariates.pc50.txt.gz    --permute     1000     --normal     --grp-best     --out
         permutation.group.txt

CONDITIONAL ANALYSIS EXAMPLE

       Conditional analysis to discover independent signals.

       1 First  we  need  to  run  a  permutation analysis (see previous section), then calculate
         nominal p-value threshold for each gene.  Here an FDR of 5% is given as an example:

         cat permutations_*.txt | gzip -c > permutations_all.txt.gz
         Rscript ./script/qtltools_runFDR_cis.R permutations_all.txt.gz 0.05 permutations_all

       2 Now you can proceed with the  actual  conditional  analysis.   Here  splitting  into  20
         chunks, and when all complete concatenate the results:

         for j in $(seq 0 20); do
             echo  "QTLtools  cis --vcf genotypes.chr22.vcf.gz --bed genes.50percent.chr22.bed.gz
             --cov genes.covariates.pc50.txt.gz --mapping permutations_all.thresholds.txt --chunk
             $j 20 --normal --out conditional_$j\_20.txt" | qsub
         done
         cat conditional_*.txt > conditional_all.txt

       3 If  you  are interested in the most significant variants per independent signal, you can
         filter the results, using the backward p-value:

         awk '{ if ($21 == 1) print $0 }' conditional_all.txt > conditional_top_variants.txt

BUGS

       o Versions up to and including 1.2, suffer from a bug  in  reading  missing  genotypes  in
         VCF/BCF files.  This bug affects variants with a DS field in their genotype's FORMAT and
         have a missing genotype (DS field is .) in one of the samples, in which  case  genotypes
         for  all  the  samples  are  set  to missing, effectively removing this variant from the
         analyses.

       Please submit bugs to <https://github.com/qtltools/qtltools>

CITATION

       Delaneau, O., Ongen, H., Brown, A. et al. A complete tool set for molecular QTL  discovery
       and analysis. Nat Commun 8, 15452 (2017).  <https://doi.org/10.1038/ncomms15452>

AUTHORS

       Olivier Delaneau (olivier.delaneau@gmail.com), Halit Ongen (halitongen@gmail.com)