Genome Analysis Wiki - User contributions [en]

Biostatistics 830: Main Page

2013-09-07T03:27:28Z

Ben Lerch: /* Target Audience */

== Objective ==

Gene mapping studies study the relationship between genetic variation and susceptibility to human disease. These studies are changing rapidly with the availability of techniques for very large scale genetic analysis, whether based on sequencing or on genotyping. Biostatistics 830 is a Ph.D. level course that dissects some recently developed methods and the principles behind their implementation. It is meant to provide students with a toolkit to facilitate development and implementation of new statistical methods.

For additional information, see also [[Biostatistics 830: Core Competencies|Core Competencies in Biostatistics Program covered by this course]].

== Target Audience ==

It is highly recommended that students registering for Biostatistics 830 should have previously completed [[Biostatistics 666]] and [[Biostatistics 615/815]], which are courses introducing methods for genetic analysis and programming principles, respectively.

== Scheduling ==

For Fall 2013, classes are scheduled for Mondays and Wednesdays, 3:00 - 4:30 pm.

The final grade will take into account your performance in problem sets and worksheets as well as your participation in class.

== Standards of Academic Conduct ==

The following is an extract from the School of Public Health's Student Code of Conduct [http://www.sph.umich.edu/academics/policies/conduct.html]:

''Student academic misconduct includes behavior involving plagiarism, cheating, fabrication, falsification of records or official documents, intentional misuse of equipment or materials, and aiding and abetting the perpetration of such acts. The preparation of reports, papers, and examinations, assigned on an individual basis, must represent each student’s own effort. Reference sources should be indicated clearly. The use of assistance from other students or aids of any kind during a written examination, except when the use of books or notes has been approved by an instructor, is a violation of the standard of academic conduct.''

In the context of this course, any work you hand-in should be your own and any material that is a transcript (or interpreted transcript) of work by others must be clearly labeled as such.

== Required Reading ==

* Albers CA, Lunter G, MacArthur DG, McVean G, Ouwehand WH, Durbin R (2010) Dindel: accurate indel calls from short-read data. Genome Res. 21:961-73

* Browning SR, Browning BL (2007) Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. ''Am J Hum Genet.'' '''81''':1084-97. PMID: 17924348

* Coventry A, Bull-Otterson LM, Liu X, Clark AG, Maxwell TJ, Crosby J, Hixson JE, Rea TJ, Muzny DM, Lewis LR, Wheeler DA, Sabo A, Lusk C, Weiss KG, Akbar H, Cree A, Hawes AC, Newsham I, Varghese RT, Villasana D, Gross S, Joshi V, Santibanez J, Morgan M, Chang K, Iv WH, Templeton AR, Boerwinkle E, Gibbs R, Sing CF (2010) Deep resequencing reveals excess rare recent variants consistent with explosive population growth. ''Nat Commun.'' '''1''':131. PMID: 21119644

* Delaneau O, Zagury JF, Marchini J (2013) Improved whole-chromosome phasing for disease and population genetic studies. ''Nat Methods.'' '''10''':5-6. PMID: 23269371

* Howie B, Fuchsberger C, Stephens M, Marchini J, Abecasis GR (2012) Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. ''Nat Genet.'' '''44''':955-9. PMID: 22820512

* Iqbal Z, Caccamo M, Turner I, Flicek P, McVean G (2012) De novo assembly and genotyping of variants using colored de Bruijn graphs. ''Nat Genet.'' '''44''':226-32. PMID: 22231483

* Jun G, Flickinger M, Hetrick KN, Romm JM, Doheny KF, Abecasis GR, Boehnke M, Kang HM (2012) Detecting and estimating contamination of human DNA samples in sequencing and array-based genotype data. ''Am J Hum Genet.'' '''91''':839-48. PMID: 23103226

* Li H, Ruan J, Durbin R (2008) Mapping short DNA sequencing reads and calling variants using mapping quality scores. ''Genome Res.'' '''18''':1851-8. PMID: 18714091

* Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. ''Bioinformatics.'' '''25''':1754-60. PMID: 19451168

* Li H, Durbin R (2011) Inference of human population history from individual whole-genome sequences. ''Nature.'' '''475''':493-6. PMID: 21753753

* Li Y, Willer CJ, Ding J, Scheet P, Abecasis GR (2010) MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes. ''Genet Epidemiol.'' '''34''':816-34. PMID: 21058334

* Lin DY, Zeng D (2010) Meta-analysis of genome-wide association studies: no efficiency gain in using individual participant data. ''Genet Epidemiol.'' '''34''':60-6. PMID: 19847795

* Liu et al (2013) http://arxiv.org/abs/1305.1318

* Wen X, Stephens M (2010) Using linear predictors to impute allele frequencies from summary or pooled genotype data. ''Ann Appl Stat.'' '''4''':1158-1182. PMID: 21479081

* Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X (2011) Rare-variant association testing for sequencing data with the sequence kernel association test. Am J Hum Genet. 89:82-93

* Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de Bruijn graphs. ''Genome Res.'' '''18''':821-9. PMID: 18349386

== Course History ==

This course is an ad-hoc course, first taught by Goncalo Abecasis in the Fall of 2013.

QPLOT

2013-08-06T16:53:23Z

Ben Lerch: /* Parameters */

= Introduction =

The qplot program calculates various summary statistics some of which are plotted in a PDF file. These statistics can be used to assess the sequencing quality of sequence reads mapped to the reference genome. The main statistics are empirical Phred scores which are calculated based on the background mismatch rate. Background mismatch rate is the rate that sequenced bases are different from the reference genome, EXCLUDING dbSNP positions. Other statistics include GC biases, insert size distribution, depth distribution, genome coverage, empirical Q20 count, and so on.

In the following sections, we will guide you through: [[#Where to Find It |how to obtain qplot]], [[#Usage |how to use qplot]], [[#Built-in example |example outputs]], [[#anchorOfInteractiveQplot |interactive diagnostic plots]], and [[#Diagnose sequencing quality |real applications]] in which qplot has helped identify sequencing problems.

= Where to Find It =

You can obtain qplot in two ways:

(1) Download the pre-compiled binary along with the source code as described in [[#Binary Download|Binary Download]].

(2) Download source code only and compile it on your own machine. Please follow the instruction in [[#Source Code Distribution|Source Code Distribution]] on fetching source code and building instructions.

== Binary Download ==

We have prepared a pre-compiled (under Ubuntu) qplot along with source code . You can download it from: [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot.20130627.tar.gz qplot.20130627.tar.gz (File Size: 1.7G)]

The executable file is under qplot/bin/qplot.

In addition, we provided the necessary input files under qplot/data/ (NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100).

You can also find an example BAM input file under qplot/example/chrom20.9M.10M.bam. It is taken from the 1000 Genome Project with sequencing reads aligned to chromosome 20 positions 8M to 9M.

== Source Code Distribution ==

We provide a source code only download in [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-source.20130627.tar.gz qplot-source.20130627.tar.gz]. Optionally, you can download example file and/or data file:

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-example.tar.gz example]: example input file, and expected outputs if you following the [[#Built-in example | direction]].

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-data.tar.gz resources data]: necessary input files for qplot, including NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100.

You can put above file(s) in the same folder and follow these steps:

* 1. Unarchive downloaded file
tar zvxf qplot-source.20130627.tar.gz

A new folder ''qplot'' will be created.

* 2. Build libStatGen
cd qplot
(cd ../libStatGen; make cloneLib)

This step will download a necessary software library [http://genome.sph.umich.edu/wiki/C%2B%2B_Library:_libStatGen libStatGen] and compile source code into a binary code library.

* 3. Build qplot
make

This step will then build qplot. Upon success, the executable qplot can be found under qplot/bin/.

* 4. (Optional) unarchive example and/or data
tar zvxf qplot-example.tar.gz

An example file, ''chrom20.9M.10M.bam'', will be extracted to qplot/example/. It contains ~1.1 million aligned Illumina sequencing reads of NA12878 from 1000 Genome Project. Example command line, ''cmd.sh'', example outputs, ''qplot.pdf'', ''qplot.stats'', and ''qplot.R'' are also provided and will be extracted qplot/example/ as well.

tar zvxf qplot-data.tar.gz

Three files will be extracted to qplot/data/: ''human.g1k.v37-bs.umfa'' is binary NCBI reference genome build 37; ''dbSNP130.UCSC.coordinates.tbl'' is dbSNP version 130; and ''human.g1k.w100.gc'' is pre-calculated GC content with windows size 100.



= Usage =

== Command line ==

After you obtain the qplot executable (either by compiling the source code or by downloading the pre-compiled binary file), you will find the executable file under qplot/bin/qplot.

Here is the qplot help page by invoking qplot without any command line arguments:

some_linux_host > qplot/bin/qplot
The following parameters are available. Ones with "[]" are in effect:

References : --reference [/net/fantasia/home/zhanxw/software/qplot/data/human.g1k.v37.fa],
--dbsnp [/net/fantasia/home/zhanxw/software/qplot/data/dbSNP130.UCSC.coordinates.tbl]
GC content file options : --winsize [100]
Region list : --regions [], --invertRegion
Flag filters : --read1_skip, --read2_skip, --paired_skip,
--unpaired_skip
Dup and QCFail : --dup_keep, --qcfail_keep
Mapping filters : --minMapQuality [0.00]
Records to process : --first_n_record [-1]
Lanes to process : --lanes []
Read group to process : --readGroup []
Input file options : --noeof
Output files : --plot [], --stats [], --Rcode [], --xml []
Plot labels : --label [], --bamLabel []
Obsoleted (DO NOT USE) : --gccontent [], --create_gc

== Input files ==

qplot runs on the input BAM/SAM file(s) specified on the command-line after all other parameters.

Additionally, three (3) precomputed files are required.

* <code>--reference</code>

The reference genome is the same as karma reference genome. If the index files do not exist, qplot will create the index files '''automatically''' using the input reference fasta file.

* <code>--dbsnp</code>

This file has two columns. First column is the chromosome name which must be consistent with the reference created above. Second column is 1-based SNP position. If you want to create your own dbSNP data from downloaded UCSC dbSNP file, one way to do it is: <code>cat dbsnp_129_b36.rod|grep "single" | awk '$4-$3==1' |cut -f2,4 > dbSNP_129_b36.tbl</code>

* <code> **OBSOLETED** --gccontent, --create_gc </code>

Although GC content can be calculated on the fly each time, it is much more efficient to load a precomputed GC content from a file.
GC content file name is automatically determined in this format: <reference_genome_base_file_name>.winsize<gc_content_window_size>.gc.
For example, if your reference genome is human.g1k.v37.fa and the window size is 100, then the GC content file name is: human.g1k.v37.winsize100.gc .

As it said, there is no need to use --gccontent to specify GC content file in each run.

* <code> input files </code>

QPLOT take SAM/BAM files.

''Note'': Before running qplot, it is critical to check how the chromosome names are coded. Some BAM/SAM files use just numbers, others use chr + numbers. '''You need to make sure that the chromosome names from the reference and dbSNP are consistent with the BAM/SAM files.'''

== Parameters ==

Some of the command line parameters are described here, but most are self explanatory.

*Flag filter

By default all reads are processed. If it is desired to check only the first read of a pair, use <code>--read2_skip</code> to ignore the second read. And so on.

*Duplication and QCFail

By default reads marked as duplication and QCFail are ignored but can be retained by
--dup_keep
or
--qcfail_keep

*Records to process

The <code>--first_n_record</code> option followed by a number, '''n''', will enable qplot to read the first '''n''' reads to test the bam files and verify it works.

* Lanes to process (only works for Illumina sequences)

If the input bam files have more than one lane and only some of them need to be checked, use something like <code>--lanes 1,3,5</code> to specify that only lanes 1, 3, and 5 need to be checked.

'''NOTE''' In order for this to work, the lane info has to be encoded in the read name such that the lane number is the second field with the delimiter ":".

* Read group to process :

Read group option can restrict qplot to process a subset of reads. For example, if BAM contain the following @RG tags:

@RG ID:UM0348_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM

If specify nothing or not using "--readGroup", QPLOT by default will process all reads;
If specify "--readGroup UM0348", then only read group UM0348_1, UM_0348_2, UM_0348_3, UM_0348_4 will be processed;
If specify "--readGroup UM0348_1", then only one read group UM0348_1 will be processed.

* Input file options :

BAM files are compress by BGZF algorithm and it should contain EOF by default. QPLOT will by default stop working when it does not found a valid EOF tag inside BAM files.
However, you can force QPLOT to continue process using --noeof. But you should be aware the input files may be corrupted.

* Mapping filters

Qplot will exclude reads with lower mapping qualities than the user specified parameter, <code>--minMapQuality</code>. By default, mapped reads with all mapping quality will be included in the analysis.

*Region list

If the interest of qplot is a list of regions, e.g. exons, this can be achieved by providing a list of regions. The regions should be in the form of "chr start end label" each line in the file (NOTE: ''start'' and ''end'' position are inclusive and they follow the convention of [http://genome.ucsc.edu/FAQ/FAQformat#format1 BED file]).
In order for this option to work, within each chromosome (contig) the regions have to be sorted by starting position, and also the input bam files have to be sorted.
For example, you can create a text file, region.txt like following:

1 100 500 region_A
1 600 800 region_B
2 100 300 region_C

Then specifying <code> --regions region.txt</code> enables qplot to calculate various statistics out of sequenced bases only within the above 3 regions.

Qplot also provides the <code>--invertRegion</code> option. Enabling this option tells qplot to operate on those sequence bases that are outside the given region.

* Plot labels

Two kinds of labels are enabled. <code>--label</code> is the label for the plot (default is empty) which is appended to the title of each subplot. <code>--bamLabels</code> followed by a column separated list of labels provides the labels for each input SAM/BAM file, e.g. sample ID (default is numbers 1, 2, ... until the number of input bam files). For example:
--label Run100 --bamLabels s1,s2,s3,s4,s5,s6,s7,s8

== Output files ==

There are three (optional) output files.
* <code>--plot ''qa.pdf''</code>

Qplot will generate a PDF file named ''qa.pdf'' containing 2 pages each with 4 figures. The plot is generated using Rscript.

* <code>--stats ''qa.stats''</code>

Qplot will generate a text file named ''qa.stats'' containing various summary statistics for each input BAM/SAM file.

* <code>--Rcode ''qa.R''</code>

Qplot will generate ''qa.R'' which is the R code used for plotting the figures in the ''qa.pdf'' file. If Rscript is not installed in the system, you can use the qa.R to generate the figures on other machines, or extract plotting data from each run and combine multiple runs together to generate more comprehensive plots (See [[#Example | Example]]).

= Example =

Qplot can generate diagnostic graphs, related R code, and summary statistics for each SAM/BAM file.

== Built-in example ==

In the pre-compiled binary download, you will find a subdirectory named examples. We provide a sample file from the 1000 Genome project, it contains aligned reads on chromosome 20 from position 8 Mbp to 9Mbp. You can invoke qplot using the following command line:

../bin/qplot --reference ../data/human.g1k.v37.umfa --dbsnp ../data/dbSNP130.UCSC.coordinates.tbl --gccontent ../data/human.g1k.w100.gc --plot qplot.pdf --stats qplot.stats --Rcode qplot.R --label "chr20:9M-10M" chrom20.9M.10M.bam

Sample outputs are listed below:

1) Figure: [[Media:qplot.pdf | qplot.pdf]]

2) Summary statistics:
Stats\BAM chrom20.9M.10M.bam
TotalReads(e6) 1.11
MappingRate(%) 97.24
MapRate_MQpass(%) 97.24
TargetMapping(%) 0.00
ZeroMapQual(%) 2.39
MapQual<10(%) 2.86
PairedReads(%) 83.76
ProperPaired(%) 71.34
MappedBases(e9) 0.04
Q20Bases(e9) 0.04
Q20BasesPct(%) 88.63
MeanDepth 42.22
GenomeCover(%) 0.03
EPS_MSE 1.81
EPS_Cycle_Mean 18.71
GCBiasMSE 0.01
ISize_mode 137
ISize_medium 184
DupRate(%) 5.90
QCFailRate(%) 0.00
BaseComp_A(%) 29.9
BaseComp_C(%) 20.1
BaseComp_G(%) 20.2
BaseComp_T(%) 29.8
BaseComp_O(%) 0.1

== Gallery of examples ==

Here we show qplot can be applied in various sequencing scenarios. Also users can customize statistics generated by qplot to their needs.

* Whole genome sequencing with 24-multiplexing

With a customized script, we aggregated 24 bar-coded samples in the same graph.
The graph will help compare sequencing quality between samples.

[[Media: qplot.Pool.9847.pdf | QPlot of 24 samples(PDF) ]]

* Interactive qplot


Qplot can be interactive. In the following example, you can use mouse scroll to zoom in and zoom out on each graph and pan to a certain part of the graph.
By presenting qplot data on a web page, users can easily identify problematic sequencing samples. Users of qplot can customize its outputs into web page format greatly easing the data exploring process.

[http://www-personal.umich.edu/~zhanxw/qplot.Pool.9847.html QPlot of 24 samples(HTML) ]

== Diagnose sequencing quality ==

Qplot is designed and implemented for the need of checking sequencing quality.
Besides the example of analyzing RNA-seq data as shown in our manuscript,
here we demonstrate two additional scenarios in which qplot can help identify problems after obtaining sequencing data.

* Base quality distributed abnormally

[[Media: WrongBaseQual.pdf | Example of qplot helping to identify wrong phred base quality]]

By checking the first graph "Empirical vs reported Phred score", we found reported base qualities are shifted to the right.
In this particular example, '33' was incorrectly added to all base qualities.
When such data used in variant calling, we may increase false positive SNP variants.

* Bar-coded samples

[[Media: WrongBarCoding.pdf | Example of qplot identifying the effect of ignoring bar-coding]]

By checking "Empirical phred score by cycle" (top right graph on the first page), we noticed the empirical qualities in the first several cycles are abnormally low. This phenomenon leads us to hypothesize that the first several bases have different properties. Further investigation confirmed that this sequencing was done using bar-coded DNA samples, but the analysis did not properly de-multiplex each sample.

= Contact =

Questions and requests should be sent to Bingshan Li ([mailto:bingshan@umich.edu bingshan@umich.edu]) or Xiaowei Zhan ([mailto:zhanxw@umich.edu zhanxw@umich.edu]) or Goncalo Abecasis ([mailto:goncalo@umich.edu goncalo@umich.edu])

EPACTS

2013-08-02T17:12:13Z

Ben Lerch: /* Creating marker group file */

'''EPACTS''' (Efficient and Parallelizable Association Container Toolbox) is a versatile software pipeline to perform various statistical tests for identifying genome-wide association from sequence data through a user-friendly interface, both to scientific analysts and to method developers.

== Join in EPACTS mailing list ==

Please join in the [http://groups.google.com/group/epacts EPACTS Google Group] to ask / discuss / comment about EPACTS.

== Lastest ChangeLog ==
* March 25th, 2013 : EPACTS v3.2.3 release
** Relaxed the checking of low-rank matrix in SKAT tests (to avoid unncessary skipping of genes)
* March 13th, 2013 : EPACTS v3.2.2 release
** Fixed an error which occasionally report mismatches in the number of samples
* March 9th, 2013 : EPACTS v3.2.1 release
**Fixed errors in loading the dynamic library
** Fixed errors in SKAT-O (thanks to Anubha Mahajan and Jason Flannick)
** Fixed bugs in emmax-CMC
** Added emmax-SKAT (contributed by Seunngeun Lee)
** And additional minor bug fixes
See [[#Full ChangeLog]] for full details

== Key Features ==

EPACTS currently provides the following set of key features
* Robust support for widely used format of sequence-based genotypes (VCF) and phenotypes with pedigree (PED)
** Efficient library for accessing VCF file to reduce computational burden to analyze large-scale sequencing data
** Support selecting markers by arbitrary combination of substring matching.
** Support for using genotype dosages instead of hard genotype calls
** Utilize PED format to perform test across multiple traits.
* Supports a large number of widely used statistical tests for single variant association and burden tests.
** See the "Currently Supported Statistical Tests" section below for more information
* Easy to Highly Parallelize Jobs
** Makefile-based partition into and ligation of multiple subtasks
** Parallel run of job is simply adding one parameter when running EPACTS
* Integrative and versatile framework that allows easy addition of additional statistical test
** Core input/output routines are implemented in C++
** Most statistical tests (except for EMMAX) are implemented in R
** Adding a simple R function to implement additional statistical test (See [[#Implementing Additional Statistical Tests]] for details)
* Useful utilities for post-association-analysis tasks
** Automatic functional annotation of associated variants
** Automatic generation of QQ and Manhattan Plot
** (TBA) Zoom plot for the significant associations

== Obtaining EPACTS ==

* The official release of EPACTS software is available at http://www.sph.umich.edu/csg/kang/epacts/ .
** From the CSG cluster, it is available at /net/fantasia/home/bin/epacts/
* Note that R (version 2.10 or higher) and gnuplot (version 4.2 or higher) must be installed in order to run EPACTS correctly.

== Currently Supported Statistical Tests ==

EPACTS supports the following sets of widely used statistical tests for single variant tests and burden tests

=== Single Variant Tests ===

<noinclude>
{|</noinclude> border="1" cellpadding="4" cellspacing="0" style="margin: 1em 1em 1em 0; background: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse; font-size: 95%; clear: center;"<noinclude>
|-
! Test Name
! Phenotypes
! Covariates
! Computational Time
! Description
| Implemented by
|-
| b.wald
| Binary
| YES (Joint)
| Slow
| Logisitic Wald Test
| Hyun Min Kang (simply used glm in R)
|-
| b.score
| Binary
| YES (Regressed Out)
| Fast
| Logistic Score Test (from Lin DY and Tang ZZ, AJHG 2011 89:354-67)
| Clement Ma & Hyun Min Kang
|-
| b.firth
| Binary
| YES (Joint)
| Slow
| Firth Bias-Corrected Logistic Likelihood Ratio Test
| Clement Ma
|-
| b.lrt
| Binary
| YES (Joint)
| Slow
| Likelihood Ratio Test
| Clement Ma
|-
| b.glrt
| Binary
| NO
| Fast
| Genotype Likelihood Ratio Test (use GL or PL field in VCF to perform case-control test)
| Hyun Min Kang
|-
| q.lm
| Quantitative
| YES (Joint)
| Slow
| Linear Wald Test
| Hyun Min Kang (as implemented in lm in R)
|-
| q.score
| Quantitative
| YES (Regressed Out)
| Fast
| Quantitative Score Test (from Lin DY and Tang ZZ, AJHG 2011 89:354-67)
| Clement Ma
|-
| q.linear
| Quantitative
| YES (Regressed Out)
| Fast
| Linear Wald Test
| Hyun Min Kang
|-
| q.reverse
| Quantitative
| YES (Joint)
| Slow
| Reverse regression of phenotypes on binary genotypes (dominant model)
| Hyun Min Kang
|-
| q.wilcox
| Quantitative
| YES (Regressed Out First)
| Slow
| Nonparametric Reverse regression of phenotypes on binary genotypes (dominant model)
| Hyun Min Kang
|-
| q.emmax
| Quantitative
| YES (Regressed Out First)
| Slow
| EMMAX ( Kang et al (2010) Nat Genet 42:348-54 )
| Hyun Min Kang
|-
|}

=== Gene-wise or group-wise tests ===

<noinclude>
{|</noinclude> border="1" cellpadding="4" cellspacing="0" style="margin: 1em 1em 1em 0; background: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse; font-size: 95%; clear: center;"<noinclude>
|-
! Test Name
! Phenotypes
! Covariates
! Computational Time
! Description
| Implemented by
|-
| b.collapse
| Binary
| YES (Joint Estimation)
| Slow
| Logistic Wald Test between binary phenotypes and 0/1 collapsed variables
| Hyun Min Kang
|-
| b.madsen
| Binary
| NO
| Slow
| Wilcoxon Rank Sum Test between binary phenotypes and weighted rare variant scores (slightly different version from the published method - it uses pooled allele frequency across cases and controls for weighting each variant)
| Hyun Min Kang
|-
| b.wcnt
| Binary
| YES (Joint Estimation)
| Slow
| Logistic Wald Test between binary phenotypes and weighted rare variant scores
| Hyun Min Kang
|-
| q.reverse
| Quantitative
| YES (Joint Estimation)
| Slow
| Reverse regression of phenotypes on binary collapsed variables
| Hyun Min Kang
|-
| q.wilcox
| Quantitative
| YES (Regressed Out First)
| Slow
| Nonparametric Reverse regression of phenotypes on collapsed variables
| Hyun Min Kang
|-
| skat
| Binary/Quantitative
| YES (Joint Estimation)
| Slow
| SKAT-O Test by Lee et al, Biostatistics (2012)
| Seunggeun Lee (adaptive by Xueling Sim and Hyun Min Kang)
|-
| VT
| Binary/Quantitative
| YES (Regressed out first)
| Slow
| Variable Threshold Test with adaptive permutation Price et al, AJHG (2010) 86:832-8
| Hyun Min Kang
|-
| emmaxCMC
| Binary/Quantitative
| YES (Regressed Out First)
| Slow
| Collapsing burden test using EMMAX
| Hyun Min Kang
|-
| emmaxVT
| Binary/Quantitative
| YES (Regressed Out First)
| Slow
| Variable-threshold burden test using EMMAX
| Hyun Min Kang
|-
| mmskat
| Quantitative
| YES (Regressed Out First)
| Slow
| SKAT test using EMMAX
| Seunggeun Lee & Hyun Min Kang
|}

== Installation Details ==

If you want to use EPACTS in an Ubuntu platform, following the step below

*Download EPACTS source distribution at http://www.sph.umich.edu/csg/kang/epacts/download/EPACTS-3.2.3.tar.gz (100MB)
*Uncompress EPACTS package, and install the package using the following set of commands

tar xzvf EPACTS-3.2.3.tar.gz
cd EPACTS-3.2.3
./configure --prefix=/path/to/install
make
make install

(Important Note: '''make sure to specify --prefix=/path/to/install''' to avoid installing to the default path /usr/local/, which you may not have the permission. /home/your_userid/epacts might be a good one, if you are not sure where to install)

* Now ${EPACTS_DIR} represents the '/path/to/install' directory

* Download the reference FASTA files from 1000 Genomes FTP automatically by running the following commands

${EPACTS_DIR}/bin/epacts download

(For advanced users, to save time for downloading the FASTA files (~900MB), you may copy a local copy of GRCh37 FASTA file and the index file to ${EPACTS_DIR}/share/EPACTS/)

*Perform a test run by running the following command

${EPACTS_DIR}/bin/test_run_epacts.sh

In order to use EPACTS in the CSG cluster, you do not need to install them. You can directly use or make a copy of the in-house release version at

/net/fantasia/home/hmkang/bin/epacts/

== Getting Started With Examples ==
If you are using EPACTS from the CSG cluster, please set the following environment variable
EPACTS_DIR=/net/fantasia/home/hmkang/bin/epacts (in bash)
setenv EPACTS_DIR /net/fantasia/home/hmkang/bin/epacts (in csh)

If you downloaded EPACTS binary and please set EPACTS_DIR to the full path of the downloaded and uncompressed directory.

=== All-in-one example ===

To get started with EPACTS, run the following command will perform an example run
${EPACTS_DIR}/bin/test_run_epacts.sh

You will find a series of lines in test_run_epacts.sh script commented out for each possible test.

The example phenotype (PED format) and genotype (VCF format) can be found at
${EPACTS_DIR}/share/EPACTS/

=== Single Variant Test ===

Or You can run EPACTS command yourself by running
${EPACTS_DIR}/epacts single \
--vcf ${EPACTS_DIR}/data/1000G_exome_chr20_example_softFiltered.calls.vcf.gz \
--ped ${EPACTS_DIR}/data/1000G_dummy_pheno.ped \
--min-maf 0.001 --chr 20 --pheno DISEASE --cov AGE --cov SEX --test b.score --anno \
--out out/test --run 2

The command above will perform single variant association test using a dummy case-control phenotype file and a subset of 1000 genomes exome VCF file (chr20) using score test statistic for all variants over 1% of higher MAF using 2 parallel runs.

You will see the 4 output files as the main outcome of the analysis

==== Output Text of All Test Statistics ====

The filename is out/test.single.b.score.epacts.gz and the contents will look like
$ zcat out/test.single.b.score.epacts.gz | head
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE N.CASE N.CTRL AF.CASE AF.CTRL
20 68303 68303 20:68303_A/G_Upstream:DEFB125 266 1 1 0.0018797 NA NA NA NA NA NA
20 68319 68319 20:68319_C/A_Upstream:DEFB125 266 1.4467e-36 1 0 NA NA NA NA NA NA
20 68396 68396 20:68396_C/T_Nonsynonymous:DEFB125 266 1 1 0.0018797 NA NA NA NA NA NA
20 76635 76635 20:76635_A/T_Intron:DEFB125 266 1.534e-37 1 0 NA NA NA NA NA NA
20 76689 76689 20:76689_T/C_Synonymous:DEFB125 266 0 1 0 NA NA NA NA NA NA
20 76690 76690 20:76690_T/C_Nonsynonymous:DEFB125 266 1 1 0.0018797 NA NA NA NA NA NA
20 76700 76700 20:76700_G/A_Nonsynonymous:DEFB125 266 0 1 0 NA NA NA NA NA NA
20 76726 76726 20:76726_C/G_Nonsynonymous:DEFB125 266 0 1 0 NA NA NA NA NA NA
20 76771 76771 20:76771_C/T_Nonsynonymous:DEFB125 266 3 1 0.0056391 0.68484 0.40587 145 121 0.013793 0.0082645

==== Output Text of Top Associations ====

Same type of file but containing top 5,000 association will be stored at out/test.epacts.top5000

$ head out/test.single.b.score.epacts.top5000
#CHROM BEGIN END MARKER_ID NS AC CALLRATE MAF PVALUE SCORE N.CASE N.CTRL AF.CASE AF.CTRL
20 1610894 1610894 20:1610894_G/A_Synonymous:SIRPG 266 138.64 1 0.26061 6.9939e-05 3.9765 145 121 0.65177 0.36476
20 4162411 4162411 20:4162411_T/C_Intron:SMOX 266 204 1 0.38346 0.00055583 -3.4523 145 121 0.62759 0.93388
20 34061918 34061918 20:34061918_T/C_Intron:CEP250 266 41.815 1 0.0786 0.00095471 3.3035 145 121 0.22543 0.075436
20 4155948 4155948 20:4155948_G/A_Intron:SMOX 266 215 1 0.40414 0.0020792 -3.0787 145 121 0.68276 0.95868
20 4680251 4680251 20:4680251_A/G_Nonsynonymous:PRNP 266 186 1 0.34962 0.0025962 3.0119 145 121 0.8069 0.57025
20 36668874 36668874 20:36668874_G/A_Synonymous:RPRD1B 266 96 1 0.18045 0.003031 2.9646 145 121 0.44828 0.2562
20 36641871 36641871 20:36641871_G/A_Synonymous:TTI1 266 10 1 0.018797 0.004308 -2.8547 145 121 0.0068966 0.07438
20 1616892 1616892 20:1616892_A/G_Synonymous:SIRPG 266 144 1 0.27068 0.0051239 2.7991 145 121 0.63449 0.42975
20 25038372 25038372 20:25038372_G/A_Intron:ACSS1 266 103.3 1 0.19418 0.005748 2.7618 145 121 0.47201 0.28813

The key columns represents:
* '''NS''' : Number of phenotyped samples with non-missing genotypes
* '''AC''' : Total Non-reference Allele Count
* '''CALLRATE''' : Fraction of non-missing genotypes.
* '''MAF''' : Minor allele frequencies
* '''PVALUE''' : P-value of single variant test
* '''AF.CASE''' : Non-reference allele frequencies for cases
* '''AF.CTRL''' : Non-reference allele frequencies for controls

==== Q-Q plot of test statistics (stratified by MAF) ====

The file out/test.b.score.epacts.qq.pdf will be generated as shown below

[[File:test_b_score_epacts_qq.png]]

==== Manhattan Plot of Test Statistics ====

The file out/test.b.score.epacts.mh.pdf will be generated for chr20 only.

[[File:test_b_score_epacts_mh.png]]

An example Genome-wide manhattan plot (from a genome-wide run) will look like below

[[File:tes_b_score_epacts_mh_gw.png]]

=== Gene-wise or group-wise burden test ===

Gene-wise or group-wise burden test requires two steps. First, 'group' file containing the list of
markers per group needs to be generated. Second, group-wise burden test needs to be run

==== Creating marker group file ====

The marker group file has the following format

[GROUP_ID] [MARKER_ID_1] [MARKER_ID_2] .... [MARKER_ID_N]

where
* [GROUP_ID] is a string representing the group (e.g. gene name)
* [MARKER_ID_K] is a marker key as a format of [CHROM]:[POS]_[REF]/[ALT] (NOTE THAT THIS IS DIFFERENT FROM TYPICAL VCF MARKER ID field)

Note that [MARKER_ID_K] has to be sorted by increasing order of genomic coordinate

In order to create gene-level group file from typically formatted VCF file, one may use the following utility

${EPACTS_DIR}/epacts make-group --vcf [input-vcf] --out [output-group-file] --format [epacts, annovar, chaos or gatk] --nonsyn

The above command create a file [output-group-file] containing a list of missense and nonsense variants per each gene. To incorporate different types of functional annotations, use --type option as follows

${EPACTS_DIR}/epacts make-group --vcf [input-vcf] --out [output-group-file] --format [epacts, annovar, chaos or gatk] --type [function_type_1] --type [function_type_2] ...

Type 'epacts makegroup -man' for the detailed documentation

==== Annotating VCF file using EPACTS ====

If the VCF is not annotated, 'epacts makegroup' cannot be used. In order to annotate VCF, one can use the example VCF using ANNOVAR as follows:

${EPACTS_DIR}/epacts anno \
--in ${EPACTS_DIR}/data/1000G_exome_chr20_example_softFiltered.calls.vcf.gz \
--out ${EPACTS_DIR}/data/1000G_exome_chr20_example_softFiltered.calls.anno.vcf.gz

The epacts anno script will add "ANNO=[function]:[genename]" entry into the INFO field based on gencodeV7 (default) or refGene database.

It is important to check whether the VCF file is already annotated or not in order to avoid no or redundant annotation.

==== Running Groupwise Test ====

To perform a groupwise burden test on the example VCF (annotated as above), run the following command

${EPACTS_DIR}/epacts group --vcf ${EPACTS_DIR}/data/1000G_exome_chr20_example_softFiltered.calls.anno.vcf.gz \
--groupf ${EPACTS_DIR}/data/1000G_exome_chr20_example_softFiltered.calls.anno.grp --out out/test.gene.skat \
--ped ${EPACTS_DIR}/data/1000G_dummy_pheno.ped --maxAF 0.05 \
--chr 20 --pheno QT --cov AGE --cov SEX --test skat --skat-o --run 2

==== Example Output ====
$ head out/test.gene.skat.epacts.top5000
#CHROM BEGIN END MARKER_ID NS FRAC_WITH_RARE NUM_ALL_VARS NUM_PASS_VARS NUM_SING_VARS PVALUE STATRHO
20 62607037 62608720 20:62607037-62608720_SAMD10 266 0.14662 9 5 1 0.0020064 1
20 2816211 2820493 20:2816211-2820493_FAM113A 266 0.011278 12 2 1 0.0032542 0
20 47245987 47361692 20:47245987-47361692_PREX1 266 0.1391 54 9 6 0.0054849 1
20 34761734 34810279 20:34761734-34810279_EPB41L1 266 0.071429 14 7 5 0.0068492 0.2
20 61340671 61391602 20:61340671-61391602_NTSR1 266 0.11278 24 9 3 0.011063 1
20 48561952 48568644 20:48561952-48568644_RNF114 266 0.011278 4 2 1 0.015175 0.2
20 60962895 60963559 20:60962895-60963559_RPS21 266 0.06015 6 3 2 0.016409 0
20 55904961 55917801 20:55904961-55917801_SPO11 266 0.011278 11 3 3 0.018031 0

The key columns represents:
* '''NS''' : Number of phenotyped samples with non-missing genotypes
* '''FRAC_WITH_RARE''' : Fraction of individual carrying rare variants below --max-maf (default : 0.05) threshold.
* '''NUM_ALL_VARS''' : Number of all variants defining the group.
* '''NUM_PASS_VARS''' : Number of variants passing the --min-maf, --min-mac, --max-maf, --min-callrate thresholds
* '''NUM_SING_VARS''' : Number of singletons among variants in NUM_PASS_VARS
* '''PVALUE''' : P-value of burden tests
* Other columns are test specific auxiliary columns. For example, in the VT test, the optimal MAF threshold is recorded as an auxiliary output column.

=== Specialized Instruction for EMMAX tests ===

EMMAX (Efficient Mixed Model Association eXpedited - Kang et al (2010) Nat Genet 42:348-54) is an efficient implementation of mixed model association accounting for sample structure including population structure and hidden relatedness. Currently EPACTS supports EMMAX association mapping in single variant test and CMC-like burden tests.

Because EMMAX is based on linear model, the method fits better to quantiative traits than binary traits. However, p-values for binary traits are expected to be valid in the spirit of Armitage trend test, although the estimated effect size may not be precise.

In order to run EMMAX analysis from sequence-based genotypes. We recommend running EPACTS multiple times using the following procedure.

==== Single Variant EMMAX Association Analysis ====

* '''Creating Kinship Matrix''' : From VCF, we recommend to set a MAF (e.g. 0.01) and call rate (e.g. 0.95) threshold to select high-quality markers to generate kinship matrix as follows.
${EPACTS_DIR}/epacts make-kin \
--vcf [input.vcf.gz] --ped [input.ped (Optional)] --min-maf 0.01 --minCallRate 0.95 \
--sepchr (if VCF is separated by chromosome) --out [outprefix.kinf] --run [# of parallel jobs]

If you provide [input.ped] file, then it will calculate the subset the individuals contained in the PED file.

The procedure above will create a file [outprefix.kinf] after splitting and merging the genomes into multiple pieces. If only a certain subset of SNPs needs to be considered due to target regions, LD-pruning, or any other reasons, a VCF containing the subset of markers must be created beforehand and should be used as input VCF file.

* '''Perform Single Variant Association''' : From VCF and PED, we recommend to use less stringent MAF threshold (e.g. 0.001) and call rate (e.g. 0.50) to perform single variant association
${EPACTS_DIR}/epacts single \
--vcf [input.vcf.gz] --ped [input.ped] --min-maf 0.001 --kin [outputprefix.kinf] \
--sepchr --pheno [PHENO_NAME] --cov [COV1] --cov [COV2] --test q.emmax \
--out [outprefix] --run [# of parallel jobs]

The procedure above will perform single variant association analysis compatible to other types of single variant association analyses implemented in EPACTS

==== Burden-style EMMAX Association Analysis ====

In order to run EMMAX analysis from sequence-based genotypes. We recommend running EPACTS multiple times using the following procedure.

* '''Creating Kinship Matrix''' : See 'Creating Kinship Matrix' section in [[#Single Variant EMMAX Association Analysis]]
* '''Create Marker Group'''
** By annotating the VCF and extracting missense and nonsense variants
*** [[#Annotating VCF file using ANNOVAR]] - This step will be required to create marker group file
*** [[#Creating marker group file]] - Assume that [group.grp] file is produced
** Or, by creating your own marker group information
*** See [[#Creating marker group file]] for details
* Run CMC-style burden test by
${EPACTS_DIR}/epacts group --groupf [group.grp] \
--vcf [input.vcf.gz] --ped [input.ped] --max-maf [max-MAF-for-rare-variants] \
--kin [outputprefix.kinf] --sepchr --pheno [PHENO_NAME] --cov [COV1] --cov [COV2] \
--test emmaxCMC --out [outprefix]
* Run Variable Threshold burden test by
${EPACTS_DIR}/epacts group --groupf [group.grp] \
--vcf [input.vcf.gz] --ped [input.ped] --max-maf [max-MAF-for-rare-variants] \
--kin [outputprefix.kinf] --sepchr --pheno [PHENO_NAME] --cov [COV1] --cov [COV2] \
--test emmaxVT --out [outprefix]

== Preparing Your Own Input Data ==

=== VCF file for Genotypes ===

EPACTS support VCF files as input for association with the following requirement
* Input VCF file must be bgzipped and tabixed before running association to allow efficient random access of the file. Below is an example command to conver plain VCF into bgzipped and tabixed VCF
bgzip input.vcf ## this command will produce input.vcf.gz
tabix -pvcf -f input.vcf.gz ## this command will produce input.vcf.gz.tbi
* If the VCF file is separated by chromosome, the VCF file must contain the string "chr1" in the chromosome 1 file, and corresponding chromosome name for other chromosomes.
* Sample IDs in the VCF file must be consistent to those from PED file
* Currently EPACTS only support bi-allelic variants, but it handles SNPs, INDELs, snd SVs.
* Currently, EPACTS only support VCF aligned with NCBI build 37 coordinates
* An example VCF file from 1000 genome project is below.
$ zcat example/1000G_integrated_phase1_chr20.vcf.gz | cut -f 1-10 | head -50
##fileformat=VCFv4.1
##INFO=<ID=LCSNP,Number=0,Type=Flag,Description="Genotype likelihood in Low coverage VCF in data integration">
##INFO=<ID=EXSNP,Number=0,Type=Flag,Description="Genotype likelihood in Exome VCF in data integration">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Genotype likelihood in INDEL VCF in data integration">
##INFO=<ID=SV,Number=0,Type=Flag,Description="Genotype likelihood in SV VCF in data integration">
##INFO=<ID=BAVGPOST,Number=1,Type=Float,Description="Average posterior probability from beagle">
##INFO=<ID=BRSQ,Number=1,Type=Float,Description="Genotype imputation quality estimate from beagle">
##INFO=<ID=LDAF,Number=1,Type=Float,Description="MLE Allele Frequency Accounting for LD">
##INFO=<ID=AVGPOST,Number=1,Type=Float,Description="Average posterior probability from MaCH/Thunder">
##INFO=<ID=RSQ,Number=1,Type=Float,Description="Genotype imputation quality from MaCH/Thunder">
##INFO=<ID=ERATE,Number=1,Type=Float,Description="Per-marker Mutation rate from MaCH/Thunder">
##INFO=<ID=THETA,Number=1,Type=Float,Description="Per-marker Transition rate from MaCH/Thunder">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=HOMLEN,Number=.,Type=Integer,Description="Length of base pair identical micro-homology at event breakpoints">
##INFO=<ID=HOMSEQ,Number=.,Type=String,Description="Sequence of base pair identical micro-homology at event breakpoints">
##INFO=<ID=SOURCE,Number=.,Type=String,Description="Source of deletion call">
##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##INFO=<ID=AC,Number=.,Type=Integer,Description="Alternate Allele Count">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total Allele Count">
##ALT=<ID=DEL,Description="Deletion">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DS,Number=1,Type=Float,Description="Genotype dosage from MaCH/Thunder">
##FORMAT=<ID=GL,Number=.,Type=Float,Description="Genotype Likelihoods">
##FORMAT=<ID=BD,Number=1,Type=Float,Description="Genotype dosage from beagle">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00096
20 60479 . C T 100 PASS LCSNP;EXSNP;BAVGPOST=1.000;BRSQ=0.894;LDAF=0.0020;AVGPOST=0.9995;RSQ=0.8779;ERATE=0.0005;THETA=0.0008;AC=4;AN=2184 GT:DS:GL:BD 0|0:0.000:-0.19,-0.46,-2.68:0.0022
20 60522 . T TC 1588 PASS INDEL;BAVGPOST=1.000;BRSQ=0.994;LDAF=0.0116;AVGPOST=0.9980;RSQ=0.9327;ERATE=0.0004;THETA=0.0167;AC=24;AN=2184 GT:DS:GL:BD 0|0:0.000:0.00,-0.90,-9.20:0
20 60571 . C A 100 PASS LCSNP;EXSNP;BAVGPOST=0.999;BRSQ=0.813;LDAF=0.0029;AVGPOST=0.9986;RSQ=0.8085;ERATE=0.0014;THETA=0.0014;AC=5;AN=2184 GT:DS:GL:BD 0|0:0.000:-0.05,-0.96,-5.00:0.0008
20 60795 . G C 100 PASS LCSNP;EXSNP;BAVGPOST=1.000;BRSQ=0.930;LDAF=0.0006;AVGPOST=0.9996;RSQ=0.7205;ERATE=0.0003;THETA=0.0041;AC=1;AN=2184 GT:DS:GL:BD 0|0:0.000:-0.03,-1.21,-5.00:0.0001
20 60810 . G GA 127 PASS INDEL;BAVGPOST=1.000;BRSQ=0.862;LDAF=0.0013;AVGPOST=0.9987;RSQ=0.5684;ERATE=0.0004;THETA=0.0061;AC=2;AN=2184 GT:DS:GL:BD 0|0:0.000:0.00,-1.80,-18.80:0

=== PED file for Phenotypes and Covariates ===

EPACTS accepts a PED format supported by MERLIN or PLINK software to represent phenotypes. For example, the example.ped file and example.dat file can represent the phenotypes and corresponding column name (from 6th column and after).

$ head example.ped
13281 NA12344 NA12347 NA12348 1 1 94.17 66.1
13281 NA12347 0 0 1 1 109.54 44.0
13281 NA12348 0 0 2 2 119.40 46.6
1328 NA06984 0 0 1 2 87.72 39.3
1328 NA06989 0 0 2 1 100.60 41.7
1328 NA12329 NA06984 NA06989 2 1 100.85 46.4
13291 NA06986 0 0 1 2 91.94 61.9
13291 NA06995 NA07435 NA07037 1 2 104.36 57.4
13291 NA06997 NA06986 NA07045 2 2 107.53 53.1

$ cat example.dat
A DISEASE
T QT
T AGE

EPACTS also accept a PED format with header information. The above file can be combined into one file as follows

$ head data/1000G_dummy_pheno.ped
#FAM_ID IND_ID FAT_ID MOT_ID SEX DISEASE QT AGE
13281 NA12344 NA12347 NA12348 1 1 94.17 66.1
13281 NA12347 0 0 1 1 109.54 44.0
13281 NA12348 0 0 2 2 119.40 46.6
1328 NA06984 0 0 1 2 87.72 39.3
1328 NA06989 0 0 2 1 100.60 41.7
1328 NA12329 NA06984 NA06989 2 1 100.85 46.4
13291 NA06986 0 0 1 2 91.94 61.9
13291 NA06995 NA07435 NA07037 1 2 104.36 57.4
13291 NA06997 NA06986 NA07045 2 2 107.53 53.1

The column names can be used to identify the names of phenotypes and covariates in the analysis.

== Frequently Asked Questions ==
=== Installation ===
# How should I install EPACTS?
#* See [[EPACTS#Installation_Details | Installation Details]]
# I am having the following error message '''configure: error: libR.{so,a} was not found. Please install it at http://www.r-project.org/ first'''. What do I have to do?
#* First, you need to find out where R was installed. Try to type "locate libR.so" and see if it returns anything
#* If "locate libR.so" returns you something, as explained [[EPACTS#Installation_Details | Installation Details]], try to add "LDFLAGS=-L/path/to/R/library" and rerun '''configure''' and '''make'''
#* If you cannot find libR.so, you make have to recompile R with --enable-R-shlib option as described in http://cran.r-project.org/doc/manuals/R-admin.html#Installation

=== Input Files ===
# What is VCF?
#* VCF refers to Variant Call Format
#* See [[http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41 1000 Genomes wiki page]] for the detailed description of VCF format
# Should input VCF be compressed into certain format?
#* Correct. EPACTS assumes that VCF file is bgzipped and tabixed already.
#* See [[#VCF file for Genotypes]] for details.
# What are the additional requirements for input VCF file?
#* Input VCF file used for association mapping must contain individual genotype information at 10-th or higher order columns.
#* GT field must be encoded as haploid or diploid
#* Bi-allelic SNPs only : Currently EPACTS may not handle multi-allelic SNPs correctly.
#* If non-GT field is used, the field is considered as dosage and should be a single numeric value.
# What are the acceptable input format to encode phenotypes and covariates?
#* See [[#PED file for Phenotypes and Covariates]] for the detailed information
# How should I encode binary phenotypes?
#* If you encode your phenotypes into two different numeric values (e.g. 0/1 or 1/2), EPACTS will automatically recognize them as binary phenotypes and encode them into 1/2 values. Higher value will be considered as cases for case-control association
# How should I encode missing genotypes?
#* The default code missing phenotypes in EPACTS are 'NA'
#* One may use --missing option to specify different types of missing values
#* The encoding of missing genotypes follows the VCF specificiation
# How do I match the relationship between VCF and PED input files?
#* EPACTS will assume that the individual IDs in each VCF and PED file are unique, and they follow the saming convention. Thus, the individual IDs overlapping between VCF and PED files will be considered in the associations
# How the individuals with missing phenotypes are handled?
#* Currently, EPACTS will automatically remove the individuals without phenotypes or covariates. If one wants to use imputed covariates to increase sample size, the PED file must contain the imputed covariate values.
#* Markers with missing genotypes won't be discarded automatically. It can be explicitly discarded by --minCallRate option when performing association

=== Output Files ===
# Which output files should I be looking at?
#* [[#Output Text of Top Associations]] is the key file to look at the individual top associations
#* [[#Q-Q plot of test statistics (stratified by MAF)]] will be important to see the global distribution of test statistics and examine if there are apparent inflation of test statistics
#* [[#Manhattan Plot of Test Statistics]] will inform us the genome-wide distribution of association signals
#* [[#Output Text of All Test Statistics]] will contain the full information of test results across all units tested
# The Q-Q and Manhattan plots cannot be found. Why?
#* It is probably because gnuplot 4.2 or higher is not installed in your system, or they are included but cannot be found in your ${PATH}. Please visit [[http://gnuplot.info/ GNUPLOT web page]] for installation.

=== More questions ===
# If you have more questions, please contact [[mailto:hmkang@umich.edu Hyun Min Kang]].

== Detailed Options ==

The detailed options can viewed by running the following commands
${EPACTS_DIR}/bin/epacts -man (for overall structure)
${EPACTS_DIR}/bin/epacts single -man (for single variant test)
${EPACTS_DIR}/bin/epacts group -man (for groupwise test)
${EPACTS_DIR}/bin/epacts anno -man (for annotation)
${EPACTS_DIR}/bin/epacts plot -man (for QQ and Manhattan plot)
${EPACTS_DIR}/bin/epacts zoom -man (for zoom plot)
${EPACTS_DIR}/bin/epacts meta -man (for meta-analysis)
${EPACTS_DIR}/bin/epacts makegroup -man (for creating gene group)

== Implementing Additional Statistical Tests ==

In order to add additional statistical test to EPACTS, the following procedure are recommended

# Create a file named 'single.[testname].R' for single variant test or 'gene.[testname].R' for gene-level test under ${EPACTS_DIR}/share/EPACTS/
# Test your implementation using --test [testname] option to perform sanity check and debugging
# If you want to add your test in the official in-house version, please send your code to Hyun

Below is an example of a single variant test implementation ( single.q.lm.R )
## Core functions of EPACTS to perform association

##################################################################
## SINGLE VARIANT TEST
## INPUT VARIABLES:
## n : total # of individuals
## NS : number of called samples
## AC : allele count
## MAF : minor allele frequency
## vids : indices from 1:nrow(NS) after AF/AC threshold
## genos : genotype matrix (after AF/AC threshold)
## EXPECTED OUTPUT : list(p, addcols, addnames) for each genos row
## p : p-value
## add : additional columns to add
## cname : column names for additional columns
##################################################################

## single.lm() : Use built-in lm() function to perform association
## KEY FEATURES : SIMPLE, BUT MAY BE SLOW
## GOOD SNIPPLET TO START A NEW FUNCTION
## TRAITS : QUANTITATIVE
## RETURNS : PVALUE, BETA, SEBETA, TSTAT
## MISSING VALUES : IGNORED
single.q.lm <- function() {
cname <- c("BETA","SEBETA","TSTAT") # column names for additional variables in the EPACTS output
m <- nrow(genos)
p <- rep(NA,m)
add <- matrix(NA,m,3) ## BETA, SEBETA, TSTAT
if ( m > 0 ) {
for(i in 1:m) {
r <- summary(lm(pheno~genos[i,]+cov-1))$coefficients[1,] # run simple linear regression
p[i] <- r[4] # store p-value to p[i]
add[i,] <- r[1:3] # store additional variables to add[i,]
}
}
return(list(p=p,add=add,cname=cname))
}

As described in the comment, you may assume that the following variables are available for use for testing association across m markers
* n (scalar) : total number of individuals
* NS (M * 1 vector) : Number of called samples for each marker
* AC (M * 1 vector) : Non-reference allele count for each marker
* MAF (M * 1 vector) : Minor allele frequency
* vids (m * 1 vector) : indices of markers passing the inclusion criteria (e.g. MAF threshold) among 1:M
* genos (m * n matrix) : genotype matrix as a input for association test

The output variables to generate is as follows
* p (m * 1 vector) : p-value matrix as output
* add (m * c matrix) : additional columns as output of test (such as SCORE, BETA, etc)
* cname (c * 1 vector) : column names of add

In the output files, the following columns will be displayed
# MARKER : Marker ID
# NS : Number of called samples
# AC : Non-ref allele count
# CALLRATE : Call rate = NS/n
# MAF : Minor allele frequency
# PVALUE : P-values
# Additional columns specified by return values 'add'

Below is an example of a gene-lvel variant test implementation ( single.q.lm.R )

##################################################################
## GENE-LEVEL BURDEN TEST
## INPUT VARIABLES:
## n : total # of individuals
## genos : genotype matrix for each gene
## NS : number of called samples for each marker
## AC : allele count for each marker
## MAC : minor allele count for each marker
## MAF : minor allele frequency
## vids : indices from 1:n after AF/AC threshold
## EXPECTED OUTPUT : list(p, addcols, addnames) for each genos row
## p : p-value
## add : additional column values
## cname : additional column names
##################################################################

## gene.q.reverse() : Reverse logistic regression
## KEY FEATURES : 0/1 collapsing variable ~ rare variants
## TRAITS : QUANTITATIVE (GAUSSIAN)
## RETURNS : PVALUE, BETA, SEBETA, ZSTAT
## MISSING VALUE : IMPUTED AS MAJOR ALLELES
gene.q.reverse <- function() {
cname <- c("BETA","SEBETA","ZSTAT")
m <- nrow(genos)
if ( m > 0 ) {
g <- as.double(colSums(genos,na.rm=T) > 0)
sg <- sum(g)
if ( ( sg > 0 ) && ( sg < n ) ) {
r <- glm(g~pheno+cov-1,family=binomial)
if ( ( r$converged ) && ( ! r$boundary ) ) {
return(list(p=summary(r)$coefficients[1,4],
add=summary(r)$coefficients[1,1:3],
cname=cname))
}
}
}
return(list(p=NA,add=rep(NA,3),cname=cname))
}

Similar to gene-level test, you may assume the following variables exist for testing A SINGLE GENE. Note that M is the number of markers spanning the gene region

* n (scalar) : total number of individuals
* NS (M * 1 vector) : Number of called samples for each marker
* AC (M * 1 vector) : Non-reference allele count for each marker
* MAC (M * 1 vector) : Minor allele count
* MAF (M * 1 vector) : Minor allele frequency
* vids (m * 1 vector) : indices of markers passing the inclusion criteria (e.g. MAF threshold) among 1:M
* genos (m * n matrix) : genotype matrix as a input for association test

The output variables to generate is as follows
* p (scalar) : p-value matrix as output
* add (c * 1 vector) : additional columns as output of test (such as SCORE, BETA, etc)
* cname (c * 1 vector) : column names of add

In the output files, the following columns will be displayed
# MARKER : Marker ID
# NS : Number of called samples
# MAF_BURDEN : MAF of 0/1 collapsing variables (existence of rare variants)
# NUM_ALL_VARS : Number of all variants within the gene
# NUM_RARE_VARS : Number of rare variants below the max-MAF threshold
# NUM_SING_VARS : Number of singleton variants
# PVALUE : P-value from the test
# Additional columns specified by return values 'add'

== Full ChangeLog ==
* March 25th, 2013 : EPACTS v3.2.3 release
** Relaxed the checking of low-rank matrix in SKAT tests (to avoid unncessary skipping of genes)
* March 13th, 2013 : EPACTS v3.2.2 release
** Fixed an error which occasionally report mismatches in the number of samples
* March 9th, 2013 : EPACTS v3.2.1 release
**Fixed errors in loading the dynamic library
** Fixed errors in SKAT-O (thanks to Anubha Mahajan and Jason Flannick)
** Fixed bugs in emmax-CMC
** Added emmax-SKAT (contributed by Seunngeun Lee)
** And additional minor bug fixes
* February 28th, 2013 : EPACTS v3.2.0 release
** R package installation bug (for some users) was fixed
** A bug in the MAF error for high frequency variants (AF>0.25) was now fixed
** SKAT version is updated to 0.81
** --bprange option is added to allow testing for small region size
** Additional minor bug fixes
* December 4th, 2012 : EPACTS v3.1.0 release
** Removed dependency on libR.so
** Additional minor bug fixes
** --bprange option is added to allow testing for small region size
** November 25th, 2012 : EPACTS v3.0.0 release
** Restructured with source code release (with autoconf / automake / libtools)
** Added zoom plot feature
** FRAC_BURDEN keyword was replace to FRAC_WITH_RARE for groupwise testing
* October 26th, 2012 : EPACTS v2.2.0-beta is released with the following updates
** Added --max-mac option
** Fixed Firth's bias-corrected test (by Clement Ma)
** Added more informative warning messages when index files do not exist
** Fixed the bug in the epacts-plot in plotting ties
** Fixed errors in the MAF estimates per case and control
** Fixed bug in --minRSQ option
* September 28, 2012 : EPACTS v2.11-beta is released with the following updates
** Counts and allele frequencies for case/control added for binary tests
** --max-maf parameter is added
** Fixed EMMAX error in MAF in the output
** More informative error messages
* September 27, 2012 : EPACTS v2.1-beta is released with the following updates
** EMMAX interface is changed. --kinOnly option is related with a new command '''make-kin'''
** SKAT-O is upgraded to version 0.77 with additional configurable parameter settings
** Some parameter names are renamed (e.g. --min-maf, --min-mac)
** Many minor bugs are fixed
* Jul 6, 2012 : EPACTS v2.01-beta is released with the following updates
** SKAT-O is upgraded to version 0.76
** Fixed minor bugs in option names (Thanks to Xueling Sim)
* Jul 3, 2012 : EPACTS v2.0-beta is released with the following updates
** Major restructuring of the software
** Annotation software is switched with built-in application
** Addition of SKAT-O and EMMAX burden test
** Minor bug fixes
* Apr 8, 2012 : EPACTS v1.2-alpha is released with the following updates, in addition to the following updates
** EMMAX bug in handling covariates was fixed
** Variable Threshold Test is added
** Variable Threshold Test with genomic score (e.g. GERP or PhyloP) is added.
* Apr 4, 2012 : EPACTS v1.1-alpha is released with the following updates, in addition to minor updates
** EMMAX burden test (Hyun Min Kang)
** Likelihood ratio test (Clement Ma)
** Updated version of Firth bias-corrected likelihood ratio test (Clement Ma)
** Updated version of EMMAX single variant test (Hyun Min Kang)
* Mar 29, 2012 : EPACTS v1.0-alpha is released

VerifyBamID

2013-08-02T14:40:03Z

Ben Lerch: /* Command Line Options */

[[Category:Software|VerifyBamID]]

'''verifyBamID''' is a software that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals), and checks whether the reads are contaminated as a mixture of two samples. '''verifyBamID''' can detect sample contamination and swaps when external genotypes are available. When external genotypes are not available, '''verifyBamID''' still robustly detects sample swaps.

== Download verifyBamID ==

To get a copy go to the [http://www.sph.umich.edu/csg/kang/verifyBamID/download VerifyBamID Download] download page.

== Join in verifyBamID mailing list ==

Please join in the [http://groups.google.com/group/verifybamid VerifyBamID Google Group] to ask / discuss / comment about verifyBamID.

== What's new ==

(2012/06/20)
* Fixed a bug of incorrect estimate of contamination when --chip-full option was used (Thanks to Richard Smith)
* Fixed a bug of incorrect per-readgroup output in --chip-* parameter

(2012/05/24)
* Fixed a bug of incorrect per-readgroup output (Thanks to Matthew Flickinger)
* '''(IMPORTANT)''' Add an option to remove either side of overlapping fragment. This option is turned on by default, and can be turned off usig --ignoreOverlapPair. If your sequence data has very short insert size, this update may increase the sensitivity of estimated contamination.
* Changes in the directory structure and Makefile

(2012/05/18) The new release of verifyBamID have undergone major change since the last version (as of 2011 April). Here are the highlights
* The genotype / allele frequency file is now based on VCF format rather than PLINK format.
* The reference sequence information is no longer required
* Uses Brent's method for precise estimation of contamination parameters
* Generate the depth distribution statistics.
* Estimated reference-bias parameters (useful mostly for ABI SOLiD sequence data)

== Build verifyBamID ==

The binary download of verifyBamID is available. You may use the version in Ubuntu 64-bit platform. To build verifyBamID, download the statgen library and run the following series of commands
tar xzvf verifyBamID.20120620.tar.gz
cd verifyBamID
make cloneLib (in the case ../libStatGen does not exist)
make
./bin/verifyBamID

Note that '''git clone''' command will create a directory ./libStatGen under your working directory, and '''make''' will create binary of verifyBamID under verifyBamID/bin/

verifyBamID is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, and later, there are likely other portability issues.

Currently we support verifyBamID only on Ubuntu 9.10 and later on 64-bit processors.

== Basic Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. ''verifyBamID'' checks whether reads in a BAM file match previous genotypes for a specific sample. In addition, it detects possible sample mixture from population allele frequency only, which can be particularly useful when the genotype data is not available.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, ''verifyBamID'' tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

verifyBamID --vcf [input.vcf] --bam [input.bam] --out [output.prefix] --verbose --ignoreRG

where
[input.bam] is a BAM (Binary Alignment Map) file of a sequence reads
[input.vcf] is input VCF file containing individual genotypes or AF or AC/AN fields in the INFO field. gzipped VCF is also allowed.
[outPrefix] is output prefix of output files - [outPrefix].{selfRG,selfSM,bestRG,bestSM,depthRG,depthSM} will be created.

More detailed description of command line input is below

== Preparing input files ==

verifyBamID requires two input files - VCF file containing external genotypes or allele frequency information, and the BAM file.

=== VCF input genotype file ===

The input VCF file contains (1) external genotype information and/or (2) allele frequency information as AF entry or AC/AN entries in the INFO field. (See [http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41 | VCF specification] for further details). If neither information is provided, verifyBamID will not work properly.

If external genotype information is provided, sequence+array method will identify contamination and sample swaps by comparing the concordance between the external genotypes and the sequence reads. Additionally, sequence-only method will provide additional contamination estimates by modeling the sequence reads as mixture of two unknown samples based on the allele frequency information in the VCF file.

Input VCF file needs to meet several additional contraints need to meet in order to properly run verifyBamID.
* The VCF is assumed to be well-formed. For example, verifyBamID does not check whether REF allele actually matches with reference sequence.
* The VCF should only contain SNPs. Current version of verifyBamID does not accept INDELs, MNPs, Structural Variations, or other complex variants.
* The individual IDs in the VCF file, must be identical with the individual identifier in the BAM file. Otherwise, --smID option can override the sample ID information of the BAM file to the ID that matches to the individual IDs in the VCF file.
* IMPORTANT : For targeted sequencing data, it is important to subselect the markers to only include on-target markers in the genotype file. Off-target markers are not likely to have multiple non-duplicated reads at the marker position, and it may create artifacts in the analysis due to overlapping fragments.
* Currently, verifyBamID takes only autosomal chromosomes as input VCF.

An example input VCF file (without external genotype) is provided below. Note that AC and AC entries exists in the INFO field for the allele frequency information.

#CHROM POS ID REF ALT QUAL FILTER INFO
20 61651 SNP20-9651 C A . PASS CR=99.86851;GentrainScore=0.7055;HW=0.077647716;AN=2180;AC=11
20 63231 SNP20-11231 T G . PASS CR=99.93036;GentrainScore=0.7837;HW=0.035481825;AN=2182;AC=275
20 63244 rs6139074 A C . PASS CR=98.893394;GentrainScore=0.8001;HW=7.327299E-7;AN=2162;AC=501
20 63799 rs1418258 C T . PASS CR=99.75217;GentrainScore=0.8170;HW=0.6653377;AN=2182;AC=881

=== Input BAM file ===

verifyBamID requires a sorted, indexed, base quality recalibrated, and duplication-marked BAM file. It also requires to contain "@RG" header lines to annotation different readGroups (sequencing runs and lanes). The SM tag in the "@RG" header should match with one of the genotyped sample. Otherwise, verifyBamID may not be able to test whether the sequenced sample matches with genotyped sample, but will try to detect sample mixture from allele frequency, and will try to detect the best-matching sample among the genotyped sample.

== What the default option does ==

The default option of '''verifyBamID''' is the recommended setting for the most sequencing studies to provide a rapid and informative response. The default option provides the following features:
* --free-mix is turned on for estimating contamination using sequence-only method
* --chip-mix is turned on for estimating contamination or swap using sequence+array method, if the external genotype file is provided in the VCF
* --self is turnd on : The default option does not try to compare the sequence reads to identify the best matching individual (which is possible with --best option). It only compares with the external genotypes from the same individual to the sequenced individual.
* --maxDepth 20 is used without --precise option : The default option is intended for whole genome low coverage sequencing. For the targeted exome sequencing, --maxDepth 1000 and --precise is recommended.
* --ignoreRG is not a default option, but a recommended option, when you want to check the contamination for the entire BAM rather than examining each read group separately. This option will increase the computational efficiency especially in the case whether the sequence reads are multiplexed across many sequencing runs.

== Interpreting output files ==

See also [[Understanding VerifyBamID output]].

=== Output files ===
When verifyBamID runs successfully, the following sets of files may be generated.
* [outPrefix].selfSM - Per-sample statistics describing how well the sample matches to the annotated sample.
* [outPrefix].depthSM - The depth distribution of the sequence reads per sample
* [outPrefix].selfRG - Per-readGroup statistics describing how well each lane matches to the annotated sample. (available only without --ignoreRG option)
* [outPrefix].depthRG - The depth distribution of the sequence reads per readGroup. (available only without --ignoreRG option)
* [outPrefix].bestSM - Per-sample best-match statistics with best-matching sample among the genotyped sample (available only with --best option)
* [outPrefix].bestRG - Per-readgroup best-match statistics with best-matching sample among the genotyped sample (available only with --best and without --ignoreRG option)

=== Column information in the output files ===
The .selfSM/.selfRG/.bestSM/.bestRG files have the following 19 columns per sample, or per readgroup (lane).

# SEQ_SM : Sample ID of the sequenced sample. Obtained from @RG header / SM tag in the BAM file
# RG : ReadGroup ID of sequenced lane. For [outPrefix].selfSM and [outPrefix].bestSM, these values are "ALL"
# CHIP_ID : Sample ID compared to in the genotype file. For [outPrefix].selfRG and [outPrefix].selfSM, these values should be identical to [SEQ_SM] or "NA" if the genotype of sequenced samples are unavailable. For [outPrefix].bestRG and [outPrefix].bestSM, these values should be the ID of best-matching sample among the genotype files compared to.
# # SNPs : # of SNPs passing the criteria from the VCF file
# # READS : Total # of reads loaded from the BAM file
# # AVG_DP : Average sequencing depth at the sites in the VCF file
# FREEMIX : Sequence-only estimate of contamination (0-1 scale)
# FREELK1 : Maximum log-likelihood of the sequence reads given estimated contamination under sequence-only method
# FREELK0 : Log-ikelihood of the sequence reads given no contamination under sequence-only method
# FREE_RH : Estimated reference bias parameter Pr(refBase|HET) (when --free-refBias or --free-full is used)
# FREE_RA : Estimated reference bias parameter Pr(refBase|HOMALT) (when --free-refBias or --free-full is used)
# CHIPMIX : Sequence+array estimate of contamination (NA if the external genotype is unavailable) (0-1 scale)
# CHIPLK1 : Maximum log-likelihood of the sequence reads given estimated contamination under sequence+array method (NA if the external genotypes are unavailable)
# CHIPLK0 : Log-likelihood of the sequence reads given no contamination under sequence+array method (NA if the external genotypes are unavailable)
# CHIP_RH : Estimated reference bias parameter Pr(refBase|HET) (when --chip-refBias or --chip-full is used)
# CHIP_RA : Estimated reference bias parameter Pr(refBase|HOMALT) (when --chip-refBias or --chip-full is used)
# DPREF : Depth (Coverage) of HomRef site (based on the genotypes of (SELF_SM/BEST_SM), passing mapQ, baseQual, maxDepth thresholds.
# RDPHET : DPHET/DPREF, Relative depth at Heterozygous site.
# RDPALT : DPHET/DPREF, Relative depth at HomAlt site.

=== A guideline to interpret output files ===

verifyBamID provides a series of information that is informative to determine whether the sample is possibly contaminated or swapped, but there is no single criteria that works for every circumstances. There are a few unmodeled factor in the estimation of [SELF-IBD]/[BEST-IBD] and [%MIX], so please note that the MLE estimation may not always exactly match to the true amount of contamination. Here we provide a guideline to flag potentially contaminated/swapped samples

* Each sample or lane can be checked in this way. When [CHIPMIX] >> 0.02 and/or [FREEMIX] >> 0.02, meaning 2% or more of non-reference bases are observed in reference sites, we recommend to examine the data more carefully for the possibility of contamination.
* We recommend to check each lane for the possibility of sample swaps. When [CHIPMIX] ~ 1 AND [FREEMIX] ~ 0, then it is possible that the sample is swapped with another sample. When [CHIPMIX] ~ 0 in .bestSM file, [CHIP_ID] might be actually the swapped sample. Otherwise, the swapped sample may not exist in the genotype data you have compared.
* When genotype data is not available but allele-frequency-based estimates of [FREEMIX] >= 0.03 and [FREELK1]-[FREELK0] is large, then it is possible that the sample is contaminated with other sample. We recommend to use per-sample data rather than per-lane data for checking this for low coverage data, because the inference will be more confident when there are large number of bases with depth 2 or higher.

== Command Line Options ==

The following parameters are available. Ones with "[]" are in effect:

Available Options
Input Files : --vcf [], --bam [], --subset [],
--smID []
VCF analysis options : --genoError [1.0e-03],
--minAF [0.01],
--minCallRate [0.50]
Individuals to compare with chip data : --site, --self, --best
Chip-free optimization options : --free-none, --free-mix [ON],
--free-refBias, --free-full
With-chip optimization options : --chip-none, --chip-mix [ON],
--chip-refBias, --chip-full
BAM analysis options : --ignoreRG, --ignoreOverlapPair,
--noEOF, --precise, --minMapQ [10],
--maxDepth [20], --minQ [13],
--maxQ [40], --grid [0.05]
Modeling Reference Bias : --refRef [1.00], --refHet [0.50],
--refAlt [0.00]
Output options : --out [], --verbose

Each option provides the following features:
* --vcf : specify required VCF file
* --bam : specify required BAM file (indexed with .bam.bai or .bai file)
* --subset : list of individual IDs to calculate the allele frequency. All individuals are used if unspecified
* --smID : If the individual ID in the BAM file and VCF file does not match, substitute the BAM file's ID into the specified argument
* --genoError : error rate of the external genotype file
* --minAF : minimum allele frequency of the markers to include
* --minAF : minimum call rate of the markers to include
* --site : If set, use only site information in the VCF and do not compare with the actual genotypes
* --self : Only compare the ID-matching individuals between the VCF and BAM file
* --best : Find the best matching individuals (.bestSM and .bestRG files will be produced). This option is substantially longer than the default option
* --free-none : Do not perform sequence-only method to estimate parameters
* --free-mix : (default) Estimate contamination using sequence-only method with Brent's single dimensional optimization.
* --free-refBias : Estimate the reference bias parameters using sequence-only method with Simplex method
* --free-full : Estimate both reference bias parameters and the contamination parameters using sequence-only method
* --chip-none : Do not perform sequence+array method to estimate parameters
* --free-mix : (default) Estimate contamination using sequence+array method with Brent's single dimensional optimization.
* --free-refBias : Estimate the refernece bias parameters using sequence+array method with Simplex method
* --free-full : Estimate both reference bias parameters and the contamination parameters using sequence+array method
* --ignoreRG : ignore the read grouup level comparison and compare samples only (recommended for an expedited run)
* --ignoreOverlapPair : ignore overlapping pair end fragment covering the same base. Disabling this option may decrease the sensitivity of the method when the insert size is short (with slight gain in the computational speed)
* --noEOF : do not check the EOF marker of the BAM file (for earlier version of BAM)
* --precise : calculate the likelihood in log-scale for high-depth data (recommended when --maxDepth is greater than 20. Can be a little bit slower)
* --minMapQ : minimum mapping quality of the sequence reads to compare
* --minQ : minimum base quality to include
* --maxQ : maximum base quality to cap
* --grid : the grid interval to search the optimum before running Brent's algorithm.
* --refRef : Initial Pr(refBase|HOMREFGeno) parameter
* --refHet : Initial Pr(refBase|HETGeno) parameter
* --refAlt : Initial Pr(refBase|HOMALTGeno) parameter
* --out : output file prefix (required)
* --verbose : print the progress of the method on the screeen

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== Reference ==

Please cite the following paper:

G. Jun, M. Flickinger, K. N. Hetrick, Kurt, J. M. Romm, K. F. Doheny, G. Abecasis, M. Boehnke,and H. M. Kang, ''Detecting and Estimating Contamination of Human DNA Samples in Sequencing and Array-Based Genotype Data'', American journal of human genetics doi:10.1016/j.ajhg.2012.09.004 (volume 91 issue 5 pp.839 - 848)

== Contamination in Array Data ==

[[VerifyIDintensity]] or [[BAFRegress]] can estimate sample contamination from Illumina genotype array data.

== Acknowledgements ==

VerifyBamID is a result from collaborative effort by Hyun Min Kang, Goo Jun, Matthew Flickinger, Mary Kate Wing, and Goncalo Abecasis. Please email to Hyun Min Kang [[mailto:hmkang@umich.edu| hmkang@umich.edu ]] for any questions.

GotCloud: Variant Calling Pipeline

2013-07-01T14:41:03Z

Ben Lerch: /* Results */

Back to parent: [[GotCloud]]

The Variant Calling Pipeline (previously called 'UMAKE') makes genotype calls from recalibrated BAM files. These genotype calls are output into [http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41 VCF (Variant Call Format) files].

= Running the GotCloud Variant Calling Pipeline =

The variant calling pipeline (umake) is run using <code>gotcloud snpcall</code> and <code>gotcloud ldrefine</code>. <code>gotcloud</code> is found under <code>gotcloud/</code>.

==Running the Automatic Test==

The automatic test runs the variant calling pipeline on a small testset and checks the results against expected results validating that GotCloud is installed correctly.

*Run variant calling pipeline test:
gotcloud snpcall --test OUTPUT_DIR
where OUTPUT_DIR is the directory where you want to store the test results

If you see "Successfully ran the test case, congratulations!", then you are ready to align samples.

= Overview of Variant Calling Pipeline Steps =
Here is an overview of the Variant Calling Pipeline:

[[File: umakeSteps.png]]

= Input Data=
*Aligned/Processed/Recalibrated BAM files
*Index file containing Sample IDs & BAM file names
*Reference files
*(Optional) Configuration file to override default options

== BAM files ==
The BAM files need to be duplicate-marked and base-quality recalibrated in order to obtain high quality SNP calls. Generating these BAM files from original FASTQs is documented elsewhere as part of the [[Alignment Pipeline]] of gotCloud.

== Index File ==
Each line of the index file represents each individual under the following format. Note that multiple BAMs per individual may be provided. Note that if all samples are from the same population, just specify "ALL" for the population label for each sample.
[SAMPLE_ID] [COMMA SEPARATED POPULATION LABELS] [BAM_FILE1] [BAM_FILE2] ...

Columns:
# sample id
# comma separated population labels
# BAM File 1 (preferable to have full paths to BAM files)
# BAM File 2 (if applicable)
:...

: # BAM File N

== Reference Files ==
The variant calling pipeline requires multiple reference files in order to work correctly.

* Reference Sequence in fasta format.
** Configuration File Setting: <code>REF = path/file.fa</code>
* Indel VCF File Prefix
** Configuration File Setting: <code>INDEL_PREFIX = path/indels.sites.hg19</code>
** <code>path/</code> contains <code>indels.sites.hg19.chr20.vcf</code> for each chromosome being processed
* DBSNP File vcf.gz file (must be indexed with tabix)
** Configuration File Setting: <code>DBSNP_VCF = path/dbsnp_135.b37.vcf.gz</code>
** <code>path/</code> contains <code>dbsnp_135_b37.rod.chr20.map</code> for each chromosome being processed
* HapMap3 polymorphic site vcf.gz file (must be indexed with tabix)
** Configuration File Setting: <code>HM3_VCF = path/hapmap_3.3.b37.sites.vcf.gz<code>
** <code>path/</code> contains <code>hapmap3.qc.poly.chr20.bim</code> & <code>hapmap3.qc.poly.chr20.frq</code> for each chromosome being processed

A set of reference files can be downloaded from: [[ftp://share.sph.umich.edu/1000genomes/umake-resources/ | FTP Download of Full Resource Files]]

Configuration File Example Reference Settings:
REF = path/file.fa
INDEL_PREFIX = path/indels.sites.hg19
DBSNP_VCF = path/dbsnp_135_b37.vcf.gz
HM3_VCF = path/hapmap_3.3.b37.sites.vcf.gz

== Configuration File ==
Configuration file contains the run-time options including the software binaries and command line arguments. A default configuration file is automatically loaded. Users must specify their own configuration file specifying just the values different than the defaults.

Comments begin with a <code>#</code>

Format: KEY = value

Where KEY is the item being set and value is its new value

===Required User Config Files Settings===
The following Config File Settings must be specified by the user:
* CHRS = space separated list of chromosomes you want
* BAM_INDEX = path to the Index File of BAMs

===Required on Command-Line or in Config File===
The following Command-Line or Config File Settings must be specified by the user:
* --outdir/OUT_DIR= path to desired output directory

===Targeted/Exome Sequencing Settings===
If you are running Targeted/Exome Sequencing, the user should specify:
* Write loci file when performing pileup
** WRITE_TARGET_LOCI = TRUE
* Specify the output sub-directory to store target information, for example: targetDir
** Should not be a full path as this will co under the OUT_DIR directory.
** TARGET_DIR = targetDir

If all individuals have the same target:
* Specify the single bed file, for example: target.bed
** UNIFORM_TARGET_BED = target.bed

If not all individuals have the same target:
* Specify the file containing the sample id -> bed map, for example: targetMap.txt
** MULTIPLE_TARGET_MAP = targetMap.txt
*** Each line of the file contains [SM_ID] [TARGET_BED]

Optional Settings:
* Extend the target region by a given number of bases, for example: 50
** OFFSET_OFF_TARGET = 50

=== Configure Reference Files ===
See [[#Reference Files| Reference Files]] for information on how to specify the reference files.

=== Chromosome X Calling ===
Making calls on the X chromosome requires the user to specifty a PED file with sex information.
* PED_INDEX = pedfile.ped

== Example Configuration File ==
Example configuration file where reference files happen to be stored in /path/reference, and bam index file in path/freeze5
CHRS = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
BAM_INDEX = /path/freeze5/freeze5.bam.index ### The BAM index file described above
OUT_DIR = /path/freeze5/output ### Directory in which to put all gotcloud output
REF = /path/reference/hs37d5.fa ### Reference sequence
INDEL_PREFIX = /path/reference/1kg.pilot_release.merged.indels.sites.hg19 ### Known indel sites
HM3_VCF = /path/reference/hapmap3_r3_b37.sites.vcf.gz ### HapMap variants (requires tabix index file in same directory)
DBSNP_VCF = /path/reference/dbsnp_135.b37.sites.vcf.gz ### dbSNP variants (requires tabix index file in same directory)

= Running =

Running variant calling is straightforward:

<code>
'''gotcloud snpcall --conf vc.conf --numjobs 2
'''gotcloud ldrefine --conf vc.conf --numjobs 2
</code>

Replace vc.conf with the approprate path/name of the user's configuration file.

If <code>OUT_DIR</code> is not defined in the configuration file, add <code>--outdir</code> followed by the path to the user's desired output directory.

Update the value following <code>--numjobs</code> to the appropriate number of jobs to be run in parallel.

== Running on a Cluster ==
To run on the Cluster, the following settings need to be added to the configuration file:

<--- Following may need revision --->
TODO: COMING SOON
SLEEP_MULT = 20
REMOTE_PREFIX = # REMOTE_PREFIX : Set if cluster node see the directory differently (e.g. /net/mymachine/[original-dir])
<--- End: Following may need revision --->

Here's the same configuration file we used above but now made to run on a cluster computer with MOSIX.
== Example Configuration File ==
CHRS = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
BAM_INDEX = /path/freeze5/freeze5.bam.index
OUT_DIR = /path/freeze5/output
REF = /path/reference/hs37d5.fa
INDEL_PREFIX = /path/reference/1kg.pilot_release.merged.indels.sites.hg19
HM3_VCF = /path/reference/hapmap3_r3_b37.sites.vcf.gz
DBSNP_VCF = /path/reference/dbsnp_135.b37.sites.vcf.gz
BATCH_TYPE = mosix ### Specify MOSIX as the batch system
BATCH_OPTS = -j10,11,12,13 ### Specify available MOSIX compute nodes

= Results =

If there is a failure, you should see a message like:
make: *** [...] Error 1
Where ... is filled in with other text indicating what step failed.

On SNP Call success, you should see the following output sub-directories under your output directory:
* glfs with a bams & samples subdirectory
* pvcfs with a subdirectory per chromosome and then per region
* split with a subdirectory per chromosome
* vcfs with a subdirectory per chromosome
* (optionally your target directory)

Under the vcf/chrXX directory, there should be:
* chrXX.filtered.sites.vcf
* chrXX.filtered.sites.vcf.norm.log
* chrXX.filtered.sites.vcf.summary
* chrXX.filtered.vcf.gz
* chrXX.filtered.vcf.gz.OK
* chrXX.filtered.vcf.gz.tbi
* chrXX.hardfiltered.sites.vcf
* chrXX.hardfiltered.sites.vcf.log
* chrXX.hardfiltered.sites.vcf.summary
* chrXX.hardfiltered.vcf.gz
* chrXX.hardfiltered.vcf.gz.OK
* chrXX.hardfiltered.vcf.gz.tbi
* chrXX.merged.sites.vcf
* chrXX.merged.stats.vcf
* chrXX.merged.vcf
* chrXX.merged.vcf.OK

The .merged.vcf is the merged together versions of the separate regions in the same chromosome.

The filtered is the merged.vcf after it has been run through filters and is marked with PASS/FAIL.

Under the split/chrXX directory, there should be:
* chrXX.filtered.PASS.split.[N].vcf.gz
* chrXX.filtered.PASS.split.err
* chrXX.filtered.PASS.split.vcflist
* chrXX.filtered.PASS.gz
* subset.OK

Tutorial: GotCloud

2013-06-26T13:52:11Z

Ben Lerch: /* STEP 2 : Run GotCloud Alignment Pipeline */

= GotCloud Tutorial =
In this tutorial, we illustrate some of the essential steps in the analysis of next generation sequence data.

For a background on GotCloud and Sequence Analysis Pipelines, see [[GotCloud]]

While GotCloud can run on a cluster of machines or instances, this tutorial is just a small test that just runs on the machine the commands are run on.

GotCloud and this basic tutorial were presented at the [http://ibg.colorado.edu/dokuwiki/doku.php?id=workshop:2013:announcement 2013 IBG Workshop]. It was presented in two sessions. On Wednesday an overview was presented with steps for running the tutorial data: [[Media:IBG2013GotCloud.pdf|IBG2013GotCloud.pdf]]. On Friday more detail on the input files and what goes into generating the input files was presented: [[Media:GotCloudIBGWorkshop2013Friday.pdf|GotCloudIBGWorkshop2013Friday.pdf]].

'''This tutorial is in the process of being updated for gotcloud version 1.06 (April 17. 2013).'''

== STEP 1 : Setup GotCloud ==

[[GotCloud]] has been developed and tested on Linux Ubuntu 12.10 and 12.04.2 LTS but has not been tested on other Linux operating systems. It is not available for Windows. If you do not have your own set of machines to run on, GotCloud is also available for Ubuntu running on the Amazon Elastic Compute Cloud, see [[Amazon_Snapshot]] for more information.

We will use 3 different directories for this tutorial:
# path to the directory where gotcloud is installed, default ~/gotcloud/
# path to the directory where the example data is installed, default ~/gotcloudExample
# path to your output directory, default ~/gotcloudTutorialOut/

If the directories specified above do not reflect the directories you would like to use, replace their occurrances in the instructions below with the appropriate paths.

=== Step 1a: Install GotCloud ===
In order to run this tutorial, you need to make sure you have GotCloud installed on your system.

If you have root and would like to install gotcloud on your system, follow: [[GotCloud#Install_GotCloud_Software| root access installation instructions]]

Otherwise, you can install it in your own directory:
# Change to the directory where you want gotcloud/ installed
# Download the gotcloud tar from the ftp site.
# Extract the tar
# Build (compile) the source
#* Note: as the source builds, many messages will scroll through your terminal. You may even see some warnings. These messages are normal and expected. As long as the build does not end with an error, you have successfully built the source.

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloud_latest.tgz # Download
tar xf gotcloud_latest.tgz # Extracts into gotcloud/
cd ~/gotcloud/src; make # Build source

GotCloud requires the following tools to be installed.
You can run ~/gotcloud/scripts/check_requirements.sh
...TBD – put in required programs/tools.
* java (java-common default-jre on ubuntu)
* make (make on ubuntu)
* libssl (libssl0.9.8 on ubuntu)
* gcc 4.4 or newer

=== Step 1b: Install Example Dataset ===
Our dataset consists of 60 individuals from Great Britain (GBR) sequenced by the 1000 Genomes Project. These individuals have been sequenced to an average depth of about 4x.

To conserve time and disk-space, our analysis will focus on a small region on chromosome 20, 42900000 - 43200000.

The tutorial will run the alignment pipeline on 2 of the individuals (HG00096, HG00100). The fastqs used for this step are reduced to reads that fall into our target region.

The tutorial will then used previously aligned/mapped reads for the full 60 individuals to generate a list of polymorphic sites and estimate accurate genotypes at each of these sites.

The example dataset we'll be using is available at: ftp://share.sph.umich.edu/gotcloud/gotcloudExample.tgz

# Change directory to where you want to install the Tutorial data
# Download the dataset tar from the ftp site
# Extract the tar

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloudExample_latest.tgz # Download
tar xvf gotcloudExample_latest.tgz # Extracts into gotcloudExample/

== STEP 2 : Run GotCloud Alignment Pipeline ==
The first step in processing next generation sequence data is mapping the reads to the reference genome, generating per sample BAM files.

The alignment pipeline has multiple built-in steps to generate BAMs:
# Align the fastqs to the reference genome
#* handles both single & paired end
# Merge the results from multiple fastqs into 1 file per sample
# Mark Duplicate Reads are marked
# Recalibrate Base Qualities

This processing results in 1 BAM file per sample.

The alignment pipeline also includes Quality Control (QC) steps:
# Visualization of various quality measures (QPLOT)
# Screen for sample contamination & swap (VerifyBamID)

Run the alignment pipeline (the example aligns 2 samples) :
~/gotcloud/gotcloud align --conf ~/gotcloudExample/[[#Alignment Configuration File|GBR2align.conf]] --outdir [[#Alignment Output Directory|~/gotcloudTutorialOut]] --baseprefix ~/gotcloudExample

Upon successful completion of the alignment pipeline (about 1-3 minutes), you will see the following message:
Processing finished in n secs with no errors reported

The final BAM files produced by the alignment pipeline are:
ls ~/gotcloudTutorialOut/bams
In this directory you will see:
* BAM (.bam) files - 1 per sample
** HG00096.recal.bam
** HG00100.recal.bam
* BAM index files (.bai) – 1 per sample
** HG00096.recal.bam.bai
** HG00100.recal.bam.bai
* Indicator files that the step completed successfully:
** HG00096.recal.bam.done
** HG00100.recal.bam.done

The Quality Control (QC) files are:
ls ~/gotcloudTutorialOut/QCFiles
In this directory you will see:
* VerifyBamID output files:
** HG00096.genoCheck.depthRG
** HG00096.genoCheck.depthSM
** HG00096.genoCheck.selfRG
** HG00096.genoCheck.selfSM
** HG00100.genoCheck.depthRG
** HG00100.genoCheck.depthSM
** HG00100.genoCheck.selfRG
** HG00100.genoCheck.selfSM

* VerifyBamID step completion files – 1 per sample
** HG00096.genoCheck.done
** HG00100.genoCheck.done

* QPLOT output files
** HG00096.qplot.R
** HG00096.qplot.stats
** HG00100.qplot.R
** HG00100.qplot.stats

* QPLOT step completion files – 1 per sample
** HG00096.qplot.done
** HG00100.qplot.done

For information on the VerifyBamID output, see: [[Understanding VerifyBamID output]]

For information on the QPLOT output, see: [[Understanding QPLOT output]]

== STEP 3 : Run GotCloud Variant Calling Pipeline ==
The next step is to analyze BAM files by calling SNPs and generating a VCF file containing the variant calls.

The variant calling pipeline has multiple built-in steps to generate BAMs:
# Filter out reads with low mapping quality
# Per Base Alignment Quality Adjustment (BAQ)
# Resolve overlapping paired end reads
# Generate genotype likelihood files
# Perform variant calling
# Extract features from variant sites
# Perform variant filtering

To speed variant calling, each chromosome is broken up into smaller regions which are processed separately. While initially split by sample, the per sample data gets merged and is processed together for each region. These regions are later merged to result in a single Variant Call File (VCF) per chromosome. For the tutorial all of the data falls within a single region.

Run the variant calling pipeline:
~/gotcloud/gotcloud snpcall --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --region 20:42900000-43200000 --baseprefix ~/gotcloudExample

Upon successful completion of the variant calling pipeline (about 3-4 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

On SNP Call success, the VCF files of interest are:
ls ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered*

This gives you the following files:
* '''chr20.filtered.vcf.gz ''' - vcf for whole chromosome after it has been run through hardfilters and SVM filters and marked with PASS/FAIL including per sample genotypes
* chr20.filtered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
* chr20.filtered.sites.vcf.norm.log - log file
* chr20.filtered.sites.vcf.summary - summary of filters applied
* chr20.filtered.vcf.gz.OK - indicator that the filtering completed successfully
* chr20.filtered.vcf.gz.tbi - index file for the vcf file

Also in the ~/gotcloudTutorialOut/vcfs/chr20 directory are intermediate files:
* the whole chromosome variant calls prior to any filtering:
** chr20.merged.sites.vcf - without per sample genotypes
** chr20.merged.stats.vcf
** chr20.merged.vcf - including per sample genotypes
** chr20.merged.vcf.OK - indicator that the step completed successfully
* the hardfiltered (pre-svm filtered) variant calls:
** chr20.filtered.vcf.gz - vcf for whole chromosome after it has been run through hard filters
** chr20.hardfiltered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
** chr20.hardfiltered.sites.vcf.log - log file
** chr20.hardfiltered.sites.vcf.summary - summary of filters applied
** chr20.hardfiltered.vcf.gz.OK - indicator that the filtering completed successfully
** chr20.hardfiltered.vcf.gz.tbi - index file for the vcf file
* 40000001.45000000 subdirectory contains the data for just that region.

The ~/gotcloudTutorialOut/split/chr20 folder contains a VCF with just the sites that pass the filters.
ls ~/gotcloudTutorialOut/split/chr20/
* '''chr20.filtered.PASS.vcf.gz ''' – vcf of just sites that pass all filters
* chr20.filtered.PASS.split.1.vcf.gz - intermediate file
* chr20.filtered.PASS.split.err - log file
* chr20.filtered.PASS.split.vcflist - list of intermediate files
* subset.OK

In addition to the vcfs subdirectory, there are additional intermediate files/directories:
* glfs – holds genotype likelihood format [[GLF]] files split by chromosome, region, and sample
* pvcfs – holds intermediate vcf files split by chromosome and region

Note: the tutorial does not produce a target directory, but if you run with targeted data, you may see that.

== STEP 4 : Run GotCloud Genotype Refinement Pipeline ==
The next step is to perform genotype refinement using linkage disequilibrium information using [http://faculty.washington.edu/browning/beagle/beagle.html Beagle] & [[ThunderVCF]].

Run the LD-aware genotype refinement pipeline:
~/gotcloud/gotcloud ldrefine --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --baseprefix ~/gotcloudExample

Upon successful completion of this pipeline (about 10 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

The output from the beagle step of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz.tbi

The output from the thunderVcf (final step) of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz.tbi

== STEP 5 : Run GotCloud Association Analysis Pipeline (EPACTS) ==

We will assume that the EPACTS are installed in the following directory
setenv EPACTS /path/to/epacts
(If you need to install EPACTS, please refer to the documentation at [[EPACTS#Installation_Details]])

$EPACTS/epacts single --vcf ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered.vcf.gz --ped ~/gotcloudExample/test.GBR60.ped \\
--out ~/gotcloudTutorialOut/epacts --test q.linear --run 1 --top 1 --chr 20

Upon successful run, you will see files starting with ~/gotcloudTutorialOut/epacts
ls ~/gotcloudTutorialOut/epacts*

To see the top associated variants, you can run
less ~/gotcloudTutorialOut/epacts.epacts.top5000

To see the locus-zoom like plot, you can type the following command (assuming GNU gnuplot 4.2 or higher version was installed)
xpdf ~/gotcloudTutorialOut/epacs.zoom.20.42987877.pdf

Click [[Media:EPACTS TEST.zoom.20.42987877.pdf | Exampe LocusZoom PDF]] to see the expected output pdf

= Frequently Asked Questions (FAQs) =

'''I ran the tutorial example successfully, how can I run it with my real sequence data?'''

Congratulations for your successful run of your [[GotCloud]] Tutorial. Please see [[#Tutorial Inputs]] section to prepare your own input files for your sequence data. You will need to specify the FASTQ files associated with its sample names as explained. In addition, you will need to download the full reference and resource file across whole genome (the Tutorial contains only chr20 portion to make it compact) See [[#Alignment Configuration File]] section for the detailed information. Also, please refer to the original documentation of [[GotCloud]] for more detailed guide on installation beyond the scope of the tutorial.

= Input Files for GotCloud Tutorial =

This section describes the input files needed for the GotCloud tutorial. You don't need to know this detail to run the tutorial, but if you're interested in understanding the structure of GotCloud pipeline and run with your own sample, this would be a good starting point

== Alignment Pipeline ==
=== List of Input Files needed for Alignment ===
The command-line inputs to the tutorial alignment pipeline are:
# [[#Alignment Configuration File|Configuration File (--conf)]]
#* Specifies the configuration file to use when running
# [[#Alignment Output Directory|Output Directory (--outdir)]]
#* Directory where the output should be placed.

Additional information required to run the alignment pipeline:
# [[#Alignment FASTQ Index File|Index file of FASTQs]]
# [[#Alignment Reference Files|Alignment Reference Files]]
For the tutorial, these values are specified in the configuration file.

=== Alignment Configuration File ===
The configuration file contains KEY = VALUE settings that override defaults and set specific values for the given run.

<pre>
INDEX_FILE = GBR60fastq.index
############
# References
REF_DIR = chr20Ref
AS = NCBI37
FA_REF = $(REF_DIR)/human_g1k_v37_chr20.fa
DBSNP_VCF = $(REF_DIR)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_DIR)/hapmap_3.3.b37.sites.chr20.vcf.gz
</pre>

This configuration file sets:
* [[#Alignment FASTQ Index File|INDEX_FILE]] - file containing the fastqs to be processed as well as the read group information for these fastqs.
* Reference Information: see [[#Alignment Reference Files|Alignment Reference Files]] for more information
** AS - assembly value to put in the BAM

The index file and chromosome 20 references used in this tutorial are included with the example data under the ~/gotcloudExample directory. The tutorial uses chromosome 20 only references in order to speed the processing time.

When running with your own data, you will need to update the:
* Index File to contain the information for your own fastq files
** See [[#Alignment FASTQ Index File|Alignment FASTQ Index File]] for more information on the contents of the index file.
* The reference files to be whole genome references
** If you are just running chromosome 20, you can use the tutorial references
** Whole genome reference files can be downloaded from [[GotCloudReference]].
** See [[#Alignment Reference Files|Alignment Reference Files]] for more information.

Note: It is recommended that you use absolute paths (full path names, like “/home/mktrost/gotcloudReference” rather than just “gotcloudReference”). This example does not use absolute paths in order to be flexible to where the data is installed, but using relative paths requires it to be run from the correct directory.

=== Reference Files ===

Reference files are required for running both the alignment and variant calling pipelines.

The configuration keys for setting these are:
* FA_REF - Genome sequence reference files (needed for both pipelines)
* DBSNP_VCF – DBSNP site VCF file (needed for both pipelines)
* HM3_VCF - HAPMAP site VCF file (needed for both pipelines)
* INDEL_PREFIX - Indel sites file (need for variant calling pipeline)

The tutorial configuration file is setup to point to the required chromosome 20 reference files which are included with the tutorial example data in ~/gotcloudExample/chr20Ref/.

If you are running more than just chromosome 20, you will need whole genome reference files which can be downloaded from [[GotCloudReference]].

The configuration settings for these files are setup in the default configuration so do not need to be specified. You just need to set REF_DIR in your configuration file to the path where you installed your reference files.

To learn more about the reference files that are required, see [[GotCloud: Reference Files]].

=== Alignment Output Directory ===
This setting tells the alignment pipeline where to write the output and intermediate files.

The output directory will be created if it doesn't already exist and will contain the following Directories/files:
* bams - directory containing the final bams/bai files
** HG00096.OK - indicates that this sample completed alignment processing
** HG00100.OK - indicates that this sample completed alignment processing
* failLogs - directory containing logs from steps that failed
** this directory is only created if an error is detected
* Makefiles - directory containing the makefiles with commands for processing each sample
** biopipe_HG00096.Makefile – commands for processing sample HG00096
** biopipe_HG00100.Makefile – commands for processing sample HG00100
** biopipe_HG00096.Makefile.log – log file from running the associated Makefile
** biopipe_HG00100.Makefile.log – log file from running the associated Makefile
* QCFiles - directory containing the QC Results
**
* tmp - directory containing temporary alignment files
** bwa.sai.t – contains temporary files for the 1st step of bwa that generates sai files from fastq files
*** fastq
**** *.done – indicator files that the step to generate the file completed
**** HG0096
** alignment.bwa – contains temporary files for the 2nd step of bwa that generates BAM files
*** fastq
**** *.done – indicator files that the step to generate the file completed
** alignment.pol - contains temporary files for the polish bam step that cleans up BAM files
** alignment.dedup – contains temporary files for the deduping step
*** *.done – indicator files that the step to generate the file completed
*** *.metrics – metrics files for the deduping step
** alignment.recal – contains temporary files for the deduping step
*** *.log – contains information about the recalibration step
*** *.qemp – contains the recalibration tables used for recalibrating each BAM file

===Alignment FASTQ Index File===

There are four fastq files in {ROOT_DIR}/test/align/fastq/Sample_1 and four fastq files in {ROOT_DIR}/test/align/fastq/Sample_2, both in paired-end format. Normally, we would need to build an index file for these files. Conveniently, an index file (indexFile.txt) already exists for the automatic test samples. It can be found in {ROOT_DIR}/test/align/, and contains the following information in tab-delimited format:

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/Sample_1/File1_R1.fastq.gz fastq/Sample_1/File1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/Sample_1/File2_R1.fastq.gz fastq/Sample_1/File2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/Sample_2/File1_R1.fastq.gz fastq/Sample_2/File1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/Sample_2/File2_R1.fastq.gz fastq/Sample_2/File2_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA

If you are in the {ROOT_DIR}/test/align directory, you can use this file as-is. If you prefer, you can create a new index file and change the MERGE_NAME, RGID, SAMPLE, LIBRARY, CENTER, or PLATFORM values. It is recommended that you do not modify existing files in {ROOT_DIR}/test/align.

If you want to run this example from a different directory, make sure the FASTQ1 and FASTQ2 paths are correct. That is, each of the FASTQ1 and FASTQ2 entry in the index file should look like the following:

{ROOT_DIR}/test/align/fastq/Sample_1/File1_R1.fastq.gz

Alternately, if you want to run this example from a different directory, but do not want to edit the index file, you can create a relative path to the test fastq files so their path agrees with that listed in the index file:

ln -s {ROOT_DIR}/test/align/fastq fastq

This will create a symbolic link to the test fastq directory from your current directory.

(More information about: [[Mapping_Pipeline#Sequence_Index_File|the index file]].)

==Analyzing a Sample==

Using UMAKE, you can analyze BAM files by calling SNPs, and generate a VCF file containing the results. Once again, we can analyze BAM files used in the automatic test. For this example, we have 60 BAM files, which can be found in {ROOT_DIR}/test/umake/bams. These contain sequence information for a targeted region in chromosome 20.

In addition to the BAM files, you will need three files to run UMAKE: an index file, a configuration file, and a bed file (needed to analyze BAM files from targeted/exome sequencing).

===Index file===

First, you need a list of all the BAM files to be analyzed. Conveniently, the a test index file (umake_test.index) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

NA12272 ALL bams/NA12272.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
NA12004 ALL bams/NA12004.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
...
NA12874 ALL bams/NA12874.mapped.LS454.ssaha2.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam

You can use this file directly if you change your current directory to {ROOT_DIR}/test/umake/.

Alternately, if you want to copy and use this index file to a different directory, you can create a symbolic link to the bams folder as follows:

ln -s {ROOT_DIR}/test/umake/bams bams

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Index_File|the index file]].)

===BED file===

This file contains a single line:

chr20 20000050 20300000

You can copy this to the current directory and use it as-is.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Targeted.2FExome_Sequencing_Settings|targeted/exome sequencing settings]].)

===Configuration file===

A configuration file (umake_test.conf) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

CHRS = 20
BAM_INDEX = GBR60bam.index
############
# References
REF_ROOT = chr20Ref
#
REF = $(REF_ROOT)/human_g1k_v37_chr20.fa
INDEL_PREFIX = $(REF_ROOT)/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_VCF = $(REF_ROOT)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_ROOT)/hapmap_3.3.b37.sites.chr20.vcf.gz

CHRS = 20
TEST_ROOT = $(UMAKE_ROOT)/test/umake
BAM_INDEX = $(TEST_ROOT)/umake_test.index
OUT_PREFIX = umake_test
REF_ROOT = $(TEST_ROOT)/ref
#
REF = $(REF_ROOT)/karma.ref/human.g1k.v37.chr20.fa
INDEL_PREFIX = $(REF_ROOT)/indels/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_PREFIX = $(REF_ROOT)/dbSNP/dbsnp_135_b37.rod
HM3_PREFIX = $(REF_ROOT)/HapMap3/hapmap3_r3_b37_fwd.consensus.qc.poly
#
RUN_INDEX = TRUE # create BAM index file
RUN_PILEUP = TRUE # create GLF file from BAM
RUN_GLFMULTIPLES = TRUE # create unfiltered SNP calls
RUN_VCFPILEUP = TRUE # create PVCF files using vcfPileup and run infoCollector
RUN_FILTER = TRUE # filter SNPs using vcfCooker
RUN_SPLIT = TRUE # split SNPs into chunks for genotype refinement
RUN_BEAGLE = FALSE # BEAGLE - MUST SET AFTER FINISHING PREVIOUS STEPS
RUN_SUBSET = FALSE # SUBSET FOR THUNDER - MAY BE SET WITH BEAGLE STEP TOGETHER
RUN_THUNDER = FALSE # THUNDER - MUST SET AFTER FINISHING PREVIOUS STEPS
###############################################################################
WRITE_TARGET_LOCI = TRUE # FOR TARGETED SEQUENCING ONLY -- Write loci file when performing pileup
UNIFORM_TARGET_BED = $(TEST_ROOT)/umake_test.bed # Targeted sequencing : When all individuals has the same target. Otherwise, comment it out
OFFSET_OFF_TARGET = 50 # Extend target by given # of bases
MULTIPLE_TARGET_MAP = # Target per individual : Each line contains [SM_ID] [TARGET_BED]
TARGET_DIR = target # Directory to store target information
SAMTOOLS_VIEW_TARGET_ONLY = TRUE # When performing samtools view, exclude off-target regions (may make command line too long)

If you are running this from a different directory, you will want to change some of the lines as follows:

BAM_INDEX = {CURRENT_DIR}/umake_test.index
UNIFORM_TARGET_BED = {CURRENT_DIR}/umake_test.bed

where {CURRENT_DIR} is the absolute path to the directory that contains the index and bed files.

An additional option can be added in the configuration file:

OUT_DIR = {OUT_DIR}

where {OUT_DIR} is the name of directory in which you want the output to be stored. If you do not specify this in the configuration file, you will need to add an extra parameter when you run UMAKE in the next step.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Configuration_File|the configuration file]], [[Variant_Calling_Pipeline_(UMAKE)#Reference_Files|reference files]].)

==Further Information==

[[Mapping_Pipeline|Mapping (Alignment) Pipeline]]

[[Variant_Calling_Pipeline_(UMAKE)|Variant Calling Pipeline (UMAKE)]]

Tutorial: GotCloud

2013-06-26T12:09:34Z

Ben Lerch: /* STEP 2 : Run GotCloud Alignment Pipeline */

= GotCloud Tutorial =
In this tutorial, we illustrate some of the essential steps in the analysis of next generation sequence data.

For a background on GotCloud and Sequence Analysis Pipelines, see [[GotCloud]]

While GotCloud can run on a cluster of machines or instances, this tutorial is just a small test that just runs on the machine the commands are run on.

GotCloud and this basic tutorial were presented at the [http://ibg.colorado.edu/dokuwiki/doku.php?id=workshop:2013:announcement 2013 IBG Workshop]. It was presented in two sessions. On Wednesday an overview was presented with steps for running the tutorial data: [[Media:IBG2013GotCloud.pdf|IBG2013GotCloud.pdf]]. On Friday more detail on the input files and what goes into generating the input files was presented: [[Media:GotCloudIBGWorkshop2013Friday.pdf|GotCloudIBGWorkshop2013Friday.pdf]].

'''This tutorial is in the process of being updated for gotcloud version 1.06 (April 17. 2013).'''

== STEP 1 : Setup GotCloud ==

[[GotCloud]] has been developed and tested on Linux Ubuntu 12.10 and 12.04.2 LTS but has not been tested on other Linux operating systems. It is not available for Windows. If you do not have your own set of machines to run on, GotCloud is also available for Ubuntu running on the Amazon Elastic Compute Cloud, see [[Amazon_Snapshot]] for more information.

We will use 3 different directories for this tutorial:
# path to the directory where gotcloud is installed, default ~/gotcloud/
# path to the directory where the example data is installed, default ~/gotcloudExample
# path to your output directory, default ~/gotcloudTutorialOut/

If the directories specified above do not reflect the directories you would like to use, replace their occurrances in the instructions below with the appropriate paths.

=== Step 1a: Install GotCloud ===
In order to run this tutorial, you need to make sure you have GotCloud installed on your system.

If you have root and would like to install gotcloud on your system, follow: [[GotCloud#Install_GotCloud_Software| root access installation instructions]]

Otherwise, you can install it in your own directory:
# Change to the directory where you want gotcloud/ installed
# Download the gotcloud tar from the ftp site.
# Extract the tar
# Build (compile) the source
#* Note: as the source builds, many messages will scroll through your terminal. You may even see some warnings. These messages are normal and expected. As long as the build does not end with an error, you have successfully built the source.

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloud_latest.tgz # Download
tar xf gotcloud_latest.tgz # Extracts into gotcloud/
cd ~/gotcloud/src; make # Build source

GotCloud requires the following tools to be installed.
You can run ~/gotcloud/scripts/check_requirements.sh
...TBD – put in required programs/tools.
* java (java-common default-jre on ubuntu)
* make (make on ubuntu)
* libssl (libssl0.9.8 on ubuntu)
* gcc 4.4 or newer

=== Step 1b: Install Example Dataset ===
Our dataset consists of 60 individuals from Great Britain (GBR) sequenced by the 1000 Genomes Project. These individuals have been sequenced to an average depth of about 4x.

To conserve time and disk-space, our analysis will focus on a small region on chromosome 20, 42900000 - 43200000.

The tutorial will run the alignment pipeline on 2 of the individuals (HG00096, HG00100). The fastqs used for this step are reduced to reads that fall into our target region.

The tutorial will then used previously aligned/mapped reads for the full 60 individuals to generate a list of polymorphic sites and estimate accurate genotypes at each of these sites.

The example dataset we'll be using is available at: ftp://share.sph.umich.edu/gotcloud/gotcloudExample.tgz

# Change directory to where you want to install the Tutorial data
# Download the dataset tar from the ftp site
# Extract the tar

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloudExample_latest.tgz # Download
tar xvf gotcloudExample_latest.tgz # Extracts into gotcloudExample/

== STEP 2 : Run GotCloud Alignment Pipeline ==
The first step in processing next generation sequence data is mapping the reads to the reference genome, generating per sample BAM files.

The alignment pipeline has multiple built-in steps to generate BAMs:
# Align the fastqs to the reference genome
#* handles both single & paired end
# Merge the results from multiple fastqs into 1 file per sample
# Mark Duplicate Reads are marked
# Recalibrate Base Qualities

This processing results in 1 BAM file per sample.

The alignment pipeline also includes Quality Control (QC) steps:
# Visualization of various quality measures (QPLOT)
# Screen for sample contamination & swap (VerifyBamID)

Run the alignment pipeline (the example aligns 2 samples) :
~/gotcloud/gotcloud align --conf ~/gotcloudExample/[[#Alignment Configuration File|GBR2align.conf]] --outdir [[#Alignment Output Directory|~/gotcloudTutorialOut]] --baseprefix ~/gotcloudExample

Upon successful completion of the alignment pipeline (about 1-3 minutes), you will see the following message:
Processing finished in n secs with no errors reported

The final BAM files produced by the alignment pipeline are:
ls ~/gotcloudTutorialOut/bams
In this directory you will see:
* BAM (.bam) files - 1 per sample
** HG00096.recal.bam
** HG00100.recal.bam
* BAM index files (.bai) – 1 per sample
** HG00096.recal.bam.bai
** HG00100.recal.bam.bai
* BAM checksum files (.md5) – 1 per sample
** HG00096.recal.bam.md5
** HG00100.recal.bam.md5
* Indicator files that the step completed successfully:
** HG00096.recal.bam.done
** HG00100.recal.bam.done

The Quality Control (QC) files are:
ls ~/gotcloudTutorialOut/QCFiles
In this directory you will see:
* VerifyBamID output files:
** HG00096.genoCheck.depthRG
** HG00096.genoCheck.depthSM
** HG00096.genoCheck.selfRG
** HG00096.genoCheck.selfSM
** HG00100.genoCheck.depthRG
** HG00100.genoCheck.depthSM
** HG00100.genoCheck.selfRG
** HG00100.genoCheck.selfSM

* VerifyBamID step completion files – 1 per sample
** HG00096.genoCheck.done
** HG00100.genoCheck.done

* QPLOT output files
** HG00096.qplot.R
** HG00096.qplot.stats
** HG00100.qplot.R
** HG00100.qplot.stats

* QPLOT step completion files – 1 per sample
** HG00096.qplot.done
** HG00100.qplot.done

For information on the VerifyBamID output, see: [[Understanding VerifyBamID output]]

For information on the QPLOT output, see: [[Understanding QPLOT output]]

== STEP 3 : Run GotCloud Variant Calling Pipeline ==
The next step is to analyze BAM files by calling SNPs and generating a VCF file containing the variant calls.

The variant calling pipeline has multiple built-in steps to generate BAMs:
# Filter out reads with low mapping quality
# Per Base Alignment Quality Adjustment (BAQ)
# Resolve overlapping paired end reads
# Generate genotype likelihood files
# Perform variant calling
# Extract features from variant sites
# Perform variant filtering

To speed variant calling, each chromosome is broken up into smaller regions which are processed separately. While initially split by sample, the per sample data gets merged and is processed together for each region. These regions are later merged to result in a single Variant Call File (VCF) per chromosome. For the tutorial all of the data falls within a single region.

Run the variant calling pipeline:
~/gotcloud/gotcloud snpcall --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --region 20:42900000-43200000 --baseprefix ~/gotcloudExample

Upon successful completion of the variant calling pipeline (about 3-4 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

On SNP Call success, the VCF files of interest are:
ls ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered*

This gives you the following files:
* '''chr20.filtered.vcf.gz ''' - vcf for whole chromosome after it has been run through hardfilters and SVM filters and marked with PASS/FAIL including per sample genotypes
* chr20.filtered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
* chr20.filtered.sites.vcf.norm.log - log file
* chr20.filtered.sites.vcf.summary - summary of filters applied
* chr20.filtered.vcf.gz.OK - indicator that the filtering completed successfully
* chr20.filtered.vcf.gz.tbi - index file for the vcf file

Also in the ~/gotcloudTutorialOut/vcfs/chr20 directory are intermediate files:
* the whole chromosome variant calls prior to any filtering:
** chr20.merged.sites.vcf - without per sample genotypes
** chr20.merged.stats.vcf
** chr20.merged.vcf - including per sample genotypes
** chr20.merged.vcf.OK - indicator that the step completed successfully
* the hardfiltered (pre-svm filtered) variant calls:
** chr20.filtered.vcf.gz - vcf for whole chromosome after it has been run through hard filters
** chr20.hardfiltered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
** chr20.hardfiltered.sites.vcf.log - log file
** chr20.hardfiltered.sites.vcf.summary - summary of filters applied
** chr20.hardfiltered.vcf.gz.OK - indicator that the filtering completed successfully
** chr20.hardfiltered.vcf.gz.tbi - index file for the vcf file
* 40000001.45000000 subdirectory contains the data for just that region.

The ~/gotcloudTutorialOut/split/chr20 folder contains a VCF with just the sites that pass the filters.
ls ~/gotcloudTutorialOut/split/chr20/
* '''chr20.filtered.PASS.vcf.gz ''' – vcf of just sites that pass all filters
* chr20.filtered.PASS.split.1.vcf.gz - intermediate file
* chr20.filtered.PASS.split.err - log file
* chr20.filtered.PASS.split.vcflist - list of intermediate files
* subset.OK

In addition to the vcfs subdirectory, there are additional intermediate files/directories:
* glfs – holds genotype likelihood format [[GLF]] files split by chromosome, region, and sample
* pvcfs – holds intermediate vcf files split by chromosome and region

Note: the tutorial does not produce a target directory, but if you run with targeted data, you may see that.

== STEP 4 : Run GotCloud Genotype Refinement Pipeline ==
The next step is to perform genotype refinement using linkage disequilibrium information using [http://faculty.washington.edu/browning/beagle/beagle.html Beagle] & [[ThunderVCF]].

Run the LD-aware genotype refinement pipeline:
~/gotcloud/gotcloud ldrefine --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --baseprefix ~/gotcloudExample

Upon successful completion of this pipeline (about 10 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

The output from the beagle step of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz.tbi

The output from the thunderVcf (final step) of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz.tbi

== STEP 5 : Run GotCloud Association Analysis Pipeline (EPACTS) ==

We will assume that the EPACTS are installed in the following directory
setenv EPACTS /path/to/epacts
(If you need to install EPACTS, please refer to the documentation at [[EPACTS#Installation_Details]])

$EPACTS/epacts single --vcf ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered.vcf.gz --ped ~/gotcloudExample/test.GBR60.ped \\
--out ~/gotcloudTutorialOut/epacts --test q.linear --run 1 --top 1 --chr 20

Upon successful run, you will see files starting with ~/gotcloudTutorialOut/epacts
ls ~/gotcloudTutorialOut/epacts*

To see the top associated variants, you can run
less ~/gotcloudTutorialOut/epacts.epacts.top5000

To see the locus-zoom like plot, you can type the following command (assuming GNU gnuplot 4.2 or higher version was installed)
xpdf ~/gotcloudTutorialOut/epacs.zoom.20.42987877.pdf

Click [[Media:EPACTS TEST.zoom.20.42987877.pdf | Exampe LocusZoom PDF]] to see the expected output pdf

= Frequently Asked Questions (FAQs) =

'''I ran the tutorial example successfully, how can I run it with my real sequence data?'''

Congratulations for your successful run of your [[GotCloud]] Tutorial. Please see [[#Tutorial Inputs]] section to prepare your own input files for your sequence data. You will need to specify the FASTQ files associated with its sample names as explained. In addition, you will need to download the full reference and resource file across whole genome (the Tutorial contains only chr20 portion to make it compact) See [[#Alignment Configuration File]] section for the detailed information. Also, please refer to the original documentation of [[GotCloud]] for more detailed guide on installation beyond the scope of the tutorial.

= Input Files for GotCloud Tutorial =

This section describes the input files needed for the GotCloud tutorial. You don't need to know this detail to run the tutorial, but if you're interested in understanding the structure of GotCloud pipeline and run with your own sample, this would be a good starting point

== Alignment Pipeline ==
=== List of Input Files needed for Alignment ===
The command-line inputs to the tutorial alignment pipeline are:
# [[#Alignment Configuration File|Configuration File (--conf)]]
#* Specifies the configuration file to use when running
# [[#Alignment Output Directory|Output Directory (--outdir)]]
#* Directory where the output should be placed.

Additional information required to run the alignment pipeline:
# [[#Alignment FASTQ Index File|Index file of FASTQs]]
# [[#Alignment Reference Files|Alignment Reference Files]]
For the tutorial, these values are specified in the configuration file.

=== Alignment Configuration File ===
The configuration file contains KEY = VALUE settings that override defaults and set specific values for the given run.

<pre>
INDEX_FILE = GBR60fastq.index
############
# References
REF_DIR = chr20Ref
AS = NCBI37
FA_REF = $(REF_DIR)/human_g1k_v37_chr20.fa
DBSNP_VCF = $(REF_DIR)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_DIR)/hapmap_3.3.b37.sites.chr20.vcf.gz
</pre>

This configuration file sets:
* [[#Alignment FASTQ Index File|INDEX_FILE]] - file containing the fastqs to be processed as well as the read group information for these fastqs.
* Reference Information: see [[#Alignment Reference Files|Alignment Reference Files]] for more information
** AS - assembly value to put in the BAM

The index file and chromosome 20 references used in this tutorial are included with the example data under the ~/gotcloudExample directory. The tutorial uses chromosome 20 only references in order to speed the processing time.

When running with your own data, you will need to update the:
* Index File to contain the information for your own fastq files
** See [[#Alignment FASTQ Index File|Alignment FASTQ Index File]] for more information on the contents of the index file.
* The reference files to be whole genome references
** If you are just running chromosome 20, you can use the tutorial references
** Whole genome reference files can be downloaded from [[GotCloudReference]].
** See [[#Alignment Reference Files|Alignment Reference Files]] for more information.

Note: It is recommended that you use absolute paths (full path names, like “/home/mktrost/gotcloudReference” rather than just “gotcloudReference”). This example does not use absolute paths in order to be flexible to where the data is installed, but using relative paths requires it to be run from the correct directory.

=== Reference Files ===

Reference files are required for running both the alignment and variant calling pipelines.

The configuration keys for setting these are:
* FA_REF - Genome sequence reference files (needed for both pipelines)
* DBSNP_VCF – DBSNP site VCF file (needed for both pipelines)
* HM3_VCF - HAPMAP site VCF file (needed for both pipelines)
* INDEL_PREFIX - Indel sites file (need for variant calling pipeline)

The tutorial configuration file is setup to point to the required chromosome 20 reference files which are included with the tutorial example data in ~/gotcloudExample/chr20Ref/.

If you are running more than just chromosome 20, you will need whole genome reference files which can be downloaded from [[GotCloudReference]].

The configuration settings for these files are setup in the default configuration so do not need to be specified. You just need to set REF_DIR in your configuration file to the path where you installed your reference files.

To learn more about the reference files that are required, see [[GotCloud: Reference Files]].

=== Alignment Output Directory ===
This setting tells the alignment pipeline where to write the output and intermediate files.

The output directory will be created if it doesn't already exist and will contain the following Directories/files:
* bams - directory containing the final bams/bai files
** HG00096.OK - indicates that this sample completed alignment processing
** HG00100.OK - indicates that this sample completed alignment processing
* failLogs - directory containing logs from steps that failed
** this directory is only created if an error is detected
* Makefiles - directory containing the makefiles with commands for processing each sample
** biopipe_HG00096.Makefile – commands for processing sample HG00096
** biopipe_HG00100.Makefile – commands for processing sample HG00100
** biopipe_HG00096.Makefile.log – log file from running the associated Makefile
** biopipe_HG00100.Makefile.log – log file from running the associated Makefile
* QCFiles - directory containing the QC Results
**
* tmp - directory containing temporary alignment files
** bwa.sai.t – contains temporary files for the 1st step of bwa that generates sai files from fastq files
*** fastq
**** *.done – indicator files that the step to generate the file completed
**** HG0096
** alignment.bwa – contains temporary files for the 2nd step of bwa that generates BAM files
*** fastq
**** *.done – indicator files that the step to generate the file completed
** alignment.pol - contains temporary files for the polish bam step that cleans up BAM files
** alignment.dedup – contains temporary files for the deduping step
*** *.done – indicator files that the step to generate the file completed
*** *.metrics – metrics files for the deduping step
** alignment.recal – contains temporary files for the deduping step
*** *.log – contains information about the recalibration step
*** *.qemp – contains the recalibration tables used for recalibrating each BAM file

===Alignment FASTQ Index File===

There are four fastq files in {ROOT_DIR}/test/align/fastq/Sample_1 and four fastq files in {ROOT_DIR}/test/align/fastq/Sample_2, both in paired-end format. Normally, we would need to build an index file for these files. Conveniently, an index file (indexFile.txt) already exists for the automatic test samples. It can be found in {ROOT_DIR}/test/align/, and contains the following information in tab-delimited format:

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/Sample_1/File1_R1.fastq.gz fastq/Sample_1/File1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/Sample_1/File2_R1.fastq.gz fastq/Sample_1/File2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/Sample_2/File1_R1.fastq.gz fastq/Sample_2/File1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/Sample_2/File2_R1.fastq.gz fastq/Sample_2/File2_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA

If you are in the {ROOT_DIR}/test/align directory, you can use this file as-is. If you prefer, you can create a new index file and change the MERGE_NAME, RGID, SAMPLE, LIBRARY, CENTER, or PLATFORM values. It is recommended that you do not modify existing files in {ROOT_DIR}/test/align.

If you want to run this example from a different directory, make sure the FASTQ1 and FASTQ2 paths are correct. That is, each of the FASTQ1 and FASTQ2 entry in the index file should look like the following:

{ROOT_DIR}/test/align/fastq/Sample_1/File1_R1.fastq.gz

Alternately, if you want to run this example from a different directory, but do not want to edit the index file, you can create a relative path to the test fastq files so their path agrees with that listed in the index file:

ln -s {ROOT_DIR}/test/align/fastq fastq

This will create a symbolic link to the test fastq directory from your current directory.

(More information about: [[Mapping_Pipeline#Sequence_Index_File|the index file]].)

==Analyzing a Sample==

Using UMAKE, you can analyze BAM files by calling SNPs, and generate a VCF file containing the results. Once again, we can analyze BAM files used in the automatic test. For this example, we have 60 BAM files, which can be found in {ROOT_DIR}/test/umake/bams. These contain sequence information for a targeted region in chromosome 20.

In addition to the BAM files, you will need three files to run UMAKE: an index file, a configuration file, and a bed file (needed to analyze BAM files from targeted/exome sequencing).

===Index file===

First, you need a list of all the BAM files to be analyzed. Conveniently, the a test index file (umake_test.index) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

NA12272 ALL bams/NA12272.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
NA12004 ALL bams/NA12004.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
...
NA12874 ALL bams/NA12874.mapped.LS454.ssaha2.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam

You can use this file directly if you change your current directory to {ROOT_DIR}/test/umake/.

Alternately, if you want to copy and use this index file to a different directory, you can create a symbolic link to the bams folder as follows:

ln -s {ROOT_DIR}/test/umake/bams bams

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Index_File|the index file]].)

===BED file===

This file contains a single line:

chr20 20000050 20300000

You can copy this to the current directory and use it as-is.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Targeted.2FExome_Sequencing_Settings|targeted/exome sequencing settings]].)

===Configuration file===

A configuration file (umake_test.conf) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

CHRS = 20
BAM_INDEX = GBR60bam.index
############
# References
REF_ROOT = chr20Ref
#
REF = $(REF_ROOT)/human_g1k_v37_chr20.fa
INDEL_PREFIX = $(REF_ROOT)/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_VCF = $(REF_ROOT)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_ROOT)/hapmap_3.3.b37.sites.chr20.vcf.gz

CHRS = 20
TEST_ROOT = $(UMAKE_ROOT)/test/umake
BAM_INDEX = $(TEST_ROOT)/umake_test.index
OUT_PREFIX = umake_test
REF_ROOT = $(TEST_ROOT)/ref
#
REF = $(REF_ROOT)/karma.ref/human.g1k.v37.chr20.fa
INDEL_PREFIX = $(REF_ROOT)/indels/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_PREFIX = $(REF_ROOT)/dbSNP/dbsnp_135_b37.rod
HM3_PREFIX = $(REF_ROOT)/HapMap3/hapmap3_r3_b37_fwd.consensus.qc.poly
#
RUN_INDEX = TRUE # create BAM index file
RUN_PILEUP = TRUE # create GLF file from BAM
RUN_GLFMULTIPLES = TRUE # create unfiltered SNP calls
RUN_VCFPILEUP = TRUE # create PVCF files using vcfPileup and run infoCollector
RUN_FILTER = TRUE # filter SNPs using vcfCooker
RUN_SPLIT = TRUE # split SNPs into chunks for genotype refinement
RUN_BEAGLE = FALSE # BEAGLE - MUST SET AFTER FINISHING PREVIOUS STEPS
RUN_SUBSET = FALSE # SUBSET FOR THUNDER - MAY BE SET WITH BEAGLE STEP TOGETHER
RUN_THUNDER = FALSE # THUNDER - MUST SET AFTER FINISHING PREVIOUS STEPS
###############################################################################
WRITE_TARGET_LOCI = TRUE # FOR TARGETED SEQUENCING ONLY -- Write loci file when performing pileup
UNIFORM_TARGET_BED = $(TEST_ROOT)/umake_test.bed # Targeted sequencing : When all individuals has the same target. Otherwise, comment it out
OFFSET_OFF_TARGET = 50 # Extend target by given # of bases
MULTIPLE_TARGET_MAP = # Target per individual : Each line contains [SM_ID] [TARGET_BED]
TARGET_DIR = target # Directory to store target information
SAMTOOLS_VIEW_TARGET_ONLY = TRUE # When performing samtools view, exclude off-target regions (may make command line too long)

If you are running this from a different directory, you will want to change some of the lines as follows:

BAM_INDEX = {CURRENT_DIR}/umake_test.index
UNIFORM_TARGET_BED = {CURRENT_DIR}/umake_test.bed

where {CURRENT_DIR} is the absolute path to the directory that contains the index and bed files.

An additional option can be added in the configuration file:

OUT_DIR = {OUT_DIR}

where {OUT_DIR} is the name of directory in which you want the output to be stored. If you do not specify this in the configuration file, you will need to add an extra parameter when you run UMAKE in the next step.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Configuration_File|the configuration file]], [[Variant_Calling_Pipeline_(UMAKE)#Reference_Files|reference files]].)

==Further Information==

[[Mapping_Pipeline|Mapping (Alignment) Pipeline]]

[[Variant_Calling_Pipeline_(UMAKE)|Variant Calling Pipeline (UMAKE)]]

GotCloud: Alignment Pipeline

2013-06-21T20:37:09Z

Ben Lerch: /* Optional Configurable Settings */

= Alignment Pipeline =

Back to parent: [[GotCloud]]

The Alignment/Mapping Pipeline takes FASTQ files and generates recalibrated BAM files from them.

== Running the GotCloud Alignment Pipeline ==

The alignment pipeline is run using the <code>align</code> option of the <code>gotcloud</code> script. This option calls <code>align.pl</code> found in the <code>bin/</code> directory under the <code>gotcloud</code> installation.

===Running the Automated Test===

The automated test runs the alignment pipeline on a small set of test data and checks that the results against expected results validating that GotCloud is installed correctly.

*Run alignment pipeline test:
gotcloud align --test OUTPUT_DIR
where OUTPUT_DIR is the directory where you want to store the test results

If you see "Successfully ran the test case, congratulations!", then you are ready to align samples.

== Overview of Alignment Pipeline Steps ==
Here is an overview of the Alignment Pipeline:

[[File:MappingSteps.png]]

== Input Data:==
*Raw Sequence (FASTQ) files
*Sequence Index file containing fastqs & RG info
*Reference files
*(Optional) Configuration file to override default options

=== Raw Sequence (FASTQ) files ===

These are the FASTQ files that need to be mapped to BAM files.

These files are specified in the [[#Sequence Index File|Sequence Index File]].

=== Sequence Index File ===
This file specifies the FASTQ files that need to be processed and the Read Group information for them.

This file is specified either via the command line parameter <code>--index_file</code> or via the configuration file setting <code>INDEX_FILE</code>.

The command-line setting takes precedence over the configuration file setting.

The Sequence Index is a tab delimited file that starts with a header line. The columns may be in any order.

Following the header line, there is one line per single-end read and one line per paired-end read (only 1 line per pair).

Required Column Names:
* MERGE_NAME - base name for the resulting BAM file for the sample (used to group multiple fastqs or fastq pairs into a single BAM)
* FASTQ1 - name of the fastq or the first in the pair if paired-end. (Only 1 line per pair)

Optional Column Names:
* FASTQ2 - name of the 2nd fastq in paired-end reads. Specify '.' if the column exists, but this line is single-ended.
* RGID - Read Group ID for this entry
* SAMPLE - Sample Name for this entry
* LIBRARY - Library for this entry
* CENTER - Center Name for this entry
* PLATFORM - Platform for this entry

The RGID, SAMPLE, LIBRARY, CENTER, and PLATFORM are used to populate the Read Group information for this entry. These fields are optional. Either leave the column header out of the file or specify '.' if the column header exists, but the data is N/A. As long as the RGID field is specified non-N/A fields are added to the BAM file.

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/S1/F1_R1.fastq.gz fastq/S1/F1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/S1/F2_R1.fastq.gz fastq/S1/F2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/S2/F1_R1.fastq.gz fastq/S2/F1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/S2/F2.fastq.gz . RGID2 SampleID2 Lib2 UM ILLUMINA

The <code>--fastq</code>/<code>FASTQ</code> setting can be used to specify a prefix to the FASTQ1/FASTQ2 file paths that should be applied before using the files.

=== Reference Files ===

The following Reference Files are required:
* Reference File fasta files
** Files required: .fa, -bs.umfa, .GCContent, .amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa
*** If you don't have the -bs.umfa file, the software will try to create it in the same directory as the reference fasta.
*** .GCContent can be generated using qplot, see: [[QPLOT#Input_files| QPLOT: Input Files: --gccontent]] and name the resulting file as <code>.fa.GCcontent</code>
*** Use <code>bin/bwa index ref.fa</code> if you need to generate the bwa reference files (.amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa)
** Configuration Name: FA_REF - specify the ref.fa/ref.fa.gz name
* DBSNP File
** tab delimited file/VCF, can be compressed
*** 1st column -> chromosome
*** 2nd column -> 1-based position
** Configuration Name: DBSNP_VCF
* PLINK-compatible binary genotype files
** Files required: .bed, .bin, .fam
** Configuration Name: PLINK

=== Configuration File ===
Configuration file contains the run-time options including the software binaries and command line arguments. A default configuration file is automatically loaded. Users may specify their own configuration file specifying just the values different than the defaults. The configuration file is not required if there are no values to override.

Comments begin with a <code>#</code>

Format: KEY = value

Where KEY is the item being set and value is its new value

See [[#Command-Line Options|Command-Line Options]] for values that can be set either via command line or via configuration.

Note: Command-line options take priority over configuration file settings

==== Required Settings ====
See [[#Reference Files|Reference Files]] for the required reference file settings.

See [[#Sequence Index File|Sequence Index File]] for how to set the index file either via command line options or via configuration.

==== Turning Off Optional Steps====
Quality Control steps can be disabled.

To Disable QPLOT, set:
RUN_QPLOT = 0

To Disable VerifyBamID, set:
RUN_VERIFY_BAM_ID = 0

==== Optional Configurable Settings ====
You may want to adjust the amount of memory/threads that are used:

There are additional configurable settings, but these are the ones most likely to be adjusted.

* BWA_THREADS = -t N
** Fill in the N with the number of threads you want BWA to run with, default is 1
* BWA_MAX_MEM = 2000000000
** Maximum amount of memory used by samtools sort after running bwa
* JAVA_MEM = -Xmx4g
** Set the maximum size of the java memory allocation pool. Default is 4g, adjust that as necessary.
*BATCH_TYPE = mosix
** Tells the cluster gateway to use mosix to send jobs to the client nodes.
*BATCH_OPTS = -j36,37,38,39,40,41,45,46,47,48,49
** Specifies which client nodes mosix should send jobs to.

== Running the Alignment Pipeline ==

=== Command-Line Options ===
* help - print usage
* test OUTPUT_DIR - run the test example placing the output in a user specified OUTPUT_DIR. No other options are required.
* out_dir OUTPUT_DIR - directory for the output
** May also be specified via OUT_DIR in the configuration file
** Required to be set either via command-line or configuration
* conf CONFIG_FILE - configuration file
* index_file INDEX_FILE_NAME - name of the index file
** May also be specified via INDEX_FILE in the configuration file
** Required to be set either via command-line or configuration
* ref_dir REFERENCE_DIR - value to set config key REF_DIR to, overriding other values, REF_DIR can then be used inside config files.
** May also be specified via REF_DIR in the configuration file
* fastq FASTQ_PATH - prefix path to the fastq files specified in the INDEX_FILE
** May also be specified via FASTQ in the configuration file
* keepTmp - Do not remove the temporary files (removed by default)
** May also be specified via KEEP_TMP in the configuration file
* numcs N - Replace N with the number of samples that should be processed in parallel
* numjobs N - Replace N with the number of targets in each makefile that should be run in parallel

Note: Command-line options take priority over configuration file settings

===Running the Alignment Pipeline===
Run <code>gotcloud align</code> with the appropriate command-line parameters.

Example:
gotcloud align --conf config.txt --outdir output

This step generates 1 Makefile per sample in the output/Makefiles/ directory and then automatically runs them. The Makefiles contain all of the information to run each sample.

If you only want to generate the makefiles and not run them, use the <code>--dryrun</code> option. It will generate the Makefiles and print instructions for running the Makefiles.

Each Makefile is independent and can be run in parallel and across a cloud.

On success, you will see:
Processing finished in nn secs with no errors reported
and should see the following subdirectories under the user specified output directory:
* bams/
* Makefiles/
* QCFiles/ (if all quality control is not disabled)
* tmp/

You should see a <code>.OK</code> for each Sample in the index file.

If you do not see these <code>.OK</code> files, then your Alignment Pipeline failed.

On success, the bams/ directory contains the final BAMs and bais.

If processing fails part way through, you can pick up where you left off by rerunning gotcloud or the make command.

GotCloud: Alignment Pipeline

2013-06-21T20:31:11Z

Ben Lerch: /* Command-Line Options */

= Alignment Pipeline =

Back to parent: [[GotCloud]]

The Alignment/Mapping Pipeline takes FASTQ files and generates recalibrated BAM files from them.

== Running the GotCloud Alignment Pipeline ==

The alignment pipeline is run using the <code>align</code> option of the <code>gotcloud</code> script. This option calls <code>align.pl</code> found in the <code>bin/</code> directory under the <code>gotcloud</code> installation.

===Running the Automated Test===

The automated test runs the alignment pipeline on a small set of test data and checks that the results against expected results validating that GotCloud is installed correctly.

*Run alignment pipeline test:
gotcloud align --test OUTPUT_DIR
where OUTPUT_DIR is the directory where you want to store the test results

If you see "Successfully ran the test case, congratulations!", then you are ready to align samples.

== Overview of Alignment Pipeline Steps ==
Here is an overview of the Alignment Pipeline:

[[File:MappingSteps.png]]

== Input Data:==
*Raw Sequence (FASTQ) files
*Sequence Index file containing fastqs & RG info
*Reference files
*(Optional) Configuration file to override default options

=== Raw Sequence (FASTQ) files ===

These are the FASTQ files that need to be mapped to BAM files.

These files are specified in the [[#Sequence Index File|Sequence Index File]].

=== Sequence Index File ===
This file specifies the FASTQ files that need to be processed and the Read Group information for them.

This file is specified either via the command line parameter <code>--index_file</code> or via the configuration file setting <code>INDEX_FILE</code>.

The command-line setting takes precedence over the configuration file setting.

The Sequence Index is a tab delimited file that starts with a header line. The columns may be in any order.

Following the header line, there is one line per single-end read and one line per paired-end read (only 1 line per pair).

Required Column Names:
* MERGE_NAME - base name for the resulting BAM file for the sample (used to group multiple fastqs or fastq pairs into a single BAM)
* FASTQ1 - name of the fastq or the first in the pair if paired-end. (Only 1 line per pair)

Optional Column Names:
* FASTQ2 - name of the 2nd fastq in paired-end reads. Specify '.' if the column exists, but this line is single-ended.
* RGID - Read Group ID for this entry
* SAMPLE - Sample Name for this entry
* LIBRARY - Library for this entry
* CENTER - Center Name for this entry
* PLATFORM - Platform for this entry

The RGID, SAMPLE, LIBRARY, CENTER, and PLATFORM are used to populate the Read Group information for this entry. These fields are optional. Either leave the column header out of the file or specify '.' if the column header exists, but the data is N/A. As long as the RGID field is specified non-N/A fields are added to the BAM file.

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/S1/F1_R1.fastq.gz fastq/S1/F1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/S1/F2_R1.fastq.gz fastq/S1/F2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/S2/F1_R1.fastq.gz fastq/S2/F1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/S2/F2.fastq.gz . RGID2 SampleID2 Lib2 UM ILLUMINA

The <code>--fastq</code>/<code>FASTQ</code> setting can be used to specify a prefix to the FASTQ1/FASTQ2 file paths that should be applied before using the files.

=== Reference Files ===

The following Reference Files are required:
* Reference File fasta files
** Files required: .fa, -bs.umfa, .GCContent, .amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa
*** If you don't have the -bs.umfa file, the software will try to create it in the same directory as the reference fasta.
*** .GCContent can be generated using qplot, see: [[QPLOT#Input_files| QPLOT: Input Files: --gccontent]] and name the resulting file as <code>.fa.GCcontent</code>
*** Use <code>bin/bwa index ref.fa</code> if you need to generate the bwa reference files (.amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa)
** Configuration Name: FA_REF - specify the ref.fa/ref.fa.gz name
* DBSNP File
** tab delimited file/VCF, can be compressed
*** 1st column -> chromosome
*** 2nd column -> 1-based position
** Configuration Name: DBSNP_VCF
* PLINK-compatible binary genotype files
** Files required: .bed, .bin, .fam
** Configuration Name: PLINK

=== Configuration File ===
Configuration file contains the run-time options including the software binaries and command line arguments. A default configuration file is automatically loaded. Users may specify their own configuration file specifying just the values different than the defaults. The configuration file is not required if there are no values to override.

Comments begin with a <code>#</code>

Format: KEY = value

Where KEY is the item being set and value is its new value

See [[#Command-Line Options|Command-Line Options]] for values that can be set either via command line or via configuration.

Note: Command-line options take priority over configuration file settings

==== Required Settings ====
See [[#Reference Files|Reference Files]] for the required reference file settings.

See [[#Sequence Index File|Sequence Index File]] for how to set the index file either via command line options or via configuration.

==== Turning Off Optional Steps====
Quality Control steps can be disabled.

To Disable QPLOT, set:
RUN_QPLOT = 0

To Disable VerifyBamID, set:
RUN_VERIFY_BAM_ID = 0

==== Optional Configurable Settings ====
You may want to adjust the amount of memory/threads that are used:

There are additional configurable settings, but these are the ones most likely to be adjusted.

* BWA_THREADS = -t N
** Fill in the N with the number of threads you want BWA to run with, default is 1
* BWA_MAX_MEM = 2000000000
** Maximum amount of memory used by samtools sort after running bwa
* JAVA_MEM = -Xmx4g
** Set the maximum size of the java memory allocation pool. Default is 4g, adjust that as necessary.

== Running the Alignment Pipeline ==

=== Command-Line Options ===
* help - print usage
* test OUTPUT_DIR - run the test example placing the output in a user specified OUTPUT_DIR. No other options are required.
* out_dir OUTPUT_DIR - directory for the output
** May also be specified via OUT_DIR in the configuration file
** Required to be set either via command-line or configuration
* conf CONFIG_FILE - configuration file
* index_file INDEX_FILE_NAME - name of the index file
** May also be specified via INDEX_FILE in the configuration file
** Required to be set either via command-line or configuration
* ref_dir REFERENCE_DIR - value to set config key REF_DIR to, overriding other values, REF_DIR can then be used inside config files.
** May also be specified via REF_DIR in the configuration file
* fastq FASTQ_PATH - prefix path to the fastq files specified in the INDEX_FILE
** May also be specified via FASTQ in the configuration file
* keepTmp - Do not remove the temporary files (removed by default)
** May also be specified via KEEP_TMP in the configuration file
* numcs N - Replace N with the number of samples that should be processed in parallel
* numjobs N - Replace N with the number of targets in each makefile that should be run in parallel

Note: Command-line options take priority over configuration file settings

===Running the Alignment Pipeline===
Run <code>gotcloud align</code> with the appropriate command-line parameters.

Example:
gotcloud align --conf config.txt --outdir output

This step generates 1 Makefile per sample in the output/Makefiles/ directory and then automatically runs them. The Makefiles contain all of the information to run each sample.

If you only want to generate the makefiles and not run them, use the <code>--dryrun</code> option. It will generate the Makefiles and print instructions for running the Makefiles.

Each Makefile is independent and can be run in parallel and across a cloud.

On success, you will see:
Processing finished in nn secs with no errors reported
and should see the following subdirectories under the user specified output directory:
* bams/
* Makefiles/
* QCFiles/ (if all quality control is not disabled)
* tmp/

You should see a <code>.OK</code> for each Sample in the index file.

If you do not see these <code>.OK</code> files, then your Alignment Pipeline failed.

On success, the bams/ directory contains the final BAMs and bais.

If processing fails part way through, you can pick up where you left off by rerunning gotcloud or the make command.

GotCloud: Alignment Pipeline

2013-06-21T20:29:52Z

Ben Lerch: /* Command-Line Options */

= Alignment Pipeline =

Back to parent: [[GotCloud]]

The Alignment/Mapping Pipeline takes FASTQ files and generates recalibrated BAM files from them.

== Running the GotCloud Alignment Pipeline ==

The alignment pipeline is run using the <code>align</code> option of the <code>gotcloud</code> script. This option calls <code>align.pl</code> found in the <code>bin/</code> directory under the <code>gotcloud</code> installation.

===Running the Automated Test===

The automated test runs the alignment pipeline on a small set of test data and checks that the results against expected results validating that GotCloud is installed correctly.

*Run alignment pipeline test:
gotcloud align --test OUTPUT_DIR
where OUTPUT_DIR is the directory where you want to store the test results

If you see "Successfully ran the test case, congratulations!", then you are ready to align samples.

== Overview of Alignment Pipeline Steps ==
Here is an overview of the Alignment Pipeline:

[[File:MappingSteps.png]]

== Input Data:==
*Raw Sequence (FASTQ) files
*Sequence Index file containing fastqs & RG info
*Reference files
*(Optional) Configuration file to override default options

=== Raw Sequence (FASTQ) files ===

These are the FASTQ files that need to be mapped to BAM files.

These files are specified in the [[#Sequence Index File|Sequence Index File]].

=== Sequence Index File ===
This file specifies the FASTQ files that need to be processed and the Read Group information for them.

This file is specified either via the command line parameter <code>--index_file</code> or via the configuration file setting <code>INDEX_FILE</code>.

The command-line setting takes precedence over the configuration file setting.

The Sequence Index is a tab delimited file that starts with a header line. The columns may be in any order.

Following the header line, there is one line per single-end read and one line per paired-end read (only 1 line per pair).

Required Column Names:
* MERGE_NAME - base name for the resulting BAM file for the sample (used to group multiple fastqs or fastq pairs into a single BAM)
* FASTQ1 - name of the fastq or the first in the pair if paired-end. (Only 1 line per pair)

Optional Column Names:
* FASTQ2 - name of the 2nd fastq in paired-end reads. Specify '.' if the column exists, but this line is single-ended.
* RGID - Read Group ID for this entry
* SAMPLE - Sample Name for this entry
* LIBRARY - Library for this entry
* CENTER - Center Name for this entry
* PLATFORM - Platform for this entry

The RGID, SAMPLE, LIBRARY, CENTER, and PLATFORM are used to populate the Read Group information for this entry. These fields are optional. Either leave the column header out of the file or specify '.' if the column header exists, but the data is N/A. As long as the RGID field is specified non-N/A fields are added to the BAM file.

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/S1/F1_R1.fastq.gz fastq/S1/F1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/S1/F2_R1.fastq.gz fastq/S1/F2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/S2/F1_R1.fastq.gz fastq/S2/F1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/S2/F2.fastq.gz . RGID2 SampleID2 Lib2 UM ILLUMINA

The <code>--fastq</code>/<code>FASTQ</code> setting can be used to specify a prefix to the FASTQ1/FASTQ2 file paths that should be applied before using the files.

=== Reference Files ===

The following Reference Files are required:
* Reference File fasta files
** Files required: .fa, -bs.umfa, .GCContent, .amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa
*** If you don't have the -bs.umfa file, the software will try to create it in the same directory as the reference fasta.
*** .GCContent can be generated using qplot, see: [[QPLOT#Input_files| QPLOT: Input Files: --gccontent]] and name the resulting file as <code>.fa.GCcontent</code>
*** Use <code>bin/bwa index ref.fa</code> if you need to generate the bwa reference files (.amb, .ann, .bwt, .pac, .rbwt, .rpac, .rsa, .sa)
** Configuration Name: FA_REF - specify the ref.fa/ref.fa.gz name
* DBSNP File
** tab delimited file/VCF, can be compressed
*** 1st column -> chromosome
*** 2nd column -> 1-based position
** Configuration Name: DBSNP_VCF
* PLINK-compatible binary genotype files
** Files required: .bed, .bin, .fam
** Configuration Name: PLINK

=== Configuration File ===
Configuration file contains the run-time options including the software binaries and command line arguments. A default configuration file is automatically loaded. Users may specify their own configuration file specifying just the values different than the defaults. The configuration file is not required if there are no values to override.

Comments begin with a <code>#</code>

Format: KEY = value

Where KEY is the item being set and value is its new value

See [[#Command-Line Options|Command-Line Options]] for values that can be set either via command line or via configuration.

Note: Command-line options take priority over configuration file settings

==== Required Settings ====
See [[#Reference Files|Reference Files]] for the required reference file settings.

See [[#Sequence Index File|Sequence Index File]] for how to set the index file either via command line options or via configuration.

==== Turning Off Optional Steps====
Quality Control steps can be disabled.

To Disable QPLOT, set:
RUN_QPLOT = 0

To Disable VerifyBamID, set:
RUN_VERIFY_BAM_ID = 0

==== Optional Configurable Settings ====
You may want to adjust the amount of memory/threads that are used:

There are additional configurable settings, but these are the ones most likely to be adjusted.

* BWA_THREADS = -t N
** Fill in the N with the number of threads you want BWA to run with, default is 1
* BWA_MAX_MEM = 2000000000
** Maximum amount of memory used by samtools sort after running bwa
* JAVA_MEM = -Xmx4g
** Set the maximum size of the java memory allocation pool. Default is 4g, adjust that as necessary.

== Running the Alignment Pipeline ==

=== Command-Line Options ===
* help - print usage
* test OUTPUT_DIR - run the test example placing the output in a user specified OUTPUT_DIR. No other options are required.
* out_dir OUTPUT_DIR - directory for the output
** May also be specified via OUT_DIR in the configuration file
** Required to be set either via command-line or configuration
* conf CONFIG_FILE - configuration file
* index_file INDEX_FILE_NAME - name of the index file
** May also be specified via INDEX_FILE in the configuration file
** Required to be set either via command-line or configuration
* ref_dir REFERENCE_DIR - value to set config key REF_DIR to, overriding other values, REF_DIR can then be used inside config files.
** May also be specified via REF_DIR in the configuration file
* fastq FASTQ_PATH - prefix path to the fastq files specified in the INDEX_FILE
** May also be specified via FASTQ in the configuration file
* keepTmp - Do not remove the temporary files (removed by default)
** May also be specified via KEEP_TMP in the configuration file
* numjobs N - Replace N with the number of jobs that should be run in parallel

Note: Command-line options take priority over configuration file settings

===Running the Alignment Pipeline===
Run <code>gotcloud align</code> with the appropriate command-line parameters.

Example:
gotcloud align --conf config.txt --outdir output

This step generates 1 Makefile per sample in the output/Makefiles/ directory and then automatically runs them. The Makefiles contain all of the information to run each sample.

If you only want to generate the makefiles and not run them, use the <code>--dryrun</code> option. It will generate the Makefiles and print instructions for running the Makefiles.

Each Makefile is independent and can be run in parallel and across a cloud.

On success, you will see:
Processing finished in nn secs with no errors reported
and should see the following subdirectories under the user specified output directory:
* bams/
* Makefiles/
* QCFiles/ (if all quality control is not disabled)
* tmp/

You should see a <code>.OK</code> for each Sample in the index file.

If you do not see these <code>.OK</code> files, then your Alignment Pipeline failed.

On success, the bams/ directory contains the final BAMs and bais.

If processing fails part way through, you can pick up where you left off by rerunning gotcloud or the make command.

QPLOT

2013-06-17T17:45:17Z

Ben Lerch: /* Input files */

= Introduction =

The qplot program calculates various summary statistics some of which are plotted in a PDF file. These statistics can be used to assess the sequencing quality of sequence reads mapped to the reference genome. The main statistics are empirical Phred scores which are calculated based on the background mismatch rate. Background mismatch rate is the rate that sequenced bases are different from the reference genome, EXCLUDING dbSNP positions. Other statistics include GC biases, insert size distribution, depth distribution, genome coverage, empirical Q20 count, and so on.

In the following sections, we will guide you through: [[#Where to Find It |how to obtain qplot]], [[#Usage |how to use qplot]], [[#Built-in example |example outputs]], [[#anchorOfInteractiveQplot |interactive diagnostic plots]], and [[#Diagnose sequencing quality |real applications]] in which qplot has helped identify sequencing problems.

= Where to Find It =

You can obtain qplot in two ways:

(1) Download the pre-compiled binary along with the source code as described in [[#Binary Download|Binary Download]].

(2) Download source code only and compile it on your own machine. Please follow the instruction in [[#Source Code Distribution|Source Code Distribution]] on fetching source code and building instructions.

== Binary Download ==

We have prepared a pre-compiled (under Ubuntu) qplot along with source code . You can download it from: [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot.20120602.tar.gz qplot.20120602.tar.gz (File Size: 1.7G)]

The executable file is under qplot/bin/qplot.

In addition, we provided the necessary input files under qplot/data/ (NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100).

You can also find an example BAM input file under qplot/example/chrom20.9M.10M.bam. It is taken from the 1000 Genome Project with sequencing reads aligned to chromosome 20 positions 8M to 9M.

== Source Code Distribution ==

We provide a source code only download in [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-source.20120602.tar.gz qplot-source.20120602.tar.gz]. Optionally, you can download example file and/or data file:

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-example.tar.gz example]: example input file, and expected outputs if you following the [[#Built-in example | direction]].

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-data.tar.gz resources data]: necessary input files for qplot, including NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100.

You can put above file(s) in the same folder and follow these steps:

* 1. Unarchive downloaded file
tar zvxf qplot-source.20120602.tar.gz

A new folder ''qplot'' will be created.

* 2. Build libStatGen
cd qplot
make libStatGen

This step will download a necessary software library [http://genome.sph.umich.edu/wiki/C%2B%2B_Library:_libStatGen libStatGen] and compile source code into a binary code library.

* 3. Build qplot
make all

This step will then build qplot. Upon success, the executable qplot can be found under qplot/bin/.

* 4. (Optional) unarchive example and/or data
tar zvxf qplot-example.tar.gz

An example file, ''chrom20.9M.10M.bam'', will be extracted to qplot/example/. It contains ~1.1 million aligned Illumina sequencing reads of NA12878 from 1000 Genome Project. Example command line, ''cmd.sh'', example outputs, ''qplot.pdf'', ''qplot.stats'', and ''qplot.R'' are also provided and will be extracted qplot/example/ as well.

tar zvxf qplot-data.tar.gz

Three files will be extracted to qplot/data/: ''human.g1k.v37-bs.umfa'' is binary NCBI reference genome build 37; ''dbSNP130.UCSC.coordinates.tbl'' is dbSNP version 130; and ''human.g1k.w100.gc'' is pre-calculated GC content with windows size 100.



= Usage =

== Command line ==

After you obtain the qplot executable (either by compiling the source code or by downloading the pre-compiled binary file), you will find the executable file under qplot/bin/qplot.

Here is the qplot help page by invoking qplot without any command line arguments:

some_linux_host > qplot/bin/qplot

References : --reference [/net/fantasia/home/zhanxw/software/qplot/data/human.g1k.v37.fa],
--dbsnp [/net/fantasia/home/zhanxw/software/qplot/data/dbSNP130.UCSC.coordinates.tbl],
--gccontent [/net/fantasia/home/zhanxw/software/qplot/data/human.g1k.w100.gc]
Create gcContent file : --create_gc [], --winsize [100]
Region list : --regions [], --invertRegion
Flag filters : --read1_skip, --read2_skip, --paired_skip,
--unpaired_skip
Dup and QCFail : --dup_keep, --qcfail_keep
Mapping filters : --minMapQuality [0.00]
Records to process : --first_n_record [-1]
Lanes to process : --lanes []
Read group to process : --readGroup []
Input file options : --noeof
Output files : --plot [], --stats [], --Rcode [], --xml []
Plot labels : --label [], --bamLabel []

== Input files ==

qplot runs on the input BAM/SAM file(s) specified on the command-line after all other parameters.

Additionally, three (3) precomputed files are required.

* <code>--reference</code>

The reference genome is the same as karma reference genome. If the index files do not exist, qplot will create the index files '''automatically''' using the input reference fasta file.

* <code>--dbsnp</code>

This file has two columns. First column is the chromosome name which must be consistent with the reference created above. Second column is 1-based SNP position. If you want to create your own dbSNP data from downloaded UCSC dbSNP file, one way to do it is: <code>cat dbsnp_129_b36.rod|grep "single" | awk '$4-$3==1' |cut -f2,4 > dbSNP_129_b36.tbl</code>

* <code>--gccontent</code>

Although GC content can be calculated on the fly each time, it is much more efficient to load a precomputed GC content from a file. To generate the file, use the following command:
qplot --reference reference.fa --windowsize winsize --create_gc reference.gc

''Note'': Before running qplot, it is critical to check how the chromosome names are coded. Some BAM/SAM files use just numbers, others use chr + numbers. '''You need to make sure that the chromosome names from the reference and dbSNP are consistent with the BAM/SAM files.'''

== Parameters ==

Some of the command line parameters are described here, but most are self explanatory.

*Flag filter

By default all reads are processed. If it is desired to check only the first read of a pair, use <code>--read2_skip</code> to ignore the second read. And so on.

*Duplication and QCFail

By default reads marked as duplication and QCFail are ignored but can be retained by
--dup_keep
or
--qcfail_keep

*Records to process

The <code>--first_n_record</code> option followed by a number, '''n''', will enable qplot to read the first '''n''' reads to test the bam files and verify it works.

* Lanes to process (only works for Illumina sequences)

If the input bam files have more than one lane and only some of them need to be checked, use something like <code>--lanes 1,3,5</code> to specify that only lanes 1, 3, and 5 need to be checked.

'''NOTE''' In order for this to work, the lane info has to be encoded in the read name such that the lane number is the second field with the delimiter ":".

* Read group to process :

Read group option can restrict qplot to process a subset of reads. For example, if BAM contain the following @RG tags:

@RG ID:UM0348_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM

If specify nothing or not using "--readGroup", QPLOT by default will process all reads;
If specify "--readGroup UM0348", then only read group UM0348_1, UM_0348_2, UM_0348_3, UM_0348_4 will be processed;
If specify "--readGroup UM0348_1", then only one read group UM0348_1 will be processed.

* Input file options :

BAM files are compress by BGZF algorithm and it should contain EOF by default. QPLOT will by default stop working when it does not found a valid EOF tag inside BAM files.
However, you can force QPLOT to continue process using --noeof. But you should be award the input files may be corrupted.

* Mapping filters

Qplot will exclude reads with lower mapping qualities than the user specified parameter, <code>--minMapQuality</code>. By default, mapped reads with all mapping quality will be included in the analysis.

*Region list

If the interest of qplot is a list of regions, e.g. exons, this can be achieved by providing a list of regions. The regions should be in the form of "chr start end label" each line in the file (NOTE: ''start'' and ''end'' position are inclusive and they follow the convention of [http://genome.ucsc.edu/FAQ/FAQformat#format1 BED file]).
In order for this option to work, within each chromosome (contig) the regions have to be sorted by starting position, and also the input bam files have to be sorted.
For example, you can create a text file, region.txt like following:

1 100 500 region_A
1 600 800 region_B
2 100 300 region_C

Then specifying <code> --regions region.txt</code> enables qplot to calculate various statistics out of sequenced bases only within the above 3 regions.

Qplot also provides the <code>--invertRegion</code> option. Enabling this option tells qplot to operate on those sequence bases that are outside the given region.

* Plot labels

Two kinds of labels are enabled. <code>--label</code> is the label for the plot (default is empty) which is appended to the title of each subplot. <code>--bamLabels</code> followed by a column separated list of labels provides the labels for each input SAM/BAM file, e.g. sample ID (default is numbers 1, 2, ... until the number of input bam files). For example:
--label Run100 --bamLabels s1,s2,s3,s4,s5,s6,s7,s8

== Output files ==

There are three (optional) output files.
* <code>--plot ''qa.pdf''</code>

Qplot will generate a PDF file named ''qa.pdf'' containing 2 pages each with 4 figures. The plot is generated using Rscript.

* <code>--stats ''qa.stats''</code>

Qplot will generate a text file named ''qa.stats'' containing various summary statistics for each input BAM/SAM file.

* <code>--Rcode ''qa.R''</code>

Qplot will generate ''qa.R'' which is the R code used for plotting the figures in the ''qa.pdf'' file. If Rscript is not installed in the system, you can use the qa.R to generate the figures on other machines, or extract plotting data from each run and combine multiple runs together to generate more comprehensive plots (See [[#Example | Example]]).

= Example =

Qplot can generate diagnostic graphs, related R code, and summary statistics for each SAM/BAM file.

== Built-in example ==

In the pre-compiled binary download, you will find a subdirectory named examples. We provide a sample file from the 1000 Genome project, it contains aligned reads on chromosome 20 from position 8 Mbp to 9Mbp. You can invoke qplot using the following command line:

../bin/qplot --reference ../data/human.g1k.v37.umfa --dbsnp ../data/dbSNP130.UCSC.coordinates.tbl --gccontent ../data/human.g1k.w100.gc --plot qplot.pdf --stats qplot.stats --Rcode qplot.R --label "chr20:9M-10M" chrom20.9M.10M.bam

Sample outputs are listed below:

1) Figure: [[Media:qplot.pdf | qplot.pdf]]

2) Summary statistics:
Stats\BAM chrom20.9M.10M.bam
TotalReads(e6) 1.11
MappingRate(%) 97.24
MapRate_MQpass(%) 97.24
TargetMapping(%) 0.00
ZeroMapQual(%) 2.39
MapQual<10(%) 2.86
PairedReads(%) 83.76
ProperPaired(%) 71.34
MappedBases(e9) 0.04
Q20Bases(e9) 0.04
Q20BasesPct(%) 88.63
MeanDepth 42.22
GenomeCover(%) 0.03
EPS_MSE 1.81
EPS_Cycle_Mean 18.71
GCBiasMSE 0.01
ISize_mode 137
ISize_medium 184
DupRate(%) 5.90
QCFailRate(%) 0.00
BaseComp_A(%) 29.9
BaseComp_C(%) 20.1
BaseComp_G(%) 20.2
BaseComp_T(%) 29.8
BaseComp_O(%) 0.1

== Gallery of examples ==

Here we show qplot can be applied in various sequencing scenarios. Also users can customize statistics generated by qplot to their needs.

* Whole genome sequencing with 24-multiplexing

With a customized script, we aggregated 24 bar-coded samples in the same graph.
The graph will help compare sequencing quality between samples.

[[Media: qplot.Pool.9847.pdf | QPlot of 24 samples(PDF) ]]

* Interactive qplot


Qplot can be interactive. In the following example, you can use mouse scroll to zoom in and zoom out on each graph and pan to a certain part of the graph.
By presenting qplot data on a web page, users can easily identify problematic sequencing samples. Users of qplot can customize its outputs into web page format greatly easing the data exploring process.

[http://www-personal.umich.edu/~zhanxw/qplot.Pool.9847.html QPlot of 24 samples(HTML) ]

== Diagnose sequencing quality ==

Qplot is designed and implemented for the need of checking sequencing quality.
Besides the example of analyzing RNA-seq data as shown in our manuscript,
here we demonstrate two additional scenarios in which qplot can help identify problems after obtaining sequencing data.

* Base quality distributed abnormally

[[Media: WrongBaseQual.pdf | Example of qplot helping to identify wrong phred base quality]]

By checking the first graph "Empirical vs reported Phred score", we found reported base qualities are shifted to the right.
In this particular example, '33' was incorrectly added to all base qualities.
When such data used in variant calling, we may increase false positive SNP variants.

* Bar-coded samples

[[Media: WrongBarCoding.pdf | Example of qplot identifying the effect of ignoring bar-coding]]

By checking "Empirical phred score by cycle" (top right graph on the first page), we noticed the empirical qualities in the first several cycles are abnormally low. This phenomenon leads us to hypothesize that the first several bases have different properties. Further investigation confirmed that this sequencing was done using bar-coded DNA samples, but the analysis did not properly de-multiplex each sample.

= Contact =

Questions and requests should be sent to Bingshan Li ([mailto:bingshan@umich.edu bingshan@umich.edu]) or Xiaowei Zhan ([mailto:zhanxw@umich.edu zhanxw@umich.edu]) or Goncalo Abecasis ([mailto:goncalo@umich.edu goncalo@umich.edu])

QPLOT

2013-06-17T17:44:58Z

Ben Lerch: /* Input files */

= Introduction =

The qplot program calculates various summary statistics some of which are plotted in a PDF file. These statistics can be used to assess the sequencing quality of sequence reads mapped to the reference genome. The main statistics are empirical Phred scores which are calculated based on the background mismatch rate. Background mismatch rate is the rate that sequenced bases are different from the reference genome, EXCLUDING dbSNP positions. Other statistics include GC biases, insert size distribution, depth distribution, genome coverage, empirical Q20 count, and so on.

In the following sections, we will guide you through: [[#Where to Find It |how to obtain qplot]], [[#Usage |how to use qplot]], [[#Built-in example |example outputs]], [[#anchorOfInteractiveQplot |interactive diagnostic plots]], and [[#Diagnose sequencing quality |real applications]] in which qplot has helped identify sequencing problems.

= Where to Find It =

You can obtain qplot in two ways:

(1) Download the pre-compiled binary along with the source code as described in [[#Binary Download|Binary Download]].

(2) Download source code only and compile it on your own machine. Please follow the instruction in [[#Source Code Distribution|Source Code Distribution]] on fetching source code and building instructions.

== Binary Download ==

We have prepared a pre-compiled (under Ubuntu) qplot along with source code . You can download it from: [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot.20120602.tar.gz qplot.20120602.tar.gz (File Size: 1.7G)]

The executable file is under qplot/bin/qplot.

In addition, we provided the necessary input files under qplot/data/ (NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100).

You can also find an example BAM input file under qplot/example/chrom20.9M.10M.bam. It is taken from the 1000 Genome Project with sequencing reads aligned to chromosome 20 positions 8M to 9M.

== Source Code Distribution ==

We provide a source code only download in [http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-source.20120602.tar.gz qplot-source.20120602.tar.gz]. Optionally, you can download example file and/or data file:

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-example.tar.gz example]: example input file, and expected outputs if you following the [[#Built-in example | direction]].

[http://www.sph.umich.edu/csg/zhanxw/software/qplot/qplot-data.tar.gz resources data]: necessary input files for qplot, including NCBI human genome build v37, dbSNP 130, and pre-computed GC file with windows size 100.

You can put above file(s) in the same folder and follow these steps:

* 1. Unarchive downloaded file
tar zvxf qplot-source.20120602.tar.gz

A new folder ''qplot'' will be created.

* 2. Build libStatGen
cd qplot
make libStatGen

This step will download a necessary software library [http://genome.sph.umich.edu/wiki/C%2B%2B_Library:_libStatGen libStatGen] and compile source code into a binary code library.

* 3. Build qplot
make all

This step will then build qplot. Upon success, the executable qplot can be found under qplot/bin/.

* 4. (Optional) unarchive example and/or data
tar zvxf qplot-example.tar.gz

An example file, ''chrom20.9M.10M.bam'', will be extracted to qplot/example/. It contains ~1.1 million aligned Illumina sequencing reads of NA12878 from 1000 Genome Project. Example command line, ''cmd.sh'', example outputs, ''qplot.pdf'', ''qplot.stats'', and ''qplot.R'' are also provided and will be extracted qplot/example/ as well.

tar zvxf qplot-data.tar.gz

Three files will be extracted to qplot/data/: ''human.g1k.v37-bs.umfa'' is binary NCBI reference genome build 37; ''dbSNP130.UCSC.coordinates.tbl'' is dbSNP version 130; and ''human.g1k.w100.gc'' is pre-calculated GC content with windows size 100.



= Usage =

== Command line ==

After you obtain the qplot executable (either by compiling the source code or by downloading the pre-compiled binary file), you will find the executable file under qplot/bin/qplot.

Here is the qplot help page by invoking qplot without any command line arguments:

some_linux_host > qplot/bin/qplot

References : --reference [/net/fantasia/home/zhanxw/software/qplot/data/human.g1k.v37.fa],
--dbsnp [/net/fantasia/home/zhanxw/software/qplot/data/dbSNP130.UCSC.coordinates.tbl],
--gccontent [/net/fantasia/home/zhanxw/software/qplot/data/human.g1k.w100.gc]
Create gcContent file : --create_gc [], --winsize [100]
Region list : --regions [], --invertRegion
Flag filters : --read1_skip, --read2_skip, --paired_skip,
--unpaired_skip
Dup and QCFail : --dup_keep, --qcfail_keep
Mapping filters : --minMapQuality [0.00]
Records to process : --first_n_record [-1]
Lanes to process : --lanes []
Read group to process : --readGroup []
Input file options : --noeof
Output files : --plot [], --stats [], --Rcode [], --xml []
Plot labels : --label [], --bamLabel []

== Input files ==

qplot runs on the input BAM/SAM file(s) specified on the command-line after all other parameters.

Additionally, three (3) precomputed files are required.

* <code>--reference</code>

The reference genome is the same as karma reference genome. If the index files do not exist, qplot will create the index files '''automatically''' using the input reference fasta file.

* <code>--dbsnp</code>

This file has two columns. First column is the chromosome name which must be consistent with the reference created above. Second column is 1-based SNP position. If you want to create your own dbSNP data from downloaded UCSC dbSNP file, one way to do it is: <code>cat dbsnp_129_b36.rod|grep "single" | awk '$4-$3==1' |cut -f2,4 > dbSNP_129_b36.tbl</code>

* <code>--gccontent</code>

Although GC content can be calculated on the fly each time, it is much more efficient to load a precomputed GC content from a file. To generate the file, use the following command:
qplot --refence reference.fa --windowsize winsize --create_gc reference.gc

''Note'': Before running qplot, it is critical to check how the chromosome names are coded. Some BAM/SAM files use just numbers, others use chr + numbers. '''You need to make sure that the chromosome names from the reference and dbSNP are consistent with the BAM/SAM files.'''

== Parameters ==

Some of the command line parameters are described here, but most are self explanatory.

*Flag filter

By default all reads are processed. If it is desired to check only the first read of a pair, use <code>--read2_skip</code> to ignore the second read. And so on.

*Duplication and QCFail

By default reads marked as duplication and QCFail are ignored but can be retained by
--dup_keep
or
--qcfail_keep

*Records to process

The <code>--first_n_record</code> option followed by a number, '''n''', will enable qplot to read the first '''n''' reads to test the bam files and verify it works.

* Lanes to process (only works for Illumina sequences)

If the input bam files have more than one lane and only some of them need to be checked, use something like <code>--lanes 1,3,5</code> to specify that only lanes 1, 3, and 5 need to be checked.

'''NOTE''' In order for this to work, the lane info has to be encoded in the read name such that the lane number is the second field with the delimiter ":".

* Read group to process :

Read group option can restrict qplot to process a subset of reads. For example, if BAM contain the following @RG tags:

@RG ID:UM0348_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0348_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_1:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_2:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_3:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM
@RG ID:UM0360_4:1 PL:ILLUMINA LB:M5390 SM:M5390 CN:UM

If specify nothing or not using "--readGroup", QPLOT by default will process all reads;
If specify "--readGroup UM0348", then only read group UM0348_1, UM_0348_2, UM_0348_3, UM_0348_4 will be processed;
If specify "--readGroup UM0348_1", then only one read group UM0348_1 will be processed.

* Input file options :

BAM files are compress by BGZF algorithm and it should contain EOF by default. QPLOT will by default stop working when it does not found a valid EOF tag inside BAM files.
However, you can force QPLOT to continue process using --noeof. But you should be award the input files may be corrupted.

* Mapping filters

Qplot will exclude reads with lower mapping qualities than the user specified parameter, <code>--minMapQuality</code>. By default, mapped reads with all mapping quality will be included in the analysis.

*Region list

If the interest of qplot is a list of regions, e.g. exons, this can be achieved by providing a list of regions. The regions should be in the form of "chr start end label" each line in the file (NOTE: ''start'' and ''end'' position are inclusive and they follow the convention of [http://genome.ucsc.edu/FAQ/FAQformat#format1 BED file]).
In order for this option to work, within each chromosome (contig) the regions have to be sorted by starting position, and also the input bam files have to be sorted.
For example, you can create a text file, region.txt like following:

1 100 500 region_A
1 600 800 region_B
2 100 300 region_C

Then specifying <code> --regions region.txt</code> enables qplot to calculate various statistics out of sequenced bases only within the above 3 regions.

Qplot also provides the <code>--invertRegion</code> option. Enabling this option tells qplot to operate on those sequence bases that are outside the given region.

* Plot labels

Two kinds of labels are enabled. <code>--label</code> is the label for the plot (default is empty) which is appended to the title of each subplot. <code>--bamLabels</code> followed by a column separated list of labels provides the labels for each input SAM/BAM file, e.g. sample ID (default is numbers 1, 2, ... until the number of input bam files). For example:
--label Run100 --bamLabels s1,s2,s3,s4,s5,s6,s7,s8

== Output files ==

There are three (optional) output files.
* <code>--plot ''qa.pdf''</code>

Qplot will generate a PDF file named ''qa.pdf'' containing 2 pages each with 4 figures. The plot is generated using Rscript.

* <code>--stats ''qa.stats''</code>

Qplot will generate a text file named ''qa.stats'' containing various summary statistics for each input BAM/SAM file.

* <code>--Rcode ''qa.R''</code>

Qplot will generate ''qa.R'' which is the R code used for plotting the figures in the ''qa.pdf'' file. If Rscript is not installed in the system, you can use the qa.R to generate the figures on other machines, or extract plotting data from each run and combine multiple runs together to generate more comprehensive plots (See [[#Example | Example]]).

= Example =

Qplot can generate diagnostic graphs, related R code, and summary statistics for each SAM/BAM file.

== Built-in example ==

In the pre-compiled binary download, you will find a subdirectory named examples. We provide a sample file from the 1000 Genome project, it contains aligned reads on chromosome 20 from position 8 Mbp to 9Mbp. You can invoke qplot using the following command line:

../bin/qplot --reference ../data/human.g1k.v37.umfa --dbsnp ../data/dbSNP130.UCSC.coordinates.tbl --gccontent ../data/human.g1k.w100.gc --plot qplot.pdf --stats qplot.stats --Rcode qplot.R --label "chr20:9M-10M" chrom20.9M.10M.bam

Sample outputs are listed below:

1) Figure: [[Media:qplot.pdf | qplot.pdf]]

2) Summary statistics:
Stats\BAM chrom20.9M.10M.bam
TotalReads(e6) 1.11
MappingRate(%) 97.24
MapRate_MQpass(%) 97.24
TargetMapping(%) 0.00
ZeroMapQual(%) 2.39
MapQual<10(%) 2.86
PairedReads(%) 83.76
ProperPaired(%) 71.34
MappedBases(e9) 0.04
Q20Bases(e9) 0.04
Q20BasesPct(%) 88.63
MeanDepth 42.22
GenomeCover(%) 0.03
EPS_MSE 1.81
EPS_Cycle_Mean 18.71
GCBiasMSE 0.01
ISize_mode 137
ISize_medium 184
DupRate(%) 5.90
QCFailRate(%) 0.00
BaseComp_A(%) 29.9
BaseComp_C(%) 20.1
BaseComp_G(%) 20.2
BaseComp_T(%) 29.8
BaseComp_O(%) 0.1

== Gallery of examples ==

Here we show qplot can be applied in various sequencing scenarios. Also users can customize statistics generated by qplot to their needs.

* Whole genome sequencing with 24-multiplexing

With a customized script, we aggregated 24 bar-coded samples in the same graph.
The graph will help compare sequencing quality between samples.

[[Media: qplot.Pool.9847.pdf | QPlot of 24 samples(PDF) ]]

* Interactive qplot


Qplot can be interactive. In the following example, you can use mouse scroll to zoom in and zoom out on each graph and pan to a certain part of the graph.
By presenting qplot data on a web page, users can easily identify problematic sequencing samples. Users of qplot can customize its outputs into web page format greatly easing the data exploring process.

[http://www-personal.umich.edu/~zhanxw/qplot.Pool.9847.html QPlot of 24 samples(HTML) ]

== Diagnose sequencing quality ==

Qplot is designed and implemented for the need of checking sequencing quality.
Besides the example of analyzing RNA-seq data as shown in our manuscript,
here we demonstrate two additional scenarios in which qplot can help identify problems after obtaining sequencing data.

* Base quality distributed abnormally

[[Media: WrongBaseQual.pdf | Example of qplot helping to identify wrong phred base quality]]

By checking the first graph "Empirical vs reported Phred score", we found reported base qualities are shifted to the right.
In this particular example, '33' was incorrectly added to all base qualities.
When such data used in variant calling, we may increase false positive SNP variants.

* Bar-coded samples

[[Media: WrongBarCoding.pdf | Example of qplot identifying the effect of ignoring bar-coding]]

By checking "Empirical phred score by cycle" (top right graph on the first page), we noticed the empirical qualities in the first several cycles are abnormally low. This phenomenon leads us to hypothesize that the first several bases have different properties. Further investigation confirmed that this sequencing was done using bar-coded DNA samples, but the analysis did not properly de-multiplex each sample.

= Contact =

Questions and requests should be sent to Bingshan Li ([mailto:bingshan@umich.edu bingshan@umich.edu]) or Xiaowei Zhan ([mailto:zhanxw@umich.edu zhanxw@umich.edu]) or Goncalo Abecasis ([mailto:goncalo@umich.edu goncalo@umich.edu])

Tutorial: GotCloud

2013-06-16T23:03:49Z

Ben Lerch: /* Frequently Asked Questions (FAQs) */

= GotCloud Tutorial =
In this tutorial, we illustrate some of the essential steps in the analysis of next generation sequence data.

For a background on GotCloud and Sequence Analysis Pipelines, see [[GotCloud]]

While GotCloud can run on a cluster of machines or instances, this tutorial is just a small test that just runs on the machine the commands are run on.

GotCloud and this basic tutorial were presented at the [http://ibg.colorado.edu/dokuwiki/doku.php?id=workshop:2013:announcement 2013 IBG Workshop]. It was presented in two sessions. On Wednesday an overview was presented with steps for running the tutorial data: [[Media:IBG2013GotCloud.pdf|IBG2013GotCloud.pdf]]. On Friday more detail on the input files and what goes into generating the input files was presented: [[Media:GotCloudIBGWorkshop2013Friday.pdf|GotCloudIBGWorkshop2013Friday.pdf]].

'''This tutorial is in the process of being updated for gotcloud version 1.06 (April 17. 2013).'''

== STEP 1 : Setup GotCloud ==

[[GotCloud]] has been developed and tested on Linux Ubuntu 12.10 and 12.04.2 LTS but has not been tested on other Linux operating systems. It is not available for Windows. If you do not have your own set of machines to run on, GotCloud is also available for Ubuntu running on the Amazon Elastic Compute Cloud, see [[Amazon_Snapshot]] for more information.

We will use 3 different directories for this tutorial:
# path to the directory where gotcloud is installed, default ~/gotcloud/
# path to the directory where the example data is installed, default ~/gotcloudExample
# path to your output directory, default ~/gotcloudTutorialOut/

If the directories specified above do not reflect the directories you would like to use, replace their occurrances in the instructions below with the appropriate paths.

=== Step 1a: Install GotCloud ===
In order to run this tutorial, you need to make sure you have GotCloud installed on your system.

If you have root and would like to install gotcloud on your system, follow: [[GotCloud#Install_GotCloud_Software| root access installation instructions]]

Otherwise, you can install it in your own directory:
# Change to the directory where you want gotcloud/ installed
# Download the gotcloud tar from the ftp site.
# Extract the tar
# Build (compile) the source
#* Note: as the source builds, many messages will scroll through your terminal. You may even see some warnings. These messages are normal and expected. As long as the build does not end with an error, you have successfully built the source.

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloud_latest.tgz # Download
tar xf gotcloud_latest.tgz # Extracts into gotcloud/
cd ~/gotcloud/src; make # Build source

GotCloud requires the following tools to be installed.
You can run ~/gotcloud/scripts/check_requirements.sh
...TBD – put in required programs/tools.
* java (java-common default-jre on ubuntu)
* make (make on ubuntu)
* libssl (libssl0.9.8 on ubuntu)
* gcc 4.4 or newer

=== Step 1b: Install Example Dataset ===
Our dataset consists of 60 individuals from Great Britain (GBR) sequenced by the 1000 Genomes Project. These individuals have been sequenced to an average depth of about 4x.

To conserve time and disk-space, our analysis will focus on a small region on chromosome 20, 42900000 - 43200000.

The tutorial will run the alignment pipeline on 2 of the individuals (HG00096, HG00100). The fastqs used for this step are reduced to reads that fall into our target region.

The tutorial will then used previously aligned/mapped reads for the full 60 individuals to generate a list of polymorphic sites and estimate accurate genotypes at each of these sites.

The example dataset we'll be using is available at: ftp://share.sph.umich.edu/gotcloud/gotcloudExample.tgz

# Change directory to where you want to install the Tutorial data
# Download the dataset tar from the ftp site
# Extract the tar

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloudExample_latest.tgz # Download
tar xvf gotcloudExample_latest.tgz # Extracts into gotcloudExample/

== STEP 2 : Run GotCloud Alignment Pipeline ==
The first step in processing next generation sequence data is mapping the reads to the reference genome, generating per sample BAM files.

The alignment pipeline has multiple built-in steps to generate BAMs:
# Align the fastqs to the reference genome
#* handles both single & paired end
# Merge the results from multiple fastqs into 1 file per sample
# Mark Duplicate Reads are marked
# Recalibrate Base Qualities

This processing results in 1 BAM file per sample.

The alignment pipeline also includes Quality Control (QC) steps:
# Visualization of various quality measures (QPLOT)
# Screen for sample contamination & swap (VerifyBamID)

Run the alignment pipeline (the example aligns 2 samples) :
~/gotcloud/gotcloud align --conf ~/gotcloudExample/[[#Alignment Configuration File|GBR2align.conf]] --outdir [[#Alignment Output Directory|~/gotcloudTutorialOut]] --baseprefix ~/gotcloudExample

Upon successful completion of the alignment pipeline (about 1-3 minutes), you will see the following message:
Processing finished in nn secs with no errors reported

The final BAM files produced by the alignment pipeline are:
ls ~/gotcloudTutorialOut/bams
In this directory you will see:
* BAM (.bam) files - 1 per sample
** HG00096.recal.bam
** HG00100.recal.bam
* BAM index files (.bai) – 1 per sample
** HG00096.recal.bam.bai
** HG00100.recal.bam.bai
* BAM checksum files (.md5) – 1 per sample
** HG00096.recal.bam.md5
** HG00100.recal.bam.md5
* Indicator files that the step completed successfully:
** HG00096.recal.bam.done
** HG00100.recal.bam.done

The Quality Control (QC) files are:
ls ~/gotcloudTutorialOut/QCFiles
In this directory you will see:
* VerifyBamID output files:
** HG00096.genoCheck.depthRG
** HG00096.genoCheck.depthSM
** HG00096.genoCheck.selfRG
** HG00096.genoCheck.selfSM
** HG00100.genoCheck.depthRG
** HG00100.genoCheck.depthSM
** HG00100.genoCheck.selfRG
** HG00100.genoCheck.selfSM

* VerifyBamID step completion files – 1 per sample
** HG00096.genoCheck.done
** HG00100.genoCheck.done

* QPLOT output files
** HG00096.qplot.R
** HG00096.qplot.stats
** HG00100.qplot.R
** HG00100.qplot.stats

* QPLOT step completion files – 1 per sample
** HG00096.qplot.done
** HG00100.qplot.done

For information on the VerifyBamID output, see: [[Understanding VerifyBamID output]]

For information on the QPLOT output, see: [[Understanding QPLOT output]]

== STEP 3 : Run GotCloud Variant Calling Pipeline ==
The next step is to analyze BAM files by calling SNPs and generating a VCF file containing the variant calls.

The variant calling pipeline has multiple built-in steps to generate BAMs:
# Filter out reads with low mapping quality
# Per Base Alignment Quality Adjustment (BAQ)
# Resolve overlapping paired end reads
# Generate genotype likelihood files
# Perform variant calling
# Extract features from variant sites
# Perform variant filtering

To speed variant calling, each chromosome is broken up into smaller regions which are processed separately. While initially split by sample, the per sample data gets merged and is processed together for each region. These regions are later merged to result in a single Variant Call File (VCF) per chromosome. For the tutorial all of the data falls within a single region.

Run the variant calling pipeline:
~/gotcloud/gotcloud snpcall --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --region 20:42900000-43200000 --baseprefix ~/gotcloudExample

Upon successful completion of the variant calling pipeline (about 3-4 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

On SNP Call success, the VCF files of interest are:
ls ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered*

This gives you the following files:
* '''chr20.filtered.vcf.gz ''' - vcf for whole chromosome after it has been run through hardfilters and SVM filters and marked with PASS/FAIL including per sample genotypes
* chr20.filtered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
* chr20.filtered.sites.vcf.norm.log - log file
* chr20.filtered.sites.vcf.summary - summary of filters applied
* chr20.filtered.vcf.gz.OK - indicator that the filtering completed successfully
* chr20.filtered.vcf.gz.tbi - index file for the vcf file

Also in the ~/gotcloudTutorialOut/vcfs/chr20 directory are intermediate files:
* the whole chromosome variant calls prior to any filtering:
** chr20.merged.sites.vcf - without per sample genotypes
** chr20.merged.stats.vcf
** chr20.merged.vcf - including per sample genotypes
** chr20.merged.vcf.OK - indicator that the step completed successfully
* the hardfiltered (pre-svm filtered) variant calls:
** chr20.filtered.vcf.gz - vcf for whole chromosome after it has been run through hard filters
** chr20.hardfiltered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
** chr20.hardfiltered.sites.vcf.log - log file
** chr20.hardfiltered.sites.vcf.summary - summary of filters applied
** chr20.hardfiltered.vcf.gz.OK - indicator that the filtering completed successfully
** chr20.hardfiltered.vcf.gz.tbi - index file for the vcf file
* 40000001.45000000 subdirectory contains the data for just that region.

The ~/gotcloudTutorialOut/split/chr20 folder contains a VCF with just the sites that pass the filters.
ls ~/gotcloudTutorialOut/split/chr20/
* '''chr20.filtered.PASS.vcf.gz ''' – vcf of just sites that pass all filters
* chr20.filtered.PASS.split.1.vcf.gz - intermediate file
* chr20.filtered.PASS.split.err - log file
* chr20.filtered.PASS.split.vcflist - list of intermediate files
* subset.OK

In addition to the vcfs subdirectory, there are additional intermediate files/directories:
* glfs – holds genotype likelihood format [[GLF]] files split by chromosome, region, and sample
* pvcfs – holds intermediate vcf files split by chromosome and region

Note: the tutorial does not produce a target directory, but if you run with targeted data, you may see that.

== STEP 4 : Run GotCloud Genotype Refinement Pipeline ==
The next step is to perform genotype refinement using linkage disequilibrium information using [http://faculty.washington.edu/browning/beagle/beagle.html Beagle] & [[ThunderVCF]].

Run the LD-aware genotype refinement pipeline:
~/gotcloud/gotcloud ldrefine --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --baseprefix ~/gotcloudExample

Upon successful completion of this pipeline (about 10 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

The output from the beagle step of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz.tbi

The output from the thunderVcf (final step) of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz.tbi

== STEP 5 : Run GotCloud Association Analysis Pipeline (EPACTS) ==

We will assume that the EPACTS are installed in the following directory
setenv EPACTS /path/to/epacts
(If you need to install EPACTS, please refer to the documentation at [[EPACTS#Installation_Details]])

$EPACTS/epacts single --vcf ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered.vcf.gz --ped ~/gotcloudExample/test.GBR60.ped \\
--out ~/gotcloudTutorialOut/epacts --test q.linear --run 1 --top 1 --chr 20

Upon successful run, you will see files starting with ~/gotcloudTutorialOut/epacts
ls ~/gotcloudTutorialOut/epacts*

To see the top associated variants, you can run
less ~/gotcloudTutorialOut/epacts.epacts.top5000

To see the locus-zoom like plot, you can type the following command (assuming GNU gnuplot 4.2 or higher version was installed)
xpdf ~/gotcloudTutorialOut/epacs.zoom.20.42987877.pdf

Click [[Media:EPACTS TEST.zoom.20.42987877.pdf | Exampe LocusZoom PDF]] to see the expected output pdf

= Frequently Asked Questions (FAQs) =

'''I ran the tutorial example successfully, how can I run it with my real sequence data?'''

Congratulations for your successful run of your [[GotCloud]] Tutorial. Please see [[#Tutorial Inputs]] section to prepare your own input files for your sequence data. You will need to specify the FASTQ files associated with its sample names as explained. In addition, you will need to download the full reference and resource file across whole genome (the Tutorial contains only chr20 portion to make it compact) See [[#Alignment Configuration File]] section for the detailed information. Also, please refer to the original documentation of [[GotCloud]] for more detailed guide on installation beyond the scope of the tutorial.

= Input Files for GotCloud Tutorial =

This section describes the input files needed for the GotCloud tutorial. You don't need to know this detail to run the tutorial, but if you're interested in understanding the structure of GotCloud pipeline and run with your own sample, this would be a good starting point

== Alignment Pipeline ==
=== List of Input Files needed for Alignment ===
The command-line inputs to the tutorial alignment pipeline are:
# [[#Alignment Configuration File|Configuration File (--conf)]]
#* Specifies the configuration file to use when running
# [[#Alignment Output Directory|Output Directory (--outdir)]]
#* Directory where the output should be placed.

Additional information required to run the alignment pipeline:
# [[#Alignment FASTQ Index File|Index file of FASTQs]]
# [[#Alignment Reference Files|Alignment Reference Files]]
For the tutorial, these values are specified in the configuration file.

=== Alignment Configuration File ===
The configuration file contains KEY = VALUE settings that override defaults and set specific values for the given run.

<pre>
INDEX_FILE = GBR60fastq.index
############
# References
REF_DIR = chr20Ref
AS = NCBI37
FA_REF = $(REF_DIR)/human_g1k_v37_chr20.fa
DBSNP_VCF = $(REF_DIR)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_DIR)/hapmap_3.3.b37.sites.chr20.vcf.gz
</pre>

This configuration file sets:
* [[#Alignment FASTQ Index File|INDEX_FILE]] - file containing the fastqs to be processed as well as the read group information for these fastqs.
* Reference Information: see [[#Alignment Reference Files|Alignment Reference Files]] for more information
** AS - assembly value to put in the BAM

The index file and chromosome 20 references used in this tutorial are included with the example data under the ~/gotcloudExample directory. The tutorial uses chromosome 20 only references in order to speed the processing time.

When running with your own data, you will need to update the:
* Index File to contain the information for your own fastq files
** See [[#Alignment FASTQ Index File|Alignment FASTQ Index File]] for more information on the contents of the index file.
* The reference files to be whole genome references
** If you are just running chromosome 20, you can use the tutorial references
** Whole genome reference files can be downloaded from [[GotCloudReference]].
** See [[#Alignment Reference Files|Alignment Reference Files]] for more information.

Note: It is recommended that you use absolute paths (full path names, like “/home/mktrost/gotcloudReference” rather than just “gotcloudReference”). This example does not use absolute paths in order to be flexible to where the data is installed, but using relative paths requires it to be run from the correct directory.

=== Reference Files ===

Reference files are required for running both the alignment and variant calling pipelines.

The configuration keys for setting these are:
* FA_REF - Genome sequence reference files (needed for both pipelines)
* DBSNP_VCF – DBSNP site VCF file (needed for both pipelines)
* HM3_VCF - HAPMAP site VCF file (needed for both pipelines)
* INDEL_PREFIX - Indel sites file (need for variant calling pipeline)

The tutorial configuration file is setup to point to the required chromosome 20 reference files which are included with the tutorial example data in ~/gotcloudExample/chr20Ref/.

If you are running more than just chromosome 20, you will need whole genome reference files which can be downloaded from [[GotCloudReference]].

The configuration settings for these files are setup in the default configuration so do not need to be specified. You just need to set REF_DIR in your configuration file to the path where you installed your reference files.

To learn more about the reference files that are required, see [[GotCloud: Reference Files]].

=== Alignment Output Directory ===
This setting tells the alignment pipeline where to write the output and intermediate files.

The output directory will be created if it doesn't already exist and will contain the following Directories/files:
* bams - directory containing the final bams/bai files
** HG00096.OK - indicates that this sample completed alignment processing
** HG00100.OK - indicates that this sample completed alignment processing
* failLogs - directory containing logs from steps that failed
** this directory is only created if an error is detected
* Makefiles - directory containing the makefiles with commands for processing each sample
** biopipe_HG00096.Makefile – commands for processing sample HG00096
** biopipe_HG00100.Makefile – commands for processing sample HG00100
** biopipe_HG00096.Makefile.log – log file from running the associated Makefile
** biopipe_HG00100.Makefile.log – log file from running the associated Makefile
* QCFiles - directory containing the QC Results
**
* tmp - directory containing temporary alignment files
** bwa.sai.t – contains temporary files for the 1st step of bwa that generates sai files from fastq files
*** fastq
**** *.done – indicator files that the step to generate the file completed
**** HG0096
** alignment.bwa – contains temporary files for the 2nd step of bwa that generates BAM files
*** fastq
**** *.done – indicator files that the step to generate the file completed
** alignment.pol - contains temporary files for the polish bam step that cleans up BAM files
** alignment.dedup – contains temporary files for the deduping step
*** *.done – indicator files that the step to generate the file completed
*** *.metrics – metrics files for the deduping step
** alignment.recal – contains temporary files for the deduping step
*** *.log – contains information about the recalibration step
*** *.qemp – contains the recalibration tables used for recalibrating each BAM file

===Alignment FASTQ Index File===

There are four fastq files in {ROOT_DIR}/test/align/fastq/Sample_1 and four fastq files in {ROOT_DIR}/test/align/fastq/Sample_2, both in paired-end format. Normally, we would need to build an index file for these files. Conveniently, an index file (indexFile.txt) already exists for the automatic test samples. It can be found in {ROOT_DIR}/test/align/, and contains the following information in tab-delimited format:

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/Sample_1/File1_R1.fastq.gz fastq/Sample_1/File1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/Sample_1/File2_R1.fastq.gz fastq/Sample_1/File2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/Sample_2/File1_R1.fastq.gz fastq/Sample_2/File1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/Sample_2/File2_R1.fastq.gz fastq/Sample_2/File2_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA

If you are in the {ROOT_DIR}/test/align directory, you can use this file as-is. If you prefer, you can create a new index file and change the MERGE_NAME, RGID, SAMPLE, LIBRARY, CENTER, or PLATFORM values. It is recommended that you do not modify existing files in {ROOT_DIR}/test/align.

If you want to run this example from a different directory, make sure the FASTQ1 and FASTQ2 paths are correct. That is, each of the FASTQ1 and FASTQ2 entry in the index file should look like the following:

{ROOT_DIR}/test/align/fastq/Sample_1/File1_R1.fastq.gz

Alternately, if you want to run this example from a different directory, but do not want to edit the index file, you can create a relative path to the test fastq files so their path agrees with that listed in the index file:

ln -s {ROOT_DIR}/test/align/fastq fastq

This will create a symbolic link to the test fastq directory from your current directory.

(More information about: [[Mapping_Pipeline#Sequence_Index_File|the index file]].)

==Analyzing a Sample==

Using UMAKE, you can analyze BAM files by calling SNPs, and generate a VCF file containing the results. Once again, we can analyze BAM files used in the automatic test. For this example, we have 60 BAM files, which can be found in {ROOT_DIR}/test/umake/bams. These contain sequence information for a targeted region in chromosome 20.

In addition to the BAM files, you will need three files to run UMAKE: an index file, a configuration file, and a bed file (needed to analyze BAM files from targeted/exome sequencing).

===Index file===

First, you need a list of all the BAM files to be analyzed. Conveniently, the a test index file (umake_test.index) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

NA12272 ALL bams/NA12272.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
NA12004 ALL bams/NA12004.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
...
NA12874 ALL bams/NA12874.mapped.LS454.ssaha2.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam

You can use this file directly if you change your current directory to {ROOT_DIR}/test/umake/.

Alternately, if you want to copy and use this index file to a different directory, you can create a symbolic link to the bams folder as follows:

ln -s {ROOT_DIR}/test/umake/bams bams

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Index_File|the index file]].)

===BED file===

This file contains a single line:

chr20 20000050 20300000

You can copy this to the current directory and use it as-is.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Targeted.2FExome_Sequencing_Settings|targeted/exome sequencing settings]].)

===Configuration file===

A configuration file (umake_test.conf) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

CHRS = 20
BAM_INDEX = GBR60bam.index
############
# References
REF_ROOT = chr20Ref
#
REF = $(REF_ROOT)/human_g1k_v37_chr20.fa
INDEL_PREFIX = $(REF_ROOT)/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_VCF = $(REF_ROOT)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_ROOT)/hapmap_3.3.b37.sites.chr20.vcf.gz

CHRS = 20
TEST_ROOT = $(UMAKE_ROOT)/test/umake
BAM_INDEX = $(TEST_ROOT)/umake_test.index
OUT_PREFIX = umake_test
REF_ROOT = $(TEST_ROOT)/ref
#
REF = $(REF_ROOT)/karma.ref/human.g1k.v37.chr20.fa
INDEL_PREFIX = $(REF_ROOT)/indels/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_PREFIX = $(REF_ROOT)/dbSNP/dbsnp_135_b37.rod
HM3_PREFIX = $(REF_ROOT)/HapMap3/hapmap3_r3_b37_fwd.consensus.qc.poly
#
RUN_INDEX = TRUE # create BAM index file
RUN_PILEUP = TRUE # create GLF file from BAM
RUN_GLFMULTIPLES = TRUE # create unfiltered SNP calls
RUN_VCFPILEUP = TRUE # create PVCF files using vcfPileup and run infoCollector
RUN_FILTER = TRUE # filter SNPs using vcfCooker
RUN_SPLIT = TRUE # split SNPs into chunks for genotype refinement
RUN_BEAGLE = FALSE # BEAGLE - MUST SET AFTER FINISHING PREVIOUS STEPS
RUN_SUBSET = FALSE # SUBSET FOR THUNDER - MAY BE SET WITH BEAGLE STEP TOGETHER
RUN_THUNDER = FALSE # THUNDER - MUST SET AFTER FINISHING PREVIOUS STEPS
###############################################################################
WRITE_TARGET_LOCI = TRUE # FOR TARGETED SEQUENCING ONLY -- Write loci file when performing pileup
UNIFORM_TARGET_BED = $(TEST_ROOT)/umake_test.bed # Targeted sequencing : When all individuals has the same target. Otherwise, comment it out
OFFSET_OFF_TARGET = 50 # Extend target by given # of bases
MULTIPLE_TARGET_MAP = # Target per individual : Each line contains [SM_ID] [TARGET_BED]
TARGET_DIR = target # Directory to store target information
SAMTOOLS_VIEW_TARGET_ONLY = TRUE # When performing samtools view, exclude off-target regions (may make command line too long)

If you are running this from a different directory, you will want to change some of the lines as follows:

BAM_INDEX = {CURRENT_DIR}/umake_test.index
UNIFORM_TARGET_BED = {CURRENT_DIR}/umake_test.bed

where {CURRENT_DIR} is the absolute path to the directory that contains the index and bed files.

An additional option can be added in the configuration file:

OUT_DIR = {OUT_DIR}

where {OUT_DIR} is the name of directory in which you want the output to be stored. If you do not specify this in the configuration file, you will need to add an extra parameter when you run UMAKE in the next step.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Configuration_File|the configuration file]], [[Variant_Calling_Pipeline_(UMAKE)#Reference_Files|reference files]].)

==Further Information==

[[Mapping_Pipeline|Mapping (Alignment) Pipeline]]

[[Variant_Calling_Pipeline_(UMAKE)|Variant Calling Pipeline (UMAKE)]]

Tutorial: GotCloud

2013-06-16T21:35:43Z

Ben Lerch: /* STEP 2 : Run GotCloud Alignment Pipeline */

= GotCloud Tutorial =
In this tutorial, we illustrate some of the essential steps in the analysis of next generation sequence data.

For a background on GotCloud and Sequence Analysis Pipelines, see [[GotCloud]]

While GotCloud can run on a cluster of machines or instances, this tutorial is just a small test that just runs on the machine the commands are run on.

GotCloud and this basic tutorial were presented at the [http://ibg.colorado.edu/dokuwiki/doku.php?id=workshop:2013:announcement 2013 IBG Workshop]. It was presented in two sessions. On Wednesday an overview was presented with steps for running the tutorial data: [[Media:IBG2013GotCloud.pdf|IBG2013GotCloud.pdf]]. On Friday more detail on the input files and what goes into generating the input files was presented: [[Media:GotCloudIBGWorkshop2013Friday.pdf|GotCloudIBGWorkshop2013Friday.pdf]].

'''This tutorial is in the process of being updated for gotcloud version 1.06 (April 17. 2013).'''

== STEP 1 : Setup GotCloud ==

[[GotCloud]] has been developed and tested on Linux Ubuntu 12.10 and 12.04.2 LTS but has not been tested on other Linux operating systems. It is not available for Windows. If you do not have your own set of machines to run on, GotCloud is also available for Ubuntu running on the Amazon Elastic Compute Cloud, see [[Amazon_Snapshot]] for more information.

We will use 3 different directories for this tutorial:
# path to the directory where gotcloud is installed, default ~/gotcloud/
# path to the directory where the example data is installed, default ~/gotcloudExample
# path to your output directory, default ~/gotcloudTutorialOut/

If the directories specified above do not reflect the directories you would like to use, replace their occurrances in the instructions below with the appropriate paths.

=== Step 1a: Install GotCloud ===
In order to run this tutorial, you need to make sure you have GotCloud installed on your system.

If you have root and would like to install gotcloud on your system, follow: [[GotCloud#Install_GotCloud_Software| root access installation instructions]]

Otherwise, you can install it in your own directory:
# Change to the directory where you want gotcloud/ installed
# Download the gotcloud tar from the ftp site.
# Extract the tar
# Build (compile) the source
#* Note: as the source builds, many messages will scroll through your terminal. You may even see some warnings. These messages are normal and expected. As long as the build does not end with an error, you have successfully built the source.

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloud_latest.tgz # Download
tar xf gotcloud_latest.tgz # Extracts into gotcloud/
cd ~/gotcloud/src; make # Build source

GotCloud requires the following tools to be installed.
You can run ~/gotcloud/scripts/check_requirements.sh
...TBD – put in required programs/tools.
* java (java-common default-jre on ubuntu)
* make (make on ubuntu)
* libssl (libssl0.9.8 on ubuntu)
* gcc 4.4 or newer

=== Step 1b: Install Example Dataset ===
Our dataset consists of 60 individuals from Great Britain (GBR) sequenced by the 1000 Genomes Project. These individuals have been sequenced to an average depth of about 4x.

To conserve time and disk-space, our analysis will focus on a small region on chromosome 20, 42900000 - 43200000.

The tutorial will run the alignment pipeline on 2 of the individuals (HG00096, HG00100). The fastqs used for this step are reduced to reads that fall into our target region.

The tutorial will then used previously aligned/mapped reads for the full 60 individuals to generate a list of polymorphic sites and estimate accurate genotypes at each of these sites.

The example dataset we'll be using is available at: ftp://share.sph.umich.edu/gotcloud/gotcloudExample.tgz

# Change directory to where you want to install the Tutorial data
# Download the dataset tar from the ftp site
# Extract the tar

cd ~
wget ftp://share.sph.umich.edu/gotcloud/gotcloudExample_latest.tgz # Download
tar xvf gotcloudExample_latest.tgz # Extracts into gotcloudExample/

== STEP 2 : Run GotCloud Alignment Pipeline ==
The first step in processing next generation sequence data is mapping the reads to the reference genome, generating per sample BAM files.

The alignment pipeline has multiple built-in steps to generate BAMs:
# Align the fastqs to the reference genome
#* handles both single & paired end
# Merge the results from multiple fastqs into 1 file per sample
# Mark Duplicate Reads are marked
# Recalibrate Base Qualities

This processing results in 1 BAM file per sample.

The alignment pipeline also includes Quality Control (QC) steps:
# Visualization of various quality measures (QPLOT)
# Screen for sample contamination & swap (VerifyBamID)

Run the alignment pipeline (the example aligns 2 samples) :
~/gotcloud/gotcloud align --conf ~/gotcloudExample/[[#Alignment Configuration File|GBR2align.conf]] --outdir [[#Alignment Output Directory|~/gotcloudTutorialOut]] --baseprefix ~/gotcloudExample

Upon successful completion of the alignment pipeline (about 1-3 minutes), you will see the following message:
Processing finished in nn secs with no errors reported

The final BAM files produced by the alignment pipeline are:
ls ~/gotcloudTutorialOut/bams
In this directory you will see:
* BAM (.bam) files - 1 per sample
** HG00096.recal.bam
** HG00100.recal.bam
* BAM index files (.bai) – 1 per sample
** HG00096.recal.bam.bai
** HG00100.recal.bam.bai
* BAM checksum files (.md5) – 1 per sample
** HG00096.recal.bam.md5
** HG00100.recal.bam.md5
* Indicator files that the step completed successfully:
** HG00096.recal.bam.done
** HG00100.recal.bam.done

The Quality Control (QC) files are:
ls ~/gotcloudTutorialOut/QCFiles
In this directory you will see:
* VerifyBamID output files:
** HG00096.genoCheck.depthRG
** HG00096.genoCheck.depthSM
** HG00096.genoCheck.selfRG
** HG00096.genoCheck.selfSM
** HG00100.genoCheck.depthRG
** HG00100.genoCheck.depthSM
** HG00100.genoCheck.selfRG
** HG00100.genoCheck.selfSM

* VerifyBamID step completion files – 1 per sample
** HG00096.genoCheck.done
** HG00100.genoCheck.done

* QPLOT output files
** HG00096.qplot.R
** HG00096.qplot.stats
** HG00100.qplot.R
** HG00100.qplot.stats

* QPLOT step completion files – 1 per sample
** HG00096.qplot.done
** HG00100.qplot.done

For information on the VerifyBamID output, see: [[Understanding VerifyBamID output]]

For information on the QPLOT output, see: [[Understanding QPLOT output]]

== STEP 3 : Run GotCloud Variant Calling Pipeline ==
The next step is to analyze BAM files by calling SNPs and generating a VCF file containing the variant calls.

The variant calling pipeline has multiple built-in steps to generate BAMs:
# Filter out reads with low mapping quality
# Per Base Alignment Quality Adjustment (BAQ)
# Resolve overlapping paired end reads
# Generate genotype likelihood files
# Perform variant calling
# Extract features from variant sites
# Perform variant filtering

To speed variant calling, each chromosome is broken up into smaller regions which are processed separately. While initially split by sample, the per sample data gets merged and is processed together for each region. These regions are later merged to result in a single Variant Call File (VCF) per chromosome. For the tutorial all of the data falls within a single region.

Run the variant calling pipeline:
~/gotcloud/gotcloud snpcall --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --region 20:42900000-43200000 --baseprefix ~/gotcloudExample

Upon successful completion of the variant calling pipeline (about 3-4 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

On SNP Call success, the VCF files of interest are:
ls ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered*

This gives you the following files:
* '''chr20.filtered.vcf.gz ''' - vcf for whole chromosome after it has been run through hardfilters and SVM filters and marked with PASS/FAIL including per sample genotypes
* chr20.filtered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
* chr20.filtered.sites.vcf.norm.log - log file
* chr20.filtered.sites.vcf.summary - summary of filters applied
* chr20.filtered.vcf.gz.OK - indicator that the filtering completed successfully
* chr20.filtered.vcf.gz.tbi - index file for the vcf file

Also in the ~/gotcloudTutorialOut/vcfs/chr20 directory are intermediate files:
* the whole chromosome variant calls prior to any filtering:
** chr20.merged.sites.vcf - without per sample genotypes
** chr20.merged.stats.vcf
** chr20.merged.vcf - including per sample genotypes
** chr20.merged.vcf.OK - indicator that the step completed successfully
* the hardfiltered (pre-svm filtered) variant calls:
** chr20.filtered.vcf.gz - vcf for whole chromosome after it has been run through hard filters
** chr20.hardfiltered.sites.vcf - vcf for whole chromosome after it has been run through filters and marked with PASS/FAIL without the per sample genotypes
** chr20.hardfiltered.sites.vcf.log - log file
** chr20.hardfiltered.sites.vcf.summary - summary of filters applied
** chr20.hardfiltered.vcf.gz.OK - indicator that the filtering completed successfully
** chr20.hardfiltered.vcf.gz.tbi - index file for the vcf file
* 40000001.45000000 subdirectory contains the data for just that region.

The ~/gotcloudTutorialOut/split/chr20 folder contains a VCF with just the sites that pass the filters.
ls ~/gotcloudTutorialOut/split/chr20/
* '''chr20.filtered.PASS.vcf.gz ''' – vcf of just sites that pass all filters
* chr20.filtered.PASS.split.1.vcf.gz - intermediate file
* chr20.filtered.PASS.split.err - log file
* chr20.filtered.PASS.split.vcflist - list of intermediate files
* subset.OK

In addition to the vcfs subdirectory, there are additional intermediate files/directories:
* glfs – holds genotype likelihood format [[GLF]] files split by chromosome, region, and sample
* pvcfs – holds intermediate vcf files split by chromosome and region

Note: the tutorial does not produce a target directory, but if you run with targeted data, you may see that.

== STEP 4 : Run GotCloud Genotype Refinement Pipeline ==
The next step is to perform genotype refinement using linkage disequilibrium information using [http://faculty.washington.edu/browning/beagle/beagle.html Beagle] & [[ThunderVCF]].

Run the LD-aware genotype refinement pipeline:
~/gotcloud/gotcloud ldrefine --conf ~/gotcloudExample/[[GBR60vc.conf]] --outdir ~/gotcloudTutorialOut --numjobs 2 --baseprefix ~/gotcloudExample

Upon successful completion of this pipeline (about 10 minutes), you will see the following message:
Commands finished in nnn secs with no errors reported

The output from the beagle step of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz ~/gotcloudTutorialOut/beagle/chr20/chr20.filtered.PASS.beagled.vcf.gz.tbi

The output from the thunderVcf (final step) of the genotype refinement pipeline is found in:
ls ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz ~/gotcloudTutorialOut/thunder/chr20/GBR/chr20.filtered.PASS.beagled.GBR.thunder.vcf.gz.tbi

== STEP 5 : Run GotCloud Association Analysis Pipeline (EPACTS) ==

We will assume that the EPACTS are installed in the following directory
setenv EPACTS /path/to/epacts
(If you need to install EPACTS, please refer to the documentation at [[EPACTS#Installation_Details]])

$EPACTS/epacts single --vcf ~/gotcloudTutorialOut/vcfs/chr20/chr20.filtered.vcf.gz --ped ~/gotcloudExample/test.GBR60.ped \\
--out ~/gotcloudTutorialOut/epacts --test q.linear --run 1 --top 1 --chr 20

Upon successful run, you will see files starting with ~/gotcloudTutorialOut/epacts
ls ~/gotcloudTutorialOut/epacts*

To see the top associated variants, you can run
less ~/gotcloudTutorialOut/epacts.epacts.top5000

To see the locus-zoom like plot, you can type the following command (assuming GNU gnuplot 4.2 or higher version was installed)
xpdf ~/gotcloudTutorialOut/epacs.zoom.20.42987877.pdf

Click [[Media:EPACTS TEST.zoom.20.42987877.pdf | Exampe LocusZoom PDF]] to see the expected output pdf

= Frequently Asked Questions (FAQs) =

'''I ran the tutorial example successfully, how can I run it with my real sequence data?'''

Congratulations for your successful run of your [[GotCloud]] Tutorial. Please see [[#Tutorial Inputs]] section to prepare your own input files for your sequence data. You will need to specify the FASTQ files associated with its sample names as explained. In addition, you will need to download the full reference and resource file across whole genome (the Tutorial contains only chr20 portion to make it compact) See [[#Alignment Configuration File]] section for the detailed information. Also, please refer to the original documentation of [[GotCloud]] for more detailed guide on installation beyond the scope of tutorial

= Input Files for GotCloud Tutorial =

This section describes the input files needed for the GotCloud tutorial. You don't need to know this detail to run the tutorial, but if you're interested in understanding the structure of GotCloud pipeline and run with your own sample, this would be a good starting point

== Alignment Pipeline ==
=== List of Input Files needed for Alignment ===
The command-line inputs to the tutorial alignment pipeline are:
# [[#Alignment Configuration File|Configuration File (--conf)]]
#* Specifies the configuration file to use when running
# [[#Alignment Output Directory|Output Directory (--outdir)]]
#* Directory where the output should be placed.

Additional information required to run the alignment pipeline:
# [[#Alignment FASTQ Index File|Index file of FASTQs]]
# [[#Alignment Reference Files|Alignment Reference Files]]
For the tutorial, these values are specified in the configuration file.

=== Alignment Configuration File ===
The configuration file contains KEY = VALUE settings that override defaults and set specific values for the given run.

<pre>
INDEX_FILE = GBR60fastq.index
############
# References
REF_DIR = chr20Ref
AS = NCBI37
FA_REF = $(REF_DIR)/human_g1k_v37_chr20.fa
DBSNP_VCF = $(REF_DIR)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_DIR)/hapmap_3.3.b37.sites.chr20.vcf.gz
</pre>

This configuration file sets:
* [[#Alignment FASTQ Index File|INDEX_FILE]] - file containing the fastqs to be processed as well as the read group information for these fastqs.
* Reference Information: see [[#Alignment Reference Files|Alignment Reference Files]] for more information
** AS - assembly value to put in the BAM

The index file and chromosome 20 references used in this tutorial are included with the example data under the ~/gotcloudExample directory. The tutorial uses chromosome 20 only references in order to speed the processing time.

When running with your own data, you will need to update the:
* Index File to contain the information for your own fastq files
** See [[#Alignment FASTQ Index File|Alignment FASTQ Index File]] for more information on the contents of the index file.
* The reference files to be whole genome references
** If you are just running chromosome 20, you can use the tutorial references
** Whole genome reference files can be downloaded from [[GotCloudReference]].
** See [[#Alignment Reference Files|Alignment Reference Files]] for more information.

Note: It is recommended that you use absolute paths (full path names, like “/home/mktrost/gotcloudReference” rather than just “gotcloudReference”). This example does not use absolute paths in order to be flexible to where the data is installed, but using relative paths requires it to be run from the correct directory.

=== Reference Files ===

Reference files are required for running both the alignment and variant calling pipelines.

The configuration keys for setting these are:
* FA_REF - Genome sequence reference files (needed for both pipelines)
* DBSNP_VCF – DBSNP site VCF file (needed for both pipelines)
* HM3_VCF - HAPMAP site VCF file (needed for both pipelines)
* INDEL_PREFIX - Indel sites file (need for variant calling pipeline)

The tutorial configuration file is setup to point to the required chromosome 20 reference files which are included with the tutorial example data in ~/gotcloudExample/chr20Ref/.

If you are running more than just chromosome 20, you will need whole genome reference files which can be downloaded from [[GotCloudReference]].

The configuration settings for these files are setup in the default configuration so do not need to be specified. You just need to set REF_DIR in your configuration file to the path where you installed your reference files.

To learn more about the reference files that are required, see [[GotCloud: Reference Files]].

=== Alignment Output Directory ===
This setting tells the alignment pipeline where to write the output and intermediate files.

The output directory will be created if it doesn't already exist and will contain the following Directories/files:
* bams - directory containing the final bams/bai files
** HG00096.OK - indicates that this sample completed alignment processing
** HG00100.OK - indicates that this sample completed alignment processing
* failLogs - directory containing logs from steps that failed
** this directory is only created if an error is detected
* Makefiles - directory containing the makefiles with commands for processing each sample
** biopipe_HG00096.Makefile – commands for processing sample HG00096
** biopipe_HG00100.Makefile – commands for processing sample HG00100
** biopipe_HG00096.Makefile.log – log file from running the associated Makefile
** biopipe_HG00100.Makefile.log – log file from running the associated Makefile
* QCFiles - directory containing the QC Results
**
* tmp - directory containing temporary alignment files
** bwa.sai.t – contains temporary files for the 1st step of bwa that generates sai files from fastq files
*** fastq
**** *.done – indicator files that the step to generate the file completed
**** HG0096
** alignment.bwa – contains temporary files for the 2nd step of bwa that generates BAM files
*** fastq
**** *.done – indicator files that the step to generate the file completed
** alignment.pol - contains temporary files for the polish bam step that cleans up BAM files
** alignment.dedup – contains temporary files for the deduping step
*** *.done – indicator files that the step to generate the file completed
*** *.metrics – metrics files for the deduping step
** alignment.recal – contains temporary files for the deduping step
*** *.log – contains information about the recalibration step
*** *.qemp – contains the recalibration tables used for recalibrating each BAM file

===Alignment FASTQ Index File===

There are four fastq files in {ROOT_DIR}/test/align/fastq/Sample_1 and four fastq files in {ROOT_DIR}/test/align/fastq/Sample_2, both in paired-end format. Normally, we would need to build an index file for these files. Conveniently, an index file (indexFile.txt) already exists for the automatic test samples. It can be found in {ROOT_DIR}/test/align/, and contains the following information in tab-delimited format:

MERGE_NAME FASTQ1 FASTQ2 RGID SAMPLE LIBRARY CENTER PLATFORM
Sample1 fastq/Sample_1/File1_R1.fastq.gz fastq/Sample_1/File1_R2.fastq.gz RGID1 SampleID1 Lib1 UM ILLUMINA
Sample1 fastq/Sample_1/File2_R1.fastq.gz fastq/Sample_1/File2_R2.fastq.gz RGID1a SampleID1 Lib1 UM ILLUMINA
Sample2 fastq/Sample_2/File1_R1.fastq.gz fastq/Sample_2/File1_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA
Sample2 fastq/Sample_2/File2_R1.fastq.gz fastq/Sample_2/File2_R2.fastq.gz RGID2 SampleID2 Lib2 UM ILLUMINA

If you are in the {ROOT_DIR}/test/align directory, you can use this file as-is. If you prefer, you can create a new index file and change the MERGE_NAME, RGID, SAMPLE, LIBRARY, CENTER, or PLATFORM values. It is recommended that you do not modify existing files in {ROOT_DIR}/test/align.

If you want to run this example from a different directory, make sure the FASTQ1 and FASTQ2 paths are correct. That is, each of the FASTQ1 and FASTQ2 entry in the index file should look like the following:

{ROOT_DIR}/test/align/fastq/Sample_1/File1_R1.fastq.gz

Alternately, if you want to run this example from a different directory, but do not want to edit the index file, you can create a relative path to the test fastq files so their path agrees with that listed in the index file:

ln -s {ROOT_DIR}/test/align/fastq fastq

This will create a symbolic link to the test fastq directory from your current directory.

(More information about: [[Mapping_Pipeline#Sequence_Index_File|the index file]].)

==Analyzing a Sample==

Using UMAKE, you can analyze BAM files by calling SNPs, and generate a VCF file containing the results. Once again, we can analyze BAM files used in the automatic test. For this example, we have 60 BAM files, which can be found in {ROOT_DIR}/test/umake/bams. These contain sequence information for a targeted region in chromosome 20.

In addition to the BAM files, you will need three files to run UMAKE: an index file, a configuration file, and a bed file (needed to analyze BAM files from targeted/exome sequencing).

===Index file===

First, you need a list of all the BAM files to be analyzed. Conveniently, the a test index file (umake_test.index) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

NA12272 ALL bams/NA12272.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
NA12004 ALL bams/NA12004.mapped.ILLUMINA.bwa.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam
...
NA12874 ALL bams/NA12874.mapped.LS454.ssaha2.CEU.low_coverage.20101123.chrom20.20000001.20300000.bam

You can use this file directly if you change your current directory to {ROOT_DIR}/test/umake/.

Alternately, if you want to copy and use this index file to a different directory, you can create a symbolic link to the bams folder as follows:

ln -s {ROOT_DIR}/test/umake/bams bams

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Index_File|the index file]].)

===BED file===

This file contains a single line:

chr20 20000050 20300000

You can copy this to the current directory and use it as-is.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Targeted.2FExome_Sequencing_Settings|targeted/exome sequencing settings]].)

===Configuration file===

A configuration file (umake_test.conf) already exists in {ROOT_DIR}/test/umake/. It contains the following information:

CHRS = 20
BAM_INDEX = GBR60bam.index
############
# References
REF_ROOT = chr20Ref
#
REF = $(REF_ROOT)/human_g1k_v37_chr20.fa
INDEL_PREFIX = $(REF_ROOT)/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_VCF = $(REF_ROOT)/dbsnp135_chr20.vcf.gz
HM3_VCF = $(REF_ROOT)/hapmap_3.3.b37.sites.chr20.vcf.gz

CHRS = 20
TEST_ROOT = $(UMAKE_ROOT)/test/umake
BAM_INDEX = $(TEST_ROOT)/umake_test.index
OUT_PREFIX = umake_test
REF_ROOT = $(TEST_ROOT)/ref
#
REF = $(REF_ROOT)/karma.ref/human.g1k.v37.chr20.fa
INDEL_PREFIX = $(REF_ROOT)/indels/1kg.pilot_release.merged.indels.sites.hg19
DBSNP_PREFIX = $(REF_ROOT)/dbSNP/dbsnp_135_b37.rod
HM3_PREFIX = $(REF_ROOT)/HapMap3/hapmap3_r3_b37_fwd.consensus.qc.poly
#
RUN_INDEX = TRUE # create BAM index file
RUN_PILEUP = TRUE # create GLF file from BAM
RUN_GLFMULTIPLES = TRUE # create unfiltered SNP calls
RUN_VCFPILEUP = TRUE # create PVCF files using vcfPileup and run infoCollector
RUN_FILTER = TRUE # filter SNPs using vcfCooker
RUN_SPLIT = TRUE # split SNPs into chunks for genotype refinement
RUN_BEAGLE = FALSE # BEAGLE - MUST SET AFTER FINISHING PREVIOUS STEPS
RUN_SUBSET = FALSE # SUBSET FOR THUNDER - MAY BE SET WITH BEAGLE STEP TOGETHER
RUN_THUNDER = FALSE # THUNDER - MUST SET AFTER FINISHING PREVIOUS STEPS
###############################################################################
WRITE_TARGET_LOCI = TRUE # FOR TARGETED SEQUENCING ONLY -- Write loci file when performing pileup
UNIFORM_TARGET_BED = $(TEST_ROOT)/umake_test.bed # Targeted sequencing : When all individuals has the same target. Otherwise, comment it out
OFFSET_OFF_TARGET = 50 # Extend target by given # of bases
MULTIPLE_TARGET_MAP = # Target per individual : Each line contains [SM_ID] [TARGET_BED]
TARGET_DIR = target # Directory to store target information
SAMTOOLS_VIEW_TARGET_ONLY = TRUE # When performing samtools view, exclude off-target regions (may make command line too long)

If you are running this from a different directory, you will want to change some of the lines as follows:

BAM_INDEX = {CURRENT_DIR}/umake_test.index
UNIFORM_TARGET_BED = {CURRENT_DIR}/umake_test.bed

where {CURRENT_DIR} is the absolute path to the directory that contains the index and bed files.

An additional option can be added in the configuration file:

OUT_DIR = {OUT_DIR}

where {OUT_DIR} is the name of directory in which you want the output to be stored. If you do not specify this in the configuration file, you will need to add an extra parameter when you run UMAKE in the next step.

(More information about: [[Variant_Calling_Pipeline_(UMAKE)#Configuration_File|the configuration file]], [[Variant_Calling_Pipeline_(UMAKE)#Reference_Files|reference files]].)

==Further Information==

[[Mapping_Pipeline|Mapping (Alignment) Pipeline]]

[[Variant_Calling_Pipeline_(UMAKE)|Variant Calling Pipeline (UMAKE)]]