Minimac: 1000 Genomes Imputation Cookbook

This page documents how to impute 1000 Genome SNPs using

• MaCH (concurrent phasing approach).
• or minimac (pre-phasing approach).

Shoud I use MaCH or minimac?

The imputation quality of both approaches is comparable.

If you are carrying out only an one-time imputation and you are focusing only on autosomes - then go with MaCH.

If you are planning to re-impute into new reference panels and your are interested into the X-chromosome - then minimac is the right choice for you.

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 [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).

The 1000 Genome pilot project genotypes use NCBI Build 37.

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, this August 2010 release is substantially better than previous 1000 Genome Project genotype call sets.

MaCH Imputation

Estimating Model Parameters

The first step for genotype imputation analyses using MaCH is to build a model that relates your dataset to a reference set of haplotypes. This model will include information on the length of haplotype stretches shared between your sample and the reference panel (in a .rec file) and information on the similarity of marker genotypes between your sample and the reference panel (in a .err file). The contents of the second file reflect the combined effects of genotyping error and differences in genotyping assays between the two samples.

In our view, MaCH's ability to estimate these parameters makes it especially robust to problems in genotyping, differences in genotyping assays and in the ability to adjust to varying degrees of similarity between the sample being imputed and a reference panel.

Selecting a Subset of Samples

Using all availalble samples to estimate model parameters can be quite time consuming and is typically not necessary. Instead, we recommend that you select a set of 200 samples to use in this step of analysis. If you are not worried about batch effects, you could just take the first 200 samples in your pedigree file, for example:

  head -200 chr1.ped > chr1.first200.ped

And then run mach as:

  mach -d chr1.dat -p chr1.first200.ped -s CEU.1000G/chr1.snps -h CEU.1000G/chr1.hap --greedy -r 100 -o chr1.parameters 

Optionally, you can add the --compact option to the end of the command line to reduce memory use.

If you'd like to play it safe and sample 200 individuals at random, you could use PedScript and a series of command similar to this one:

   pedscript << EOF
   READ DAT chr1.dat
   READ PED chr1.ped
   
   NUKE
   TRIM

   WRITE DAT chr1.rand200.dat
   SAMPLE 200 PERSONS TO chr1.rand200.ped

   QUIT
   EOF

Then, run mach just as before:

  mach -d chr1.rand200.dat -p chr1.rand200.ped -s CEU.1000G/chr1.snps -h CEU.1000G/chr1.hap --greedy -r 100 -o chr1.parameters 

You can decide on the exact number of samples to use for this analysis ... A number of between 100-500 individuals should be plenty. On a compute cluster, you should be able to run this step for all chromosomes in parallel; the exact syntax will depend on your cluster. In ours, we'd try something like:

   foreach chr (`seq 1 22`)

     pedscript << EOF
        READ DAT chr$chr.dat
        READ PED chr$chr.ped
        NUKE
        TRIM
        WRITE DAT chr$chr.rand200.dat
        SAMPLE 200 PERSONS TO chr$chr.rand200.ped
        QUIT
EOF

     runon -m 4096 mach -d chr$chr.rand200.dat -p chr$chr.rand200.ped -s CEU.1000G/chr$chr.snps -h CEU.1000G/chr$chr.hap \
                        --greedy -r 100 -o chr$chr.parameters --compact >& chr$chr.model.log &

   end

Imputation

Once model parameters have been estimated, you are ready to impute. This step should be comparatively fast, at least on a per individual basis.

Your command line should look like this:

   mach -d chr1.dat -p chr1.ped -s CEU.1000G/chr1.snps -h CEU.1000G/chr1.hap \
        --errorMap chr1.parameters.erate --crossoverMap chr1.parameters.rec \
        --compact --greedy --autoFlip --mle --mldetails >& chr$chr.impute.log 

Again, 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 -s CEU.1000G/chr$chr.snps -h CEU.1000G/chr$chr.hap \
           --errorMap chr$chr.parameters.erate --crossoverMap chr$chr.parameters.rec \
           --compact --greedy --autoFlip --mle --mldetails >& chr$chr.impute.log 

   end

minimac Imputation

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

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 100 - 400 is typical.
--rounds 20	Iterations of the Markov sampler to use for haplotyping. Typically, using 20 - 100 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.
--phase	Tell MaCH to estimate phased haplotypes for each individual.
--compact	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 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.

Running Minimac

A typical minimac command line might look like this:

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

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 >

FAM1003 ID1234 0 0 M A/0 A/0 C/0

FAM1004 ID5678 0 0 M 0/0 C/0 G/0

...

<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 drops 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 is probably flipped.
• empRSQ - this is the actual R2 value, comparing imputed and true genotypes.

These statistics can be found in the *.info file

Quality Filtering

For 1000 Genome SNPs, we currently recommend that any markers with estimated r² of <0.5 should be treated with caution. This is a bit more conservative than the threshold of 0.3 we recommend for HapMap; but 1000 Genome SNP genotypes are also (as of June 2010) of lower quality.

Questions and Comments

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