|
|
| Line 27: |
Line 27: |
| | | | |
| | = Performance = | | = Performance = |
| − |
| |
| − | == Pre-phasing ==
| |
| − | For the pre-phasing step, cost increases quadratically with the number of states and linearly with the number of rounds. The following table provides a simple example.
| |
| − |
| |
| − | {| class="wikitable" border="1" cellpadding="2"
| |
| − | |- bgcolor="lightgray"
| |
| − | ! States
| |
| − | ! Cost per round
| |
| − | |-
| |
| − | | 100 states
| |
| − | | 3 min
| |
| − | |-
| |
| − | | 200 states
| |
| − | | 12 min = (3 min * 2<sup>2</sup>)
| |
| − | |-
| |
| − | | 400 states
| |
| − | | 48 min = (3 min * 4<sup>2</sup>)
| |
| − | |-
| |
| − | | 500 states
| |
| − | | 75 min = (3 min * 5<sup>2</sup>)
| |
| − | |}
| |
| − |
| |
| − | So, in this case running haplotyping with 500 states and 10 rounds would require 75 min * 10 = 750 min.
| |
| − |
| |
| − | Typically, haplotype quality improves rapidly with the number of states but only slowly with the number of rounds. We recommend running ~20 rounds of the MaCH haplotyper and selecting a number of states as high as your patience will allow (but ideally greater than 200).
| |
| − |
| |
| | = Getting Started = | | = Getting Started = |
| − |
| |
| − | Using minimac for genotype imputation involves two steps. First, you will have to estimate haplotypes for your entire sample -- this will be the more computationally demanding step. Once that is done, you will be ready to quickly impute missing genotypes using the reference panel of your choice.
| |
| − |
| |
| − | == Estimating Haplotypes for Your Sample ==
| |
| − |
| |
| − | For the haplotyping step, we current recommend using [[MaCH]] with the --phase command line option. As input [[MaCH]] will need [[Merlin]] format pedigree and data files. All markers should be ordered according to their physical position and alleles should be labeled on the forward strand.
| |
| − |
| |
| − | === Preparing Your Data ===
| |
| − |
| |
| − | To get started, you will need to store your data in [[Merlin]] format pedigree and data files, one per chromosome. For details, of the Merlin file format, see the [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin Tutorial].
| |
| − |
| |
| − | Within each file, markers should be stored by chromosome position. Alleles should be stored in the forward strand and can be encoded as 'A', 'C', 'G' or 'T' (there is no need to use numeric identifiers for each allele).
| |
| − |
| |
| − | We recommend that, if at all possible, you should phase your chromosomes according to NCBI build 37. Future releases of the 1000 Genomes Reference panel and other public sets of reference haplotypes are expected to use this genome build.
| |
| − |
| |
| − | If figuring out position and strand for each marker seems like hard work, don't despair. For you, this should be the hardest bit of the entire process! For the computer, the fun is about to start.
| |
| − |
| |
| − | ==== NCBI build 36 / NCBI build 37 ====
| |
| − | The 1000G June reference panel is on build 36, future 1000G reference panels will be on build 37. This has also some impact on the pre-phasing:
| |
| − | some SNPs will ordered differently in each genome build (we expect order will be more accurate in more recent builds!). Purists will claim that imputation using NCBI build 36 and NCBI build 37 reference panels requires phasing using the exact same reference panel -- others will claim this makes little difference.
| |
| − |
| |
| − | === Running MaCH ===
| |
| − |
| |
| − | A typical MaCH command line to estimate phased haplotypes might look like this:
| |
| − |
| |
| − | mach1 -d sample.dat -p sample.ped --rounds 20 --states 200 --phase --interim 5 --sample 5 --compact
| |
| − |
| |
| − | This will request that MaCH estimate haplotypes for your sample, using 20 iterations of its Markov sampler and conditioning each update on up to 200 haplotypes. A summary description of these parameters follows (but for a more complete description, you should go to the MaCH website):
| |
| − |
| |
| − | {| class="wikitable" border="1" cellpadding="2"
| |
| − | |- bgcolor="lightgray"
| |
| − | ! Parameter
| |
| − | ! Description
| |
| − | |-
| |
| − | |style=white-space:nowrap|<code>-d sample.dat</code>
| |
| − | | Data file in [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin format]. Markers should be listed according to their order along the chromosome.
| |
| − | |-
| |
| − | | <code>-p sample.ped</code>
| |
| − | | Pedigree file in [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html Merlin format]. Alleles should be labeled on the forward strand.
| |
| − | |-
| |
| − | | <code>--states 200</code>
| |
| − | | Number of haplotypes to consider during each update. Increasing this value will typically lead to better haplotypes, but can dramatically increase computing time and memory use. A value of 200 - 400 is typical.
| |
| − | |-
| |
| − | | <code>--rounds 20</code>
| |
| − | | Iterations of the Markov sampler to use for haplotyping. Typically, using 20 - 30 rounds should give good results. To obtain better results, it is usually better to increase the <code>--states</code> parameter.
| |
| − | |-
| |
| − | | <code>--interim 5</code>
| |
| − | | Request that intermediate results should be saved to disk periodically. These will facilitate analyses in case a run doesn't complete.
| |
| − | |-
| |
| − | | <code>--sample 5</code>
| |
| − | | Request that random (but plausible) sets of haplotypes for each individual should be drawn every 5 iterations. This parameter is optional, but for some rare variant analyses, these alternative haplotypes can be very useful.
| |
| − | |-
| |
| − | | <code>--phase</code>
| |
| − | | Tell [[MaCH]] to estimate phased haplotypes for each individual.
| |
| − | |-
| |
| − | | <code>--compact</code>
| |
| − | | Reduce memory use at the cost of approximately doubling runtime. This option is recommended for most GWAS scale datasets and computing platforms.
| |
| − | |}
| |
| − |
| |
| − | == Imputation into Phased Haplotypes ==
| |
| − |
| |
| − | Imputing genotypes using '''minimac''' is an easy and straightforward process: after selecting a set of reference haplotypes, plugging-in the target haplotypes from the previous step and setting the number of rounds to use for the model parameter estimation, imputation should proceed rapidly.
| |
| − |
| |
| − | === Creating SNP List File ===
| |
| − |
| |
| − | '''Minimac''' requires a file listing markers in the haplotype file. This file can be easily generated by extracting the second column from the .dat file. In a standard Unix system, a command like this should do:
| |
| − |
| |
| − | cut -f 2 -d " " sample.dat > target.snps
| |
| − |
| |
| − | === Running Minimac ===
| |
| − |
| |
| − | A typical minimac command line might look like this:
| |
| − |
| |
| − | ==== using a VCF reference panel ====
| |
| − | minimac --vcfReference --refHaps ref.vcf.gz --haps target.hap.gz --snps target.snps.gz --rounds 5 --states 200 --prefix results
| |
| − |
| |
| − | '''Note''': GWAS SNPs (file --snps target.snps.gz) are by default expected to be in the chr:pos format e.g. 1:1000 and on build37/hg19; otherwise, please set the --rs flag
| |
| − |
| |
| − | ==== using a MaCH reference panel ====
| |
| − |
| |
| − | minimac --refHaps ref.hap.gz --refSnps ref.snps.gz --haps target.hap.gz --snps target.snps.gz --rounds 5 --states 200 --prefix results
| |
| − |
| |
| − | A detailed description of all minimac options is available [[Minimac Command Reference|elsewhere]]. Here is a brief description of the above parameters:
| |
| − |
| |
| − | {| class="wikitable" border="1" cellpadding="2"
| |
| − | |- bgcolor="lightgray"
| |
| − | ! Parameter
| |
| − | ! Description
| |
| − | |-
| |
| − | | <code>--refSnps ref.snps.gz </code>
| |
| − | | List of SNPs in the reference panel
| |
| − | |-
| |
| − | | <code>--refHaps ref.hap.gz </code>
| |
| − | | Reference haplotypes (e.g. from [http://www.sph.umich.edu/csg/abecasis/MACH/download/ MaCH download page])
| |
| − | |-
| |
| − | | <code>--vcfReference </code>
| |
| − | | This option specifies that the provided --refHaps file is provided in VCF format , no --refSNPs file needed.
| |
| − | |-
| |
| − | | <code>--rs </code>
| |
| − | | In combination with --vcfReference, allows to use rs GWAS SNP identifiers
| |
| − | |-
| |
| − | | <code>--snps target.snps.gz </code>
| |
| − | | SNPs in phased haplotypes. These should largely be a subset of the SNPs in the reference panel.
| |
| − | |-
| |
| − | | <code>--haps target.hap.gz </code>
| |
| − | | Phased haplotypes where missing genotypes will be imputed.
| |
| − | |-
| |
| − | | <code>--sample target.sample </code>
| |
| − | | Sample list in ShapeIT format.
| |
| − | |-
| |
| − | | <code>--shape_haps target.hap.gz </code>
| |
| − | | ShapeIT phased haplotypes where missing genotypes will be imputed.
| |
| − | |-
| |
| − | | <code>--chr 22</code>
| |
| − | | Chromosome for which we will carry out imputation (needed to run ShapeIT with chr:pos identifiers - default setting).
| |
| − | |-
| |
| − | | <code>--rounds 5</code>
| |
| − | | Rounds of optimization for model parameters, which describe population recombination rates and per SNP error rates.
| |
| − | |-
| |
| − | | <code>--states 200</code>
| |
| − | | Maximum number of reference (or target) haplotypes to be examined during parameter optimization.
| |
| − | |-
| |
| − | | <code>--prefix imputed</code>
| |
| − | | Optionally, a string that is used to help generate output file names.
| |
| − | |}
| |
| − |
| |
| − | === Reference Haplotypes ===
| |
| − |
| |
| − | Reference haplotypes generated by the 1000 Genomes project and formatted so that they are ready for analysis are available from the [http://www.sph.umich.edu/csg/abecasis/MACH/download/ MaCH download page]. As of this writing, the most recent set of haplotypes are based on genotype calls were generated in May 2011 and are an interim analysis of Project's Phase I data.
| |
| − |
| |
| − | === Imputation quality evaluation ===
| |
| − | To evaluate imputation quality, Minimac hides data for each genotyped SNP in turn and calculates 3 statistics:
| |
| − | * looRSQ - this is the estimated rsq for that SNP (as if SNP weren't typed).
| |
| − | * empR - this is the empirical correlation between true and imputed genotypes for the SNP. If this is negative, the SNP is probably flipped.
| |
| − | * empRSQ - this is the actual R2 value, comparing imputed and true genotypes.
| |
| − |
| |
| − | These statistics can be found in the .info file
| |
| − |
| |
| − | === Additional Sources of Information ===
| |
| − |
| |
| − | If the combination of MaCH and Minimac still runs too slowly for you, and you have access to a multi-processor compute cluster, you can look at [[ChunkChromosome]] page to learn how to conveniently split each chromosome into multiple segments that can be analyzed in parallel.
| |
| − |
| |
| − | If you are especially interested in 1000 Genomes Imputation, then you should look at the [[Minimac: 1000 Genomes Imputation Cookbook]].
| |
| − |
| |
| − | == X Chromosome Imputation ==
| |
| − | minimac supports the imputation of genotypes on the X chromosome (non-pseudo-autosomal part).
| |
| − |
| |
| − | # Split the X chromosome pedigree file by sex.
| |
| − | ## For females: follow the same protocol as for autosomes (phase with MaCH and impute with minimac).
| |
| − | ## For males
| |
| − | ### Remove the pseudo-autosomal part (for build hg18: chrX:1-2709520 and chrX:154584238-154913754 ; for build hg19 chrX:60001-2699520 and chrX:154931044-155260560)
| |
| − | ### Convert the pedigree file into a MaCH haplotype file (missing genotypes should be encoded as: "0" or "." or "N" ) and impute using minimac as described above.
| |
| − |
| |
| − |
| |
| − |
| |
| − | :::: '''<Example of a male only pedigree file>'''
| |
| − | :::: FAM1003 ID1234 0 0 M A/A A/A C/C
| |
| − | :::: FAM1004 ID5678 0 0 M 0/0 C/0 G/G
| |
| − | :::: ...
| |
| − | :::: '''<End of pedigree file>'''
| |
| − |
| |
| − | :::: ''Note that, consistent with the Merlin convention, hemizygous males are listed as if they were homozygous.''
| |
| − |
| |
| − | :::: '''<Example of the corresponding haplotype file>'''
| |
| − | :::: FAM1003->ID1234 HAPLO1 AAC
| |
| − | :::: FAM1003->ID1234 HAPLO2 AAC
| |
| − | :::: FAM1004->ID5678 HAPLO1 0CG
| |
| − | :::: FAM1004->ID5678 HAPLO2 0CG
| |
| − | :::: ...
| |
| − | :::: '''<End of the corresponding haplotype file>'''
| |
| − |
| |
| − | = Post-imputation Association Analysis =
| |
| − | == Quantitative Traits ==
| |
| − | Please use [http://www.sph.umich.edu/csg/yli/mach/download/mach2qtl.source.V108.tgz mach2qtl].
| |
| − |
| |
| − | == Binary Traits ==
| |
| − | Please use [http://www.sph.umich.edu/csg/yli/mach/download/mach2dat.source.1.0.18.tgz mach2dat]. Versions 1.0.18 and above accommodate to minimac output.
| |
| − |
| |
| − | = Reference =
| |
| − |
| |
| − | If you use minimac, please cite:
| |
| − |
| |
| − | Howie B, Fuchsberger C, Stephens M, Marchini J, and Abecasis GR.
| |
| − | Fast and accurate genotype imputation in genome-wide association studies
| |
| − | through pre-phasing. Nature Genetics 2012 [http://www.nature.com/ng/journal/vaop/ncurrent/full/ng.2354.html]
| |
| | | | |
| | = Questions and Comments = | | = Questions and Comments = |
