Test zCall Pull Request

.gitattributes

@@ -0,0 +1,2 @@
+*.png -crlf -diff
+*.pdf -crlf -diff

.gitignore

@@ -0,0 +1,7 @@
+# patterns for git version control to ignore
+*~
+*.pyc
+*.html
+*.aux
+*.dvi
+*.log

Makefile

@@ -0,0 +1,19 @@
+# makefile for zCall by Iain Bancarz, ib5@sanger.ac.uk
+
+PREFIX = 	PREFIX_DIRECTORY    # dummy value as default
+DEST =		$(PREFIX)/zCall
+SCRIPTS =	src/zcall
+ETC =       src/etc
+
+usage:
+	@echo -e "Usage: make install PREFIX=<destination directory>\nWill install to the zCall subdirectory of PREFIX.\nPREFIX must exist, zCall subdirectory will be created if necessary."
+
+install: $(PREFIX)
+	@echo -e "Installing scripts..."
+	install -d $(DEST) $(DEST)/zcall $(DEST)/etc $(DEST)/doc
+	@rm -f $(DEST)/zcall/*.pyc    # force recompiling of any .pyc files
+	install $(SCRIPTS)/*.py $(SCRIPTS)/*.r $(DEST)/zcall
+	install $(ETC)/*.ini $(DEST)/etc
+	@echo -e "Writing documentation..."
+	$(SCRIPTS)/createDocs.py --out  $(DEST)/doc
+	@echo -e "Installation complete. See $(DEST)/doc/zcall.html for class documentation."

README.md

@@ -0,0 +1,59 @@
+
+zCall: A Rare Variant Caller for Array-Based Genotyping
+=======================================================
+
+1. Overview
+-----------
+
+zCall is a variant caller specifically designed for calling rare single 
+nucleotide polymorphisms (SNPs) from array-based technology. This caller is 
+implemented as a post-processing step after a default calling algorithm has 
+been applied such as Illumina's GenCall algorithm. zCall uses the intensity 
+profile of the common allele homozygote cluster to define the location of the 
+other two genotype clusters. 
+
+The zCall code includes three prototype versions, and an extended version; 
+these are documented respectively in README_prototypes and README_extended.md. 
+The prototypes are run in a series of steps with some manual intervention; 
+extended zCall can be run with a single command and has other added 
+capabilities.
+
+2. Publication and downloads
+-----------------------------
+
+The paper describing zCall is:
+Goldstein JI, Crenshaw A, Carey J, Grant GB, Maguire J, Fromer M, 
+O'Dushlaine C, Moran JL, Chambert K, Stevens C; Swedish Schizophrenia 
+Consortium; ARRA Autism Sequencing Consortium, Sklar P, Hultman CM, Purcell S, 
+McCarroll SA, Sullivan PF, Daly MJ, Neale BM. zCall: a rare variant caller 
+for array-based genotyping: Genetics and population analysis. Bioinformatics. 
+2012 Oct 1;28(19):2543-2545. Epub 2012 Jul 27. PubMed PMID: 22843986.
+
+zCall is hosted on Github (https://github.com/jigold/zCall). Prototype versions 
+are available as .zip files, while the extended version is in the src 
+directory. The entire zCall repository can be cloned using Git or downloaded 
+from Github as a .zip file (approximately 0.8 MB). Extended zCall can be 
+installed from a download of the full repository, using the included Makefile.
+
+3. Disclaimer
+---------------
+
+The prototype and extended editions of zCall include code provided by Illumina. 
+The Illumina provided Code was provided as-is and with no warranty as to 
+performance and no warranty against it infringing any other party's 
+intellectual property rights.
+
+4. Contacts
+------------
+
+The original zCall method and prototype implementations were developed by 
+Jackie Goldstein et al.  For questions about prototype zCall or reporting 
+problems with the code, please send an email to Jackie Goldstein 
+(jigold@broadinstitute.org). For all other inquiries, please send an email to 
+both Ben Neale (bneale@broadinstitute.org) and Jackie Goldstein 
+(jigold@broadinstitute.org).
+
+Any queries concerning the extended version of zCall in the src directory 
+should be directed to Iain Bancarz (ib5@sanger.ac.uk).
+
+This document was written by Iain Bancarz.

README_extended.md

@@ -0,0 +1,132 @@
+
+zCall: Extended Version
+=======================
+
+1. Overview
+-----------
+
+This document covers the extended version of zCall, located in the `src` 
+directory on github. For other implementations, see Section 9 (History).
+
+2. Installation
+---------------
+
+Use GNU Make, with the target `install`. For example, `make install 
+/home/jsmith/foo` will create a directory `/home/jsmith/foo/zCall` and 
+install the relevant scripts and documentation. 
+
+The etc/config.ini file should be edited to reflect the user's local 
+computing environment. 
+
+Software prerequisites are:
+  * Python 2.7.x (or PyPy >= 2.1)
+  * R 2.x
+  * Plink should be in the PATH environment variable
+    See: http://pngu.mgh.harvard.edu/~purcell/plink/
+  * Plinktools should be in the PYTHONPATH environment variable
+    See: https://github.com/wtsi-npg/plinktools
+
+3. The zCall method
+-------------------
+
+Applying zCall consists of four steps:
+  1. Generate candidate zscore threshold files
+  2. Evaluate metrics on the input data for each candidate threshold
+  3. Merge evaluation results and choose an optimal threshold for calling
+  4. Apply zCall using the chosen threshold to re-call any 'no-calls' in 
+the input data.
+
+Steps 2 and 3 may be omitted in favour of simply using a default threshold; 
+JI Goldstein et al. suggest a default of z=7.
+
+4. Usage
+--------
+
+The user's PATH variable must include the path to the 'zcall' subdirectory of 
+the installation (or alternatively to src/zcall in the source code).
+
+Two approaches are supported:
+  1. A self-contained script to run zCall from start to finish as a single 
+process:  zCallComplete.py
+  2. Scripts to run each step independently, enabling parallelization.  For 
+steps 1-4 respectively: prepareThresholds.py, evaluateThresholds.py, 
+mergeEvaluation.py, runZCall.py
+
+Any of these scripts can be run with -h or --help for detailed help and 
+usage information. Note that if calling with a single default threshold is 
+desired, this can be achieved by simply using zCallComplete.py with a 
+threshold "range" of only one value.
+
+The original authors of zCall recommend that samples which fail quality 
+control should be excluded from zCall calibration and re-calling.
+
+5. Other scripts and modules
+----------------------------
+
+  * Additional command-line scripts, run with -h or --help for usage: 
+appendPEDline.py, createDocs.py, textToSampleJson.py
+  * Other Python files are modules which contain relevant classes but are not intended to be run as scripts.
+  * R script used to find thresholds: findBetas.r
+
+6. Data formats
+---------------
+
+Input data consists of binary .gtc files, one for each sample. In addition, 
+zCall requires an appropriate .egt cluster file and .bpm.csv manifest file. 
+The .egt and .bpm.csv files are proprietary Illumina formats and can be 
+downloaded from: http://support.illumina.com/downloads.ilmn
+
+Final genotype output is in Plink format (see 
+http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml).  Intermediate 
+metadata files are written in JSON, a simple text-based format for storing 
+data structures (see http://www.json.org/).
+
+7. Further reading
+------------------
+
+The document src/doc/zCallExtended.tex describes how the extended version 
+expands on the prototypes, including automated evaluation of candidate zscore 
+thresholds. HTML documentation for Python modules can be generated by running 
+src/zcall/createDocs.py. The zCall method is discussed in the paper by JI 
+Goldstein et al. (see Section 8 below).
+
+8. Strand normalization
+-----------------------
+
+Normalization of the beadpool manifest to the Illumina TOP strand is not 
+implemented by zCall, in either the prototype or extended software. 
+
+* If normalization is required, the simtools package developed at the 
+Wellcome Trust Sanger Institute can be used: https://github.com/wtsi-npg/simtools
+
+* For details of the normalization procedure, see: http://res.illumina.com/documents/products/technotes/technote_topbot.pdf‎
+
+9. History
+----------
+
+The zCall method was originally developed by Jackie Goldstein 
+(jigold@broadinstitute.org) et al. and published in the following paper:
+Goldstein JI, Crenshaw A, Carey J, Grant GB, Maguire J, Fromer M, 
+O'Dushlaine C, Moran JL, Chambert K, Stevens C; Swedish Schizophrenia 
+Consortium; ARRA Autism Sequencing Consortium, Sklar P, Hultman CM, Purcell S, 
+McCarroll SA, Sullivan PF, Daly MJ, Neale BM. zCall: a rare variant caller 
+for array-based genotyping: Genetics and population analysis. Bioinformatics. 
+2012 Oct 1;28(19):2543-2545. Epub 2012 Jul 27. PubMed PMID: 22843986.
+
+zCall has been substantially extended by Iain Bancarz (ib5@sanger.ac.uk), to 
+allow incorporation into the WTSI Genotyping Pipeline (WTSI-GP, see 
+https://github.com/wtsi-npg/genotyping). The extension includes metrics for 
+evaluation of the 'zscore' threshold parameter, automated implementation of 
+the complete zcall workflow, and support for Plink binary output.  WTSI-GP will 
+support zCall as part of a fully automated pipeline workflow, with parallel 
+processing on LSF.
+
+The prototype implementation of zCall consisted of several self-contained 
+versions for different input formats. These versions were committed to Github 
+as .zip files. Extended zCall is based on 'zCall_Version1.3_AutoCall.zip'. 
+Individual files in extended zCall are under version control on Github, and 
+are located in the src directory.
+
+Illumina provided code was provided as-is and with no warranty as to 
+performance and no warranty against it infringing any other party's 
+intellectual property rights.

README → README_prototypes

@@ -1,7 +1,7 @@
 --------------------------------------------------------------
 zCall: A Rare Variant Caller for Array-based Genotyping
 
-For questions about implementing zCall or reporting problems with the code, please send an email to Jackie Goldstein (jigold@broadinstitute.org). For all other inquiries, please send an email to both Ben Neale (bneale@broadinstitute.org) and Jackie Goldstein (jigold@broadinstitute.org).
+For questions about implementing zCall or reporting problems with the code, please send an email to Jackie Goldstein (jigold@broadinstitute.org). For all other inquiries, please send an email to both Ben Neale (bneale@broadinstitute.org) and Jackie Goldstein (jigold@broadinstitute.org).
 
 *** The Illumina provided Code was provided as-is and with no warranty as to performance and no warranty against it infringing any other party's intellectual property rights.
 
@@ -11,15 +11,21 @@ Goldstein JI, Crenshaw A, Carey J, Grant GB, Maguire J, Fromer M, O'Dushlaine C,
 
 
 I. Overview
-zCall is a variant caller specifically designed for calling rare single nucleotide polymorphisms (SNPs) from array-based technology. This caller is implemented as a post-processing step after a default calling algorithm has been applied such as Illumina�s GenCall algorithm. zCall uses the intensity profile of the common allele homozygote cluster to define the location of the other two genotype clusters. The outputs of zCall are either a PLINK PED and MAP file or a TPED and TFAM file depending on what the input files to zCall are. See http://pngu.mgh.harvard.edu/~purcell/plink/ for more details about PLINK.
+
+zCall is a variant caller specifically designed for calling rare single nucleotide polymorphisms (SNPs) from array-based technology. This caller is implemented as a post-processing step after a default calling algorithm has been applied such as IlluminaÕs GenCall algorithm. zCall uses the intensity profile of the common allele homozygote cluster to define the location of the other two genotype clusters. The outputs of zCall are either a PLINK PED and MAP file or a TPED and TFAM file depending on what the input files to zCall are. See http://pngu.mgh.harvard.edu/~purcell/plink/ for more details about PLINK.
 
 Three different versions of zCall have been provided depending on whether AutoCall or GenomeStudio has been run. Each version includes a set of Python scripts and an R script as well as line-by line instructions on how to run zCall located in the README file. Scripts are available at https://github.com/jigold/zCall. Examples of input file format are available at https://github.com/jigold/zCall; however, these files are subsets of actual inputs and are not usable for testing out the code.
 
 
 II. Implementation
-zCall is implemented as a set of Python scripts and one R script that are run on the command line. The scripts have been written to work with Illumina�s software output files such as an EGT file, GTC files, a .bpm.csv file, and a GenomeStudio Report. Examples of input files are located in the examples directory. We have written three versions of zCall that work with different input files based on whether AutoCall or GenomeStudio has been used for the initial calling:Version 1 (AutoCall): 
-Input files are Illumina�s binary GTC files, a .bpm.csv file, and a binary EGT file and the output files are a PLINK .ped and .map file. This version should be used after AutoCall has been run. Because the inputs are GTC files, calling of samples can be parallelized and the need to upload IDAT files into GenomeStudio can be avoided. In addition, this version avoids making large intermediate files with genotypes and intensity data.Version 2 (GenomeStudio - Thresholds derived from EGT file): 
-Input files are Illumina�s binary EGT file, an Illumina probe manifest file, and an Illumina GenomeStudio final report with the following fields: genotype, normalized X intensity, and normalized Y intensity where the normalized intensities are between 0 and 2. The output files are PLINK .tped and .tfam files where the alleles are A/B. PLINK can later be used to transform the A/B alleles into A,T,G,C using the --update-alleles flag.
+
+zCall is implemented as a set of Python scripts and one R script that are run on the command line. The scripts have been written to work with IlluminaÕs software output files such as an EGT file, GTC files, a .bpm.csv file, and a GenomeStudio Report. Examples of input files are located in the examples directory. We have written three versions of zCall that work with different input files based on whether AutoCall or GenomeStudio has been used for the initial calling:
+
+Version 1 (AutoCall): 
+Input files are IlluminaÕs binary GTC files, a .bpm.csv file, and a binary EGT file and the output files are a PLINK .ped and .map file. This version should be used after AutoCall has been run. Because the inputs are GTC files, calling of samples can be parallelized and the need to upload IDAT files into GenomeStudio can be avoided. In addition, this version avoids making large intermediate files with genotypes and intensity data.
+
+Version 2 (GenomeStudio - Thresholds derived from EGT file): 
+Input files are IlluminaÕs binary EGT file, an Illumina probe manifest file, and an Illumina GenomeStudio final report with the following fields: genotype, normalized X intensity, and normalized Y intensity where the normalized intensities are between 0 and 2. The output files are PLINK .tped and .tfam files where the alleles are A/B. PLINK can later be used to transform the A/B alleles into A,T,G,C using the --update-alleles flag.
 
 This version requires using GenomeStudio to generate a final report, which can be on the order of 40 GB for 5000 samples. However, it makes use of the fact that an EGT file already has the calculated means and standard deviations of each cluster when the initial clustering was done. Therefore, if a cluster file has already been generated with a large number of samples, then using the existing cluster file to generate the thresholds is much faster and more accurate than calculating the mean and standard deviation of each cluster from the GenomeStudio report.
 
@@ -48,11 +54,11 @@ zCall works by using the intensity profile of the common allele homozygote clust
 
 5. Genotypes are assigned based on where points are located with respect to the two thresholds.
 
-6. Analysts should only consider sites and samples that meet QC criteria from the original calls.
+6. Analysts should only consider sites and samples that meet QC criteria from the original calls.
 
 
 IV. Miscellaneous
 
 The run time and memory requirements are dependent on the number of probes and samples. For the Illumina HumanExome BeadChip, the run time to call one sample from a GTC file is 15 seconds and requires about 3 MB of memory (Max Memory Swap = 39 MB). For 971 samples calling from GTC files and generating one large PED file is 145 MB of memory (Max Memory Swap = 438 MB). For 384 samples from an Illumina Genome Studio report, it takes about 13 minutes to generate a new TPED and TFAM file with predetermined thresholds and 129 MB of memory (Max Memory Swap = 289 MB). It is assumed that the number, order, and names of sites match between the EGT, GTC, GenomeStudio report, and .bpm.csv files.
 
-For best results, we recommend using at least 1000 samples to derive the thresholds. This can either be an EGT file clustered with 1000 samples or from a GenomeStudio report with 1000 samples. We found that we got the best results using a z-score threshold equal to 7; however we recommend using the calibration step to verify that 7 will work well for your data. When using GenCall to generate the original genotype calls, we got the best results when using a smaller cluster file (90 samples) rather than a larger cluster file (9,479 samples). Please see Goldstein,J.I. et al. (2012) zCall: A Rare Variant Caller for Array-based Genotyping. Bioinformatics: 28(19):2543-2545 for more information.
+For best results, we recommend using at least 1000 samples to derive the thresholds. This can either be an EGT file clustered with 1000 samples or from a GenomeStudio report with 1000 samples. We found that we got the best results using a z-score threshold equal to 7; however we recommend using the calibration step to verify that 7 will work well for your data. When using GenCall to generate the original genotype calls, we got the best results when using a smaller cluster file (90 samples) rather than a larger cluster file (9,479 samples). Please see Goldstein,J.I. et al. (2012) zCall: A Rare Variant Caller for Array-based Genotyping. Bioinformatics: 28(19):2543-2545 for more information.

src/data/.gitattributes

@@ -0,0 +1,2 @@
+# do not diff the .csv manifest (treat as binary data)
+*.bpm.csv -crlf -diff

src/data/.gitignore

@@ -0,0 +1,6 @@
+# ignore output from tests
+thresholds_*.txt
+*.egt
+*.bpm.csv
+thresholds.json
+tmp*

src/data/gtc.json

@@ -0,0 +1,26 @@
+[
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R01C01.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R01C02.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R02C01.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R02C02.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R03C01.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R03C02.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R04C01.gtc"
+    },
+    {
+        "result": "/nfs/gapi/data/genotype/zcall_test/7197037167_R04C02.gtc"
+    }
+]