Low budget analysis of Direct-To-Consumer genomic testing familial data

Family SNP chip analysis

Whole family genotypes provide an additional genetic context that individual data analyses cannot offer, enabling enhanced error correction and inheritance state analysis¹⁰. We found that data downloaded from 23andMe at different times may vary, probably as a consequence of changes in genotype-calling algorithms. To ensure consistency in our analyses, we downloaded the most recent data, for all family members, on May 30^th 2012.

In order to facilitate comparison between different genotypes, we excluded non-autosomal SNP data (chromosomes X, Y and MT). All v3 chips had a total of 930,342 autosomal SNPs; the v2 SNP chip had 556,694. The reported ‘no call’ rate (shown as ‘- -’ in the downloaded data) varied slightly for each individual (Table 1), but overall, when genotypes are expressed as percentages of the total, some differences are observed for Son (v2) when compared to all other v3 individuals (Figure 3). V3 individuals show a very similar distribution of genotypes.

Figure 3. A distribution of all the different occurring genotypes as a percentage of the total for all individuals is shown.

For the purposes of unbiased comparison, only autosome data is included. I and D indicate insertion and deletion, respectively. Son’s percentages (v2) show slight differences to all other v3 individuals whose genotype proportions are more similar.

Calculation of error rates

23andMe reports a 98% or greater call rate¹¹, meaning that the chip can provide accurate data for more than 98% of those variants in any particular person. When an allele variant present in heterozygous state is "undercalled" (not observed), the locus may be reported as being homozygous for the other variant, leading to missed heterozygosity. Such sites may significantly impact the disease risks predicted for the individual. Under the simplifying assumption of uniform undercall probability, we estimated the number of heterozygous sites mistakenly reported as homozygous (Table 2). This means that for Son, 1 in every 400 sites is mistakenly called. For Father, 1 in every 200. Next, we analyzed Mendelian Inheritance Errors (MIEs). If at one site a reported genotype is ‘CC’ but the genotypes for both parents is ‘TT’, one possible explanation is that one of the parents is actually heterozygous ‘CT’ but was undercalled as ‘CC’, and likewise the son is heterozygous ‘CT’ but was undercalled ‘TT’. Given 5 people, there are 10 possible pairwise relations. Four of these represent direct parent/offspring relations, for which discrepancies can be counted as MIEs (Table 3).

One SNP in the Daughter (chr4:7, 9957, 622) is in disagreement with both parents: Father=CC, Mother=CC, Daughter=TT.

For the remaining 6 relationships, only part of the genome is expected to be identical by descent (IBD). Fully incompatible sites can be numbered as "pseudo-MIEs" (Table 4).

The reduced numbers of MIEs and pseudo-MIEs between the Son and all other family members were due to the lower total number of SNPs assayed for the Son. The Daughter has more MIEs relative to the Mother than relative to the Father. This may be due to having inherited a few deleted segments leading to a hemizygous state: hemizygous sites are reported as homozygous for the allele present, leading to an accumulation of apparent MIE sites. For example, at chr2:41093584 (rs12465519) (Table 5). The ISCA analysis explains how we inferred a deletion from observed discordant genotypes.

SNP similarity analyses

We compared the 23andMe genotypes of the family members at three levels: 1) pairwise all-against-all comparison, 2) inheritance state analysis in family quartets, and 3) analysis of population admixture.

1. Pairwise comparison. We performed an all-against-all genotype comparison among all family members. As an external point of reference, we included in our analysis a male individual of Indian ethnic background (denoted as non-CEU, chip version 3). Table 6 shows a summary of the SNP similarity patterns found among all family members and the non-CEU individual. When comparing the similarity of tested SNPs between family members (for information on how similarity scores are calculated see Methods section) we find that Son is most similar to Daughter (85.7%), then to Father (83.8%) and Mother (83.7%), and least to Aunt (78.7%) and the non-CEU individual (75.4%). Among individuals sharing the same platform, we find that siblings, Aunt and Mother share the greatest number of identical genotypes (84.5%). When comparing whose SNPs are the most similar to Father, neither Mother nor Aunt are found significantly different (p-value = 0.5561). The non-CEU individual has, on average, a similarity of 75.3% to all family members tested on platform v3. Contrary to our expectation that both Sister and Son should be equally similar to Mother and Father, we found that Daughter is 84.8% identical to Mother and 85.0% to Father, potentially explained by the inheritance of a deletion (see previous section). This difference Father/Daughter vs. Mather/Daughter similarity was found significant (p-value = 1.731e-05). Son however, exhibited 83.7% identity to Mother and 83.8% to Father and his tested SNPs were not found to be significantly different to either of them (p-value = 0.2751).

2. Inheritance state analysis. The availability of nuclear families with two or more offspring enables the identification of inheritance states¹⁰. These represent whether the offspring inherited the same alleles from both parents ("identical" state), the same allele from one parent but different alleles from the other ("haploidentical" state, maternal or paternal according to the parent from which the same allele was inherited), or different alleles from both parents ("nonidentical" state). Some family genotypes are consistent with all inheritance states (e.g. when all family members are homozygous) and are thus non-informative. Some family genotypes are consistent only with a subset of the inheritance states. For example if both parents are heterozygous A/G, and both offspring are homozygous A/A, clearly the offspring inherited the same alleles from both parents, which is consistent only with the "identical" state. Some family genotypes are consistent with two inheritance states. By combining the evidence from individual SNPs along a chromosome, it is possible to identify contiguous blocks of consistent inheritance, bounded by recombination events in either parent (Figure 3). The overall fraction of the genome present in each inheritance state (Figure 3 inset) deviates little from the expected 25%.

We found that 23andMe genotypes are sufficient for performing this analysis and identifying well-defined inheritance state blocks, despite covering a very small fraction of the genome (approx. half a million shared SNPs, limited by the lower density version used for Son). This is due to the largely uniform sampling of SNPs along the genome.

3. Admixture analysis. Visualization of admixture for Mother and Father was done with ADMIXTURE¹² in the context of other similar Southern European individuals. Figure 2 showed an admixture mapping for a selection of individuals from the Eurogenes Genetic Ancestry project¹³. Mother was denoted as ES7 (red arrow) and Father as ES8 (blue arrow). Mother and Father seemed to be markedly different yet in and around the Portuguese and Spanish cluster of individuals.

Combining SNP chip data and exome SNPs in Son for SNPedia annotation

To leverage all genotype data contained in all different sources, we combined the SNP chip data v2 (574,406 SNPs) with those found in Son’s exome data. 10,203 genotypes were annotated in SNPedia when exome SNPs were pooled with the SNP chip version 2 (generated on 11^th June, 2012). Processing just the exome, only 925 genotypes were found annotated in SNPedia. It is not surprising that so few additional SNPs were found, as the exome comprises a very small percentage of the total genome. Two SNPs were discovered to have conflicting genotypes between the two platforms: rs12344615 reported as ‘AG’ and ‘GG’ and rs2290272 reported as ‘CT’ and ‘TT’ respectively. The most likely informative SNPs from the exome data are summarized as judged by their observed frequency in HapMap⁸.

Exome data summary statistics

For the analysis of similarities between genotypes only 23andMe data was analyzed. Although exome sequencing is not widely marketed yet within the DTC providers as an option, this is likely to change in the near future. With our current budget constraints, we were able to sequence Son’s exome and relevant SNP and variation data was added for further analysis. A BAM file was created out of 4 FASTQ files downloaded from a server from the Beijing Genomics Institute. A compressed VCF file was also created. A total of 2.54 Gigabases of sequence was aligned at high quality. Summary metrics of the exome were calculated using Picard¹⁴ and showed a minimum of 61.42% of the on-target regions were covered with a depth at least 20x. Genotyping with GATK¹⁵ identified 37,702 variations relative to the reference genome (GRCh37). This was noted to be lower than expected if additional samples had been genotyped concurrently. 97% of the variants identified were within a known gene, 58% of them overlapped a protein domain and 5,565 were non-synonymous (15%) with serious predicted consequences on the protein product (as determined by SIFT¹⁶, PolyPhen¹⁷ or Condel¹⁸). Of these 5,565 potentially pathogenic SNPs, 413 had not been previously identified (verified against dbSNP release 132). This represents a normal figure for the private, novel, non-synonymous changes carried by most individuals. No more serious, novel changes were identified (such as stop codons gained or lost).

Visualization with SNPedia tools

We used SNPedia’s Hilbert curve visualization tool to compare chromosomes between different individuals. Figure 5 shows an example of a Hilbert curve comparison between chromosome 1 of Mother and all family members and the non-CEU individual. Each pixel corresponds to a SNP, colored according to four categories: match (light blue), half-match (dark blue), mismatch (red) and no data (grey). Two patches of light blue are apparent in the Mother-Aunt comparison of chromosome 1, and more appear in other chromosomes (not shown). These correspond to ‘identical’ segments in which the Mother and the Aunt inherited identical haplotypes from both their parents.

Figure 4. Inheritance State Consistency Analysis (ISCA) plot for the Father-Mother-Son-Daughter quartet, depicting for each autosome the number of informative SNPs supporting each of the four possible inheritance states: "identical" ("id", red), "haploidentical maternal" ("hm", green), "haploidentical paternal" ("hp", yellow) and "nonidentical" ("ni", blue).

SNPs consistent with two inheritance states contribute 0.5 weight to each. SNP counts are binned in non-overlapping 1 Mb windows; within each window, the four inheritance states are sorted by decreasing level of support. Regions without support typically overlap centromeric repeats and heterochromatic regions. Pie chart inset: fraction of the genome observed in each inheritance state.

Figure 5. A graphical representation of the SNPs from chromosome 1 for Mother com pared with herself, Father, Son, Sister, Aunt and non-CEU.

Each pixel represents a SNP. Light blue represents match, dark blue half-match and red conflict. SNPs in Son that are not present in the genotypes of the other individuals are represented in grey.

The next most similar graph can be seen to be the Mother/Daughter comparison. Much of the graph comparing Mother and Son is grey, representing SNPs present only on 23andMe v3 but absent from v2.

Inference of phenotypes using SNP data

We inferred phenotypes from the observed genotypes, by comparing to all available SNPedia SNP annotations. We analyzed family genotypes using the Promethease tool¹⁹, which allows annotation of observed genotypes from DTC analyses with SNPedia-annotated phenotype associations. These are collected directly from the scientific literature. All family members were analyzed using Promethease and associated phenotypes were collected for further interpretation. SNPedia’s annotation Magnitude, denoting a subjective degree of importance for an observed trait as judged by the curator, was used to discriminate between traits that should be considered further in our analyses. A magnitude of 0 denotes a common genotype for which no associated phenotypic data is known. Magnitude >3 is defined as ‘probably’ interesting. The maximum number is 10. In order to compare family phenotype annotations, all SNPs or ‘genosets’ (groups of SNPs) with equal or greater than magnitude 3 were extracted and summarized in Table 5. Results for Son are not directly comparable as he has fewer SNPs analyzed.

Based on observed results, the maternal hereditary line seems to carry greater risks related to diabetes and cardiovascular/thrombosis related conditions. In addition, both males (Son and Father) have greater risk of baldness as well as mixed results in terms of their ability to metabolize drugs. All the family shares a common trait of substantially increased odds of developing V617F-positive Myeloproliferative neoplasms.

SNP chip data

Five 23andMe genome analysis kits were purchased at two time points. Son’s kit was bought in May 2009 and was tested with 23andMe version 2 (~576,000 SNPs). The other 4 family members were analyzed in one batch, with kits bought in December 2010 and results returned in February 2011. The second batch used 23andMe’s chip version 3 with ~967,000 SNPs per genome analyzed. After discussion of results, consent was given by all family members of the family to publish their genotypes.

Calculation of error rates

Assume N SNPs were tested, and of these, fN are truly heterozygous. We wished to compute f (heterozygous fraction) from the observed numbers of homozygous, heterozygous and failed SNPs.

For each SNP with a dbSNP ‘rs’ identifier, we assumed that 1) the SNP is biallelic, 2) it is present in diploid state. We further assumed a probability x of not observing a given allele (the "undercall rate"), and assumed this probability was equal for both alleles, at all sites. Finally, we assumed that the probability of observing a wrong allele was zero ("overcall rate").

If the true state of the SNP was heterozygous, the following could have happened. 1) If neither allele was observed (double undercall), the SNP was called "NULL" (with conditional probability = x²). 2) If one allele was not observed (single undercall), the SNP was called "HOM" (conditional probability = 2x). 3) If both alleles were observed, the SNP was called "HET" (conditional probability = 1-2x-x2).

If the true state of the SNP was homozygous, there was only one type of allele to be observed. Thus, the following could have happened. 1) If the allele was not observed, the SNP was called "NULL" (conditional probability = x). Otherwise, 2) the SNP was called "HOM" (conditional probability = 1-x).

The expected frequencies for NULL, HET and HOM were:

Since HOM + HET + NULL = 1, the undercall rate x was given by solving:

x³ + (1+HOM)x² + (3HET+2HOM-3)x + NULL = 0

The heterozygous fraction f was then given by: f = HET/(1-2x-x²), and the number of "missing" heterozygous sites could be computed by: missing = (f-HET)N. Finally, the number of heterozygous sites reported as homozygous was given by:

het2hom = missing - N*NULL*HET/(HET+HOM)

Extraction of SNPs from exome data

Raw reads were aligned to the reference GRCh37 using bwa 0.61²¹. Local realignment was performed around indels with the Genome Analysis Toolkit (GATK v1.4)¹⁵ framework for variation discovery and genotyping using next-generation DNA sequencing data. Optical and PCR duplicates were marked in BAM files using Picard 1.62¹⁴. Original HiSeq base quality scores were recalibrated using GATK TableRecalibration and variants called with GATK UnifiedGenotyper. Indels and SNPs were hard-filtered according to Broad Institute best-practice guidelines²² to eliminate false positive calls and produce the final VCF.

Phenotype inference using promethease

Promethease, a SNPedia tool for phenotype inference, was used for assignment of SNPedia annotations to observed SNPs. SNPedia annotations contain manually curated SNPs that summarize phenotype associations observed in a particular population. Phenotype associations were inferred using SNPedia’s SNP ids, which correspond to dbSNP²³. In our analysis only SNPedia annotations of >= 3 Magnitude were examined. Magnitude is a subjective score that helps prioritize SNPs according to their expected importance and the phenotypic annotation itself in the form of free text. Magnitude is assigned by SNPedia entry curators.

Calculation of similarity scores

We calculated similarity scores using our own Perl scripts and MySQL. Similarity between any given two individuals was calculated as the total number of matching SNPs plus half-matches divided by 2. Negative results were counted as the total number of conflicts plus the number of half-matches divided by 2. As the platforms used are different for Son as compared to the rest, Son can only be compared relative to himself and not against other individuals. Similarities between individuals were statistically tested using the R package for Pearson’s Chi-squared with Yates’ correction.