SCORE-SeqTDS: Score Tests for Sequencing Studies with Trait-Dependent Sampling

General information

SCORE-SeqTDS is a command-line program written in the C language to implement the methods described in Lin et al. (2013) for analyzing primary and secondary quantitative traits under trait-dependent sampling. The primary trait is the trait that is used to select subjects for sequencing, and all other traits are treated as secondary. Each quantitative trait is related to a genetic variable and possibly covariates through a linear regression model. Both the maximum likelihood estimation (MLE) and standard least-squares (LS) methods are available. The MLE method properly accounts for trait-dependent sampling whereas the LS method does not. The LS method is the ideal choice for random sampling and is approximately correct for analyzing secondary quantitative traits in case-control or case-only studies with rare diseases. SCORE-SeqTDS performs the LS analysis on secondary quantitative traits for random sampling, case-control and case-only sampling. For random sampling, all traits are treated as secondary (because the sampling does not depend on any particular trait.) The sampling scheme is specified through the option -sampling. A table in the OPTIONS section below summarizes the available analysis options for different sampling schemes.

The genetic variable pertains to the genotype in the single-variant analysis and to the burden score in the rare-variant analysis. The burden score may be determined externally (default) or internally. Use the option -gfile to specify the file that contains the external genetic variables. Otherwise, use the options -gfile and -mfile to specify the genotype file and the mapping file, respectively, for the internal creation of the burden scores. Use the option -test to request one of the six tests (T1, T5, MB, VT, SKAT, and customized test) under the additive genetic model (default). The T1, T5, VT, and SKAT tests under the dominant (recessive) genetic model can be obtained by using the option -dominant (-recessive). There are options for the minor allele frequency (MAF) upper bound, the minor allele count (MAC) lower bound and the call rate (CR) lower bound. Under the additive or dominant (recessive) model, MAC is defined as the number of subjects with at least one (two) observed mutation. A genetic variable is excluded from analysis if its MAC is smaller than the MAC lower bound or its CR is smaller than the CR lower bound. When the burden score is created internally, a variant is deleted if its MAF is greater than the MAF upper bound or its CR is smaller than the CR lower bound. By default, the MAF upper bound is 0.05, the MAC lower bound is 1 and the CR lower bound is 0. The MAFs may be determined internally (i.e., calculated from the genotype file) or externally (i.e., input in the mapping file). The MAF thresholds in the VT test can be determined internally (i.e., based on the unique MAFs among the aggregating variants) or externally (i.e., input in the mapping file).

If the quantitative trait of interest is measured in multiple studies, either as primary or secondary trait, then SCORE-SeqTDS can be used to analyze the trait of interest for one study at a time, and the results can be combined through meta-analysis. For T1, T5 and MB, SCORE-SeqTDS outputs the score statistics and their variance estimates. For the VT, SKAT and customized tests, SCORE-SeqTDS outputs the score vector and the information matrix. The results from multiple studies can be combined through the accompanying program MASS.

SYNOPSIS

SCORE-SeqTDS [-sampling sampling] [-method method] [-nsec nsec] [-pfile_seq phenofile.seq] [-cov cov] [-pfile_nonseq phenofile.nonseq] [-rounding NMAX] [-gfile genofile] [-mfile mapfile] [-wfile wfile] [-test test] [-dominant ] [-recessive ][-ofile outfile] [-msglog msglog] [-MAF MAF_UB] [-MAC MAC_LB] [-CR CR_LB] [-log ] [-INT ] [-R-INT SD_TYPE]

OPTIONS

1. MLE method

The MLE method is based on the maximum likelihood estimation of an observed-data likelihood that reflects trait-dependent sampling and thus is valid and efficient. Let Y₁ be the primary trait and Y₂ be the secondary trait. Also, let G₁ and G₂ denote the corresponding genetic variables and Z₁ and Z₂ denote the corresponding sets of covariates. The joint distribution of Y₁ and Y₂ is formulated through a bivariate linear regression model in which (G₁, Z₁) and (G₂, Z₂) are treated as independent variables. Let β₁ and β₂ denote the regression coefficients of G₁ for Y₁ and G₂ for Y₂, respectively. Then the null hypothesis of no genetic effect on the primary trait corresponds to H₀⁽¹⁾ : β₁ = 0 while the null hypothesis of no genetic effect on the secondary trait corresponds to H₀⁽²⁾ : β₂ = 0. An EM algorithm is used to maximize the observed-data likelihood and to calculate the score statistics for testing H₀⁽¹⁾ : β₁ = 0 and H₀⁽²⁾ : β₂ = 0. The estimation results for β₁ and β₂ are also provided.

2. LS method

The LS method is based on the standard least-squares estimation. For studying the genetic effect on the primary trait, the LS method is applied to the linear regression model relating Y₁ to (G₁, Z₁). For studying the genetic effect on the secondary trait, there are two versions of the LS method: LS-M is based on the linear regression of Y₂ on (G₂, Z₂) and LS-C is based on the linear regression of Y₂ on (Y₁, G₂, Z₂). All hypothesis tests are based on the score statistics.

INPUT FILES

The input files consist of a phenotype file for the sequenced subjects, a phenotype file for the non-sequenced subjects, a genotype file for the sequenced subjects and a mapping file. In all input files, the field separator character should be tab-delimiter and missing values are denoted by “-999” or “NA”.

PHENOTYPE FILE FOR SEQUENCED SUBJECTS

The 1^st row contains the header. The remaining rows represent sequenced subjects. Under trait dependent-sampling or case-control sampling, the 1^st column pertains to the primary trait. For trait-dependent sampling, the primary trait is quantitative. For case-control sampling, the primary trait takes the value 0 or 1. Under case-only or random sampling, the 1^st column is empty. If there are nsec secondary quantitative traits, then the 2^nd to the (nsec + 1)^th columns pertain to the secondary traits. The remaining columns are covariates if there are any. The order of subjects in the rows of this file should be identical to the order of subjects in the columns of the genotype file. No missing value is allowed on the primary trait and all variables should be numeric.

FILE FOR SELECTING COVARIATES

If the covariates are different among traits, a file with the selection indicators for covariates should be specified. This file contains one matrix with entries 0 or 1. The value of 1 indicates the covariate on the column is used by the trait on the row, the value of 0 indicates otherwise. The order of the traits on the rows and covariates on the columns in this file should be the same as the order in the phenotype file for sequenced subjects.

PHENOTYPE FILE FOR NON-SEQUENCED SUBJECTS

The first column contains the values of the primary trait for the non-sequenced subjects and the second column contains the frequencies of the trait values. If only the first column is provided, the trait values are all assumed to have frequencies of 1. No missing value is allowed. This file is not needed if the LS method is specified.

Note: To speed up computation, it is suggested that the values of the primary trait (1^st column) in both phenotype files are rounded such that there is a total of at most 500 distinct values. If that is not the case, the program will automatically round the values to provide at most 500 unique values.

GENOTYPE FILE, MAPPING FILE, AND WEIGHT FILE

The table below shows the option setup, genotype file, mapping file, and weight file for different scenarios.

OUTPUT

OUTPUT FILE

1. Tests with a single genetic variable (single-variant, T1, T5, or MB test)
   The format is shown as follows:The 1^st column contains the gene ID and the remaining columns under each main header <Trait_Method> contain the test statistic, the p-value, the genetic effect estimate, its standard error estimate, and the warning message for each “Trait” and “Method” combination. For example, if the user wants to analyze two secondary traits, say BMI and BP, using the MLE method (option: “-nsec 2″), there will be 10 columns following the ID column in the output file, and the two main headers will be “BMI_MLE” and “BP_MLE”. If the user wants to analyze these two secondary traits using the LS method (option: “-method LS -nsec 2″), there will be 20 columns following the ID column since there are two versions of the LS method for the secondary trait(s), and the 4 main headers will be “BMI_LS-M”, “BMI_LS-C”, “BP_LS-M”, and “BP_LS-C”. Invalid values are denoted by “NA” in the output file.

   Different types of warning message are as follows:

   Singular_covariates: Covariate matrix becomes singular after removing subjects with missing genetic variable or missing “Trait”.

   Non-convergence: EM algorithm does not converge.

   Low_MAC: MAC is smaller than MAC_LB after removing subjects with missing “Trait”.

   Singular_score: The genetic variable contains only one distinct value after removing subjects with missing “Trait”.

   Undefined_V: Variance estimate of the score statistic is undefined (numerically 0 or negative).

   Pvalue_Fail: p-value calculation fails to obtain the desired level of accuracy.

   Var_Null: No valid genetic variable for the test.

2. Tests with multiple genetic variables (VT, SKAT, or customized test)The 1^st column contains the gene ID and the remaining columns under each main header <Trait_Method> contain the test statistic, the p-value, and the warning message for each “Trait” and “Method” combination.

LOG FILE

Log files include a set of files for each “Trait” and “Method” combination. These files contains the MAC, the score vector, the information matrix, and some other related information. The name of the file is “log_trait_method”: log is the name specified by the option -msglog; method is the method used to analyze the trait. For example, if the user wants to analyze two secondary traits (option “-nsec 2”) under trait dependent sampling (option “-sampling TDS”), say BMI and BP, using the MLE method (option “-method MLE”), there will be 2 log files with names “log_BMI_MLE” and “log_BP_MLE”. If the user wants to analyze these two secondary traits using the LS method (option “-method LS”), then there will be 4 log files with names “log_BMI_LS-M”, “log_BMI_LS-C”, “log_BP_LS-M”, and “log_BP_LS-C”. Invalid values in the output files are denoted by “NA”. The log files for multiples studies can be combined through the accompanying program MASS. We describe below the log file for each test.

1. Tests with a single genetic variable (single-variant, T1, T5, or MB test)
   The 1^st column contains the gene ID and the 2^nd column contains the MAC. The remaining columns contain the score statistic and its variance estimate.
2. VT test
   For each gene, “log_trait_method” contains the MAF thresholds, the MACs, the (vector-valued) score statistic, and the information matrix. The following is an example for two genes:
3. SKAT test
   For a set of SNPs within a gene, “log_trait_method” contains the SNP IDs, the MAFs, the MACs, the number of non-missing genotypes, the counts of homozygous reference, heterozygous, and homozygous alternative, the (vector-valued) score statistic, and the information matrix. The following is an example for two genes:
4. Customized test
   For a set of genetic scores based on the customized weight matrix within a gene, “log_trait_method” contains the weight identifiers, the MACs, the score statistic and the information matrix. The following is an example for two genes:

EXAMPLE

The example includes a phenotype file for the sequenced subjects “ldl_seq.txt” , a phenotype file for the non-sequenced subjects “ldl_nonseq.txt” and a genotype file “ldl_geno100.txt“. We want to apply the MLE method to two secondary traits BMI and BP (-nsec 2), apply the MAC lower bound of 2 (-MAC 2) and perform the log-transformation on the primary trait (-log). Enter the command

$ SCORE-SeqTDS -nsec 2 -pfile_seq ldl_seq.txt -pfile_nonseq ldl_nonseq.txt -gfile ldl_geno100.txt -ofile ldl_output.txt -MAC 2 -log

to obtain the results given in the output file “ldl_output.txt” , and the log files “log_BMI_MLE” and “log_BP_MLE” .

Alternatively, we can use “ldl_nonseq2.txt” as the phenotype file for the non-sequenced subjects and obtain the same results by entering the command

$ SCORE-SeqTDS -nsec 2 -pfile_seq ldl_seq.txt -pfile_nonseq ldl_nonseq2.txt -gfile ldl_geno100.txt -ofile ldl_output.txt -MAC 2 -log

MULTIPLE CPUs

If there are a large number of genetic variables, it is preferable to use multiple CPUs. To this end, the user can divide the genotype file into several files, each of which contains a subset of the original genetic variables. Then the software can be run separately for each genotype file on a different CPU. In this way, it takes approximately 1/K (K is the number of jobs that are split into) of the original time to complete the whole set of analysis. We provide a R script to illustrate this approach. In that example, 100 genotype files are created from the original genotype file and 100 jobs are submitted simultaneously to multiple CPUs.

REFERENCE

Dan-Yu Lin, Donglin Zeng, and Zheng-Zheng Tang (2013). Quantitative Trait Analysis in Sequencing Studies Under Trait-Dependent Sampling. Proceedings of the National Academy of Sciences of the United States of America, 110, 12247-12252.

DOWNLOAD

SCORE-SeqTDS for 64-bit x86 based Linux [updated November 18, 2015]

Example files [updated November 27, 2013]

VERSION HISTORY