We implemented GATK UnifiedGenotyper followed by imputation to maximize the amount of genetic information for downstream analyses and since genotype likelihood scores can be generated for imputation. We opted to impute diploid genotypes and missing sites for the individuals in our data set (using the 1000 Genomes Project Consortium reference panel).

The 1000 Genomes Phase 3 genetic variants reference panel was used for genotyping and imputation, as provided by BEAGLE (ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/). After removing variants that are not SNPs, multiallelic SNPs and X- and Y-chromosomes, 77,818,345 variants remained. UnifiedGenotyper was used to obtain genotypes and likelihood scores with the following parameters: --genotyping_mode --alleles <1000 Genomes reference panel>, --GENOTYPE_GIVEN_ALLELES, --output_mode EMIT_ALL_SITES, --AllSitesPLs -R <hg19 reference FASTA> (github.com/smmarciniak/aDNA_osteo_height). Due to the potential for post-mortem damage impacting C>T and G>A allele changes, the per chromosome VCF files of the called genotypes were filtered for potential post-mortem damage (modifying https://github.com/ryhui/imputation-pipeline/pmd_filter.py for multi-individual VCF files and removing deamination in both directions). Potential deamination signals C>T (T>C) and G>A (A>G) genotypes were replaced in ancient individuals heterozygous for these genotypes with ‘./.’ 

Genotype likelihoods were then estimated using the per-chromosome VCF files, followed by imputation of missing SNPs based on the genotype probability score using the 1000 Genomes phase 3 haplotypes (http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/) and GRCh37 genomic maps (http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/). Parameters for estimating genotype likelihoods were: gprobs=true, gl=<input genotypes from UnifiedGenotyper> ref=<Beagle imputation reference panel>, map=<hg19 recombination map>. Imputation parameters were: gt=<GL_output_VCF>, gprobs=true, impute=true, ref=<Beagle imputation reference panel>, map=<GRCh37 recombination map> (github.com/smmarciniak/aDNA_osteo_height). This resulted in 30,761,499 markers imputed/genotyped across 167 individuals. Prior to downstream analyses, the imputed VCF was filtered for a minimum genotype probability of 0.99 to maximize confident genotype calls post-imputation and filter out less confident calls. We repeated the above pipeline of genotyping calling and imputation without filtering for potential deamination signals and the results were consistent with the deamination filtered data.