Genome Analysis Wiki - User contributions [en]

Rvtests

2016-05-12T04:31:25Z

Zhanxw:

<span style="color: red">'''(Updated: January 2015; latest documentation can always be found at: [http://zhanxw.github.io/rvtests/ http://zhanxw.github.io/rvtests/])'''</span>

= Introduction =
<span style="color: red">'''(Updated: January 2015; latest documentation can always be found at: [http://zhanxw.github.io/rvtests/ http://zhanxw.github.io/rvtests/])'''</span>

Rvtests, which stands for Rare Variant tests, is a flexible software package for genetic association studies. It is designed to support unrealted individual or related (family-based) individuals. Both quantitative trait and binary trait are supported. It includes a variety of association tests (e.g. single variant score test, burden test, variable threshold test, SKAT test, fast linear mixed model score test). It takes [http://www.1000genomes.com/ VCF] format as genotype input file and takes PLINK format phenotype file and covariate file. From our practice, it is capable to analyze 8,000 related individuals using less than 400 Mb memory.

= Download =

Source files can be downloaded from [https://github.com/zhanxw/rvtests/archive/master.zip github] or [https://github.com/zhanxw/rvtests github page].<br />Executable binary files (for Linux 64bit) can be downloaded from [https://github.com/zhanxw/rvtests/releases/download/v1.8.6/rvtests-20150104.tar.gz here].

= Quick Tutorial =

Here is a quick example of how to use ''rvtests'' software in typical use cases.

== Single variant tests ==

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --single wald,score</pre>
This specifies single variant Wald and score test for association<br />tests for every variant in the <code>input.vcf</code> file. The 6th column of the phenotype file, <code>phenotype.ped</code>, which is in PLINK format, is used. Rvtests will automatically check whether the phenotype is binary trait or quantitative trait.<br />For binary trait, the recommended way of coding is to code controls as 1, cases as 2, missing phenotypes as -9 or 0.

For other types of association tests, you can refer to [[#models|Models]]

== Groupwise tests ==

Groupwise tests includes three major kinds of tests.

* Burden tests: group variants, which are usually less than 1% or 5% rare variants, for association tests. The category includes: CMC test, Zeggini test, Madsen-Browning test, CMAT test, and rare-cover test.<br />
* Variable threshold tests: group variants under different frequency thresholds.<br />
* Kernel methods: suitable to tests rare variants having different directions of effects. These includes SKAT test and KBAC test.

All above tests requires to group variants into a unit. The simplist case is to use gene as grouping unit. For different grouping method, see [[#Grouping|Grouping]].

To perform rare variant tests by gene, you need to use <code>--geneFile</code> to specify the gene range in a refFlat format. We provided different gene definitions in the [[#Resources|Resources]] section. You can use <code>--gene</code> to specify which gene(s) to test. For example, specify <code>--gene CFH,ARMS2</code> will perform association tests on CFH and ARMS2 genes. If there is no providing <code>--gene</code> option, all genes will be tests.

The following command line demonstrate how to use CMC method, variable threshold method(proposed by Price) and kernel based method (SKAT by Shawn Lee and KBAC by<br />Dajiang Liu) to test every gene listed in ''refFlat_hg19_uniq_gene.txt.gz''.

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --geneFile refFlat_hg19_uniq_gene.txt.gz --burden cmc --vt price --kernel skat,kbac</pre>
== Related individual tests ==

To test related individuals, you will need to first create a kinship matrix:

<pre>vcf2kinship --inVcf input.vcf --bn --out output</pre>
The option <code>--bn</code> means calculating empirical kinship using Balding-Nicols method. You can specifiy <code>--ibs</code> to obtain IBS kinship or use <code>--pedigree input.ped</code> to calculate kinship from known pedigree information.

Then you can use linear mixed model based association tests such as Fast-LMM score test, Fast-LMM LRT test and Grammar-gamma tests. An exemplar command is shown:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --kinship output.kinship --single famScore,famLRT,famGrammarGamma</pre>
== Meta-analysis tests ==

The meta-analysis models outputs association test results and genotype covariance matrix. These statistics can be used in rare variant association analysis.<br />We provide single variant score test and generate genotype covariance matrix.<br />You can use command:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output</pre>
In a more realistic scenario, you may want to adjust for covariates and want to inverse normalized residuals obtained in null model ([http://www.nature.com/ng/journal/vaop/ncurrent/full/ng.2852.html link] to our methodology paper), then this command will work:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out output</pre>
Here the <code>--covar</code> specify a covariate file, and <code>--covar-name</code> specify which covariates can used in the analysis. Covariate file format can be found [[#Covariate%20file|here]]. <code>--inverseNormal --useResidualAsPhenotype</code> specifies trait transformation method. That means first fit a regression model of the phenotype on covariates (intercept automatically added), then the residuals are inverse normalized. Trait transformation details can be found [[#Trait%20transformation|here]].

We support both unrelated individuals and related indivudlas (e.g. family data). You need to append <code>--kinship input.kinship</code> to the command line:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output --kinship input.kinship</pre>
The file <code>input.kinship</code> is calculated by <code>vcf2kinship</code> program, and usage to this program is described in [[#related-individual-tests|Related individual tests]].

=== Dominant models and recessive models ===

Dominant and recessive disease models are supported by appending "dominant" and/or "recessive" after "--meta" option. For example, use "--meta dominant,recessive" will<br />generate two sets of files. For dominant model, they are "prefix.MetaDominant.assoc" and "prefix.MetaDominantCov.assoc.gz"; for recessive model,<br />they are "prefix.MetaRecessive.assoc" and "prefix.MetaRecessiveCov.assoc.gz". Internally, in dominant models, genotypes 0/1/2 are coded as 0/1/1; in recessive models, genotypes 0/1/2 are<br />coded as 0/0/1. Missing genotypes will be imputed to the mean.

= Input files =

== Genotype file (VCF) ==

Rvtests supports VCF (Variant Call Format) files. Files in both plain txt format or gzipped format are supported. To use group-based rare variant tests, indexed the VCF files using [http://samtools.sourceforge.net/tabix.shtml tabix] are required.

Here are the commands to convert plain text format to bgzipped VCF format:

<pre>(grep ^"#" $your_old_vcf; grep -v ^"#" $your_old_vcf | sed 's:^chr::ig' | sort -k1,1n -k2,2n) | bgzip -c > $your_vcf_file
tabix -f -p vcf $your_vcf_file</pre>
The above commands will (1) remove the <code>chr</code> prefix from chromosome names; (2) sort VCF files by chromosome first, then by chromosomal positions; (3) compress using bgzip; (4) create tabix index.

Rvtests support genotype dosages. Use <code>--dosage DosageTag</code> to specify the dosage tag. For example, if VCF format field is "GT:EC" and individual genotype fields is "0/0:0.02", you can use <code>--dosage EC</code>, and rvtests will use the dosage 0.02 in the regression models.

== Phenotype file ==

You can use <code>--mpheno $phenoypeColumnNumber</code> or <code>--pheno-name</code> to specify a given phenotype.

An example phenotype file, (<code>example.pheno</code>), has the following format:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.7642934435605 -0.733862638327895 -0.980843608339726 2
P2 P2 0 0 0 0.457111744989746 0.623297281416372 -2.24266162284447 1
P3 P3 0 0 0 0.566689682543218 1.44136462889459 -1.6490100777089 1
P4 P4 0 0 0 0.350528353203767 -1.79533911725537 -1.11916876241804 1
P5 P5 0 0 1 2.72675074738545 -1.05487747371158 -0.33586430010589 2</pre>
Phenotype file is specified by the option <code>--pheno example.pheno</code> . The default phenotype column header is “<code>y1</code>”. If you want to use alternative columns as phenotype for association analysis (e.g the column with header y2), you may specify the phenotype by column or by name using either

* --mpheno 2<br />
* --pheno-name y2

'''NOTE:''' to use “<code>--pheno-name</code>”, the header line must starts with “<code>fid iid</code>” as PLINK requires.

In phenotype file, missing values can be denoted by NA or any non-numeric values. Individuals with missing phenotypes will be automatically dropped from subsequent association analysis. For each missing phenotype value, a warning will be generated and recorded in the log file.

When the phenotype values are only 0, 1 and 2, rvtests will automatically treat it as binary traits. However, if you want to treat it as continuous trait, please use "<code>--qtl</code>" option.

== Covariate file ==

You can use <code>--covar</code> and <code>--covar-name</code> to specify covariates that will be used for single variant association analysis. This is an optional parameter. If you do not have covariate in the data, this option can be ignored.

The covariate file, (e.g. <code>example.covar</code>) has a similar format as the phenotype file:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.911 -1.465 -0.817 1
P2 P2 0 0 0 2.146 -2.451 -0.178 2
P3 P3 0 0 0 1.086 -1.194 -0.899 1
P4 P4 0 0 0 0.704 -1.052 -0.237 1
P5 P5 0 0 1 2.512 -3.085 -2.579 1</pre>
The covariate file is specified by the <code>--covar</code> option (e.g. <code>--covar example.covar</code>). To specify covariates that will be used in the association analysis, the option <code>--covar-name</code> can be used. For example, when age, bmi and 3 PCs are used for association analysis, the following option can be specified for the rvtest program, i.e.<br /><code>--covar example.covar --covar-name age,bmi,pc1,pc2,pc3</code>.

Note: Missing data in the covariate file can be labeled by any non-numeric value (e.g. NA). They will be automatically imputed to the mean value in the data file.

== Trait transformation ==

In this meta-analysis, we use inversed normal transformed residuals in the association analysis, which is achieved by using a combination of <code>--inverseNormal</code> and <code>--useResidualAsPhenotype</code>. Specifically, we first fit the null model by regressing phenotype on covariates. The residuals are then inverse normal transformed (see Appendix A more detailed formulae for transformation). Transformed residuals will be used to obtain score statistics.

In meta analysis, an exemplar command for using rvtest looks like the following:

<pre>./rvtest --inVcf $vcf --pheno $example.pheno --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out $output_prefix </pre>
= Models =

Rvtests support various association models.

== Single variant tests ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only null model is used to performed the test
|-
|Wald test
|align="center"|wald
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only fit alternative model, and effect size will be estimated
|-
|Exact test
|align="center"|exact
|align="center"|B
|align="center"|N
|align="center"|U
|Fisher's test
|-
|Fam LRT
|align="center"|famLRT
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model
|-
|Fam Score
|align="center"|famScore
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model style likelihood ratio test
|-
|Grammar-gamma
|align="center"|famGrammarGamma
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Grammar-gamma method
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--single score</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Burden tests ==

{|
!Burden tests
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants by Bingshan Li.
|-
|Zeggini
|align="center"|zeggini
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Aggregate counts of rare variants by Morris Zeggini.
|-
|Madsen-Browning
|align="center"|mb
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Madsen.
|-
|Fp
|align="center"|fp
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Danyu Lin.
|-
|Exact CMC
|align="center"|exactCMC
|align="center"|B
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Fisher's exact test.
|-
|CMC Wald
|align="center"|cmcWald
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Wald test.
|-
|RareCover
|align="center"|rarecover
|align="center"|B
|align="center"|N
|align="center"|U
|Find optimal grouping unit for rare variant tests by Thomas Hoffman.
|-
|CMAT
|align="center"|cmat
|align="center"|B
|align="center"|N
|align="center"|U
|Test non-coding variants by Matt Z.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--burden cmc</code> will apply CMC test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Variable threshold models ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Variable threshold model
|align="center"|vt
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Every rare-variant frequency cutoffs are tests by Alkes Price.
|-
|Variable threshold CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|This models is natiive so that it output CMC test statistics under all possible frequency cutoffs.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--vt price</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Kernel models ==

{|
!Kernel
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|SKAT
|align="center"|skat
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Sequencing kernel association test by Shawn Lee.
|-
|KBAC
|align="center"|kbac
|align="center"|B
|align="center"|N
|align="center"|U
|Kernel-based adaptive clustering model by Dajiang Liu.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--kernel skat</code> will apply SKAT test.<br />To further customize SKAT test, you can use ''--kernel skat[nPerm=100:alpha=0.001:beta1=1:beta2=20]'' to specify permutation counts, type-1 error,<br />beta distribution parameters for upweighting rare variants. Rvtests will output a message showing:

<pre>[INFO] SKAT test significance will be evaluated using 10000 permutations at alpha = 0.001 (beta1 = 1.00, beta2 = 20.00)</pre>
(##) In trait column, B and Q stand for binary, quantitiave trait.

== Meta-analysis models ==

{|
!Type
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|standard score tests
|-
|Dominant model
|align="center"|dominant
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under dominant disease model
|-
|Recessive model
|align="center"|recessive
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under recessive disease model
|-
|Covariance
|align="center"|cov
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|covariance matrix
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--meta score,cov</code> will generate score statistics and covariance matrix for meta-analysis.<br />(##) In trait column, B and Q stand for (b)inary, (q)uantitiave trait.

== Utility models ==

Rvtests has an usually option <code>--outputRaw</code>. When specify this, rvtests can output genotypes, phenotype, covariates(if any) and collapsed genotype to tabular files. These files can be imported into other software (e.g. R) for further analysis.

= Association test options =

== Sample inclusion/exclusion ==

Rvtests can flexibly specify which sample(s) to include or exclude:

<pre> --peopleIncludeID : give IDs of people that will be included in study
--peopleIncludeFile : from given file, set IDs of people that will be included in study
--peopleExcludeID : give IDs of people that will be included in study
--peopleExcludeFile : from given file, set IDs of people that will be included in study</pre>
<code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> are used to include/exclude samples from command line.<br />For example, specify <code>--peopleIncludeID A,B,C</code> will include A, B and C sample from the VCF files if they exists.<br /><code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> followed by a file name will include or exclude the IDs in the file.<br />So to include sample A, B and C, you can provide a file, <code>people.txt</code>, looks like:

<pre>A
B
C</pre>
Then use <code>--peopleIncludeFile people.txt</code> to include them in the analysis.

== Variant site filters ==

It is common that different frequency cutoffs are applied in rare-variant analysis.<br />Therefore, rvtests specify frequency cutoffs.

Frequency Cutoff

<pre> --freqUpper : Specify upper minor allele frequency bound to be included in analysis
--freqLower : Specify lower minor allele frequency bound to be included in analysis</pre>
If you specify <code>--freqLower 0.01 --freqUpper 0.05</code>, only the variants with minor allele ferquncy between 0.01 and 0.05 (boundary inclusive) will be analyzed.

Similar to sample inclusion/exclusion options, you can specify a range of variants to be included by<br />specifying <code>--rangeList</code> option. For example <code>--rangeList 1:100-200</code> will include the chromosome 1 position 100bp to 200bp region.<br />Alternatively, use a separate file, <code>range.txt</code>, and <code>--rangeFile range.txt</code> to speicify association tests range.

<pre> --rangeList : Specify some ranges to use, please use chr:begin-end format.
--rangeFile : Specify the file containing ranges, please use chr:begin-end format.
--siteFile : Specify the file containing sites to include, please use "chr pos" format.</pre>
It is supported to filter variant site by site depth, minor allele count or annotation (annotated VCF file is needed).

<pre> --siteDepthMin : Specify minimum depth(inclusive) to be incluced in analysis
--siteDepthMax : Specify maximum depth(inclusive) to be incluced in analysis
--siteMACMin : Specify minimum Minor Allele Count(inclusive) to be incluced in analysis
--annoType : Specify annotation type that is follwed by ANNO= in the VCF INFO field, regular expression is allowed</pre>
''NOTE'': <code>--annoType Nonsynonymous</code> will only analyze nonsynonymous variants where they have <code>ANNO=Nonsynonymous</code> in the INFO field.<br />VCF with annotatino information are called annotated VCF here. And to annotate<br />a VCF file, you can use [https://github.com/zhanxw/anno ANNO], a fast and accurate annotation software.

== Genotype filters ==

Genotype with low depth or low quality can be filtered out by these options:

<pre> --indvDepthMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis
--indvDepthMax : Specify maximum depth(inclusive) of a sample to be incluced in analysis
--indvQualMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis</pre>
When genotypes are filtered, they are marked as missing genotypes.<br />Consequently, samples with missing genotype may or may not be included in the analysis.<br />That means samples with genotypes may be dropped (<code>--impute drop</code>)<br />or may still be included (<code>--impute mean</code> or <code>--impute hwe</code>).<br />By default, genotypes are imputed to its means.<br />See next section about how you like to handle missing genotypes.

== Handle missing genotypes and phenotypes ==

When genotypes are missing (e.g. genotype = "./.") or gentoypes are filtered out,<br />there are three options to handle them: (1) impute to its mean(default option); (2) impute by HWE equilibrium; (3) remove from the model.<br />Use <code>--impute [mean|hwe|drop]</code> to specify which option to use.

When quantitative phenotypes are missing, for example, some samples have gneotype files, but not phenotypes,<br />rvtests can impute missing phenotype to its mean.

''NOTE:'' Do not use <code>--imputePheno</code> for binary trait.

In summary, the following two options can be used:

<pre> --impute : Specify either of mean, hwe, and drop
--imputePheno : Impute phenotype to mean by those have genotypes but no
phenotpyes

</pre>
== Specify groups (e.g burden unit) ==

Rare variants association tests are usually performed in gruops of variants.<br />The natural grouping unit is gene. Rvtests can read gene definition file in <code>refFlat</code> format,<br />and perform association for each gene. Use <code>--geneFile</code> option to specify the gene file name.<br />For example, <code>--geneFile refFlat_hg19.txt.gz</code> will use <code>refFlat_hg19.txt.gz</code> as gene definition file,<br />and then perform association tests for every gene. Use <code>--gene</code> to specify a subset of genes to test.<br />For example, <code>--gene CFH</code> will only test CFH gene.

Alternative grouping unit can be specified as ''set''.<br />These ''sets'' are treated similar to gene.<br />You can thus use <code>--setFile</code> to define sets (similar to <code>--geneFile</code> option),<br />and use <code>--set</code> to define a specific set (similar to <code>--gene</code> option).<br />Additionally, use <code>--setList</code> can speicify a set to test from command line.

The format of a set file is: (1) set names; (2) ranges (e.g. chrom:begin-end);<br />For example, you have a set file, <code>example.set</code>, like this:

<pre>set1 1:100-200,1:250-300
set2 2:500-600</pre>
You can specify <code>--setFile example.set --set set2</code> to group variants<br />within chromosome 2, position 500 to 600bp.<br />If you want to test a particular region, for example, chromosome 2, position 500 to 550bp,<br />but do not want to make another file, you can use <code>--setList 2:500-600</code>.

In summary, options related to ''Grouping Unit'' are listed below:

<pre> --geneFile : specify a gene file (for burden tests)
--gene : specify which genes to test
--setList : specify a list to test (for burden tests)
--setFile : specify a list file (for burden tests, first two columns:
setName chr:beg-end)
--set : specify which set to test (1st column)</pre>
= Sex chromosome analysis =

Rvtests suppport X chromosome analysis. In human X chromosome, there is PAR (pseudoautosomal region) and non-PAR region.<br />For males, there are two X allele in PAR region and one allele in non-PAR region.<br />While the PAR region is treated in the same way as autosomes, rvtests treate non-PAR region differently.<br />Below we will describe the details about how rvtests handles non-PAR region.

''Prepare data''. According to VCF standard, male genotype needs to coded as 0 or 1. For compatibility, rvtests also support 0/0 or 1/1 coding.<br />In VCF files, male genotypes can be written as "0", "1", "0|0", "0/0", "1|1", "1/1". All other genotypes will be treated as missing.

''Genotype in the regression model''. For consistencmaine, male genotypes are converted to 0 or 2.

''MetaScore results''. If specify <code>--meta score</code>, the output file <code>prefix.MetaScore.assoc</code> includes PAR-region and non-PAR region analysis.<br />But in the non-PAR region, the difference is that Hardy-Weinberg P-value are calculated using female samples.

''Related individuals''. Just append <code>--xHemi</code> to the <code>vcf2kinship</code> (more details in [[#kinship-generation|Kinship generation]]) and <code>rvtest</code> command lines. Rvtests can recognize non-PAR region kinship and use it in the analysis.

''PAR region''. PAR region is defined as two regions X:60001-2699520 and X:154931044-155270560. Use <code>--xLabel</code> can speicify which chromosome has PAR region (default: 23|X)<br />and use <code>--xParRegion</code> to specify PAR region (default: hg19, meaning '60001-2699520,154931044-155260560' in the UCSC build hg19, specify "hg18" will use PAR region definition in the UCSC build hg18).

= Kinship generation =

Analysis of related individual usually requires estimation of kinship. You can a separate tool, <code>vcf2kinship</code>.<br /><code>vcf2kinship</code> is usually included in rvtests binary distribution or can be built from software source codes.

<code>vcf2kinship</code> can calcualte pedigree kinship using a pedigree input file (PED format, see [[#phenotype-file|Phenotype file]], use option <code>--ped</code>).<br />The output file name is specified by <code>--prefix</code> option. If you use <code>--prefix output</code> then the output files will include <code>output.kinship</code>.

It can also calculate empirical kinship using genotyp input file (VCF format, see [[#genotype-file-vcf|Genotype file (VCF)]], use option <code>--inVcf</code>).<br />For empiricial kinship, you also need to specify the kinship model, either Balding-Nicols model (ue option <code>--bn</code>) or Identity-by-state model (use option <code>--ibs</code>).

In sex chromosome analysis, it is often required to generate kinship on X chromsoome regions, then you need to speicfy <code>--xHemi</code>. If your input VCF file has different X chromosome label (e.g. chromosome name is '23' instead of 'X'), you can use <code>--xLabel 23</code>.

If principal component decomposition (PCA) results are needed, you can use option <code>--pca</code>. Then output files with suffix '.pca' include PCA results.

When dealing with large input files, it is often preferred to use multiple CPU to speed up calculation using the option <code>--thread N</code> in which N is the number of CPUs.

For example, to generate pedigree-based kinship (<code>--ped</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region, the command line is:

<pre>vcf2kinship --ped input.ped --xHemi --out output</pre>
To generate empirical kinship (<code>--inVcf</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region using Balding-Nicols model, the command line is:

<pre>vcf2kinship --inVcf input.vcf.gz --ped input.ped --bn --xHemi --out output</pre>
NOTE: you need to provide a pedigree file (PED) in the above case, as <code>vcf2kinship</code> need the sex information of samples.

= Frequently Asked Questions (FAQ) =

* Does rvtests suppport binary traits of related-individuals?

Not yet. It's a complex scenario and we have not found good solutions.

* Can you provide a list of command line options?

Rvtests have build help taht can be found by executing <code>rvtest --help</code>.<br />We also put all available options in this [https://github.com/zhanxw/rvtests/wiki/Command-Line-Options link].

* Can you provide standard error (SE) or confidence interval (CI) for the estimated Beta in the score model?

In the output of MetaScore model (--meta score), the standard error is the inverse of SQRT_V_STAT.<br />For example, if SQRT_V_STAT = 2, that means the standard error of estimated beta is 1/2 = 0.5.

* Why the INFORMATIVE_ALT_AC, N_REF and N_ALT columns have zero counts for certain chromosome X regions in meta-analysis models?

These counts are calculated from female individuals. If your study only has male samples, rvtests cannot report these counts. Because if a male carries a non-reference allele, we cannot conclude that this is heterozygous (0/1) site or homozygous alternatives (1/1) site.

= Feedback/Contact =

Questions and requests can be sent to Xiaowei Zhan<br />([mailto:zhanxw@umich.edu [mailto:zhanxw@umich.edu zhanxw@umich.edu]])<br />or Goncalo Abecasis<br />([mailto:goncalo@umich.edu [mailto:goncalo@umich.edu goncalo@umich.edu]])

Rvtests is a collaborative effort by Youna Hu, Bingshan Li, Dajiang Liu.

LiftOver

2015-07-15T22:33:41Z

Zhanxw: /* Resources */

LiftOver is a necesary step to bring all genetical analysis to the same reference build.
LiftOver can have three use cases:

(1) [[#Lift genome positions | Convert genome position from one genome assembly to another genome assembly]]

In most scenarios, we have known genome positions in NCBI build 36 (UCSC hg 18) and hope to lift them over to NCBI build 37 (UCSC hg19).

(2) [[#Lift dbSNP rs numbers | Convert dbSNP rs number from one build to another]]

(3) [[#Lift Merlin/PLINK format | Convert both genome position and dbSNP rs number over different versions]]

It is likely to see such type of data in Merlin/PLINK format.

We will explain the work flow for the above three cases. In the rest of this article,
our example is to lift over from lower/older build to newer/higher build, as it is the common practice.

Using different tools, liftOver can be easy.
For example, UCSC liftOver tool is able to lift [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format] file between builds.
With our customized scripts, we can also lift rsNumber and Merlin/PLINK data files.

== Lift genome positions ==
Genome positions are best represented in [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format]. UCSC provides tools to convert BED file from one genome assembly to another.

=== Binary liftOver tool ===
We need [http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver liftOver binary] from UCSC and [http://hgdownload.cse.ucsc.edu/goldenPath/hg18/liftOver/hg18ToHg19.over.chain.gz hg18 to hg 19 chain file].

Provide BED format file (e.g. input.bed)

NOTE: Use the 'chr' before each chromosome name
<pre>
chr1 743267 743268 rs3115860
chr1 766408 766409 rs12124819
chr1 773885 773886 rs17160939
</pre>
Run liftOver:

liftOver input.bed hg18ToHg19.over.chain.gz output.bed unlifted.bed

unlifted.bed file will contain all genome positions that cannot be lifted. The reason for that varies. See [[#Various reasons that lift over could fail | Various reasons that lift over could fail]]

=== Web interface ===
Alternatively, you can lift over BED file in web interface
at: [http://genome.ucsc.edu/cgi-bin/hgLiftOver Link]
Web interface can tell you why some genome position cannot
be lifted if you click "Explain failure messages"

== Lift dbSNP rs numbers ==
rs number is release by dbSNP. UCSC also make their own copy from each dbSNP version. Be aware that the same version of dbSNP from these two centers are not the same.
When we convert rs number from lower version to higher version, there are practically two ways.

=== Use RsMergeArch and SNPHistory ===
It is necessary to quickly summarize how dbSNP merge/re-activate rs number:

# when different rs number are found to refer to the same SNP, then higher rs number will be merged to lower rs number, and the merging will be recorded in RsMergeArch.bcp.gz.
# when rs number have to be retracted, rs number will be recorded in SNPHistory.bcp.gz
# a retracted SNP can be [http://www.ncbi.nlm.nih.gov/books/NBK44496/#Schema.rs4823903_which_has_merged_into re-activated] in SNPHistory.bcp.gz by adding comment

With the above in mind, we are able to combine these two tables to obtain the relationship between older rs number and new rs number.
We have developed a script (for internal use), named [http://genome.sph.umich.edu/wiki/LiftRsNumber.py liftRsNumber.py] for lift rs numbers between builds.
This scripts require RsMergeArch.bcp.gz and SNPHistory.bcp.gz, those can be found in [[#Resources | Resources]].

Example input:
<pre>
3000
3001
3002
</pre>

Command:
<pre>
python liftRsNumber.py input.rs > output.rs
</pre>

Example output:
<pre>
unchanged 3000
lifted 2032
unchanged 3002
</pre>

== Lift Merlin/PLINK format ==
In Merlin/PLINK .map files, each line contains both genome position and dbSNP rs number. Our goal here is to use both information to liftOver as many position as possible.
There are 3 methods to liftOver and we recommend the first 2 method. The first method is common and applicable in most cases, and in our observations it lifts the most genome positions, however, it does not reflect the rs number change between different dbSNP builds. The second method is more robust in the sense that each lifted rs number has valid genome position, as it lift over old rs number as the first step by using dbSNP data. The third method is not straigtforward, and we just briefly mention it.

=== Lift Merlin format ===

PLINK format and [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin format are nearly identical].
The difference is that Merlin .map file have 4 columns. We will show
the lift over procedure for PLINK format, then you can use:

awk '{print $1,$2,"\t",$3;}' PLINK.map > Merlin.map

to obtain Merlin .map file.

=== Lift PLINK format ===
[http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml PLINK] format usually referrs to .ped and .map files.

==== Method 1 ====
We mainly use UCSC LiftOver binary tools to help lift over.
We have a script [[#Resources | liftMap.py]], however, it is recommended to understand the job step by step:

(1) Convert .map to .bed file

By rearrange columns of .map file, we obtain a standard BED format file.

(2) LiftOver .bed file

Use method mentioned [[#Lift genome positions | above]] to convert .bed file from one build to another.

(3) Convert lifted .bed file back to .map file

Rearrange column of .map file to obtain .bed file in the new build.

(4) Modify .ped file

.ped file have many column files. By convention, the first six columns are family_id, person_id, father_id, mother_id, sex, and phenotype.
From the 7th column, there are two letters/digits representing a genotype at the certain marker. In step (2), as some genome positions cannot
be lifted to the new version, we need to drop their corresponding columns from .ped file to keep consistency. You can use PLINK --exclude those snps,
see [http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml#exclude Remove a subset of SNPs].

(5) (optionally) change the rs number in the .map file

Similar to the human reference build, dbSNP also have different versions. You may consider change rs number from the old dbSNP version to new dbSNP version
depending on your needs. Such steps are described in [[#Lift dbSNP rs numbers | Lift dbSNP rs numbers]].

==== Method 2 ====
The idea is to use [[#Resources |LiftRsNumber.py]] to convert old rs number to new rs number, use the data file b132_SNPChrPosOnRef_37_1.bcp.gz (a data file containing each dbSNP and its positions in NCBI build 37), and adjust .map and .ped files accordingly.

(1) Extract and lift rs numbers

Use the tools [[#Use RsMergeArch and SNPHistory | LiftRsNumber.py]] to lift the rs number in the map file from old build to new build.

(2) Lookup SNP positions from rs number

dbSNP provides a file [[#Resources | b132_SNPChrPosOnRef_37_1.bcp.gz]] which contains rsNumber, chromosome and its position.
Use this file along with the new rsNumber obtained in the first step.
In practice, some rs numbers do not exist in build 132, or not suitable to be considered ( e.g. they do not reside on human reference, or they are mapped to multiple locations, these scenarios are noted by the chromosome column with values like "AltOnly", "Multi", "NotOn", "PAR", "Un"), we can drop them in the liftover procedure.
We will obtain the rs number and its position in the new build after this step.

(3) Lift .map file and .ped file

To lift over .map files, we can scan its content line by line, and skip those not lifted rs number.
Accordingly, it is necessary to drop the un-lifted SNP genotypes from .ped file.

==== Method 3 ====
NCBI dbSNP team has provided a [[ #Resources | provisional map ]] for converting the genome position of a larget set dbSNP from NCBI build 36 to NCBI build 37.
In the second step, we have obtained unlifted genome positions, so we can try to use the table to convert those unlfted dbSNPs.
After this step, there are still some SNPs that cannot be lifted, as they are mostly located on non-reference chromosome.
Note: due to the limitation of the provisional map, some SNP can have multiple locations.
For example, we cannot convert rs10000199 to chromosome 4, 7, 12.
<pre>
10000199 A/G 4 166142415 166142415 2 3 G + 4 165922965 165922965 2 3 G +
10000199 A/G 7 4589694 4589694 2 3 C - 7 4623168 4623168 2 3 C -
10000199 A/G 12 57008620 57008620 2 3 C - 12 58722353 58722353 2 3 C -
10000199 A/G 5 156018406 156018406 2 3 C - 5 156085828 156085828 2 3 C -
</pre>

We can dissect this method into steps:

(1) Remove invalid record in dbSNP provisional map.

Provisional map have duplicated rs number or the chromsome in the new build can be "Unable to map"(UN), we need to clean this table.

(2) Use provisional map to update .map file

By joining .map file and this provisional map, we can obtain the new genome position in the new build.
Note: provisional map uses 1-based chromosomal index. Things will get tricker if we want to lift non-single site SNP e.g. AA/GG
Since provisional map provides a range in this case, it is necessary to know the genome position of that single base provided in the .map file,
and then we can look up the table, so it is not straigtforward.

(3) Adjust .map and .ped file

For those lifted dbSNP, we need to keep them in the .map files, otherwise, we need to delete them.
Accordingly, we need to deleted SNP genotypes for those cannot be lifted.

== Various reasons that lift over can fail ==

=== Genome position cannot be lifted ===
When a SNP resides in a contig that only exists in older reference build, liftOver cannot give it new genome.

You can try the following SNP (in BED format) in UCSC online liftOver site:
20 56737667 56737668 rs1073519
The error message will be: "Sequence intersects no chains"

=== SNP in higher build are located in non-referernce assembly ===
Some SNP are not in autosomes or sex chromosomes in NCBI build 37. dbSNP does not include them.
You cannot use dbSNP database to lookup its genome position by rs number.

Take rs1006094 as an example:
In NCBI dbSNP webpage, this SNP is reported as "Mapped unambiguously on non-reference assembly only"
Thus it is probably not very useful to lift this SNP.

=== rs number changed in newer dbSNP build ===
It is possible that new dbSNP build does not have certain rs numbers.
When dbSNp release new build, higher rs number may be merged to lower rs number because of those rs numbers are actually the same SNP.
This merge process can be complicate. For short description, see [[ #Use RsMergeArch and SNPHistory | Use RsMergeArch and SNPHistory ]].
For detail, see:

[http://www.ncbi.nlm.nih.gov/books/NBK44395/#FTP.do_you_have_a_table_of_merged_snps_s Finding Specific Data in dbSNP’s FTP Files]

[http://www.ncbi.nlm.nih.gov/books/NBK44468/#Build.can_two_id_numbers_correspond_to_t Merging RefSNP Numbers and RefSNP Clusters]

For example:

rs3001 has merged to rs2032.

=== Different dbSNP build ===
NCBI released dbSNP132 (VCF format), and UCSC also have their version of dbSNP132 (plain txt).
The two database files differ not only in file format, but in content.

For NCBI release, its [[#Resources| release]] will not contain:
* SNPs listed as microsatellites or named variations
* SNPs with multibyte alleles and unknown (N) adjacent base pairs
* SNPs that are not mapped on the reference genome (GRCh37)

For UCSC release, see [[ #Resources | UCSC dbSNP track note ]]

Use rs1054140 as an example:

NCBI dbSNP website gives 1 location:
[http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1054140 Link]

NCBI dbSNP VCF file has NO record.

UCSC genome browser website gives 2 locations:
[http://genome.ucsc.edu/cgi-bin/hgTracks?clade=mammal&org=Human&db=hg19&position=rs1054140&hgt.suggest=&hgt.suggestTrack=knownGene&pix=800&Submit=submit&hgsid=205770459&hgt.newJQuery=1 Link]

UCSC dbSNP file give 2 locations:
<pre>
721 chr10 17842693 17842694 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
723 chr10 18089681 18089682 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
</pre>

== Resources ==
* liftRsNumber.py [[liftRsNumber.py]] and its interal location: /net/fantasia/home/zhanxw/amd/analyze/verifyBamID/liftRsNumber.py
* liftMap.py [[liftMap.py]]
* NCBI provisional map [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.txt.gz file] and [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.info info]
* NCBI RgMergeArch [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/RsMergeArch.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=RsMergeArch schema]
* NCBI SNPHistory [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/SNPHistory.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPHistory schema]
* NCBI SNPChrPosOnRef build 132 [http://qbrc.swmed.edu/zhanxw/software/liftOver/b132_SNPChrPosOnRef_37_1.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPChrPosOnRef schema]. If this link becomes unavailable, please consider using this updated file ([ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/database/organism_data/b144_SNPChrPosOnRef_107.bcp.gz link]).

* How UCSC dbSNP differs from NCBI dbSNP [http://genomewiki.ucsc.edu/index.php/DbSNP_Track_Notes UCSC dbSNP track note]
* The dbSNP mapping process [http://www.ncbi.nlm.nih.gov/books/NBK44455/ link]
* NCBI dbSNP release 132 [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/VCF/v4.0/ByChromosomeNoGeno/00-All.vcf.gz 00-All.vcf.gz]
* UCSC dbSNP release 132 [http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz snp132.txt.gz]

== Third party contribution ==

* liftOver for BEDPE format by [mailto://doug.phanstiel@gmail.com Doug] - [http://www.sharedproteomics.com/forum/showthread.php?2930-liftOverBedpe.py-Convert-genome-coordinates-for-bedpe-files&p=4854#post4854 liftOverBedpe.py].

== Acknowledge ==

* Hyun: provides sample liftOver tool: [/net/wonderland/home/hmkang/prj/Sardinia/MetaboChip/scripts/j01-liftover-metabochip-positions.pl]
* Alex: careful examines of 0-based index in UCSC data file
* Adrian: explaination of SNPs omitted in NCBI dbSNP file
* Goncalo: all other supports

== Questions and Comments ==

Please contact [mailto:zhanxw@umich.edu Xiaowei Zhan].

RareMETALS

2015-03-25T02:37:08Z

Zhanxw: /* A Simple Tutorial for Using the rareMETALS.range function */

rareMETALS is an R-package for performing single or gene-level tests for detecting rare variant associations. For questions regarding the use of this package, please contact Dajiang Liu (dajiang at umich dot edu) or Gonçalo Abecasis (goncalo at umich dot edu). The same methodology is also implemented in command line tools. Please see [http://genome.sph.umich.edu/wiki/Rare-Metal]

== Change Log ==
* 01/24/2015 Version 5.8 is released, which fixed a serious bug for single variant unconditional association tests with group file. If you happen to run the analyses using rareMETALS.single.group() in version 5.7, the results are likely to be incorrect. Please rerun using version 5.8. Please note only rareMETALS.single.group function is affected. All other functions should not be affected by this error.
* 01/04/2015 Version 5.7 is released, which added metrics for heterogeneity of genetic effects, including I2 and Q for single variant association statistics
* 12/09/2014 Version 5.6 is released, which added function conditional.rareMETALS.range.group, and fixed a minor issue for estimating sample sizes.
* 11/19/2014 Version 5.5 is released, which fixes a few bugs on the version 5.4.
* 11/09/2014 Version 5.4 is posted with the following change 1.) Allowing for performing conditional analysis for multiple candidate variants 2.) add option correctFlip to rareMETALS.single.group, rareMETALS.range.group allowing for options to discard sites with non-matching ref or alt alleles. Default is TRUE
* 09/08/2014 Version 5.2 is posted. One change in version 5.0 and 5.1 is reverted, which could lead to undesirable effect. It improves on some border line cases as compared to Versions 4.7 - 4.9. But in general, version 5.2 and 4.7-4.9 should give very comparable results. Please update to the latest version. I would expect that version 5.2 should run stably for all models under all circumstances.
* 08/21/2014 Version 4.9 is posted. A bug is fixed for VT test. While the p-values and statistics were correct, the number of sites and the beta estimate could sometimes be incorrect in version 4.8. Now it is fixed. Please download the newest version. Thanks!
* 08/18/2014 Version 4.8 is posted. A bug for recessive model analysis is fixed. Additive and dominant models should remain unaffected. Thanks!
* 08/06/2014 Version 4.7 is posted, where a few minor bugs were fixed. Thanks to Heather Highland and Xueling Sim for careful testing!! Please update. Thanks!
* 07/15/2014 Fixed a bug in conditional.rareMETALS.single and conditional.rareMETALS.range; Please update. Thanks!
* 06/27/2014 Updated to version 4.0: Many updates are implemented, including support for group files in both single variant and gene-level association test; checks for allele flips based upon variant frequency, the detection of possible allele flips using a novel statistic based upon variations of allele frequency between studies;

== Where to download ==

The R package can be downloaded from [[Media:rareMETALS_5.8.tar.gz | rareMETALS_5.8.tar.gz]]. It will be eventually released on the Comprehensive R-archive Network. If you want to perform gene-level association test using automatically generated annotations, you will also need [[Media:refFlat_hg19.txt.gz | refFlat_hg19.txt.gz]], which is the gene definition modified from refFlat.

== How to install ==

To install the package, please use "R CMD INSTALL rareMETALS_XXX.tar.gz" command, where XXX is the version number for rareMETALS

== Supported Functionalities ==
* Marginal meta-analysis of single variant or gene-level association test
* Conditional analysis of single variant or gene-level association, for variants (gene) where there are covariance information available between candidate variants and known variants.
* Estimates of genetic effects and locus genetic variance
* Estimate measures of genetic effect heterogeneities between studies

== Exemplar Dataset==

Four datasets are useful to get you started on how to use rareMETALS R package for meta-analyses of gene-level association test

[[Media:study1.MetaScore.assoc.gz]] [[Media:study2.MetaScore.assoc.gz]] [[Media:study1.MetaCov.assoc.gz]] [[Media:study2.MetaCov.assoc.gz]]

== How to Generate Summary Association Statistics and Prepare Them for Meta-analysis ==

Meta-analysis summary association statistics can be generated by both RVTESTS and RAREMETALWORKER. Please refer to their documentations for generating summary association statistics

Once you have generated summary association statistics, you need to compress them with bgzip, and index them with tabix. If you use RAREMETALWORKER, the command should be like

'''NOTE: Tabix 1.X does not seem to support the indexing for generic tab-delimited files. To index the file, please use tabix 0.2.5 or earlier versions.

If you use RVTESTS, your command should be

bgzip study1.MetaScore.assoc

tabix -s 1 -b 2 -e 2 -S 1 study1.MetaScore.assoc.gz

tabix -s 1 -b 2 -e 2 -S 1 study1.MetaCov.assoc.gz

== A Simple Tutorial for Using the rareMETALS.single function ==

rareMETALS.single function allow you to perform meta-analyses for single variant association tests. The summary association statistics are combined using Mantel Haenszel test statistic. The details are described in our method paper Liu et al, Nat Genet, 2014.

Assume that you have a set of single variant score statistics and their covariance matrices.

Example:

library(rareMETALS)
cov.file <- c("study1.MetaCov.assoc.gz","study2.MetaCov.assoc.gz");
score.stat.file <- c("study1.MetaScore.assoc.gz","study2.MetaScore.assoc.gz")
res <- rareMETALS.single(score.stat.file,cov.file=NULL,range="19:11200093-11201275",alternative="two.sided",ix.gold=1,callrate.cutoff=0,hwe.cutoff=0);

== A Simple Tutorial for Using the rareMETALS.range function ==

res <- rareMETALS.range(score.stat.file,cov.file,range="19:11200093-11201275",range.name="LDLR",test = "GRANVIL",maf.cutoff = 0.05,alternative = c("two.sided"),ix.gold = 1,out.digits = 4,callrate.cutoff = 0,hwe.cutoff = 0,max.VT = NULL)
print(res$res.out)

<pre>
gene.name.out p.value.out statistic.out no.site.out beta1.est.out
[1,] "LDLR" "0.6064" "0.2654" "25" "-0.01729"
beta1.sd.out maf.cutoff.out direction.burden.by.study.out
[1,] "0.03357" "0.05" "--"
direction.meta.single.var.out top.singlevar.pos top.singlevar.refalt
[1,] "---++-+--+-+++++--+++++-+" "19:11200431" "C/T"
top.singlevar.pval top.singlevar.af
[1,] "0.004709" "0.01038"
pos.ref.alt.out

[1,] "19:11200093/T/C,19:11200213/G/A,19:11200235/G/A,19:11200272/C/A,19:11200282/G/A,19:11200309/C/A,19:11200412/C/T,19:11200419/C/T,19:11200431/C/T,19:1120\
0442/G/A,19:11200475/C/G,19:11200508/G/A,19:11200514/C/T,19:11200557/G/A,19:11200579/C/T,19:11200728/C/T,19:11200753/T/C,19:11200754/G/A,19:11200806/C/T,19:1\
1200839/T/A,19:11200840/C/A,19:11200896/C/T,19:11201259/G/C,19:11201274/C/T,19:11201275/A/T"
</pre>

More detailed results can be found in a list res$res.list

RareMETALS

2015-03-25T02:16:24Z

Zhanxw:

rareMETALS is an R-package for performing single or gene-level tests for detecting rare variant associations. For questions regarding the use of this package, please contact Dajiang Liu (dajiang at umich dot edu) or Gonçalo Abecasis (goncalo at umich dot edu). The same methodology is also implemented in command line tools. Please see [http://genome.sph.umich.edu/wiki/Rare-Metal]

== Change Log ==
* 01/24/2015 Version 5.8 is released, which fixed a serious bug for single variant unconditional association tests with group file. If you happen to run the analyses using rareMETALS.single.group() in version 5.7, the results are likely to be incorrect. Please rerun using version 5.8. Please note only rareMETALS.single.group function is affected. All other functions should not be affected by this error.
* 01/04/2015 Version 5.7 is released, which added metrics for heterogeneity of genetic effects, including I2 and Q for single variant association statistics
* 12/09/2014 Version 5.6 is released, which added function conditional.rareMETALS.range.group, and fixed a minor issue for estimating sample sizes.
* 11/19/2014 Version 5.5 is released, which fixes a few bugs on the version 5.4.
* 11/09/2014 Version 5.4 is posted with the following change 1.) Allowing for performing conditional analysis for multiple candidate variants 2.) add option correctFlip to rareMETALS.single.group, rareMETALS.range.group allowing for options to discard sites with non-matching ref or alt alleles. Default is TRUE
* 09/08/2014 Version 5.2 is posted. One change in version 5.0 and 5.1 is reverted, which could lead to undesirable effect. It improves on some border line cases as compared to Versions 4.7 - 4.9. But in general, version 5.2 and 4.7-4.9 should give very comparable results. Please update to the latest version. I would expect that version 5.2 should run stably for all models under all circumstances.
* 08/21/2014 Version 4.9 is posted. A bug is fixed for VT test. While the p-values and statistics were correct, the number of sites and the beta estimate could sometimes be incorrect in version 4.8. Now it is fixed. Please download the newest version. Thanks!
* 08/18/2014 Version 4.8 is posted. A bug for recessive model analysis is fixed. Additive and dominant models should remain unaffected. Thanks!
* 08/06/2014 Version 4.7 is posted, where a few minor bugs were fixed. Thanks to Heather Highland and Xueling Sim for careful testing!! Please update. Thanks!
* 07/15/2014 Fixed a bug in conditional.rareMETALS.single and conditional.rareMETALS.range; Please update. Thanks!
* 06/27/2014 Updated to version 4.0: Many updates are implemented, including support for group files in both single variant and gene-level association test; checks for allele flips based upon variant frequency, the detection of possible allele flips using a novel statistic based upon variations of allele frequency between studies;

== Where to download ==

The R package can be downloaded from [[Media:rareMETALS_5.8.tar.gz | rareMETALS_5.8.tar.gz]]. It will be eventually released on the Comprehensive R-archive Network. If you want to perform gene-level association test using automatically generated annotations, you will also need [[Media:refFlat_hg19.txt.gz | refFlat_hg19.txt.gz]], which is the gene definition modified from refFlat.

== How to install ==

To install the package, please use "R CMD INSTALL rareMETALS_XXX.tar.gz" command, where XXX is the version number for rareMETALS

== Supported Functionalities ==
* Marginal meta-analysis of single variant or gene-level association test
* Conditional analysis of single variant or gene-level association, for variants (gene) where there are covariance information available between candidate variants and known variants.
* Estimates of genetic effects and locus genetic variance
* Estimate measures of genetic effect heterogeneities between studies

== Exemplar Dataset==

Four datasets are useful to get you started on how to use rareMETALS R package for meta-analyses of gene-level association test

[[Media:study1.MetaScore.assoc.gz]] [[Media:study2.MetaScore.assoc.gz]] [[Media:study1.MetaCov.assoc.gz]] [[Media:study2.MetaCov.assoc.gz]]

== How to Generate Summary Association Statistics and Prepare Them for Meta-analysis ==

Meta-analysis summary association statistics can be generated by both RVTESTS and RAREMETALWORKER. Please refer to their documentations for generating summary association statistics

Once you have generated summary association statistics, you need to compress them with bgzip, and index them with tabix. If you use RAREMETALWORKER, the command should be like

'''NOTE: Tabix 1.X does not seem to support the indexing for generic tab-delimited files. To index the file, please use tabix 0.2.5 or earlier versions.

If you use RVTESTS, your command should be

bgzip study1.MetaScore.assoc

tabix -s 1 -b 2 -e 2 -S 1 study1.MetaScore.assoc.gz

tabix -s 1 -b 2 -e 2 -S 1 study1.MetaCov.assoc.gz

== A Simple Tutorial for Using the rareMETALS.single function ==

rareMETALS.single function allow you to perform meta-analyses for single variant association tests. The summary association statistics are combined using Mantel Haenszel test statistic. The details are described in our method paper Liu et al, Nat Genet, 2014.

Assume that you have a set of single variant score statistics and their covariance matrices.

Example:

library(rareMETALS)
cov.file <- c("study1.MetaCov.assoc.gz","study2.MetaCov.assoc.gz");
score.stat.file <- c("study1.MetaScore.assoc.gz","study2.MetaScore.assoc.gz")
res <- rareMETALS.single(score.stat.file,cov.file=NULL,range="19:11200093-11201275",alternative="two.sided",ix.gold=1,callrate.cutoff=0,hwe.cutoff=0);

== A Simple Tutorial for Using the rareMETALS.range function ==

res <- rareMETALS.range(score.stat.file,cov.file,range="19:11200093-11201275",range.name="LDLR",test = "GRANVIL",maf.cutoff = 0.05,alternative = c("two.sided"),ix.gold = 1,out.digits = 4,callrate.cutoff = 0,hwe.cutoff = 0,max.VT = NULL)

gene.name.out p.value.out statistic.out no.site.out beta1.est.out
[1,] "LDLR" "0.6064" "0.2654" "25" "-0.01729"
beta1.sd.out maf.cutoff.out direction.burden.by.study.out
[1,] "0.03357" "0.05" "--"
direction.meta.single.var.out top.singlevar.pos top.singlevar.refalt
[1,] "---++-+--+-+++++--+++++-+" "19:11200431" "C/T"
top.singlevar.pval top.singlevar.af
[1,] "0.004709" "0.01038"
pos.ref.alt.out \
\
[1,] "19:11200093/T/C,19:11200213/G/A,19:11200235/G/A,19:11200272/C/A,19:11200282/G/A,19:11200309/C/A,19:11200412/C/T,19:11200419/C/T,19:11200431/C/T,19:1120\
0442/G/A,19:11200475/C/G,19:11200508/G/A,19:11200514/C/T,19:11200557/G/A,19:11200579/C/T,19:11200728/C/T,19:11200753/T/C,19:11200754/G/A,19:11200806/C/T,19:1\
1200839/T/A,19:11200840/C/A,19:11200896/C/T,19:11201259/G/C,19:11201274/C/T,19:11201275/A/T"

More detailed results can be found in a list res$res.list

RvTests

2015-01-14T02:59:25Z

Zhanxw:

[[Category:Software]]

= Overview =

'''(See [[http://zhanxw.github.io/rvtests rvtests]] for more powerful rare-variant association test software and/or preparation for meta-analysis)'''

A few rare variants tests (Li-Leal's CMC and Madsen-Browning's weighted method) are implemented in the logisitc regression framework using C++. Please contact Youna Hu (youna@umich.edu) for comments, suggestions or questions.

The source code is located at [[File:RV3Tests.v1.tar]]

You can just download the tar file, extract it and go to the RV3Test.v1 folder and type make all to compile the code, the binary file will then be in the exectuables folder.

= Example =

See a detailed [[example]] here.

= Syntax =

This software uses command line interface as follows

RARE VARIANT ANALYSIS OPTIONS:
GENOTYPE : --genofile [pos.012],
--geneList [outGeneSorted.txt], --cutoff [0.010],
--collapseChoice [or]
PHENOTYPE : --phenofile [LDL.y.ID]
COVARIATES : --covConsider, --covfile [covFile.ID.2.txt]
PERMUTATION : --nPermute [10], --PermutationSeed [1]
GENE LEVEL TEST RESULT : --geneGlobalTestOut [globalPermuteSummary.txt],
--geneTestpvalueFile [geneTestPvalues.txt]

;GENOTYPE

;--genofile: A genotype 012 matrix (.012 is the file) This file can be prepared by using the prepare012s
source code [[File:vcfReader.v1.tar]]
Again, extract the tar file and then go into the directory to type make all to compile the code,
you then will find binary file in the executables folder.

Note: You should use Yanming's vcf annotation [http://genome.sph.umich.edu/wiki/VcfCodingSnps] on your vcf file first
to output a annotated vcf file. You SHOULD keep the log file from the annotation, which will be used to create the gene list.

Data File PREPARATION
Input files : --vcf [LDL.test.vcf], --log [], --IDfile []
Subsetting choices : --All
Output files : --outputPrefix [subsetGeno],
--outputGeneList [LDL.geneList.txt]
--vcf: Input vcf file
--log: This is the log file from Yanming's annotation output, we use this log to obtain the gene list
--IDfile specifies a file with one column of subject IDs to subject from the vcf file.
If it is not specified, then all subjects are included for the format conversion.
--All: specifies 1 to include all variants and 0 to include only nonsyn and stop annotated variants.
-- outputPrefix: Specify the prefix for the four output files which will be used in rvTests
*.012: A genotype matrix with subjects as rows and variant sites as columns.
*.012.pos: Chromosome and position numbers.
*.012.indv: Subject IDs.
*.012.frq: The frequency of the included variants.
--outputGeneList: Specify a file to store the gene list which will be used in rvTest.
The list file looks like this
1 OR4F5 69090 70008
1 SAMD11 860529 871276
1 NOC2L 879583 893918
1 KLHL17 895966 901095
1 PLEKHN1 901876 910482
1 C1orf170 910578 912021

;--geneList: This file is an output from prepare012s using the option --outputGeneList with columns as chromosome number, gene Name, start position, end position. There should be no header for this file.

THE CHROMOSOME NUMBERS SHOULD BE NUMERICS!!!! 1 - chromosome 1, DO NOT USE chr1.

;--cutoff: This is the minor allele frequency, you can specify it as 0.01, 0.05 or etc.
;--collapseChoice: Specify one of {or,sum,wt}. or: Li-Leal's CMC test, sum: Use the number of rare variants for each subject as the score, wt: Madeson-Browning's weighted rare variant score.

;PHENOTYPE
;--phenofile: A file where the first column is subject ID and the second column is phenotype (0 or 1).

;COVARIATES
;--covConsider: Default = 0, no covariate is considered. 1. covariate is considered.
;--covfile: Covariate file with the first column as subject ID and the other columns are covariates needed to be considered in the model.

;PERMUTATION
;--nPermute: Number of permutation for the evaluation of p values.
;-- PermutationSeed: Default = 1. Can be changed to other numbers too.

;GENE LEVEL TEST RESULT:
;--geneGlobalTestOut: This file stores the 5% and 95% quantiles of the p values for all the genes at each permutation
;--geneTestPvalueFile: This file gives you the gene name, number of rare variants, count of variants in case/control and p values from the RV test specified by collapseChoice.

Rvtests

2015-01-14T02:57:33Z

Zhanxw:

'''(Updated: January 2015; latest documentation can always be found at: [http://zhanxw.github.io/rvtests/ http://zhanxw.github.io/rvtests/])'''

= Introduction =

Rvtests, which stands for Rare Variant tests, is a flexible software package for genetic association studies. It is designed to support unrealted individual or related (family-based) individuals. Both quantitative trait and binary trait are supported. It includes a variety of association tests (e.g. single variant score test, burden test, variable threshold test, SKAT test, fast linear mixed model score test). It takes [http://www.1000genomes.com/ VCF] format as genotype input file and takes PLINK format phenotype file and covariate file. From our practice, it is capable to analyze 8,000 related individuals using less than 400 Mb memory.

= Download =

Source files can be downloaded from [https://github.com/zhanxw/rvtests/archive/master.zip github] or [https://github.com/zhanxw/rvtests github page].<br />Executable binary files (for Linux 64bit) can be downloaded from [https://github.com/zhanxw/rvtests/releases/download/v1.8.6/rvtests-20150104.tar.gz here].

= Quick Tutorial =

Here is a quick example of how to use ''rvtests'' software in typical use cases.

== Single variant tests ==

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --single wald,score</pre>
This specifies single variant Wald and score test for association<br />tests for every variant in the <code>input.vcf</code> file. The 6th column of the phenotype file, <code>phenotype.ped</code>, which is in PLINK format, is used. Rvtests will automatically check whether the phenotype is binary trait or quantitative trait.<br />For binary trait, the recommended way of coding is to code controls as 1, cases as 2, missing phenotypes as -9 or 0.

For other types of association tests, you can refer to [[#models|Models]]

== Groupwise tests ==

Groupwise tests includes three major kinds of tests.

* Burden tests: group variants, which are usually less than 1% or 5% rare variants, for association tests. The category includes: CMC test, Zeggini test, Madsen-Browning test, CMAT test, and rare-cover test.<br />
* Variable threshold tests: group variants under different frequency thresholds.<br />
* Kernel methods: suitable to tests rare variants having different directions of effects. These includes SKAT test and KBAC test.

All above tests requires to group variants into a unit. The simplist case is to use gene as grouping unit. For different grouping method, see [[#Grouping|Grouping]].

To perform rare variant tests by gene, you need to use <code>--geneFile</code> to specify the gene range in a refFlat format. We provided different gene definitions in the [[#Resources|Resources]] section. You can use <code>--gene</code> to specify which gene(s) to test. For example, specify <code>--gene CFH,ARMS2</code> will perform association tests on CFH and ARMS2 genes. If there is no providing <code>--gene</code> option, all genes will be tests.

The following command line demonstrate how to use CMC method, variable threshold method(proposed by Price) and kernel based method (SKAT by Shawn Lee and KBAC by<br />Dajiang Liu) to test every gene listed in ''refFlat_hg19_uniq_gene.txt.gz''.

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --geneFile refFlat_hg19_uniq_gene.txt.gz --burden cmc --vt price --kernel skat,kbac</pre>
== Related individual tests ==

To test related individuals, you will need to first create a kinship matrix:

<pre>vcf2kinship --inVcf input.vcf --bn --out output</pre>
The option <code>--bn</code> means calculating empirical kinship using Balding-Nicols method. You can specifiy <code>--ibs</code> to obtain IBS kinship or use <code>--pedigree input.ped</code> to calculate kinship from known pedigree information.

Then you can use linear mixed model based association tests such as Fast-LMM score test, Fast-LMM LRT test and Grammar-gamma tests. An exemplar command is shown:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --kinship output.kinship --single famScore,famLRT,famGrammarGamma</pre>
== Meta-analysis tests ==

The meta-analysis models outputs association test results and genotype covariance matrix. These statistics can be used in rare variant association analysis.<br />We provide single variant score test and generate genotype covariance matrix.<br />You can use command:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output</pre>
In a more realistic scenario, you may want to adjust for covariates and want to inverse normalized residuals obtained in null model ([http://www.nature.com/ng/journal/vaop/ncurrent/full/ng.2852.html link] to our methodology paper), then this command will work:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out output</pre>
Here the <code>--covar</code> specify a covariate file, and <code>--covar-name</code> specify which covariates can used in the analysis. Covariate file format can be found [[#Covariate%20file|here]]. <code>--inverseNormal --useResidualAsPhenotype</code> specifies trait transformation method. That means first fit a regression model of the phenotype on covariates (intercept automatically added), then the residuals are inverse normalized. Trait transformation details can be found [[#Trait%20transformation|here]].

We support both unrelated individuals and related indivudlas (e.g. family data). You need to append <code>--kinship input.kinship</code> to the command line:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output --kinship input.kinship</pre>
The file <code>input.kinship</code> is calculated by <code>vcf2kinship</code> program, and usage to this program is described in [[#related-individual-tests|Related individual tests]].

=== Dominant models and recessive models ===

Dominant and recessive disease models are supported by appending "dominant" and/or "recessive" after "--meta" option. For example, use "--meta dominant,recessive" will<br />generate two sets of files. For dominant model, they are "prefix.MetaDominant.assoc" and "prefix.MetaDominantCov.assoc.gz"; for recessive model,<br />they are "prefix.MetaRecessive.assoc" and "prefix.MetaRecessiveCov.assoc.gz". Internally, in dominant models, genotypes 0/1/2 are coded as 0/1/1; in recessive models, genotypes 0/1/2 are<br />coded as 0/0/1. Missing genotypes will be imputed to the mean.

= Input files =

== Genotype file (VCF) ==

Rvtests supports VCF (Variant Call Format) files. Files in both plain txt format or gzipped format are supported. To use group-based rare variant tests, indexed the VCF files using [http://samtools.sourceforge.net/tabix.shtml tabix] are required.

Here are the commands to convert plain text format to bgzipped VCF format:

<pre>(grep ^"#" $your_old_vcf; grep -v ^"#" $your_old_vcf | sed 's:^chr::ig' | sort -k1,1n -k2,2n) | bgzip -c > $your_vcf_file
tabix -f -p vcf $your_vcf_file</pre>
The above commands will (1) remove the <code>chr</code> prefix from chromosome names; (2) sort VCF files by chromosome first, then by chromosomal positions; (3) compress using bgzip; (4) create tabix index.

Rvtests support genotype dosages. Use <code>--dosage DosageTag</code> to specify the dosage tag. For example, if VCF format field is "GT:EC" and individual genotype fields is "0/0:0.02", you can use <code>--dosage EC</code>, and rvtests will use the dosage 0.02 in the regression models.

== Phenotype file ==

You can use <code>--mpheno $phenoypeColumnNumber</code> or <code>--pheno-name</code> to specify a given phenotype.

An example phenotype file, (<code>example.pheno</code>), has the following format:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.7642934435605 -0.733862638327895 -0.980843608339726 2
P2 P2 0 0 0 0.457111744989746 0.623297281416372 -2.24266162284447 1
P3 P3 0 0 0 0.566689682543218 1.44136462889459 -1.6490100777089 1
P4 P4 0 0 0 0.350528353203767 -1.79533911725537 -1.11916876241804 1
P5 P5 0 0 1 2.72675074738545 -1.05487747371158 -0.33586430010589 2</pre>
Phenotype file is specified by the option <code>--pheno example.pheno</code> . The default phenotype column header is “<code>y1</code>”. If you want to use alternative columns as phenotype for association analysis (e.g the column with header y2), you may specify the phenotype by column or by name using either

* --mpheno 2<br />
* --pheno-name y2

'''NOTE:''' to use “<code>--pheno-name</code>”, the header line must starts with “<code>fid iid</code>” as PLINK requires.

In phenotype file, missing values can be denoted by NA or any non-numeric values. Individuals with missing phenotypes will be automatically dropped from subsequent association analysis. For each missing phenotype value, a warning will be generated and recorded in the log file.

When the phenotype values are only 0, 1 and 2, rvtests will automatically treat it as binary traits. However, if you want to treat it as continuous trait, please use "<code>--qtl</code>" option.

== Covariate file ==

You can use <code>--covar</code> and <code>--covar-name</code> to specify covariates that will be used for single variant association analysis. This is an optional parameter. If you do not have covariate in the data, this option can be ignored.

The covariate file, (e.g. <code>example.covar</code>) has a similar format as the phenotype file:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.911 -1.465 -0.817 1
P2 P2 0 0 0 2.146 -2.451 -0.178 2
P3 P3 0 0 0 1.086 -1.194 -0.899 1
P4 P4 0 0 0 0.704 -1.052 -0.237 1
P5 P5 0 0 1 2.512 -3.085 -2.579 1</pre>
The covariate file is specified by the <code>--covar</code> option (e.g. <code>--covar example.covar</code>). To specify covariates that will be used in the association analysis, the option <code>--covar-name</code> can be used. For example, when age, bmi and 3 PCs are used for association analysis, the following option can be specified for the rvtest program, i.e.<br /><code>--covar example.covar --covar-name age,bmi,pc1,pc2,pc3</code>.

Note: Missing data in the covariate file can be labeled by any non-numeric value (e.g. NA). They will be automatically imputed to the mean value in the data file.

== Trait transformation ==

In this meta-analysis, we use inversed normal transformed residuals in the association analysis, which is achieved by using a combination of <code>--inverseNormal</code> and <code>--useResidualAsPhenotype</code>. Specifically, we first fit the null model by regressing phenotype on covariates. The residuals are then inverse normal transformed (see Appendix A more detailed formulae for transformation). Transformed residuals will be used to obtain score statistics.

In meta analysis, an exemplar command for using rvtest looks like the following:

<pre>./rvtest --inVcf $vcf --pheno $example.pheno --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out $output_prefix </pre>
= Models =

Rvtests support various association models.

== Single variant tests ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only null model is used to performed the test
|-
|Wald test
|align="center"|wald
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only fit alternative model, and effect size will be estimated
|-
|Exact test
|align="center"|exact
|align="center"|B
|align="center"|N
|align="center"|U
|Fisher's test
|-
|Fam LRT
|align="center"|famLRT
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model
|-
|Fam Score
|align="center"|famScore
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model style likelihood ratio test
|-
|Grammar-gamma
|align="center"|famGrammarGamma
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Grammar-gamma method
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--single score</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Burden tests ==

{|
!Burden tests
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants by Bingshan Li.
|-
|Zeggini
|align="center"|zeggini
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Aggregate counts of rare variants by Morris Zeggini.
|-
|Madsen-Browning
|align="center"|mb
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Madsen.
|-
|Fp
|align="center"|fp
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Danyu Lin.
|-
|Exact CMC
|align="center"|exactCMC
|align="center"|B
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Fisher's exact test.
|-
|CMC Wald
|align="center"|cmcWald
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Wald test.
|-
|RareCover
|align="center"|rarecover
|align="center"|B
|align="center"|N
|align="center"|U
|Find optimal grouping unit for rare variant tests by Thomas Hoffman.
|-
|CMAT
|align="center"|cmat
|align="center"|B
|align="center"|N
|align="center"|U
|Test non-coding variants by Matt Z.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--burden cmc</code> will apply CMC test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Variable threshold models ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Variable threshold model
|align="center"|vt
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Every rare-variant frequency cutoffs are tests by Alkes Price.
|-
|Variable threshold CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|This models is natiive so that it output CMC test statistics under all possible frequency cutoffs.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--vt price</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Kernel models ==

{|
!Kernel
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|SKAT
|align="center"|skat
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Sequencing kernel association test by Shawn Lee.
|-
|KBAC
|align="center"|kbac
|align="center"|B
|align="center"|N
|align="center"|U
|Kernel-based adaptive clustering model by Dajiang Liu.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--kernel skat</code> will apply SKAT test.<br />To further customize SKAT test, you can use ''--kernel skat[nPerm=100:alpha=0.001:beta1=1:beta2=20]'' to specify permutation counts, type-1 error,<br />beta distribution parameters for upweighting rare variants. Rvtests will output a message showing:

<pre>[INFO] SKAT test significance will be evaluated using 10000 permutations at alpha = 0.001 (beta1 = 1.00, beta2 = 20.00)</pre>
(##) In trait column, B and Q stand for binary, quantitiave trait.

== Meta-analysis models ==

{|
!Type
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|standard score tests
|-
|Dominant model
|align="center"|dominant
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under dominant disease model
|-
|Recessive model
|align="center"|recessive
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under recessive disease model
|-
|Covariance
|align="center"|cov
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|covariance matrix
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--meta score,cov</code> will generate score statistics and covariance matrix for meta-analysis.<br />(##) In trait column, B and Q stand for (b)inary, (q)uantitiave trait.

== Utility models ==

Rvtests has an usually option <code>--outputRaw</code>. When specify this, rvtests can output genotypes, phenotype, covariates(if any) and collapsed genotype to tabular files. These files can be imported into other software (e.g. R) for further analysis.

= Association test options =

== Sample inclusion/exclusion ==

Rvtests can flexibly specify which sample(s) to include or exclude:

<pre> --peopleIncludeID : give IDs of people that will be included in study
--peopleIncludeFile : from given file, set IDs of people that will be included in study
--peopleExcludeID : give IDs of people that will be included in study
--peopleExcludeFile : from given file, set IDs of people that will be included in study</pre>
<code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> are used to include/exclude samples from command line.<br />For example, specify <code>--peopleIncludeID A,B,C</code> will include A, B and C sample from the VCF files if they exists.<br /><code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> followed by a file name will include or exclude the IDs in the file.<br />So to include sample A, B and C, you can provide a file, <code>people.txt</code>, looks like:

<pre>A
B
C</pre>
Then use <code>--peopleIncludeFile people.txt</code> to include them in the analysis.

== Variant site filters ==

It is common that different frequency cutoffs are applied in rare-variant analysis.<br />Therefore, rvtests specify frequency cutoffs.

Frequency Cutoff

<pre> --freqUpper : Specify upper minor allele frequency bound to be included in analysis
--freqLower : Specify lower minor allele frequency bound to be included in analysis</pre>
If you specify <code>--freqLower 0.01 --freqUpper 0.05</code>, only the variants with minor allele ferquncy between 0.01 and 0.05 (boundary inclusive) will be analyzed.

Similar to sample inclusion/exclusion options, you can specify a range of variants to be included by<br />specifying <code>--rangeList</code> option. For example <code>--rangeList 1:100-200</code> will include the chromosome 1 position 100bp to 200bp region.<br />Alternatively, use a separate file, <code>range.txt</code>, and <code>--rangeFile range.txt</code> to speicify association tests range.

<pre> --rangeList : Specify some ranges to use, please use chr:begin-end format.
--rangeFile : Specify the file containing ranges, please use chr:begin-end format.
--siteFile : Specify the file containing sites to include, please use "chr pos" format.</pre>
It is supported to filter variant site by site depth, minor allele count or annotation (annotated VCF file is needed).

<pre> --siteDepthMin : Specify minimum depth(inclusive) to be incluced in analysis
--siteDepthMax : Specify maximum depth(inclusive) to be incluced in analysis
--siteMACMin : Specify minimum Minor Allele Count(inclusive) to be incluced in analysis
--annoType : Specify annotation type that is follwed by ANNO= in the VCF INFO field, regular expression is allowed</pre>
''NOTE'': <code>--annoType Nonsynonymous</code> will only analyze nonsynonymous variants where they have <code>ANNO=Nonsynonymous</code> in the INFO field.<br />VCF with annotatino information are called annotated VCF here. And to annotate<br />a VCF file, you can use [https://github.com/zhanxw/anno ANNO], a fast and accurate annotation software.

== Genotype filters ==

Genotype with low depth or low quality can be filtered out by these options:

<pre> --indvDepthMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis
--indvDepthMax : Specify maximum depth(inclusive) of a sample to be incluced in analysis
--indvQualMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis</pre>
When genotypes are filtered, they are marked as missing genotypes.<br />Consequently, samples with missing genotype may or may not be included in the analysis.<br />That means samples with genotypes may be dropped (<code>--impute drop</code>)<br />or may still be included (<code>--impute mean</code> or <code>--impute hwe</code>).<br />By default, genotypes are imputed to its means.<br />See next section about how you like to handle missing genotypes.

== Handle missing genotypes and phenotypes ==

When genotypes are missing (e.g. genotype = "./.") or gentoypes are filtered out,<br />there are three options to handle them: (1) impute to its mean(default option); (2) impute by HWE equilibrium; (3) remove from the model.<br />Use <code>--impute [mean|hwe|drop]</code> to specify which option to use.

When quantitative phenotypes are missing, for example, some samples have gneotype files, but not phenotypes,<br />rvtests can impute missing phenotype to its mean.

''NOTE:'' Do not use <code>--imputePheno</code> for binary trait.

In summary, the following two options can be used:

<pre> --impute : Specify either of mean, hwe, and drop
--imputePheno : Impute phenotype to mean by those have genotypes but no
phenotpyes

</pre>
== Specify groups (e.g burden unit) ==

Rare variants association tests are usually performed in gruops of variants.<br />The natural grouping unit is gene. Rvtests can read gene definition file in <code>refFlat</code> format,<br />and perform association for each gene. Use <code>--geneFile</code> option to specify the gene file name.<br />For example, <code>--geneFile refFlat_hg19.txt.gz</code> will use <code>refFlat_hg19.txt.gz</code> as gene definition file,<br />and then perform association tests for every gene. Use <code>--gene</code> to specify a subset of genes to test.<br />For example, <code>--gene CFH</code> will only test CFH gene.

Alternative grouping unit can be specified as ''set''.<br />These ''sets'' are treated similar to gene.<br />You can thus use <code>--setFile</code> to define sets (similar to <code>--geneFile</code> option),<br />and use <code>--set</code> to define a specific set (similar to <code>--gene</code> option).<br />Additionally, use <code>--setList</code> can speicify a set to test from command line.

The format of a set file is: (1) set names; (2) ranges (e.g. chrom:begin-end);<br />For example, you have a set file, <code>example.set</code>, like this:

<pre>set1 1:100-200,1:250-300
set2 2:500-600</pre>
You can specify <code>--setFile example.set --set set2</code> to group variants<br />within chromosome 2, position 500 to 600bp.<br />If you want to test a particular region, for example, chromosome 2, position 500 to 550bp,<br />but do not want to make another file, you can use <code>--setList 2:500-600</code>.

In summary, options related to ''Grouping Unit'' are listed below:

<pre> --geneFile : specify a gene file (for burden tests)
--gene : specify which genes to test
--setList : specify a list to test (for burden tests)
--setFile : specify a list file (for burden tests, first two columns:
setName chr:beg-end)
--set : specify which set to test (1st column)</pre>
= Sex chromosome analysis =

Rvtests suppport X chromosome analysis. In human X chromosome, there is PAR (pseudoautosomal region) and non-PAR region.<br />For males, there are two X allele in PAR region and one allele in non-PAR region.<br />While the PAR region is treated in the same way as autosomes, rvtests treate non-PAR region differently.<br />Below we will describe the details about how rvtests handles non-PAR region.

''Prepare data''. According to VCF standard, male genotype needs to coded as 0 or 1. For compatibility, rvtests also support 0/0 or 1/1 coding.<br />In VCF files, male genotypes can be written as "0", "1", "0|0", "0/0", "1|1", "1/1". All other genotypes will be treated as missing.

''Genotype in the regression model''. For consistencmaine, male genotypes are converted to 0 or 2.

''MetaScore results''. If specify <code>--meta score</code>, the output file <code>prefix.MetaScore.assoc</code> includes PAR-region and non-PAR region analysis.<br />But in the non-PAR region, the difference is that Hardy-Weinberg P-value are calculated using female samples.

''Related individuals''. Just append <code>--xHemi</code> to the <code>vcf2kinship</code> (more details in [[#kinship-generation|Kinship generation]]) and <code>rvtest</code> command lines. Rvtests can recognize non-PAR region kinship and use it in the analysis.

''PAR region''. PAR region is defined as two regions X:60001-2699520 and X:154931044-155270560. Use <code>--xLabel</code> can speicify which chromosome has PAR region (default: 23|X)<br />and use <code>--xParRegion</code> to specify PAR region (default: hg19, meaning '60001-2699520,154931044-155260560' in the UCSC build hg19, specify "hg18" will use PAR region definition in the UCSC build hg18).

= Kinship generation =

Analysis of related individual usually requires estimation of kinship. You can a separate tool, <code>vcf2kinship</code>.<br /><code>vcf2kinship</code> is usually included in rvtests binary distribution or can be built from software source codes.

<code>vcf2kinship</code> can calcualte pedigree kinship using a pedigree input file (PED format, see [[#phenotype-file|Phenotype file]], use option <code>--ped</code>).<br />The output file name is specified by <code>--prefix</code> option. If you use <code>--prefix output</code> then the output files will include <code>output.kinship</code>.

It can also calculate empirical kinship using genotyp input file (VCF format, see [[#genotype-file-vcf|Genotype file (VCF)]], use option <code>--inVcf</code>).<br />For empiricial kinship, you also need to specify the kinship model, either Balding-Nicols model (ue option <code>--bn</code>) or Identity-by-state model (use option <code>--ibs</code>).

In sex chromosome analysis, it is often required to generate kinship on X chromsoome regions, then you need to speicfy <code>--xHemi</code>. If your input VCF file has different X chromosome label (e.g. chromosome name is '23' instead of 'X'), you can use <code>--xLabel 23</code>.

If principal component decomposition (PCA) results are needed, you can use option <code>--pca</code>. Then output files with suffix '.pca' include PCA results.

When dealing with large input files, it is often preferred to use multiple CPU to speed up calculation using the option <code>--thread N</code> in which N is the number of CPUs.

For example, to generate pedigree-based kinship (<code>--ped</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region, the command line is:

<pre>vcf2kinship --ped input.ped --xHemi --out output</pre>
To generate empirical kinship (<code>--inVcf</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region using Balding-Nicols model, the command line is:

<pre>vcf2kinship --inVcf input.vcf.gz --ped input.ped --bn --xHemi --out output</pre>
NOTE: you need to provide a pedigree file (PED) in the above case, as <code>vcf2kinship</code> need the sex information of samples.

= Frequently Asked Questions (FAQ) =

* Does rvtests suppport binary traits of related-individuals?

Not yet. It's a complex scenario and we have not found good solutions.

* Can you provide a list of command line options?

Rvtests have build help taht can be found by executing <code>rvtest --help</code>.<br />We also put all available options in this [https://github.com/zhanxw/rvtests/wiki/Command-Line-Options link].

* Can you provide standard error (SE) or confidence interval (CI) for the estimated Beta in the score model?

In the output of MetaScore model (--meta score), the standard error is the inverse of SQRT_V_STAT.<br />For example, if SQRT_V_STAT = 2, that means the standard error of estimated beta is 1/2 = 0.5.

* Why the INFORMATIVE_ALT_AC, N_REF and N_ALT columns have zero counts for certain chromosome X regions in meta-analysis models?

These counts are calculated from female individuals. If your study only has male samples, rvtests cannot report these counts. Because if a male carries a non-reference allele, we cannot conclude that this is heterozygous (0/1) site or homozygous alternatives (1/1) site.

= Feedback/Contact =

Questions and requests can be sent to Xiaowei Zhan<br />([mailto:zhanxw@umich.edu [mailto:zhanxw@umich.edu zhanxw@umich.edu]])<br />or Goncalo Abecasis<br />([mailto:goncalo@umich.edu [mailto:goncalo@umich.edu goncalo@umich.edu]])

Rvtests is a collaborative effort by Youna Hu, Bingshan Li, Dajiang Liu.

Rvtests

2015-01-12T19:43:08Z

Zhanxw:

'''Table of Contents'''

* [[#introduction|Introduction]]<br />
* [[#download|Download]]<br />
* [[#quick-tutorial|Quick Tutorial]]
** [[#single-variant-tests|Single variant tests]]<br />
** [[#groupwise-tests|Groupwise tests]]<br />
** [[#related-individual-tests|Related individual tests]]<br />
** [[#meta-analysis-tests|Meta-analysis tests]]
*** [[#dominant-models-and-recessive-models|Dominant models and recessive models]]<br />
* [[#input-files|Input files]]
** [[#genotype-file-vcf|Genotype file (VCF)]]<br />
** [[#phenotype-file|Phenotype file]]<br />
** [[#covariate-file|Covariate file]]<br />
** [[#trait-transformation|Trait transformation]]<br />
* [[#models|Models]]
** [[#single-variant-tests|Single variant tests]]<br />
** [[#burden-tests|Burden tests]]<br />
** [[#variable-threshold-models|Variable threshold models]]<br />
** [[#kernel-models|Kernel models]]<br />
** [[#meta-analysis-models|Meta-analysis models]]<br />
** [[#utility-models|Utility models]]<br />
* [[#association-test-options|Association test options]]
** [[#sample-inclusionexclusion|Sample inclusion/exclusion]]<br />
** [[#variant-site-filters|Variant site filters]]<br />
** [[#genotype-filters|Genotype filters]]<br />
** [[#handle-missing-genotypes-and-phenotypes|Handle missing genotypes and phenotypes]]<br />
** [[#specify-groups-eg-burden-unit|Specify groups (e.g burden unit)]]<br />
* [[#sex-chromosome-analysis|Sex chromosome analysis]]<br />
* [[#kinship-generation|Kinship generation]]<br />
* [[#frequently-asked-questions-faq|Frequently Asked Questions (FAQ)]]<br />
* [[#feedbackcontact|Feedback/Contact]]



[https://travis-ci.org/zhanxw/rvtests [[Image:https://travis-ci.org/zhanxw/rvtests.png?branch=master|Build Status]]]

(Updated: Janurary 2015)

= Introduction =

Rvtests, which stands for Rare Variant tests, is a flexible software package for genetic association studies. It is designed to support unrealted individual or related (family-based) individuals. Both quantitative trait and binary trait are supported. It includes a variety of association tests (e.g. single variant score test, burden test, variable threshold test, SKAT test, fast linear mixed model score test). It takes [http://www.1000genomes.com/ VCF] format as genotype input file and takes PLINK format phenotype file and covariate file. From our practice, it is capable to analyze 8,000 related individuals using less than 400 Mb memory.

= Download =

Source files can be downloaded from [https://github.com/zhanxw/rvtests/archive/master.zip github] or [https://github.com/zhanxw/rvtests github page].<br />Executable binary files (for Linux 64bit) can be downloaded from [https://github.com/zhanxw/rvtests/releases/download/v1.8.6/rvtests-20150104.tar.gz here].

= Quick Tutorial =

Here is a quick example of how to use ''rvtests'' software in typical use cases.

== Single variant tests ==

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --single wald,score</pre>
This specifies single variant Wald and score test for association<br />tests for every variant in the <code>input.vcf</code> file. The 6th column of the phenotype file, <code>phenotype.ped</code>, which is in PLINK format, is used. Rvtests will automatically check whether the phenotype is binary trait or quantitative trait.<br />For binary trait, the recommended way of coding is to code controls as 1, cases as 2, missing phenotypes as -9 or 0.

For other types of association tests, you can refer to [[#models|Models]]

== Groupwise tests ==

Groupwise tests includes three major kinds of tests.

* Burden tests: group variants, which are usually less than 1% or 5% rare variants, for association tests. The category includes: CMC test, Zeggini test, Madsen-Browning test, CMAT test, and rare-cover test.<br />
* Variable threshold tests: group variants under different frequency thresholds.<br />
* Kernel methods: suitable to tests rare variants having different directions of effects. These includes SKAT test and KBAC test.

All above tests requires to group variants into a unit. The simplist case is to use gene as grouping unit. For different grouping method, see [[#Grouping|Grouping]].

To perform rare variant tests by gene, you need to use <code>--geneFile</code> to specify the gene range in a refFlat format. We provided different gene definitions in the [[#Resources|Resources]] section. You can use <code>--gene</code> to specify which gene(s) to test. For example, specify <code>--gene CFH,ARMS2</code> will perform association tests on CFH and ARMS2 genes. If there is no providing <code>--gene</code> option, all genes will be tests.

The following command line demonstrate how to use CMC method, variable threshold method(proposed by Price) and kernel based method (SKAT by Shawn Lee and KBAC by<br />Dajiang Liu) to test every gene listed in ''refFlat_hg19_uniq_gene.txt.gz''.

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --geneFile refFlat_hg19_uniq_gene.txt.gz --burden cmc --vt price --kernel skat,kbac</pre>
== Related individual tests ==

To test related individuals, you will need to first create a kinship matrix:

<pre>vcf2kinship --inVcf input.vcf --bn --out output</pre>
The option <code>--bn</code> means calculating empirical kinship using Balding-Nicols method. You can specifiy <code>--ibs</code> to obtain IBS kinship or use <code>--pedigree input.ped</code> to calculate kinship from known pedigree information.

Then you can use linear mixed model based association tests such as Fast-LMM score test, Fast-LMM LRT test and Grammar-gamma tests. An exemplar command is shown:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --out output --kinship output.kinship --single famScore,famLRT,famGrammarGamma</pre>
== Meta-analysis tests ==

The meta-analysis models outputs association test results and genotype covariance matrix. These statistics can be used in rare variant association analysis.<br />We provide single variant score test and generate genotype covariance matrix.<br />You can use command:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output</pre>
In a more realistic scenario, you may want to adjust for covariates and want to inverse normalized residuals obtained in null model ([http://www.nature.com/ng/journal/vaop/ncurrent/full/ng.2852.html link] to our methodology paper), then this command will work:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out output</pre>
Here the <code>--covar</code> specify a covariate file, and <code>--covar-name</code> specify which covariates can used in the analysis. Covariate file format can be found [[#Covariate%20file|here]]. <code>--inverseNormal --useResidualAsPhenotype</code> specifies trait transformation method. That means first fit a regression model of the phenotype on covariates (intercept automatically added), then the residuals are inverse normalized. Trait transformation details can be found [[#Trait%20transformation|here]].

We support both unrelated individuals and related indivudlas (e.g. family data). You need to append <code>--kinship input.kinship</code> to the command line:

<pre>rvtests --inVcf input.vcf --pheno phenotype.ped --meta score,cov --out output --kinship input.kinship</pre>
The file <code>input.kinship</code> is calculated by <code>vcf2kinship</code> program, and usage to this program is described in [[#related-individual-tests|Related individual tests]].

=== Dominant models and recessive models ===

Dominant and recessive disease models are supported by appending "dominant" and/or "recessive" after "--meta" option. For example, use "--meta dominant,recessive" will<br />generate two sets of files. For dominant model, they are "prefix.MetaDominant.assoc" and "prefix.MetaDominantCov.assoc.gz"; for recessive model,<br />they are "prefix.MetaRecessive.assoc" and "prefix.MetaRecessiveCov.assoc.gz". Internally, in dominant models, genotypes 0/1/2 are coded as 0/1/1; in recessive models, genotypes 0/1/2 are<br />coded as 0/0/1. Missing genotypes will be imputed to the mean.

= Input files =

== Genotype file (VCF) ==

Rvtests supports VCF (Variant Call Format) files. Files in both plain txt format or gzipped format are supported. To use group-based rare variant tests, indexed the VCF files using [http://samtools.sourceforge.net/tabix.shtml tabix] are required.

Here are the commands to convert plain text format to bgzipped VCF format:

<pre>(grep ^"#" $your_old_vcf; grep -v ^"#" $your_old_vcf | sed 's:^chr::ig' | sort -k1,1n -k2,2n) | bgzip -c > $your_vcf_file
tabix -f -p vcf $your_vcf_file</pre>
The above commands will (1) remove the <code>chr</code> prefix from chromosome names; (2) sort VCF files by chromosome first, then by chromosomal positions; (3) compress using bgzip; (4) create tabix index.

Rvtests support genotype dosages. Use <code>--dosage DosageTag</code> to specify the dosage tag. For example, if VCF format field is "GT:EC" and individual genotype fields is "0/0:0.02", you can use <code>--dosage EC</code>, and rvtests will use the dosage 0.02 in the regression models.

== Phenotype file ==

You can use <code>--mpheno $phenoypeColumnNumber</code> or <code>--pheno-name</code> to specify a given phenotype.

An example phenotype file, (<code>example.pheno</code>), has the following format:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.7642934435605 -0.733862638327895 -0.980843608339726 2
P2 P2 0 0 0 0.457111744989746 0.623297281416372 -2.24266162284447 1
P3 P3 0 0 0 0.566689682543218 1.44136462889459 -1.6490100777089 1
P4 P4 0 0 0 0.350528353203767 -1.79533911725537 -1.11916876241804 1
P5 P5 0 0 1 2.72675074738545 -1.05487747371158 -0.33586430010589 2</pre>
Phenotype file is specified by the option <code>--pheno example.pheno</code> . The default phenotype column header is “<code>y1</code>”. If you want to use alternative columns as phenotype for association analysis (e.g the column with header y2), you may specify the phenotype by column or by name using either

* --mpheno 2<br />
* --pheno-name y2

'''NOTE:''' to use “<code>--pheno-name</code>”, the header line must starts with “<code>fid iid</code>” as PLINK requires.

In phenotype file, missing values can be denoted by NA or any non-numeric values. Individuals with missing phenotypes will be automatically dropped from subsequent association analysis. For each missing phenotype value, a warning will be generated and recorded in the log file.

When the phenotype values are only 0, 1 and 2, rvtests will automatically treat it as binary traits. However, if you want to treat it as continuous trait, please use "<code>--qtl</code>" option.

== Covariate file ==

You can use <code>--covar</code> and <code>--covar-name</code> to specify covariates that will be used for single variant association analysis. This is an optional parameter. If you do not have covariate in the data, this option can be ignored.

The covariate file, (e.g. <code>example.covar</code>) has a similar format as the phenotype file:

<pre>fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.911 -1.465 -0.817 1
P2 P2 0 0 0 2.146 -2.451 -0.178 2
P3 P3 0 0 0 1.086 -1.194 -0.899 1
P4 P4 0 0 0 0.704 -1.052 -0.237 1
P5 P5 0 0 1 2.512 -3.085 -2.579 1</pre>
The covariate file is specified by the <code>--covar</code> option (e.g. <code>--covar example.covar</code>). To specify covariates that will be used in the association analysis, the option <code>--covar-name</code> can be used. For example, when age, bmi and 3 PCs are used for association analysis, the following option can be specified for the rvtest program, i.e.<br /><code>--covar example.covar --covar-name age,bmi,pc1,pc2,pc3</code>.

Note: Missing data in the covariate file can be labeled by any non-numeric value (e.g. NA). They will be automatically imputed to the mean value in the data file.

== Trait transformation ==

In this meta-analysis, we use inversed normal transformed residuals in the association analysis, which is achieved by using a combination of <code>--inverseNormal</code> and <code>--useResidualAsPhenotype</code>. Specifically, we first fit the null model by regressing phenotype on covariates. The residuals are then inverse normal transformed (see Appendix A more detailed formulae for transformation). Transformed residuals will be used to obtain score statistics.

In meta analysis, an exemplar command for using rvtest looks like the following:

<pre>./rvtest --inVcf $vcf --pheno $example.pheno --covar example.covar --covar-name age,bmi --inverseNormal --useResidualAsPhenotype --meta score,cov --out $output_prefix </pre>
= Models =

Rvtests support various association models.

== Single variant tests ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only null model is used to performed the test
|-
|Wald test
|align="center"|wald
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Only fit alternative model, and effect size will be estimated
|-
|Exact test
|align="center"|exact
|align="center"|B
|align="center"|N
|align="center"|U
|Fisher's test
|-
|Fam LRT
|align="center"|famLRT
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model
|-
|Fam Score
|align="center"|famScore
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Fast-LMM model style likelihood ratio test
|-
|Grammar-gamma
|align="center"|famGrammarGamma
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|Grammar-gamma method
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--single score</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Burden tests ==

{|
!Burden tests
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants by Bingshan Li.
|-
|Zeggini
|align="center"|zeggini
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Aggregate counts of rare variants by Morris Zeggini.
|-
|Madsen-Browning
|align="center"|mb
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Madsen.
|-
|Fp
|align="center"|fp
|align="center"|B
|align="center"|N
|align="center"|U
|Up-weight rare variant using inverse frequency from controls by Danyu Lin.
|-
|Exact CMC
|align="center"|exactCMC
|align="center"|B
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Fisher's exact test.
|-
|CMC Wald
|align="center"|cmcWald
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Collapsing and combine rare variants, then pefore Wald test.
|-
|RareCover
|align="center"|rarecover
|align="center"|B
|align="center"|N
|align="center"|U
|Find optimal grouping unit for rare variant tests by Thomas Hoffman.
|-
|CMAT
|align="center"|cmat
|align="center"|B
|align="center"|N
|align="center"|U
|Test non-coding variants by Matt Z.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--burden cmc</code> will apply CMC test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Variable threshold models ==

{|
!Single variant
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Variable threshold model
|align="center"|vt
|align="center"|B, Q
|align="center"|N
|align="center"|U
|Every rare-variant frequency cutoffs are tests by Alkes Price.
|-
|Variable threshold CMC
|align="center"|cmc
|align="center"|B, Q
|align="center"|N
|align="center"|U
|This models is natiive so that it output CMC test statistics under all possible frequency cutoffs.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--vt price</code> will apply score test.

(##) In trait column, B and Q stand for binary, quantitiave trait.

== Kernel models ==

{|
!Kernel
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|SKAT
|align="center"|skat
|align="center"|B, Q
|align="center"|Y
|align="center"|U
|Sequencing kernel association test by Shawn Lee.
|-
|KBAC
|align="center"|kbac
|align="center"|B
|align="center"|N
|align="center"|U
|Kernel-based adaptive clustering model by Dajiang Liu.
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--kernel skat</code> will apply SKAT test.<br />To further customize SKAT test, you can use ''--kernel skat[nPerm=100:alpha=0.001:beta1=1:beta2=20]'' to specify permutation counts, type-1 error,<br />beta distribution parameters for upweighting rare variants. Rvtests will output a message showing:

<pre>[INFO] SKAT test significance will be evaluated using 10000 permutations at alpha = 0.001 (beta1 = 1.00, beta2 = 20.00)</pre>
(##) In trait column, B and Q stand for binary, quantitiave trait.

== Meta-analysis models ==

{|
!Type
!align="center"|Model(#)
!align="center"|Traits(##)
!align="center"|Covariates
!align="center"|Related / unrelated
!Description
|-
|Score test
|align="center"|score
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|standard score tests
|-
|Dominant model
|align="center"|dominant
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under dominant disease model
|-
|Recessive model
|align="center"|recessive
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|score tests and covariance matrix under recessive disease model
|-
|Covariance
|align="center"|cov
|align="center"|Q
|align="center"|Y
|align="center"|R, U
|covariance matrix
|}

(#) Model columns list the regconized names in rvtests. For example, use <code>--meta score,cov</code> will generate score statistics and covariance matrix for meta-analysis.<br />(##) In trait column, B and Q stand for (b)inary, (q)uantitiave trait.

== Utility models ==

Rvtests has an usually option <code>--outputRaw</code>. When specify this, rvtests can output genotypes, phenotype, covariates(if any) and collapsed genotype to tabular files. These files can be imported into other software (e.g. R) for further analysis.

= Association test options =

== Sample inclusion/exclusion ==

Rvtests can flexibly specify which sample(s) to include or exclude:

<pre> --peopleIncludeID : give IDs of people that will be included in study
--peopleIncludeFile : from given file, set IDs of people that will be included in study
--peopleExcludeID : give IDs of people that will be included in study
--peopleExcludeFile : from given file, set IDs of people that will be included in study</pre>
<code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> are used to include/exclude samples from command line.<br />For example, specify <code>--peopleIncludeID A,B,C</code> will include A, B and C sample from the VCF files if they exists.<br /><code>--peopleIncludeID</code> and <code>--peopleExcludeID</code> followed by a file name will include or exclude the IDs in the file.<br />So to include sample A, B and C, you can provide a file, <code>people.txt</code>, looks like:

<pre>A
B
C</pre>
Then use <code>--peopleIncludeFile people.txt</code> to include them in the analysis.

== Variant site filters ==

It is common that different frequency cutoffs are applied in rare-variant analysis.<br />Therefore, rvtests specify frequency cutoffs.

Frequency Cutoff

<pre> --freqUpper : Specify upper minor allele frequency bound to be included in analysis
--freqLower : Specify lower minor allele frequency bound to be included in analysis</pre>
If you specify <code>--freqLower 0.01 --freqUpper 0.05</code>, only the variants with minor allele ferquncy between 0.01 and 0.05 (boundary inclusive) will be analyzed.

Similar to sample inclusion/exclusion options, you can specify a range of variants to be included by<br />specifying <code>--rangeList</code> option. For example <code>--rangeList 1:100-200</code> will include the chromosome 1 position 100bp to 200bp region.<br />Alternatively, use a separate file, <code>range.txt</code>, and <code>--rangeFile range.txt</code> to speicify association tests range.

<pre> --rangeList : Specify some ranges to use, please use chr:begin-end format.
--rangeFile : Specify the file containing ranges, please use chr:begin-end format.
--siteFile : Specify the file containing sites to include, please use "chr pos" format.</pre>
It is supported to filter variant site by site depth, minor allele count or annotation (annotated VCF file is needed).

<pre> --siteDepthMin : Specify minimum depth(inclusive) to be incluced in analysis
--siteDepthMax : Specify maximum depth(inclusive) to be incluced in analysis
--siteMACMin : Specify minimum Minor Allele Count(inclusive) to be incluced in analysis
--annoType : Specify annotation type that is follwed by ANNO= in the VCF INFO field, regular expression is allowed</pre>
''NOTE'': <code>--annoType Nonsynonymous</code> will only analyze nonsynonymous variants where they have <code>ANNO=Nonsynonymous</code> in the INFO field.<br />VCF with annotatino information are called annotated VCF here. And to annotate<br />a VCF file, you can use [https://github.com/zhanxw/anno ANNO], a fast and accurate annotation software.

== Genotype filters ==

Genotype with low depth or low quality can be filtered out by these options:

<pre> --indvDepthMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis
--indvDepthMax : Specify maximum depth(inclusive) of a sample to be incluced in analysis
--indvQualMin : Specify minimum depth(inclusive) of a sample to be incluced in analysis</pre>
When genotypes are filtered, they are marked as missing genotypes.<br />Consequently, samples with missing genotype may or may not be included in the analysis.<br />That means samples with genotypes may be dropped (<code>--impute drop</code>)<br />or may still be included (<code>--impute mean</code> or <code>--impute hwe</code>).<br />By default, genotypes are imputed to its means.<br />See next section about how you like to handle missing genotypes.

== Handle missing genotypes and phenotypes ==

When genotypes are missing (e.g. genotype = "./.") or gentoypes are filtered out,<br />there are three options to handle them: (1) impute to its mean(default option); (2) impute by HWE equilibrium; (3) remove from the model.<br />Use <code>--impute [mean|hwe|drop]</code> to specify which option to use.

When quantitative phenotypes are missing, for example, some samples have gneotype files, but not phenotypes,<br />rvtests can impute missing phenotype to its mean.

''NOTE:'' Do not use <code>--imputePheno</code> for binary trait.

In summary, the following two options can be used:

<pre> --impute : Specify either of mean, hwe, and drop
--imputePheno : Impute phenotype to mean by those have genotypes but no
phenotpyes

</pre>
== Specify groups (e.g burden unit) ==

Rare variants association tests are usually performed in gruops of variants.<br />The natural grouping unit is gene. Rvtests can read gene definition file in <code>refFlat</code> format,<br />and perform association for each gene. Use <code>--geneFile</code> option to specify the gene file name.<br />For example, <code>--geneFile refFlat_hg19.txt.gz</code> will use <code>refFlat_hg19.txt.gz</code> as gene definition file,<br />and then perform association tests for every gene. Use <code>--gene</code> to specify a subset of genes to test.<br />For example, <code>--gene CFH</code> will only test CFH gene.

Alternative grouping unit can be specified as ''set''.<br />These ''sets'' are treated similar to gene.<br />You can thus use <code>--setFile</code> to define sets (similar to <code>--geneFile</code> option),<br />and use <code>--set</code> to define a specific set (similar to <code>--gene</code> option).<br />Additionally, use <code>--setList</code> can speicify a set to test from command line.

The format of a set file is: (1) set names; (2) ranges (e.g. chrom:begin-end);<br />For example, you have a set file, <code>example.set</code>, like this:

<pre>set1 1:100-200,1:250-300
set2 2:500-600</pre>
You can specify <code>--setFile example.set --set set2</code> to group variants<br />within chromosome 2, position 500 to 600bp.<br />If you want to test a particular region, for example, chromosome 2, position 500 to 550bp,<br />but do not want to make another file, you can use <code>--setList 2:500-600</code>.

In summary, options related to ''Grouping Unit'' are listed below:

<pre> --geneFile : specify a gene file (for burden tests)
--gene : specify which genes to test
--setList : specify a list to test (for burden tests)
--setFile : specify a list file (for burden tests, first two columns:
setName chr:beg-end)
--set : specify which set to test (1st column)</pre>
= Sex chromosome analysis =

Rvtests suppport X chromosome analysis. In human X chromosome, there is PAR (pseudoautosomal region) and non-PAR region.<br />For males, there are two X allele in PAR region and one allele in non-PAR region.<br />While the PAR region is treated in the same way as autosomes, rvtests treate non-PAR region differently.<br />Below we will describe the details about how rvtests handles non-PAR region.

''Prepare data''. According to VCF standard, male genotype needs to coded as 0 or 1. For compatibility, rvtests also support 0/0 or 1/1 coding.<br />In VCF files, male genotypes can be written as "0", "1", "0|0", "0/0", "1|1", "1/1". All other genotypes will be treated as missing.

''Genotype in the regression model''. For consistencmaine, male genotypes are converted to 0 or 2.

''MetaScore results''. If specify <code>--meta score</code>, the output file <code>prefix.MetaScore.assoc</code> includes PAR-region and non-PAR region analysis.<br />But in the non-PAR region, the difference is that Hardy-Weinberg P-value are calculated using female samples.

''Related individuals''. Just append <code>--xHemi</code> to the <code>vcf2kinship</code> (more details in [[#kinship-generation|Kinship generation]]) and <code>rvtest</code> command lines. Rvtests can recognize non-PAR region kinship and use it in the analysis.

''PAR region''. PAR region is defined as two regions X:60001-2699520 and X:154931044-155270560. Use <code>--xLabel</code> can speicify which chromosome has PAR region (default: 23|X)<br />and use <code>--xParRegion</code> to specify PAR region (default: hg19, meaning '60001-2699520,154931044-155260560' in the UCSC build hg19, specify "hg18" will use PAR region definition in the UCSC build hg18).

= Kinship generation =

Analysis of related individual usually requires estimation of kinship. You can a separate tool, <code>vcf2kinship</code>.<br /><code>vcf2kinship</code> is usually included in rvtests binary distribution or can be built from software source codes.

<code>vcf2kinship</code> can calcualte pedigree kinship using a pedigree input file (PED format, see [[#phenotype-file|Phenotype file]], use option <code>--ped</code>).<br />The output file name is specified by <code>--prefix</code> option. If you use <code>--prefix output</code> then the output files will include <code>output.kinship</code>.

It can also calculate empirical kinship using genotyp input file (VCF format, see [[#genotype-file-vcf|Genotype file (VCF)]], use option <code>--inVcf</code>).<br />For empiricial kinship, you also need to specify the kinship model, either Balding-Nicols model (ue option <code>--bn</code>) or Identity-by-state model (use option <code>--ibs</code>).

In sex chromosome analysis, it is often required to generate kinship on X chromsoome regions, then you need to speicfy <code>--xHemi</code>. If your input VCF file has different X chromosome label (e.g. chromosome name is '23' instead of 'X'), you can use <code>--xLabel 23</code>.

If principal component decomposition (PCA) results are needed, you can use option <code>--pca</code>. Then output files with suffix '.pca' include PCA results.

When dealing with large input files, it is often preferred to use multiple CPU to speed up calculation using the option <code>--thread N</code> in which N is the number of CPUs.

For example, to generate pedigree-based kinship (<code>--ped</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region, the command line is:

<pre>vcf2kinship --ped input.ped --xHemi --out output</pre>
To generate empirical kinship (<code>--inVcf</code>) on both autosomal region and X chromosome (<code>--xHemi</code>) region using Balding-Nicols model, the command line is:

<pre>vcf2kinship --inVcf input.vcf.gz --ped input.ped --bn --xHemi --out output</pre>
NOTE: you need to provide a pedigree file (PED) in the above case, as <code>vcf2kinship</code> need the sex information of samples.

= Frequently Asked Questions (FAQ) =

* Does rvtests suppport binary traits of related-individuals?

Not yet. It's a complex scenario and we have not found good solutions.

* Can you provide a list of command line options?

Rvtests have build help taht can be found by executing <code>rvtest --help</code>.<br />We also put all available options in this [https://github.com/zhanxw/rvtests/wiki/Command-Line-Options link].

* Can you provide standard error (SE) or confidence interval (CI) for the estimated Beta in the score model?

In the output of MetaScore model (--meta score), the standard error is the inverse of SQRT_V_STAT.<br />For example, if SQRT_V_STAT = 2, that means the standard error of estimated beta is 1/2 = 0.5.

* Why the INFORMATIVE_ALT_AC, N_REF and N_ALT columns have zero counts for certain chromosome X regions in meta-analysis models?

These counts are calculated from female individuals. If your study only has male samples, rvtests cannot report these counts. Because if a male carries a non-reference allele, we cannot conclude that this is heterozygous (0/1) site or homozygous alternatives (1/1) site.

= Feedback/Contact =

Questions and requests can be sent to Xiaowei Zhan<br />([mailto:zhanxw@umich.edu [mailto:zhanxw@umich.edu zhanxw@umich.edu]])<br />or Goncalo Abecasis<br />([mailto:goncalo@umich.edu [mailto:goncalo@umich.edu goncalo@umich.edu]])

Rvtests is a collaborative effort by Youna Hu, Bingshan Li, Dajiang Liu.

[https://bitdeli.com/free [[Image:https://d2weczhvl823v0.cloudfront.net/zhanxw/rvtests/trend.png|Bitdeli Badge]]]

LiftOver

2014-10-14T05:40:36Z

Zhanxw: /* Resources */

LiftOver is a necesary step to bring all genetical analysis to the same reference build.
LiftOver can have three use cases:

(1) [[#Lift genome positions | Convert genome position from one genome assembly to another genome assembly]]

In most scenarios, we have known genome positions in NCBI build 36 (UCSC hg 18) and hope to lift them over to NCBI build 37 (UCSC hg19).

(2) [[#Lift dbSNP rs numbers | Convert dbSNP rs number from one build to another]]

(3) [[#Lift Merlin/PLINK format | Convert both genome position and dbSNP rs number over different versions]]

It is likely to see such type of data in Merlin/PLINK format.

We will explain the work flow for the above three cases. In the rest of this article,
our example is to lift over from lower/older build to newer/higher build, as it is the common practice.

Using different tools, liftOver can be easy.
For example, UCSC liftOver tool is able to lift [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format] file between builds.
With our customized scripts, we can also lift rsNumber and Merlin/PLINK data files.

== Lift genome positions ==
Genome positions are best represented in [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format]. UCSC provides tools to convert BED file from one genome assembly to another.

=== Binary liftOver tool ===
We need [http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver liftOver binary] from UCSC and [http://hgdownload.cse.ucsc.edu/goldenPath/hg18/liftOver/hg18ToHg19.over.chain.gz hg18 to hg 19 chain file].

Provide BED format file (e.g. input.bed)

NOTE: Use the 'chr' before each chromosome name
<pre>
chr1 743267 743268 rs3115860
chr1 766408 766409 rs12124819
chr1 773885 773886 rs17160939
</pre>
Run liftOver:

liftOver input.bed hg18ToHg19.over.chain.gz output.bed unlifted.bed

unlifted.bed file will contain all genome positions that cannot be lifted. The reason for that varies. See [[#Various reasons that lift over could fail | Various reasons that lift over could fail]]

=== Web interface ===
Alternatively, you can lift over BED file in web interface
at: [http://genome.ucsc.edu/cgi-bin/hgLiftOver Link]
Web interface can tell you why some genome position cannot
be lifted if you click "Explain failure messages"

== Lift dbSNP rs numbers ==
rs number is release by dbSNP. UCSC also make their own copy from each dbSNP version. Be aware that the same version of dbSNP from these two centers are not the same.
When we convert rs number from lower version to higher version, there are practically two ways.

=== Use RsMergeArch and SNPHistory ===
It is necessary to quickly summarize how dbSNP merge/re-activate rs number:

# when different rs number are found to refer to the same SNP, then higher rs number will be merged to lower rs number, and the merging will be recorded in RsMergeArch.bcp.gz.
# when rs number have to be retracted, rs number will be recorded in SNPHistory.bcp.gz
# a retracted SNP can be [http://www.ncbi.nlm.nih.gov/books/NBK44496/#Schema.rs4823903_which_has_merged_into re-activated] in SNPHistory.bcp.gz by adding comment

With the above in mind, we are able to combine these two tables to obtain the relationship between older rs number and new rs number.
We have developed a script (for internal use), named [http://genome.sph.umich.edu/wiki/LiftRsNumber.py liftRsNumber.py] for lift rs numbers between builds.
This scripts require RsMergeArch.bcp.gz and SNPHistory.bcp.gz, those can be found in [[#Resources | Resources]].

Example input:
<pre>
3000
3001
3002
</pre>

Command:
<pre>
python liftRsNumber.py input.rs > output.rs
</pre>

Example output:
<pre>
unchanged 3000
lifted 2032
unchanged 3002
</pre>

== Lift Merlin/PLINK format ==
In Merlin/PLINK .map files, each line contains both genome position and dbSNP rs number. Our goal here is to use both information to liftOver as many position as possible.
There are 3 methods to liftOver and we recommend the first 2 method. The first method is common and applicable in most cases, and in our observations it lifts the most genome positions, however, it does not reflect the rs number change between different dbSNP builds. The second method is more robust in the sense that each lifted rs number has valid genome position, as it lift over old rs number as the first step by using dbSNP data. The third method is not straigtforward, and we just briefly mention it.

=== Lift Merlin format ===

PLINK format and [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin format are nearly identical].
The difference is that Merlin .map file have 4 columns. We will show
the lift over procedure for PLINK format, then you can use:

awk '{print $1,$2,"\t",$3;}' PLINK.map > Merlin.map

to obtain Merlin .map file.

=== Lift PLINK format ===
[http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml PLINK] format usually referrs to .ped and .map files.

==== Method 1 ====
We mainly use UCSC LiftOver binary tools to help lift over.
We have a script [[#Resources | liftMap.py]], however, it is recommended to understand the job step by step:

(1) Convert .map to .bed file

By rearrange columns of .map file, we obtain a standard BED format file.

(2) LiftOver .bed file

Use method mentioned [[#Lift genome positions | above]] to convert .bed file from one build to another.

(3) Convert lifted .bed file back to .map file

Rearrange column of .map file to obtain .bed file in the new build.

(4) Modify .ped file

.ped file have many column files. By convention, the first six columns are family_id, person_id, father_id, mother_id, sex, and phenotype.
From the 7th column, there are two letters/digits representing a genotype at the certain marker. In step (2), as some genome positions cannot
be lifted to the new version, we need to drop their corresponding columns from .ped file to keep consistency. You can use PLINK --exclude those snps,
see [http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml#exclude Remove a subset of SNPs].

(5) (optionally) change the rs number in the .map file

Similar to the human reference build, dbSNP also have different versions. You may consider change rs number from the old dbSNP version to new dbSNP version
depending on your needs. Such steps are described in [[#Lift dbSNP rs numbers | Lift dbSNP rs numbers]].

==== Method 2 ====
The idea is to use [[#Resources |LiftRsNumber.py]] to convert old rs number to new rs number, use the data file b132_SNPChrPosOnRef_37_1.bcp.gz (a data file containing each dbSNP and its positions in NCBI build 37), and adjust .map and .ped files accordingly.

(1) Extract and lift rs numbers

Use the tools [[#Use RsMergeArch and SNPHistory | LiftRsNumber.py]] to lift the rs number in the map file from old build to new build.

(2) Lookup SNP positions from rs number

dbSNP provides a file [[#Resources | b132_SNPChrPosOnRef_37_1.bcp.gz]] which contains rsNumber, chromosome and its position.
Use this file along with the new rsNumber obtained in the first step.
In practice, some rs numbers do not exist in build 132, or not suitable to be considered ( e.g. they do not reside on human reference, or they are mapped to multiple locations, these scenarios are noted by the chromosome column with values like "AltOnly", "Multi", "NotOn", "PAR", "Un"), we can drop them in the liftover procedure.
We will obtain the rs number and its position in the new build after this step.

(3) Lift .map file and .ped file

To lift over .map files, we can scan its content line by line, and skip those not lifted rs number.
Accordingly, it is necessary to drop the un-lifted SNP genotypes from .ped file.

==== Method 3 ====
NCBI dbSNP team has provided a [[ #Resources | provisional map ]] for converting the genome position of a larget set dbSNP from NCBI build 36 to NCBI build 37.
In the second step, we have obtained unlifted genome positions, so we can try to use the table to convert those unlfted dbSNPs.
After this step, there are still some SNPs that cannot be lifted, as they are mostly located on non-reference chromosome.
Note: due to the limitation of the provisional map, some SNP can have multiple locations.
For example, we cannot convert rs10000199 to chromosome 4, 7, 12.
<pre>
10000199 A/G 4 166142415 166142415 2 3 G + 4 165922965 165922965 2 3 G +
10000199 A/G 7 4589694 4589694 2 3 C - 7 4623168 4623168 2 3 C -
10000199 A/G 12 57008620 57008620 2 3 C - 12 58722353 58722353 2 3 C -
10000199 A/G 5 156018406 156018406 2 3 C - 5 156085828 156085828 2 3 C -
</pre>

We can dissect this method into steps:

(1) Remove invalid record in dbSNP provisional map.

Provisional map have duplicated rs number or the chromsome in the new build can be "Unable to map"(UN), we need to clean this table.

(2) Use provisional map to update .map file

By joining .map file and this provisional map, we can obtain the new genome position in the new build.
Note: provisional map uses 1-based chromosomal index. Things will get tricker if we want to lift non-single site SNP e.g. AA/GG
Since provisional map provides a range in this case, it is necessary to know the genome position of that single base provided in the .map file,
and then we can look up the table, so it is not straigtforward.

(3) Adjust .map and .ped file

For those lifted dbSNP, we need to keep them in the .map files, otherwise, we need to delete them.
Accordingly, we need to deleted SNP genotypes for those cannot be lifted.

== Various reasons that lift over can fail ==

=== Genome position cannot be lifted ===
When a SNP resides in a contig that only exists in older reference build, liftOver cannot give it new genome.

You can try the following SNP (in BED format) in UCSC online liftOver site:
20 56737667 56737668 rs1073519
The error message will be: "Sequence intersects no chains"

=== SNP in higher build are located in non-referernce assembly ===
Some SNP are not in autosomes or sex chromosomes in NCBI build 37. dbSNP does not include them.
You cannot use dbSNP database to lookup its genome position by rs number.

Take rs1006094 as an example:
In NCBI dbSNP webpage, this SNP is reported as "Mapped unambiguously on non-reference assembly only"
Thus it is probably not very useful to lift this SNP.

=== rs number changed in newer dbSNP build ===
It is possible that new dbSNP build does not have certain rs numbers.
When dbSNp release new build, higher rs number may be merged to lower rs number because of those rs numbers are actually the same SNP.
This merge process can be complicate. For short description, see [[ #Use RsMergeArch and SNPHistory | Use RsMergeArch and SNPHistory ]].
For detail, see:

[http://www.ncbi.nlm.nih.gov/books/NBK44395/#FTP.do_you_have_a_table_of_merged_snps_s Finding Specific Data in dbSNP’s FTP Files]

[http://www.ncbi.nlm.nih.gov/books/NBK44468/#Build.can_two_id_numbers_correspond_to_t Merging RefSNP Numbers and RefSNP Clusters]

For example:

rs3001 has merged to rs2032.

=== Different dbSNP build ===
NCBI released dbSNP132 (VCF format), and UCSC also have their version of dbSNP132 (plain txt).
The two database files differ not only in file format, but in content.

For NCBI release, its [[#Resources| release]] will not contain:
* SNPs listed as microsatellites or named variations
* SNPs with multibyte alleles and unknown (N) adjacent base pairs
* SNPs that are not mapped on the reference genome (GRCh37)

For UCSC release, see [[ #Resources | UCSC dbSNP track note ]]

Use rs1054140 as an example:

NCBI dbSNP website gives 1 location:
[http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1054140 Link]

NCBI dbSNP VCF file has NO record.

UCSC genome browser website gives 2 locations:
[http://genome.ucsc.edu/cgi-bin/hgTracks?clade=mammal&org=Human&db=hg19&position=rs1054140&hgt.suggest=&hgt.suggestTrack=knownGene&pix=800&Submit=submit&hgsid=205770459&hgt.newJQuery=1 Link]

UCSC dbSNP file give 2 locations:
<pre>
721 chr10 17842693 17842694 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
723 chr10 18089681 18089682 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
</pre>

== Resources ==
* liftRsNumber.py [[liftRsNumber.py]] and its interal location: /net/dumbo/net/dumbo/home/zhanxw/amd/analyze/verifyBamID/liftRsNumber.py
* liftMap.py [[liftMap.py]]
* NCBI provisional map [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.txt.gz file] and [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.info info]
* NCBI RgMergeArch [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/RsMergeArch.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=RsMergeArch schema]
* NCBI SNPHistory [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/SNPHistory.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPHistory schema]
* NCBI SNPChrPosOnRef build 132 [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/b132_archive/organism_data/b132_SNPChrPosOnRef_37_1.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPChrPosOnRef schema]. If this link becomes unavailable, please consider using the updated file ([ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/database/organism_data/b141_SNPChrPosOnRef.bcp.gz link]).

* How UCSC dbSNP differs from NCBI dbSNP [http://genomewiki.ucsc.edu/index.php/DbSNP_Track_Notes UCSC dbSNP track note]
* The dbSNP mapping process [http://www.ncbi.nlm.nih.gov/books/NBK44455/ link]
* NCBI dbSNP release 132 [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/VCF/v4.0/ByChromosomeNoGeno/00-All.vcf.gz 00-All.vcf.gz]
* UCSC dbSNP release 132 [http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz snp132.txt.gz]

== Third party contribution ==

* liftOver for BEDPE format by [mailto://doug.phanstiel@gmail.com Doug] - [http://www.sharedproteomics.com/forum/showthread.php?2930-liftOverBedpe.py-Convert-genome-coordinates-for-bedpe-files&p=4854#post4854 liftOverBedpe.py].

== Acknowledge ==

* Hyun: provides sample liftOver tool: [/net/wonderland/home/hmkang/prj/Sardinia/MetaboChip/scripts/j01-liftover-metabochip-positions.pl]
* Alex: careful examines of 0-based index in UCSC data file
* Adrian: explaination of SNPs omitted in NCBI dbSNP file
* Goncalo: all other supports

== Questions and Comments ==

Please contact [mailto:zhanxw@umich.edu Xiaowei Zhan].

Ancestry

2014-08-12T18:53:31Z

Zhanxw: Created page with "= Introduction = Ancestry can infer ancestry using sequence reads in the Principal Component (PC) space. It is suited for targeted/exome sequencing or whole genome sequenci..."

= Introduction =

Ancestry can infer ancestry using sequence reads in the Principal Component (PC) space.
It is suited for targeted/exome sequencing or whole genome sequencing experiments.
Ancestry is implemented in C++ for fast computations.

= Download =

To get a copy of the software, please contact: zhanxw@umich.edu .

For CSG users, the binary executable is located at: /net/fantasia/home/zhanxw/spa/cpp/executable/ancestrySeq

We plan to open source this program shortly.

= Command Line Options =

== Input sequence data(--inSeq) ==

Sequence data (BAM files) need to be preprocessed in .seq format. This procedure is described [http://genome.sph.umich.edu/wiki/LASER#Seq_file here].

Seq file is generated from pileup files. It contains sequencing information and organize it in a LASER readable format.
The first two columns represent population id and individual id.
Subsequent columns are total read depths and reference base counts.
For example, column 3 and 4 are 0, 0 in the following example. That means at first marker, the sequence read depth is 0 and thus none of the reads has reference base.
We enforce tab delimiters between markers and space delimiters between each read depths and reference base counts.
On line of seq file looks like below:

NA12878.chrom22 NA12878.chrom22 0 0 0 0 0 0 0 0 0

== Input pileup sites (--inSite) ==

Site file is equivalent to BED file and it is used here to represent marker positions.

The preprocessing procedure is described [http://genome.sph.umich.edu/wiki/LASER#Site_file here].

An example site file looks like below:

CHR POS ID REF ALT
1 752566 rs3094315 G A
1 768448 rs12562034 G A
1 1005806 rs3934834 C T
1 1018704 rs9442372 A G
1 1021415 rs3737728 A G

The site file has header line, and it contains chromosome, position(1-based), id (usually marker name), ref (reference allele) and alt (alternative allele).

== Input parameter (--inModel) ==

This parameter specifies the SNP gradients and offsets. It is the output of program [http://genetics.cs.ucla.edu/spa/ spa]

== Output prefix (--out) ==

The parameter specifies the output prefix.
The main results will be stored in ''PREFIX''.loc file.

An example output file looks like below:

PopId IndvId Loc1 Loc2
MPaS3287 MPaS3287 13.5669 176.051
MPaS3287 MPaS3287.ConfInt95 6.43922,20.724 169.103,182.928

Note, when ''--ci'' option is used, the outputted IndvID column will append ".ConfInt95" indicates that the inference results is a 95% confidence interval.

== Inference option (--bootstrap) ==

When this option is specified, the program will infer ancestral locations using a bootstrap procedure.
Essentially, we resample input sequence reads, and recalculate ancestral locations after each shuffle.
The output will include ancestral locations from every resampling.

== Inference option (--ci) ==

When this option is specified, the program will infer ancestral locations based on likelihood calculations.
On a two-dimensional space, this option infer the top, bottom, left and right boundaries of an ellipse region,
and the probability that the true ancestral location fells in this region is 95%.

= Example =

A basic command looks like:

ancestrySeq --inSeq 1108.amd.to.hgdp.seq --inSite HGDP_938.site --inModel spa.model.out --out test

The result file, test.loc, includes inferred location for each sample listed in 1108.amd.to.hgdp.seq.

= Resources =

[http://genome.sph.umich.edu/wiki/LASER LASER] is a related project that starts earlier than ancestry.
Both software can perform ancestral inference.
But ancestral has computational advantages and does not sacrifice accuracies.

= Contact =

Comments on this wiki page or questions related to preparing input files for LASER can be sent to [mailto:zhanxw@umich.edu Xiaowei Zhan].
This project was helped by Chaolong Wang and was directed by Gonçalo Abecasis.

LASER

2014-07-10T14:41:35Z

Zhanxw: /* Introduction */

= Introduction =

LASER, which stands for Locating Ancestry using SEquencing Reads, is a C++ software package that can estimate individual ancestry directly from genome-wide shortgun sequencing reads without calling genotypes. The method relies on the availability of a set of reference individuals whose genome-wide SNP genotypes and ancestral information are known. We first construct a reference coordinate system by applying principal components analysis (PCA) to the genotype data of the reference individuals. Then, for each sequencing sample, use the genome-wide sequencing reads to place the sample into the reference PCA space. With an appropriate reference panel, the estimated coordinates of the sequencing samples identify their ancestral background and can be directly used to correct for population structure in association studies or to ensure adequate matching of cases and controls.

The goal of this wiki page is to help you get start using LASER.
You can also find workshop lecture notes from [http://genome.sph.umich.edu/wiki/SeqShop:_Estimates_of_Genetic_Ancestry_Practical 2014 UM Sequencing Workshop].
We encourage you to read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual] for more details.

= Download =

To get a copy of the software and manual, go to the [http://www.sph.umich.edu/csg/chaolong/LASER/ LASER Download] page.

= Workflow =

LASER generates the coordinates from both reference individuals and sequence samples. It requires essentially two input files:

[[File:LASER-Workflow.png|thumb|center|alt=LASER workflow|400px|LASER Workflow]]

*Seq file: a text file processed from BAM (alignment) files. (See [[#Process sequencing file (BAM)|Processing sequencing file]] for how to prepare seq file)
*Geno file: genotypes of reference individuals. (See [[#Geno file|Geno file]] to understand geno file format)

LASER typically outputs two coord files: (1) in reference individuals' coord file(Reference.coord), LASER outputs the reference coordinates in the PCA space; (2) in sequence samples' coord files(AllSamples.coord), LASER infers their ancestries by placing their ancestry coordinates onto reference samples' PCA space.

An example result of the coord file of sequence samples is shown below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

In the header line, popID means "population ID", indivID means "individual ID", L1 means number of loci that has been covered by at least one read, Ci means "average coverage", t means Procrustes similarity. PC1 and PC2 mean coordinates of the first and second principal components.

= Tutorial =

In this tutorial, we will show you how to prepare data and run LASER.

== Process sequencing file (BAM) ==

We illustrate how to obtain .seq file from BAM files in this section.
In this example, we use HGDP data set as a reference, which contains 938 individuals and 632,958 markers.
[[File:LASER-DataProcessing.png|thumb|center|alt=LASER workflow|400px|LASER Data Processing Procedure]]

1. Obtain pileup files from BAM files

The first step is to generate a BED file:

cat ../resource/HGDP/HGDP_938.site |awk '{if (NR > 1) {print $1, $2-1, $2;}}' > HGDP_938.bed

This BED file contains the positions of all the reference markers.

Then we use ''samtools'' to extract the sequence bases overlapping these 632,958 reference markers.
Assuming your BAM file name is ''NA12878.chrom22.recal.bam'' (our example BAM file), you can use this:

samtools mpileup -q 30 -Q 20 -f ../../LASER-resource/reference/hs37d5.fa -l HGDP_938.bed exampleBAM/NA12878.chrom22.recal.bam > NA12878.chrom22.pileup

to obtain a pileup file named ''NA12878.chrom22.pileup''. It is required to keep the ''.pileup' suffix.

2. Obtain a seq file from pileup files.

After obtaining pileup files from each BAM file, you can convert them into a single seq file before running LASER.
Use the same site file and all generated pileup files from step 1 to generate a seq file:

python pileup2seq.py -m ../resource/HGDP/HGDP_938.site -o test NA12878.chrom22.pileup

You should obtain test.seq file after this step.

== Estimate ancestries of sequence samples ==

The easiest way to perform LASER using its exemplar data is:

./laser -s pileup2seq/test.seq -g resource/HGDP/HGDP_938.geno -c resource/HGDP/HGDP_938.RefPC.coord -o test -k 2

Upon successful calculation, you will find a result file "test.SeqPC.coord".

<br>

== Interpret LASER outputs ==

Upon successfully launching LASER command line as above, the output messages should be similar to below:

===================================================================
==== LASER: Locating Ancestry from SEquencing Reads ====
==== Version 1.0 | (c) Chaolong Wang 2013 ====
====================================================================
Started at: Fri Nov 15 01:05:48 2013

938 individuals are detected in the GENO_FILE.
632958 loci are detected in the GENO_FILE.
1 individuals are detected in the SEQ_FILE.
632958 loci are detected in the SEQ_FILE.
938 individuals are detected in the COORD_FILE.
100 PCs are detected in the COORD_FILE.

Parameter values used in execution:
-------------------------------------------------
GENO_FILE (-g)resource/HGDP/HGDP_938.geno
SEQ_FILE (-s)pileup2seq/test.seq
COORD_FILE (-c)resource/HGDP/HGDP_938.RefPC.coord
OUT_PREFIX (-o)test
DIM (-k)2
MIN_LOCI (-l)100
SEQ_ERR (-e)0.01
FIRST_IND (-x)1
LAST_IND (-y)1
REPS (-r)1
OUTPUT_REPS (-R)0
CHECK_FORMAT (-fmt)10
CHECK_COVERAGE (-cov)0
PCA_MODE (-pca)0
-------------------------------------------------

Fri Nov 15 01:05:50 2013
Checking data format ...
GENO_FILE: OK.
SEQ_FILE: OK.
COORD_FILE: OK.

Fri Nov 15 01:06:01 2013
Reading reference genotypes ...

Fri Nov 15 01:09:15 2013
Reading reference PCA coordinates ...

Fri Nov 15 01:09:15 2013
Analyzing sequence samples ...
Results for the sequence samples are output to 'test.SeqPC.coord'.

Finished at: Fri Nov 15 01:09:21 2013
====================================================================

The ancestry of input samples are store in the file '''test.SeqPC.coord''', which content is shown below:

popID indivID L1 Ci t PC1 PC2
NA12878.chrom22 NA12878.chrom22 1601 0.00858193 0.977243 31.522 224.098

The ancestry coordinates for NA12878 samples are given in PC1 (31.522) and PC2 (224.098).

It is recommended to visualize this results with HGDP reference samples whose coordinates are given in file: resource/HGDP/HGDP_938.RefPC.coord

In our manuscript, an example figure is shown:

[[File:LASER paper Figure 2.png|thumb|center|alt=LASER example outputs as in Figure 2|400px|LASER Outputs]]

In this figure, 238 individuals were randomly selected from the total 938 HGDP samples as the testing set (colored symbols),
and the remaining 700 HGDP individuals were used as the reference panel (gray symbols).

= File format =

== Geno file ==

Geno file are from reference samples. LASER use genotype of these samples as a reference panel. You can obtain geno file from VCF files using [https://github.com/zhanxw/vcf2geno vcf2geno].

In our resource folder, we provide an example geno file for the HGDP data set (resource/HGDP/HGDP_938.geno):

Brahui HGDP00001 1 2 1 1 0 2 0 2 1 2 2 2 1 1 2 1 0
Brahui HGDP00003 0 0 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00005 0 2 2 0 0 1 0 2 1 2 2 2 2 1 2 2 1
Brahui HGDP00007 0 2 2 0 0 2 0 2 0 2 2 2 1 1 2 2 1
Brahui HGDP00009 0 1 0 1 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00011 1 1 2 1 1 2 1 1 1 2 2 2 1 1 2 2 0
Brahui HGDP00013 1 2 2 1 1 2 1 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00015 1 1 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00017 1 1 2 0 0 1 0 0 0 2 0 1 1 2 2 2 0
Brahui HGDP00019 0 2 2 0 0 1 0 1 0 2 1 2 2 1 2 2 0

The first and second columns represent the population id and individual id.
From the third column, each number represents a genotype.
To be consistent with the sequence data, genotypes should be given on the '''forward strand'''. Genotypes are coded by 0, 1, or 2, representing copies of the
reference allele at a locus in one individual.

In this geno file, we have 632,960 columns which contains 632,958 markers from column 3 to the last column.

== Seq file ==
Seq file is generated from pileup files. It contains sequencing information and organize it in a LASER readable format.
The first two columns represent population id and individual id.
Subsequent columns are total read depths and reference base counts.
For example, column 3 and 4 are 0, 0 in the following example. That means at first marker, the sequence read depth is 0 and thus none of the reads has reference base.
We enforce tab delimiters between markers and space delimiters between each read depths and reference base counts.
On line of seq file looks like below:

NA12878.chrom22 NA12878.chrom22 0 0 0 0 0 0 0 0 0

== Pileup file ==

Pileup file are generated using samtools. An example pileup file is shown below:

22 17094749 A 1 c D
22 17202602 T 1 . D
22 17411899 A 1 . C
22 17450515 G 2 ., 9<
22 17452966 T 1 c 5
22 17470779 C 1 , A
22 17492203 G 1 , B
22 17504945 C 3 ,.. BCA
22 17529814 T 3 .., CCC

The columns are chromosome, position (1-based), reference base, depth, bases and base qualities.

== BED file ==
BED file represents genomic regions and it follows [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 UCSC conventions]:

1 752565 752566
1 768447 768448
1 1005805 1005806
1 1018703 1018704
1 1021414 1021415

The columns are: chromosome, start position (0-based) and end position (1-based).

== Coord file ==
Coord files represent the ancestries of both reference samples and sequence samples.
An example coord file looks like below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

The columns are: popID means "population ID", indivID means "individual ID", L1 means number of loci has been covered, Ci means "average coverage", t means Procrustes similarity.
PC1, PC2 means coordinates of first and second principal components. You may notice L1, Ci, and t are omitted in the coord files of reference samples. The reason is that reference samples use genotypes and do not have coverage information.

== Site file ==
Site file is equivalent to BED file and it is used here to represent marker positions. An example site file looks like below:
CHR POS ID REF ALT
1 752566 rs3094315 G A
1 768448 rs12562034 G A
1 1005806 rs3934834 C T
1 1018704 rs9442372 A G
1 1021415 rs3737728 A G

The site file has header line, and it contains chromosome, position(1-based), id (usually marker name), ref (reference allele) and alt (alternative allele).

= Advanced options =

LASER has advanced options including (1) parallel computing; (2) increase ancestry inference accuracy using repeated runs; (3) generate PCA coordiates using genotypes.
See [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf LASER Manual] for detailed information.

= Contact =
Comments on this wiki page or questions related to preparing input files for LASER can be sent to [mailto:zhanxw@umich.edu Xiaowei Zhan].
Comments on the LASER software or the user's manual can be sent to [mailto:chaolong@umich.edu Chaolong Wang].
This project was directed by Gonçalo Abecasis and Sebastian Zöllner at the University of Michigan.

LiftOver

2014-02-14T15:12:53Z

Zhanxw: /* Resources */

LiftOver is a necesary step to bring all genetical analysis to the same reference build.
LiftOver can have three use cases:

(1) [[#Lift genome positions | Convert genome position from one genome assembly to another genome assembly]]

In most scenarios, we have known genome positions in NCBI build 36 (UCSC hg 18) and hope to lift them over to NCBI build 37 (UCSC hg19).

(2) [[#Lift dbSNP rs numbers | Convert dbSNP rs number from one build to another]]

(3) [[#Lift Merlin/PLINK format | Convert both genome position and dbSNP rs number over different versions]]

It is likely to see such type of data in Merlin/PLINK format.

We will explain the work flow for the above three cases. In the rest of this article,
our example is to lift over from lower/older build to newer/higher build, as it is the common practice.

Using different tools, liftOver can be easy.
For example, UCSC liftOver tool is able to lift [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format] file between builds.
With our customized scripts, we can also lift rsNumber and Merlin/PLINK data files.

== Lift genome positions ==
Genome positions are best represented in [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 BED format]. UCSC provides tools to convert BED file from one genome assembly to another.

=== Binary liftOver tool ===
We need [http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver liftOver binary] from UCSC and [http://hgdownload.cse.ucsc.edu/goldenPath/hg18/liftOver/hg18ToHg19.over.chain.gz hg18 to hg 19 chain file].

Provide BED format file (e.g. input.bed)

NOTE: Use the 'chr' before each chromosome name
<pre>
chr1 743267 743268 rs3115860
chr1 766408 766409 rs12124819
chr1 773885 773886 rs17160939
</pre>
Run liftOver:

liftOver input.bed hg18ToHg19.over.chain.gz output.bed unlifted.bed

unlifted.bed file will contain all genome positions that cannot be lifted. The reason for that varies. See [[#Various reasons that lift over could fail | Various reasons that lift over could fail]]

=== Web interface ===
Alternatively, you can lift over BED file in web interface
at: [http://genome.ucsc.edu/cgi-bin/hgLiftOver Link]
Web interface can tell you why some genome position cannot
be lifted if you click "Explain failure messages"

== Lift dbSNP rs numbers ==
rs number is release by dbSNP. UCSC also make their own copy from each dbSNP version. Be aware that the same version of dbSNP from these two centers are not the same.
When we convert rs number from lower version to higher version, there are practically two ways.

=== Use RsMergeArch and SNPHistory ===
It is necessary to quickly summarize how dbSNP merge/re-activate rs number:

# when different rs number are found to refer to the same SNP, then higher rs number will be merged to lower rs number, and the merging will be recorded in RsMergeArch.bcp.gz.
# when rs number have to be retracted, rs number will be recorded in SNPHistory.bcp.gz
# a retracted SNP can be [http://www.ncbi.nlm.nih.gov/books/NBK44496/#Schema.rs4823903_which_has_merged_into re-activated] in SNPHistory.bcp.gz by adding comment

With the above in mind, we are able to combine these two tables to obtain the relationship between older rs number and new rs number.
We have developed a script (for internal use), named [http://genome.sph.umich.edu/wiki/LiftRsNumber.py liftRsNumber.py] for lift rs numbers between builds.
This scripts require RsMergeArch.bcp.gz and SNPHistory.bcp.gz, those can be found in [[#Resources | Resources]].

Example input:
<pre>
3000
3001
3002
</pre>

Command:
<pre>
python liftRsNumber.py input.rs > output.rs
</pre>

Example output:
<pre>
unchanged 3000
lifted 2032
unchanged 3002
</pre>

== Lift Merlin/PLINK format ==
In Merlin/PLINK .map files, each line contains both genome position and dbSNP rs number. Our goal here is to use both information to liftOver as many position as possible.
There are 3 methods to liftOver and we recommend the first 2 method. The first method is common and applicable in most cases, and in our observations it lifts the most genome positions, however, it does not reflect the rs number change between different dbSNP builds. The second method is more robust in the sense that each lifted rs number has valid genome position, as it lift over old rs number as the first step by using dbSNP data. The third method is not straigtforward, and we just briefly mention it.

=== Lift Merlin format ===

PLINK format and [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin format are nearly identical].
The difference is that Merlin .map file have 4 columns. We will show
the lift over procedure for PLINK format, then you can use:

awk '{print $1,$2,"\t",$3;}' PLINK.map > Merlin.map

to obtain Merlin .map file.

=== Lift PLINK format ===
[http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml PLINK] format usually referrs to .ped and .map files.

==== Method 1 ====
We mainly use UCSC LiftOver binary tools to help lift over.
We have a script [[#Resources | liftMap.py]], however, it is recommended to understand the job step by step:

(1) Convert .map to .bed file

By rearrange columns of .map file, we obtain a standard BED format file.

(2) LiftOver .bed file

Use method mentioned [[#Lift genome positions | above]] to convert .bed file from one build to another.

(3) Convert lifted .bed file back to .map file

Rearrange column of .map file to obtain .bed file in the new build.

(4) Modify .ped file

.ped file have many column files. By convention, the first six columns are family_id, person_id, father_id, mother_id, sex, and phenotype.
From the 7th column, there are two letters/digits representing a genotype at the certain marker. In step (2), as some genome positions cannot
be lifted to the new version, we need to drop their corresponding columns from .ped file to keep consistency. You can use PLINK --exclude those snps,
see [http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml#exclude Remove a subset of SNPs].

(5) (optionally) change the rs number in the .map file

Similar to the human reference build, dbSNP also have different versions. You may consider change rs number from the old dbSNP version to new dbSNP version
depending on your needs. Such steps are described in [[#Lift dbSNP rs numbers | Lift dbSNP rs numbers]].

==== Method 2 ====
The idea is to use [[#Resources |LiftRsNumber.py]] to convert old rs number to new rs number, use the data file b132_SNPChrPosOnRef_37_1.bcp.gz (a data file containing each dbSNP and its positions in NCBI build 37), and adjust .map and .ped files accordingly.

(1) Extract and lift rs numbers

Use the tools [[#Use RsMergeArch and SNPHistory | LiftRsNumber.py]] to lift the rs number in the map file from old build to new build.

(2) Lookup SNP positions from rs number

dbSNP provides a file [[#Resources | b132_SNPChrPosOnRef_37_1.bcp.gz]] which contains rsNumber, chromosome and its position.
Use this file along with the new rsNumber obtained in the first step.
In practice, some rs numbers do not exist in build 132, or not suitable to be considered ( e.g. they do not reside on human reference, or they are mapped to multiple locations, these scenarios are noted by the chromosome column with values like "AltOnly", "Multi", "NotOn", "PAR", "Un"), we can drop them in the liftover procedure.
We will obtain the rs number and its position in the new build after this step.

(3) Lift .map file and .ped file

To lift over .map files, we can scan its content line by line, and skip those not lifted rs number.
Accordingly, it is necessary to drop the un-lifted SNP genotypes from .ped file.

==== Method 3 ====
NCBI dbSNP team has provided a [[ #Resources | provisional map ]] for converting the genome position of a larget set dbSNP from NCBI build 36 to NCBI build 37.
In the second step, we have obtained unlifted genome positions, so we can try to use the table to convert those unlfted dbSNPs.
After this step, there are still some SNPs that cannot be lifted, as they are mostly located on non-reference chromosome.
Note: due to the limitation of the provisional map, some SNP can have multiple locations.
For example, we cannot convert rs10000199 to chromosome 4, 7, 12.
<pre>
10000199 A/G 4 166142415 166142415 2 3 G + 4 165922965 165922965 2 3 G +
10000199 A/G 7 4589694 4589694 2 3 C - 7 4623168 4623168 2 3 C -
10000199 A/G 12 57008620 57008620 2 3 C - 12 58722353 58722353 2 3 C -
10000199 A/G 5 156018406 156018406 2 3 C - 5 156085828 156085828 2 3 C -
</pre>

We can dissect this method into steps:

(1) Remove invalid record in dbSNP provisional map.

Provisional map have duplicated rs number or the chromsome in the new build can be "Unable to map"(UN), we need to clean this table.

(2) Use provisional map to update .map file

By joining .map file and this provisional map, we can obtain the new genome position in the new build.
Note: provisional map uses 1-based chromosomal index. Things will get tricker if we want to lift non-single site SNP e.g. AA/GG
Since provisional map provides a range in this case, it is necessary to know the genome position of that single base provided in the .map file,
and then we can look up the table, so it is not straigtforward.

(3) Adjust .map and .ped file

For those lifted dbSNP, we need to keep them in the .map files, otherwise, we need to delete them.
Accordingly, we need to deleted SNP genotypes for those cannot be lifted.

== Various reasons that lift over can fail ==

=== Genome position cannot be lifted ===
When a SNP resides in a contig that only exists in older reference build, liftOver cannot give it new genome.

You can try the following SNP (in BED format) in UCSC online liftOver site:
20 56737667 56737668 rs1073519
The error message will be: "Sequence intersects no chains"

=== SNP in higher build are located in non-referernce assembly ===
Some SNP are not in autosomes or sex chromosomes in NCBI build 37. dbSNP does not include them.
You cannot use dbSNP database to lookup its genome position by rs number.

Take rs1006094 as an example:
In NCBI dbSNP webpage, this SNP is reported as "Mapped unambiguously on non-reference assembly only"
Thus it is probably not very useful to lift this SNP.

=== rs number changed in newer dbSNP build ===
It is possible that new dbSNP build does not have certain rs numbers.
When dbSNp release new build, higher rs number may be merged to lower rs number because of those rs numbers are actually the same SNP.
This merge process can be complicate. For short description, see [[ #Use RsMergeArch and SNPHistory | Use RsMergeArch and SNPHistory ]].
For detail, see:

[http://www.ncbi.nlm.nih.gov/books/NBK44395/#FTP.do_you_have_a_table_of_merged_snps_s Finding Specific Data in dbSNP’s FTP Files]

[http://www.ncbi.nlm.nih.gov/books/NBK44468/#Build.can_two_id_numbers_correspond_to_t Merging RefSNP Numbers and RefSNP Clusters]

For example:

rs3001 has merged to rs2032.

=== Different dbSNP build ===
NCBI released dbSNP132 (VCF format), and UCSC also have their version of dbSNP132 (plain txt).
The two database files differ not only in file format, but in content.

For NCBI release, its [[#Resources| release]] will not contain:
* SNPs listed as microsatellites or named variations
* SNPs with multibyte alleles and unknown (N) adjacent base pairs
* SNPs that are not mapped on the reference genome (GRCh37)

For UCSC release, see [[ #Resources | UCSC dbSNP track note ]]

Use rs1054140 as an example:

NCBI dbSNP website gives 1 location:
[http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1054140 Link]

NCBI dbSNP VCF file has NO record.

UCSC genome browser website gives 2 locations:
[http://genome.ucsc.edu/cgi-bin/hgTracks?clade=mammal&org=Human&db=hg19&position=rs1054140&hgt.suggest=&hgt.suggestTrack=knownGene&pix=800&Submit=submit&hgsid=205770459&hgt.newJQuery=1 Link]

UCSC dbSNP file give 2 locations:
<pre>
721 chr10 17842693 17842694 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
723 chr10 18089681 18089682 rs1054140 0 + T T A/T genomic single by-cluster,by-submitter ...
</pre>

== Resources ==
* liftRsNumber.py [[liftRsNumber.py]] and its interal location: /net/dumbo/net/dumbo/home/zhanxw/amd/analyze/verifyBamID/liftRsNumber.py
* liftMap.py [[liftMap.py]]
* NCBI provisional map [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.txt.gz file] and [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/misc/exchange/Remap_36_3_37_1.info info]
* NCBI RgMergeArch [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/RsMergeArch.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=RsMergeArch schema]
* NCBI SNPHistory [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/organism_data/SNPHistory.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPHistory schema]
* NCBI SNPChrPosOnRef build 132 [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/database/b132_archive/organism_data/b132_SNPChrPosOnRef_37_1.bcp.gz file] and [http://www.ncbi.nlm.nih.gov/SNP/snp_db_table_description.cgi?t=SNPChrPosOnRef schema]
* How UCSC dbSNP differs from NCBI dbSNP [http://genomewiki.ucsc.edu/index.php/DbSNP_Track_Notes UCSC dbSNP track note]
* The dbSNP mapping process [http://www.ncbi.nlm.nih.gov/books/NBK44455/ link]
* NCBI dbSNP release 132 [ftp://ftp.ncbi.nih.gov:/snp/organisms/human_9606/VCF/v4.0/ByChromosomeNoGeno/00-All.vcf.gz 00-All.vcf.gz]
* UCSC dbSNP release 132 [http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz snp132.txt.gz]

== Third party contribution ==

* liftOver for BEDPE format by [mailto://doug.phanstiel@gmail.com Doug] - [http://www.sharedproteomics.com/forum/showthread.php?2930-liftOverBedpe.py-Convert-genome-coordinates-for-bedpe-files&p=4854#post4854 liftOverBedpe.py].

== Acknowledge ==

* Hyun: provides sample liftOver tool: [/net/wonderland/home/hmkang/prj/Sardinia/MetaboChip/scripts/j01-liftover-metabochip-positions.pl]
* Alex: careful examines of 0-based index in UCSC data file
* Adrian: explaination of SNPs omitted in NCBI dbSNP file
* Goncalo: all other supports

== Questions and Comments ==

Please contact [mailto:zhanxw@umich.edu Xiaowei Zhan].

RvTests

2014-01-29T22:21:00Z

Zhanxw:

[[Category:Software]]

= Overview =
(See [[rvtests]] for more rare variants association tests or preparing meta-analysis inputs.)

A few rare variants tests (Li-Leal's CMC and Madsen-Browning's weighted method) are implemented in the logisitc regression framework using C++. Please contact Youna Hu (youna@umich.edu) for comments, suggestions or questions.

The source code is located at [[File:RV3Tests.v1.tar]]

You can just download the tar file, extract it and go to the RV3Test.v1 folder and type make all to compile the code, the binary file will then be in the exectuables folder.

= Example =

See a detailed [[example]] here.

= Syntax =

This software uses command line interface as follows

RARE VARIANT ANALYSIS OPTIONS:
GENOTYPE : --genofile [pos.012],
--geneList [outGeneSorted.txt], --cutoff [0.010],
--collapseChoice [or]
PHENOTYPE : --phenofile [LDL.y.ID]
COVARIATES : --covConsider, --covfile [covFile.ID.2.txt]
PERMUTATION : --nPermute [10], --PermutationSeed [1]
GENE LEVEL TEST RESULT : --geneGlobalTestOut [globalPermuteSummary.txt],
--geneTestpvalueFile [geneTestPvalues.txt]

;GENOTYPE

;--genofile: A genotype 012 matrix (.012 is the file) This file can be prepared by using the prepare012s
source code [[File:vcfReader.v1.tar]]
Again, extract the tar file and then go into the directory to type make all to compile the code,
you then will find binary file in the executables folder.

Note: You should use Yanming's vcf annotation [http://genome.sph.umich.edu/wiki/VcfCodingSnps] on your vcf file first
to output a annotated vcf file. You SHOULD keep the log file from the annotation, which will be used to create the gene list.

Data File PREPARATION
Input files : --vcf [LDL.test.vcf], --log [], --IDfile []
Subsetting choices : --All
Output files : --outputPrefix [subsetGeno],
--outputGeneList [LDL.geneList.txt]
--vcf: Input vcf file
--log: This is the log file from Yanming's annotation output, we use this log to obtain the gene list
--IDfile specifies a file with one column of subject IDs to subject from the vcf file.
If it is not specified, then all subjects are included for the format conversion.
--All: specifies 1 to include all variants and 0 to include only nonsyn and stop annotated variants.
-- outputPrefix: Specify the prefix for the four output files which will be used in rvTests
*.012: A genotype matrix with subjects as rows and variant sites as columns.
*.012.pos: Chromosome and position numbers.
*.012.indv: Subject IDs.
*.012.frq: The frequency of the included variants.
--outputGeneList: Specify a file to store the gene list which will be used in rvTest.
The list file looks like this
1 OR4F5 69090 70008
1 SAMD11 860529 871276
1 NOC2L 879583 893918
1 KLHL17 895966 901095
1 PLEKHN1 901876 910482
1 C1orf170 910578 912021

;--geneList: This file is an output from prepare012s using the option --outputGeneList with columns as chromosome number, gene Name, start position, end position. There should be no header for this file.

THE CHROMOSOME NUMBERS SHOULD BE NUMERICS!!!! 1 - chromosome 1, DO NOT USE chr1.

;--cutoff: This is the minor allele frequency, you can specify it as 0.01, 0.05 or etc.
;--collapseChoice: Specify one of {or,sum,wt}. or: Li-Leal's CMC test, sum: Use the number of rare variants for each subject as the score, wt: Madeson-Browning's weighted rare variant score.

;PHENOTYPE
;--phenofile: A file where the first column is subject ID and the second column is phenotype (0 or 1).

;COVARIATES
;--covConsider: Default = 0, no covariate is considered. 1. covariate is considered.
;--covfile: Covariate file with the first column as subject ID and the other columns are covariates needed to be considered in the model.

;PERMUTATION
;--nPermute: Number of permutation for the evaluation of p values.
;-- PermutationSeed: Default = 1. Can be changed to other numbers too.

;GENE LEVEL TEST RESULT:
;--geneGlobalTestOut: This file stores the 5% and 95% quantiles of the p values for all the genes at each permutation
;--geneTestPvalueFile: This file gives you the gene name, number of rare variants, count of variants in case/control and p values from the RV test specified by collapseChoice.

Template:Other uses

2014-01-29T22:11:42Z