MaCH

MaCH (MArkov Chain Haplotyping), mostly known as a software for genotype imputation, is a Hidden Markov Model (HMM) based haplotyper can that reconstruct haplotypes from genotypes of unrelated individuals. Three primary uses of MaCH are (1) to resolve haplotypes from diploid genotypes; (2) impute missing genotypes; and (3) perform disease mapping analysis.

Input Files

MaCH requires a pedigree file and a data file as input. The pedigree file stores genotypes and phenotypes for each individual. The data file describes the contents of the pedigre file.Both files should be in Merlin format, which is described at http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html

Typically, MaCH expects markers to be ordered by position in the data and pedigree files.

Pedigree File (mandatory)

Each person is listed in a separate line of the pedigree file. The lines should start with 5 canonical fields, which are:

• Family ID
• Individual ID
• Father ID (typically zero for unrelated individuals)
• Mother ID (typically zero for unrelated individuals)
• Sex (encoded as M or 1 for males and F or 2 for females)

These 5 canonical fields, which are present in all pedigree files, will typically be followed by a series of marker genotypes. Here is an example of a complete pedigree file:

<sample.ped>

 fam1 indiv1 0 0 1 ./. C/G ./. 
 fam2 indiv2 0 0 2 A/C C/C A/T 

This sample pedigree contains 2 individuals.

The first individual is a male from family fam1 and is named indiv1. Parental codes are set to zero; since there is no individual named 0 in fam1, these simply indicate the individual is a founder. The sex code is 1, indicating the individual is a male. Two of the genotypes for this individual are missing and have been set to ./.. Missing genotypes can also be coded as 0/0.

The second individual is a female founder from family fam2 with ID indiv2. There are no missing genotypes for this individual.

In a file that includes only unrelated individuals, you could set family and individual ids to be identical for every individual. Alternatively, you could simply set the individual id or family id to be 1 for all individuals.

Marker Information File (mandatory)

An example data file that matches the pedigree above might be:

<sample.dat>

 M SNP1
 M SNP2
 M SNP3

This file tells us that the three files that follow the sex code in the pedigree file store genotypes for 3 markers, named SNP1, SNP2 and SNP3. MaCH expects SNPs to be ordered by position.

Optional Input files

Reference Haplotypes

Reference Haplotypes (typically from HapMap or 1000 Genomes) are optional. Mach 1.0 accepts two different formats: MACH format or HapMap format.

MACH format SNP File

One line per SNP and one field (marker name) only.

For example:

marker1
marker2
...

MACH format Haplotype File

One line per haplotype. Each haplotype can be preceded by a series of annotation columns.

For example:

 H_0001->H_0001 HAPLO1 CGGCGCGCTTGGC
 H_0001->H_0001 HAPLO2 CGGCGCGTCCAGC
 H_0002->H_0002 HAPLO1 GGGCGCGCTTGGC
 H_0002->H_0002 HAPLO2 GGAAGCACTCGGC
 ...

HapMap format reference files

HapMap format SNP File: legend file downloaded from HapMap website

HapMap format Haplotype File: phase fileddownloaded from HapMap website

HapMap format files can be downloaded from http://hapmap.org/downloads/phasing/2006-07_phaseII/phased/ or http://hapmap.org/downloads/phasing/2007-08_rel22/phased/

When using HapMap format files, turn on --hapmapFormat option. For example:

 mach1 -d sample.dat -p sample.ped -s genotypes_chr14_CEU_r21_nr_fwd_legend.txt -h genotypes_chr14_CEU_r21_nr_fwd_phased.gz --hapmapFormat ...

Options

Input Files: --datfile Marker information file for subjects under study. --pedfile Pedigree file for subjects under study.

FAQ

Why and how to perform a 2-step imputation?

When one has a large number of individuals (>1000), we recommend a 2-step imputation to speed up.

     A 2-step imputation contains the following 2 steps:

    (step 1) a representative subset of >= 200 unrelated individuals are used to calibrate model parameters; and
    (step 2) actual genotype imputation is performed for every person using parameters inferred in step 1.

      Example command lines for a 2-step imputation:

# step 1:
mach1 -d sample.dat -p subset.ped -s chr20.snps -h chr20.hap --compact --greedy --autoFlip -r 100 -o par_infer > mach.infer.log

# step 2:
mach1 -d sample.dat -p sample.ped -s chr20.snps -h chr20.hap --compact --greedy --autoFlip --errorMap par_infer.erate --crossoverMap par_infer.rec --mle --mldetails > mach.imp.log

Where can I find combined HapMap reference files?

A: http://www.sph.umich.edu/csg/yli/mach/download/HapMap-r21.html

Where can I find HapMap III / 1000 Genomes reference files?

A: http://www.sph.umich.edu/csg/yli/mach/download/

Does --mle overwrite input genotypes?

Yes, but not often. The --mle option outputs the most likely genotype configuration taking into account observed genotypes and integration over the most similar reference haplotypes. The original genotypes will be changed only if the underlying reference haplotypes strongly contradict the input genotype.

How do I get imputation quality estimates?

A simple approach is to use --mask option. For example, --mask 0.02 masks 2% of the genotypes at random, impute them and compare with the masked original to estimate genotypic and allelic error rates. Messages like the following will be generated to stdout:

 Comparing 948352 masked genotypes with MLE estimates ...
 Estimated per genotype error rate is 0.0568
 Estimated per allele error rate is 0.0293 

A better approach is to mask a small proportion of SNPs (vs. genotypes in the above simple approach). One can generate a mask.dat from the original .dat file by simply changing the flag of a subset of markers from M to S2 without duplicating the .ped file. Post-imputation, one can use   CalcMatch and doseR2.pl to estimate genotypic/allelic error rate and correlation respectively. Both programs can be downloaded from http://www.sph.umich.edu/csg/ylwtx/software.html.

Warning: Imputation involving masked datasets should be performed separately for imputation quality estimation. For production, one should use all available information.

Shall I apply QC before or after imputation? If so, how?

We strongly recommend QC both before and after imputation. Before imputation, we recommend the standard battery of QC filters including HWE, MAF (recommended cutoff is 1% for genotyping-based GWAS), completeness, Mendelian inconsistency etc. Post-imputation, we recommend Rsq 0.3 (which removes >70% of poorly-imputed SNPs at the cost of <0.5% well-imputed SNPs) and MAF of 1%.

How do I get reference files for an region of interest?

1. For HapMapII format, download haplotypes from http://www.sph.umich.edu/csg/ylwtx/HapMapForMach.tgz 2. For MACH format, you can do the following:

• First, find the first and last SNP in the region you are interested in. Say "rsFIRST" and "rsLAST", defined according to position.
• Then:

 @ first = `grep -n rsFIRST orig.snps | cut -f1 -d ':'`
 @ last = `grep -n rsLAST orig.snps | cut -f1 -d ':'`

• Finally (assuming the third field contains the actual haplotypes, where alleles are separated by whitespace):

 awk '{print $3}' orig.hap | cut -c${first}-${last} > region.hap

Do I have to sort the pedigree file by physical positions?

If you use external reference, you do not have to as long as the external reference is in correct order. **HOWEVER**, you will probably avoid problems by including markers in the pedigree file sorted in chromosome order.

What if I specified --states R where R exceeds the maximum possible (2*number diploid individuals - 2 + number_haplotypes)?

Mach will automatically change this to the maximum possible value.

How is AL1 defined? Which allele dosage is .dose/.mldose counting?

AL1 is an arbitrary allele. Typically, it is the first allele read in the reference haplotypes (file fed to -h or --haps). The earliest versions (prior to April 2007) of mach counted the expected number copies of AL2 and more recent versions count the number of AL1. One can find out which allele is counted following the steps below.

1. . First, find the two alleles for one of the markers in your data

 prompt> head -2 mlinfo/chr21.mlinfo 
 SNP      Al1 Al2 Freq1   MAF    Quality  Rsq 
 rs885550 2   4   0.9840  0.0160  0.9682  0.992

1. . Second, check the dosage for a few individuals at this SNP.

 prompt> head -3 mldose/chr21.mldose | cut -f3 -d ' ' 
 1.962 
 1.000
 0.078

1. . Finally, compare these dosages to genotypes.

 prompt> head -1 mlgeno/chr21.mlgeno | cut -f3 -d ' ' 
 2/2 
 2/4
 4/4

In this example, you can see that the first individual has a high dosage count (1.962) and most likely genotype 2/2. The last individual has a low dosage count and most likely genotype 4/4. Thus, the output corresponds to version of Mach released after April 2007, which should tally allele 1 counts.

Note that, on the example above, .mldose could be replaced with .dose and .mlgeno could be replaced with .geno.

Based on the three files above, we've confirmed that dosage is the number of AL1 copies: you will only to check for one informative case (i.e, dosage values close to 0 or 2) since it's consistent across all individuals and all SNPs.

Can I used unphased reference?

A: Yes. You could create pedigree (.ped) and data files (.dat) that include both reference panel and sample genotypes or request that MaCH merge apppropriate files on the fly.

For example, if you have:

reference.dat

M SNP1
M SNP2
M SNP3
M SNP4
M SNP5

reference.ped

REF1 REF1 0 0 1 A/C C/C G/G G/A A/A

sample.dat

M SNP1
M SNP4 

sample.ped

1 1 0 0 1 A/A G/G

Your could create a combined data set as:

comb.dat

M SNP1
M SNP2
M SNP3
M SNP4
M SNP5

comb.ped

REF1 REF1 0 0 1 A/C C/C G/G G/A A/A

  1    1 0 0 1 A/A ./. ./. G/G ./. 

Equivalently, you could write -d reference.dat,sample.dat -p reference.ped,sample.ped on the command line and MACH would merge both files on-the-fly.

How long does imputation take?

The following factors/parameters affect computational time:

1. m, # of genotyped markers (number of markers in .dat file)
   
2. n, # of individuals
   
3. h, # of reference haplotypes (determined by --greedy or states, by default, h = 2*number diploid individuals - 2 + number_haplotypes)
   
4. r, # of rounds (-r or --rounds, --mle corresponds to 1-2 rounds)

Computational time increases linearly with m, n, r and quadratically with h. On our Xeon 3.0GHz machine, imputation with m=25K, n=250, h=120, and r=100 takes ~20 hours.

More questions?

Email Yun Li or Goncalo Abecasis.

Examples

Imputation

 mach1 -d sample.dat -p sample.ped -s hapmap.snps -h hapmap.hap -r 100 -o phase