Minimac: 1000 Genomes Imputation Cookbook

From Genome Analysis Wiki
Revision as of 11:17, 10 August 2011 by Zhanxw (talk | contribs) (Your Own Data)
Jump to: navigation, search

This page documents how to impute 1000 Genome SNPs using Minimac, which is typically the preferred approach for imputation using large reference panels such as the 1000 Genomes data. For information on how to carry out 1000 Genomes Imputation using MaCH, see MaCH: 1000 Genomes Imputation Cookbook. For pointers on how to carry out 1000 Genomes Imputation using IMPUTE2, see IMPUTE2: 1000 Genomes Imputation Cookbook.

Before reading this tutorial, you might find it useful to spend a few minutes reading through the main Minimac documentation.

Should I use MaCH or minimac?

The imputation quality of both approaches is comparable.

If you are carrying out only an one-time imputation, using a relatively small reference panel (with <200 haplotypes) and focusing only on autosomes - then go with MaCH. If you are using a large reference panel, planning to re-impute your data as new reference panels become available or interested in X-chromosome imputation - then minimac is the right choice for you.

As of June 2011, the 1000 Genomes Panel set of reference haplotypes includes >2000 haplotypes, and you are likely to be better off using minimac.

Getting Started

Your Own 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 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).

The latest reference panel generated by the 1000 Genomes project uses NCBI Build 37 (HG 19). Make sure that your data is on Build 37 (or Minimac may ignore genotyped markers whose names have changed in build 37). If you are trying to convert your data from an earlier genome build to build 37, you'll probably find the dbSNP merge table (table description on the NCBI website), which logs rs# changes between dbSNP builds, and the UCSC online liftOver tool, which converts genome positions between different genome builds, to be quite useful. We have aslo documented LiftOver link a general procedure to convert genome positions and rs number between builds.

If you are planning to use imputation with the MetaboChip, you might find a list of SNPs whose order varies between NCBI genome build 36 and 37 convenient. Here it is: List of Metabochip SNPs Whose Order Changes With Build

Quality Control

You should apply standard quality control filters to the set of SNPs that you use as input to the imputation procedure. These filters are typically study specific but usually a series of per marker and per individual quality checks.

Per marker quality checks typically evaluate:

  • Genotype completeness (markers with relatively high missingness rates are typically excluded)
  • Hardy Weinberg equilibrium (markers with clear deviations from Hardy Weinberg proportions are typically excluded)
  • Duplicate concordance (markers with high discordance rates are excluded)
  • Mendelian inconsistencies (when parent-offspring pairs are available, markers with an excess of Mendelian inconsistencies are typically excluded).
  • Polymorphism check (monomorphic markers can represent assay failures and are typically excluded).

Per individual quality checks typically include:

  • Check for per sample completeness (samples with high missingness rates are typically excluded)
  • Check for per sample heterozygosity (samples with unusual heterozygosity are typically excluded)
  • Inspection of principal component or multi-dimensional scaling plots (genetic ancestry outliers are usually excluded)
  • Check for related individuals (if related samples are identified, the analysis plan must be adjusted or they must be removed)

With older genotyping platforms, low frequency SNPs are also often excluded because they are hard to genotype accurately. With more modern genotyping arrays, the accuracy of genotype calls for low frequency SNPs is less of a concern.

Two good ways to verify that you have used appropriate quality control steps are to generate a Q-Q plot for your dataset and to calculate a genomic control parameter. If these are not satisfactory, it may be useful to investigate before proceeding to imputation.

Reference Haplotypes

Reference haplotypes generated by the 1000 Genomes project and formatted so that they are ready for analysis are available from the MaCH download page. In our hands, it is ideal to always use the most recent release since generation of additional sequence data, improvements in variant discovery, genotyping and haplotyping strategies typically result in noticeable improvements in data quality. When this page was written, the most recent release of 1000 Genomes project haplotypes was the Interim Phase I release, including 1094 individuals genotyped at >37 million autosomal sites.

Minimac Imputation

Minimac relies on a two step approach. First, the samples that are to be analyzed must be phased into a series of estimated haplotypes. Second, imputation is carried out directly into these phased haplotypes. As newer reference panels become available, only the second step must be repeated.

Pre-phasing - MaCH

A convenient way to haplotype your sample is to use MaCH. A typical MaCH command line to estimate phased haplotypes might look like this:

 mach1 -d chr1.dat -p chr1.ped --rounds 20 --states 200 --phase --interim 5 --sample 5

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):

Parameter Description
-d sample.dat Data file in Merlin format. Markers should be listed according to their order along the chromosome.
-p sample.ped Pedigree file in Merlin format. Alleles should be labeled on the forward strand.
--states 200 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. If computing time is a constraints, you may consider decreasing this parameter to 100. If you have substantial computing resources, consider increasing this value to 600 or even 800.
--rounds 20 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 --states parameter.
--interim 5 Request that intermediate results should be saved to disk periodically. These will facilitate analyses in case a run doesn't complete.
--sample 5 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. If you are not short of disk space, you should consider enabling this parameter.
--phase Tell MaCH to estimate phased haplotypes for each individual.
--compact Reduce memory use at the cost of approximately doubling runtime.

You should be able to run this step in parallel and in our cluster we'd use:

   foreach chr (`seq 1 22`)

     runon -m 4096 mach -d chr$chr.dat -p chr$chr.ped --rounds 20 --states 200 --phase --interim 5 --sample 5


Imputation into Phased Haplotypes - minimac

Imputing genotypes using minimac is a 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 estimating model parameters (which describe the length and conservation of haplotype stretches shared between the reference panel and your study samples), imputation should proceed rapidly.

A typical minimac command line might look like this:

 minimac --refHaps ref.hap.$chr.gz --refSnps ref.snps.$chr.gz --haps target.hap.$chr.gz --snps target.snps.$chr.gz --rounds 5 --states 200 --prefix results

A detailed description of all minimac options is available elsewhere. Here is a brief description of the above parameters:

Parameter Description
--refSnps ref.snps.gz List of SNPs in the reference panel
--refHaps ref.hap.gz Reference haplotypes (e.g. from MaCH download page)
--snps target.snps.gz SNPs in phased haplotypes. These should largely be a subset of the SNPs in the reference panel.
--haps target.hap.gz Phased haplotypes where missing genotypes will be imputed.
--rounds 5 Rounds of optimization for model parameters, which describe population recombination rates and per SNP error rates.
--states 200 Maximum number of reference (or target) haplotypes to be examined during parameter optimization.
--prefix imputed Optionally, a string that is used to help generate output file names.

Again, you can speed-up things by running this step in parallel and in our cluster we'd use:

   foreach chr (`seq 1 22`)

     runon -m 1024 minimac --refHaps ref.hap.$chr.gz --refSnps ref.snps.$chr.gz --haps target.hap.$chr.gz --snps target.snps.$chr.gz --rounds 5 --states 200 --prefix results


Further Time Savings

The recipe above imputes whole chromosomes, one a time. A further time savings is possible by imputing chromosomes in chunks, a process that can be facilitated using the ChunkChromosome tool. This tool automates some of the manual editing steps that would be required to divide each chromosome into more manageable chunks. It also interfaces with minimac to ensure SNPs that overlap between chunks are only imputed once.

Recall that the original analysis script might have looked like this:


# Build SNP lists for each chromosome
foreach chr (`seq 1 22`)
   awk '{ if ($1 == "M") print $2; }' < chr${chr}.dat > chr${chr}.snps

# Estimate haplotypes for all individuals, one chromosome at a time
foreach chr (`seq 1 22`)
   mach -d chr$chr.dat -p chr$chr.ped --rounds 20 --states 200 --phase --sample 5 --prefix chr$chr.haps >& chr$chr-mach.log &

# Impute into phased haplotypes
foreach chr (`seq 1 22`)
   minimac --refHaps ref.hap.$chr.gz --refSnps ref.snps.$chr.gz --haps chr$chr.haps.gz --snps chr$chr.snps --rounds 5 --states 200 --prefix chr$chr.imputed >& chr$chr-minimac.log &

A modified script, using ChunkChromosome would look like this:


@ length = 2500
@ overlap = 500

# Estimate haplotypes for all individuals, in 2500 marker chunks, with 500 marker overhang
foreach chr (`seq 1 22`)

   ChunkChromosome -d chr$chr.dat -n $length -o $overlap

   foreach chunk (chunk*-chr$chr.dat)

      mach -d $chunk -p chr$chr.ped --prefix ${chunk:r} \
           --rounds 20 --states 200 --phase --sample 5 >& ${chunk:r}-mach.log &



# Impute into phased haplotypes
foreach chr (`seq 1 22`)

   foreach chunk (chunk*-chr$chr.dat)

      set haps = /data/1000g/hap/all/20101123.chr$chr.hap.gz
      set snps = /data/1000g/snps/chr$chr.snps

      minimac --refHaps $haps --refSnps $snps --rounds 5 --states 200 \
              --haps ${chunk:r}.gz --snps ${chunk}.snps  --autoClip autoChunk-chr$chr.dat  \
              --prefix ${chunk:r}.imputed >& ${chunk:r}-minimac.log &



Although the loss in accuracy will typically be small, analyses of small chromosome chunks will necessarily result in slightly less accurate results. The differences will be large in populations where haplotype sharing is expected to extend further (for example, in Finland and Sardinia we find that haplotypes can be estimated accurately over several megabases and focusing on small chunks of chromosome will make these matches harder to identify).

X Chromosome Imputation

minimac supports the imputation of genotypes on the X chromosome (non-pseudo-autosomal part).

  1. Split the X chromosome pedigree file by sex.
    1. For females: follow the same protocol as for autosomes (phase with MaCH and impute with minimac).
    2. For males
      1. 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)
      2. 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; hemizygous males appear as homozygotes >
FAM1003 ID1234 0 0 M A/A A/A C/C
FAM1004 ID5678 0 0 M 0/0 C/C G/G
<End of pedigree file>

<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>

Imputation quality evaluation

Minimac hides each of the genotyped SNPs in turn and then 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 alleles are probably flipped.
  • empRSQ - this is the actual R2 value, comparing imputed and true genotypes.

These statistics can be found in the *.info file

Be aware that, unfortunately, imputation quality statistics are not directly comparable between different imputation programs (MaCH/minimac vs. Impute vs. Beagle etc.).

Questions and Comments

Please contact Goncalo Abecasis or Christian Fuchsberger (minimac) or Yun Li (MaCH).