Difference between revisions of "Tutorial: RAREMETAL"

From Genome Analysis Wiki
Jump to navigationJump to search
 
(7 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 +
[[Category:RAREMETAL]]
 +
[[Category:RAREMETALWORKER]]
 
== Useful Wiki Pages ==
 
== Useful Wiki Pages ==
  
Line 38: Line 40:
 
* To analyse the first study, we execute the following command:  
 
* To analyse the first study, we execute the following command:  
  
   raremetalworker  --ped example1.ped --dat example1.dat --vcf example1.vcf.gz --traitName QT1  
+
   raremetalworker  --ped example1.ped --dat example1.dat --vcf example1.vcf.gz --traitName QT1 \
 
                   --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY1  
 
                   --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY1  
  
This command transforms phenotypes to normality (--inverseNormal), calculates trait residuals after adjusting for covariates (--makeResiduals), estimates relatedness between individuals and saves this for later use (--kinSave and --kinGeno), and even generates some annotated PDF files summarizing results (--labelHits).
+
This command transforms phenotypes to normality (--inverseNormal), calculates trait residuals after adjusting for covariates (--makeResiduals), estimates relatedness between individuals and saves this for later use (--kinSave and --kinGeno), and even generates some PDF files summarizing results.
  
 
* After RAREMETALWORKER finishes running, you will see several output files:
 
* After RAREMETALWORKER finishes running, you will see several output files:
Line 53: Line 55:
 
* We next analyse the second study using a similar command:
 
* We next analyse the second study using a similar command:
  
   raremetalworker  --ped example2.ped --dat example2.dat --vcf example2.vcf.gz --traitName QT1  
+
   raremetalworker  --ped example2.ped --dat example2.dat --vcf example2.vcf.gz --traitName QT1 \
 
                     --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY2
 
                     --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY2
  
Line 74: Line 76:
 
   tabix -c "#" -s 1 -b 2 -e 2 STUDY2.QT1.singlevar.cov.txt.gz
 
   tabix -c "#" -s 1 -b 2 -e 2 STUDY2.QT1.singlevar.cov.txt.gz
  
* We next create two text files that will drive the meta-analysis. The first file lists the input files with summary statistics. Let's call it summaryfiles. In most Linux workstations, you can use the command pico or nano to create this file. These should be the contents of "summaryfiles":
+
* We next create two text files that will drive the meta-analysis. The first file lists the input files with summary statistics. Let's call it summaryfiles. In most Linux workstations, you can use the command pico or nano to create this file. These should be the contents of "summaryfiles":
  
 
   STUDY1.QT1.singlevar.score.txt.gz
 
   STUDY1.QT1.singlevar.score.txt.gz
Line 86: Line 88:
 
* Now, we are ready for meta-analysis. To perform single variant and gene-level meta-analyses all at once, use the following command:
 
* Now, we are ready for meta-analysis. To perform single variant and gene-level meta-analyses all at once, use the following command:
  
   raremetal --summaryFiles summaryfiles --covFiles covfiles --groupFile group.file --SKAT --burden --MB --VT  
+
   raremetal --summaryFiles summaryfiles --covFiles covfiles --groupFile group.file \
            --hwe 1.0e-05 --callRate 0.95 --longOutput --tabulateHits --hitsCutoff 1e-05 --prefix COMBINED.QT1 --labelHits
+
            --SKAT --burden --MB --VT --longOutput --tabulateHits --hitsCutoff 1e-05 \
 +
            --prefix COMBINED.QT1 --hwe 1.0e-05 --callRate 0.95 -
  
This command filters summary statistics based on HWE p-value and variant call rate, generates single variant meta-analysis results, generates gene-level meta-analysis results using simple burden test, variable threshold test, Madson-Browning weighted burden test, and SKAT, tabulates significant genes with detailed single variant results included, and even generates some annotated PDF files summarizing results (--labelHits).
+
This command filters summary statistics based on HWE p-value and variant call rate, generates single variant meta-analysis results, generates gene-level meta-analysis results using simple burden test, variable threshold test, Madson-Browning weighted burden test, and SKAT, tabulates significant genes with detailed single variant results included, and even generates some PDF files summarizing results.
  
 
* The following output will be generated:
 
* The following output will be generated:

Latest revision as of 17:51, 16 March 2018

Useful Wiki Pages

There are several pages in this Wiki that may be useful to RAREMETAL users. Here are links to a few:

Introduction

In this tutorial, we will use RAREMETAL to perform single variant and gene-level meta-analysis using summary statistics of two small studies. We will use RAREMETALWORKER to generate summary statistics for each study.

STEP 1: Install Software and Download Example Data Sets

  • To unpack the example dataset, use the following two Unix commands:
 tar xvzf raremetal_tutorial.tar.gz 
 cd raremetal_tutorial

STEP 2: Analyze individual samples using RAREMETALWORKER

  • Our first dataset has 743 unrelated individuals. Their phenotypes and relationships are described in a pair of pedigree and data files; you can get a summary of their contents using pedstats.
  • These individuals have been genotyped at ~1000 markers, and these genotypes are stored in a VCF file. The VCF is a text format, so you can try to peek at the contents if you like.
  • Although these individuals are "unrelated", there is the possibility that some of them are more closely related than others. In our analysis, we will calculate an empirical (i.e. data-driven) kinship (i.e. relatedness) matrix to describe the similarity between individuals.
  • To analyse the first study, we execute the following command:
 raremetalworker  --ped example1.ped --dat example1.dat --vcf example1.vcf.gz --traitName QT1 \
                  --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY1 

This command transforms phenotypes to normality (--inverseNormal), calculates trait residuals after adjusting for covariates (--makeResiduals), estimates relatedness between individuals and saves this for later use (--kinSave and --kinGeno), and even generates some PDF files summarizing results.

  • After RAREMETALWORKER finishes running, you will see several output files:
STUDY1.QT1.singlevar.score.txt     ## single variant statistics
STUDY1.QT1.singlevar.cov.txt       ## covariance matrices between score statistics
STUDY1.plots.pdf                   ## QQ plots and Manhattan plots
STUDY1.Empirical.Kinship.gz        ## Relatedness matrix
STUDY1.singlevar.log               ## Log file
  • We next analyse the second study using a similar command:
  raremetalworker  --ped example2.ped --dat example2.dat --vcf example2.vcf.gz --traitName QT1 \
                   --inverseNormal --makeResiduals --kinSave --kinGeno --prefix STUDY2

STEP 3: Run RAREMETAL for Meta-Analysis

  • In this step, we will use RAREMETAL to combine the summary statistics we just generated.
  • We will first index the result files generated by RAREMETAL. This step relies on bgzip and tabix, two tools that allow rapid indexing and retrieval of results from compressed text files.
 bgzip STUDY1.QT1.singlevar.score.txt
 tabix -c "#" -s 1 -b 2 -e 2 STUDY1.QT1.singlevar.score.txt.gz
 bgzip STUDY1.QT1.singlevar.cov.txt
 tabix -c "#" -s 1 -b 2 -e 2 STUDY1.QT1.singlevar.cov.txt.gz
  • And again for the second study:
 bgzip STUDY2.QT1.singlevar.score.txt
 tabix -c "#" -s 1 -b 2 -e 2 STUDY2.QT1.singlevar.score.txt.gz
 bgzip STUDY2.QT1.singlevar.cov.txt
 tabix -c "#" -s 1 -b 2 -e 2 STUDY2.QT1.singlevar.cov.txt.gz
  • We next create two text files that will drive the meta-analysis. The first file lists the input files with summary statistics. Let's call it summaryfiles. In most Linux workstations, you can use the command pico or nano to create this file. These should be the contents of "summaryfiles":
 STUDY1.QT1.singlevar.score.txt.gz
 STUDY2.QT1.singlevar.score.txt.gz
  • The second file lists the input files with variance-covariance information between markers. Let's call it covfiles. These should be its contents:
 STUDY1.QT1.singlevar.cov.txt.gz
 STUDY2.QT1.singlevar.cov.txt.gz
  • Now, we are ready for meta-analysis. To perform single variant and gene-level meta-analyses all at once, use the following command:
 raremetal --summaryFiles summaryfiles --covFiles covfiles --groupFile group.file \
           --SKAT --burden --MB --VT --longOutput --tabulateHits --hitsCutoff 1e-05 \
           --prefix COMBINED.QT1 --hwe 1.0e-05 --callRate 0.95 -

This command filters summary statistics based on HWE p-value and variant call rate, generates single variant meta-analysis results, generates gene-level meta-analysis results using simple burden test, variable threshold test, Madson-Browning weighted burden test, and SKAT, tabulates significant genes with detailed single variant results included, and even generates some PDF files summarizing results.

  • The following output will be generated:
 COMBINED.QT1.meta.plots.pdf (## QQ plots and manhattan plots)
 COMBINED.QT1.meta.singlevar.results 
 COMBINED.QT1.meta.burden.results
 COMBINED.QT1.meta.SKAT.results
 COMBINED.QT1.meta.VT.results
 COMBINED.QT1.meta.MB.results
 COMBINED.QT1.meta.tophits.SKAT.tbl
 COMBINED.QT1.meta.tophits.VT.tbl
 COMBINED.QT1.meta.tophits.burden.tbl
 COMBINED.QT1.meta.tophits.MB.tbl
 COMBINED.QT1.raremetal.log
  • It is probably a good idea to spend a few minutes reviewing these files - they have the key results for your analysis! A detailed description of output files is available elsewhere.
  • If you are feeling adventurous and are not yet totally confused, you can continue to explore advanced features of RAREMETAL.
  • For example, RAREMETAL can carry out conditional analyses using the summary statistics in summaryfiles and covfiles.
  • To carry out a conditional analysis, create a text file that specifies the variant that you want to condition on (lets call it "conditioningfile") and add "--condition conditioningfile" to the command line. The "conditioningfile" might look like this:
 9:505484545:C:T
  • When you now run RAREMETAL with the extra options --condition conditioningfile, results will be adjusted for the variant 9:505484545:C:T.
  • An alternative to generating a group.file is to provide RAREMETAL with an annotated VCF as input. For details, see documentation. This option is especially useful when you first use RAREMETAL's --writeVCF option to create a file listing a superset of all available variants.