This page documents how to perform variant calling from low-coverage sequencing data using glfmultiples and thunder. The pipeline was originally developed by [mailto:yunli@med.unc.edu Yun Li] for the 1000 Genomes Low Coverage Pilot Project.

== Input Data ==

To get started, you will need glf files in the standard format [http://samtools.sourceforge.net/SAM1.pdf glf format]. Sample files are available at [ftp://share.sph.umich.edu/1000genomes/pilot1/examples/glf.tgz sample glf files].

If you do not have glf files, you can generate them from bam files (bam format also specified in [http://samtools.sourceforge.net/SAM1.pdf glf format bam format]) using the following command line:

samtools pileup -g -T 1 -f ref.fa my.bam > my.glf

Note: you will need the reference fasta file ref.fa to create glf file from bam file.

== How to Run ==

This variant calling pipeline has two steps. (step 1) promotion of a set of potential polymorphisms; and (step 2) genotype/haplotype calling using LD information.

=== (step 1) Site promotion using software glfMultiples [https://www.sph.umich.edu/csg/yli/GPT_Freq.011.source.tgz GPT_Freq] ===

GPT_Freq -b my.out -p 0.9 --minDepth 10 --maxDepth 1000 *.glf

minDepth and maxDepth are the cutoffs on total depth (across all individuals). We have found it useful to exclude sites with extremely low and high total depth. Please see Important Filters below.

=== (step 2) Genotype/haplotype calling using thunder [https://www.sph.umich.edu/csg/yli/thunder/thunder.V010.source.tgz thunder_glf_freq] ===

thunder_glf_freq --shotgun my.out.$chr -r 100 --states 200 --dosage --phase --interim 25 -o my.final.out

Notes:

(1) The program thunder used in step 2 is an extension of MaCH, the genotype imputation software we have previously developed. For details regarding the shared options, please check out [http://www.sph.umich.edu/csg/yli/mach/index.html MaCH website] and [http://genome.sph.umich.edu/wiki/Mach MaCH wiki].

(2) Check out example files and command lines under examples/thunder/ in the thunder package [https://www.sph.umich.edu/csg/yli/thunder/thunder.V009.source.tgz thunder_glf_freq].

== Important Filters ==

We have found that the following filters are helpful.

=== total depth filter ===
For the 1000 Genomes Project (average depth per individual ~4X), we have found it useful to exclude sites with average total depth per individual < 0.5X or > 10X.

=== coverage filter ===
We recommend the filter of >50% individuals with coverage.

=== flanking sequence filter ===
We recommend excluding sites with >0.1% flanking 10-mer frequency among candidate sites.

The rationale is ....

=== indel filter ===
We recommend distance to known indels >= 5bp. A catalog of known indels can be found at [ indel catalog].

=== site promotion filter ===
We recommend setting parameter -p at least >= 0.9 in step 1 (running glfMultiples).

== Questions and Comments? ==

Email [mailto:yunli@med.unc.edu Yun Li].

Mach DAC

2010-12-07T13:29:00Z

Ylwtx:

This is the MaCH Divide and Conquer page, documenting how to break the genome into smaller pieces before imputation/phasing and how to ligate after imputation/phasing.

== Phasing without External Reference ==
=== 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 Merlin tutorial [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html]. 

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

=== Split Your Data ===
You can split your data using [http://www.sph.umich.edu/csg/yli/splitPed/ splitPed].

== Phasing/Imputation with External Reference ==
When you phase/impute with external reference panel, you will only need to break the reference files into parts containing subsets of markers because SNPs in your own data (pedigree files) but not in reference files will be automatically ignored by MaCH and minimac. 

You can split the reference data using [http://www.sph.umich.edu/csg/yli/splitRef/ splitRef].

== Post Phasing/Imputation Ligation ==
You can use [http://www.sph.umich.edu/csg/yli/ligateHap.V002.tgz LigateHaplotypes ] to ligate the parts.

== Questions and Comments? ==

Email [mailto:yunli@med.unc.edu Yun Li].

Mach DAC

2010-12-07T13:27:56Z

Ylwtx: /* Post Phasing/Imputation Ligation */

Mach DAC

2010-12-07T13:27:14Z

Ylwtx:

SplitRef

2010-12-07T13:19:27Z

Ylwtx: Created page with 'This page documents the splitRef program, which splits a reference haplotype file into smaller files with subsets of markers. == Input Files == === Required Input Files === ====…'

This page documents the splitRef program, which splits a reference haplotype file into smaller files with subsets of markers.

== Input Files ==
=== Required Input Files ===
==== Haplotype file (.hap) ====
File fed to -hap option. One line for one haplotype, with last field containing the actual alleles with no separators between alleles. 

==== Marker list (.snps) file ====
File fed to -snps option. One line for each marker: marker name only. 

=== Optional Input Files ===
==== Map file ====
File fed to -map option, containing chromosome, marker name, and marker coordinate (in base pairs) information for each marker. Markers should be stored in the same order as in the marker information file. 

== Options ==
=== Required options ===
==== window size ====
Window size can be specified by one of the following three options: (1) -nWindows (2) -windowSize and (3) -windowLength. 
-nWindows specifies the number of windows to split into and the program splits markers evenly into output windows. 
-windowSize specifies the number of markers in one output window. The remainder goes to the last window. 
-windowLength specifies the length (in base pairs) of one output window. The remainder goes to the last window. Note that this option is only allowed when map input file is specified. 

==== flanking region ====
Size of flanking region on each side can be specified by one of the following two options: (1) -overlapSize and (2) -overlapLength. 
-overlapSize specifies the number of markers in each flanking region (so that the total number of flanking markers for each window is twice the number specified except for the first and last window). 
-overlapLength specifies the length (in base pairs) of each flanking region (so that the total length of the flanking regions is twice the number specified except for the first and last window). 

==== Output prefix ====
Specified by -o option. 

=== Additional options ===
==== Estimate window size only ====
This is controlled by -extimateWindowOnly option. By default, splitting is performed. 
But if one only wishes to peek into how the markers are allocated into output windows, use "-extimateWindowOnly 1". 

== Example Commands ==
splitRef.pl -hap example.hap.gz -snps example.snps -map example.map -windowLength 10000000 -overlapLength 1000000
splitRef.pl -hap example.hap.gz -snps example.snps -windowSize 10000 -overlapSize 1000
splitRef.pl -hap example.hap.gz -snps example.snps -nWindows 12 -overlapSize 1000

== Download ==
You can download splitPed at [http://www.sph.umich.edu/csg/yli/splitRef/download/ splitRef Download Page].

== Questions and Comments? ==

Email [mailto:yunli@med.unc.edu Yun Li].

MaCH FAQ

2010-12-05T13:18:47Z

Ylwtx: /* How do I get reference files for an region of interest? */

== How to speed up? ==

=== minimac ===
See [http://genome.sph.umich.edu/wiki/Minimac minimac] for details.

=== Divide and Conquer ===
See [http://genome.sph.umich.edu/wiki/Mach_DAC MaCH Divide and Conquer] for details.

=== 2-step imputation ===
See [http://genome.sph.umich.edu/wiki/MaCH_FAQ#Why_and_how_to_perform_a_2-step_imputation.3F 2-step imputation] for details.

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

In step1, one can use --greedy in combination with --states XX in MaCH versions 16.b and above. We have found that using 1/3 of the reference haplotypes (with 1/9 computational time) results in almost no power loss for the current HapMap and 1000G reference panels. 

In step2, each individual is imputed independently and can therefore be split into as many as n (sample size) jobs for each chromosome for parallelism.

For other approaches to speed up, see [how to speed up].

== Can MaCH perform imputation for chromosome X? ==
Yes. See [http://genome.sph.umich.edu/wiki/MaCH:_machX MaCH X Chromosome] for details.

== Where can I find combined HapMap reference files? ==

You can find them at http://www.sph.umich.edu/csg/yli/mach/download/HapMap-r21.html or on the HapMap Project website.

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

You can find these at the MaCH download page, which is at 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 (in the second step alone if using two-step imputation). 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   [http://www.sph.umich.edu/csg/ylwtx/CalcMatch.1.0.5.tgz CalcMatch ]and [http://www.sph.umich.edu/csg/ylwtx/doseR2.tgz 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 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.

== How do I interpret the imputation quality estimates? ==
In the simple approach, you will only get concordance/error estimates. There are two aspects to check. (1) the ratio between the genotypic error and allelic error. We expect that only a small proportion of errors where one homozygote is imputed as the other homozygote. Therefore, a ~2:1 ratio is expected. (2) the absolute error rate. There are several factors influencing imputation quality including the population to be imputed, the reference population and the genotyping panel used. Typically, we expect <2% allelic error rate among Caucasians and East Asians; 3-5% among Africans and African Americans. Figure below show imputation quality from the Human Genome Diversity Project (HGDP) for 52 populations across the world and by different HapMap reference panel.

http://www.sph.umich.edu/csg/yli/figure3.gif

Table 3 in the MaCH 1.0 paper tabulates imputation quality by commercial panel in CEU, YRI, and CHB+JPT.

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

The created reference files are in MaCH format. You do NOT need to turn on --hapmapFormat option.

== Do I always have to sort the pedigree file by marker position? ==

If you use a reference set of haplotypes, 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 specify ''--states R'' where ''R'' exceeds the maximum possible (2*number diploid individuals - 2 + number_haplotypes)? ==

Mach caps the number of states at 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. 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.

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

<source lang="text">
prompt> head -2 mlinfo/chr21.mlinfo
SNP Al1 Al2 Freq1 MAF Quality Rsq
rs885550 2 4 0.9840 0.0160 0.9682 0.992
</source>

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

<source lang="text">
prompt> head -3 mldose/chr21.mldose | cut -f3 -d ' '
1.962
1.000
0.078
</source>

#. Finally, compare these dosages to genotypes.

<source lang="text">
prompt> head -1 mlgeno/chr21.mlgeno | cut -f3 -d ' '
2/2
2/4
4/4
</source>

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 an unphased reference? ==

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 big are the imputation output file? ==
For 1,000 individuals with 8 million SNPs, gz compressed geno/dose/prob files take ~5Gb/10Gb/15Gb.

== How long does imputation take? ==

The following factors/parameters affect computational time:

#m, # of genotyped markers (number of markers in .dat file) 
#n, # of individuals 
#h, # of reference haplotypes (determined by --greedy or states, by default, h = 2*number diploid individuals - 2 + number_haplotypes) 
#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 (25000*250*120^2*100/4.5/10^11).

If you have a larger number of individuals to impute (e.g., > 1,000), we recommend a 2-step imputation manner http://genome.sph.umich.edu/wiki/MaCH_FAQ#Why_and_how_to_perform_a_2-step_imputation.3F.

== undefined symbol: gzopen64 ==
If you see this message, you will need to re-compile the program. Type the following commands:

make clear
make all

New executables mach1 and thunder will then be generated under folder executables/

== More questions? ==

Email [mailto:yunli@med.unc.edu Yun Li] or [mailto:goncalo@umich.edu Goncalo Abecasis].

Mach DAC

2010-11-29T20:18:24Z

Ylwtx:

Mach DAC

2010-11-29T19:33:53Z

Ylwtx:

Mach DAC

2010-11-29T18:52:23Z

Ylwtx:

This is the MaCH Divide and Conquer page, documenting how to break the genome into smaller pieces before imputation/phasing and how to ligate after imputation/phasing.

== Phase without External Reference ==
=== 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 Merlin tutorial [http://www.sph.umich.edu/csg/abecasis/Merlin/tour/input_files.html]. 

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

=== Split Your Data ===
You can split your data using [http://www.sph.umich.edu/csg/yli/splitPed/ splitPed].

== Phase/Imputation with External Reference ==

== Questions and Comments? ==

Email [mailto:yunli@med.unc.edu Yun Li].