A computational framework to assess genome-wide distribution of polymorphic human endogenous retrovirus-K In human populations

Human Endogenous Retrovirus type K (HERV-K) is the only HERV known to be insertionally polymorphic; not all individuals have a retrovirus at a specific genomic location. It is possible that HERV-Ks contribute to human disease because people differ in both number and genomic location of these retroviruses. Indeed viral transcripts, proteins, and antibody against HERV-K are detected in cancers, auto-immune, and neurodegenerative diseases. However, attempts to link a polymorphic HERV-K with any disease have been frustrated in part because population prevalence of HERV-K provirus at each polymorphic site is lacking and it is challenging to identify closely related elements such as HERV-K from short read sequence data. We present an integrated and computationally robust approach that uses whole genome short read data to determine the occupation status at all sites reported to contain a HERV-K provirus. Our method estimates the proportion of fixed length genomic sequence (k-mers) from whole genome sequence data matching a reference set of k-mers unique to each HERV-K locus and applies mixture model-based clustering of these values to account for low depth sequence data. Our analysis of 1000 Genomes Project Data (KGP) reveals numerous differences among the five KGP super-populations in the prevalence of individual and co-occurring HERV-K proviruses; we provide a visualization tool to easily depict the proportion of the KGP populations with any combination of polymorphic HERV-K provirus. Further, because HERV-K is insertionally polymorphic, the genome burden of known polymorphic HERV-K is variable in humans; this burden is lowest in East Asian (EAS) individuals. Our study identifies population-specific sequence variation for HERV-K proviruses at several loci. We expect these resources will advance research on HERV-K contributions to human diseases.


Introduction
Endogenous retroviruses (ERVs) are derived from infectious retroviruses that integrated into a host germ cell at some time in the evolutionary history of a species [1][2][3][4][5]. ERVs in humans (HERVs) comprise up to 8% of the genome and have contributed important functions to their host [6][7][8]. The infection events that resulted in the contemporary profile of HERVs occurred prior to emergence of modern humans so most HERVs are fixed in human populations and those of closely related primates. However some HERVs are still transcriptionally active and capable of causing new germline insertions so that individuals differ in the number and genomic location occupied by an ERV, a situation termed insertional polymorphism [9][10][11]. Among all families of HERVs, HERV-K is the only one known to be insertionally polymorphic in humans. However, HERV-K genomes are closely related and as with many repetitive elements, they are difficult to accurately assign to a genomic location using standard mapping approaches [12,13].
The DNA form of a retrovirus is called a provirus and minimally encodes the structural gag and env gene, and genes for a protease and polymerase, termed pol. Viral genes are flanked by long terminal repeats (5' or 3' LTR). While there are several HERV-K that are full length, none are infectious and most contain mutations or deletions that affect the open reading frames or truncate the virus. Further, the LTRs are substrates for homologous recombination, which deletes virus genes while retaining a single, or solo, LTR at the integration site [14][15][16]. Insertional polymorphism typically refers to the presence or absence of a retrovirus at a specific locus [17,18]. However an occupied site can contain a provirus in some individuals and a solo LTR in others and hence still display polymorphism. Thus HERV-K and other HERVs have contributed to genomic diversity in the global human population in several ways [19].
The presence of antibodies to HERV proteins or HERV transcripts has spurred a quest to determine if HERVs from multiple families have a role in either proliferative or degenerative diseases in humans [20][21][22][23][24][25][26]. Although there are known mechanisms by which a HERV can cause disease; for example, by inducing genome structural variation through recombination [27][28][29][30][31], affecting host gene expression [32], and inappropriate activation of an immune response by viral RNA or proteins [23], it has been difficult to establish an etiological role of a HERV in any disease. HERV-K specifically has been associated with breast and other cancers [3,[33][34][35][36][37], and autoimmune diseases, such as rheumatoid arthritis [38,39], multiple sclerosis [22,40] and systemic lupus erythematosus [8,22,41] without definitive evidence of causality or of specific loci involved. Recently, a HERV-K envelope protein was shown to recapitulate the clinical and histological lesions characterizing Amyotropic Lateral Sclerosis [42,43], providing an important mechanistic advance of a role for a HERV-K protein in a disease. Despite growing evidence for a contribution of HERV-K transcripts or proteins to the pathogenesis of human disease, it is difficult to distinguish among HERV-K loci to investigate potential roles and, in particular, to determine if a loci that is polymorphic for presence or absence of a provirus could be involved.
In this paper, we focus on characterizing the genomic distribution of known insertionally polymorphic HERV-K proviruses in the 1000 Genomes Project (KGP) data. We present a data-mining tool and a statistical framework that accommodates low depth whole genome sequence data characteristic of the KGP-and often patient-data to estimate the presence or absence of a provirus at all loci currently known to contain a HERV-K provirus. Using these data, we determine the number of known polymorphic HERV-K proviruses per genome because HERV-Ks can affect genomic stability [44] contributing to the pathogenesis of a disease. We also provide a tool to visualize HERV-K co-occurrence in global populations to facilitate exploration of synergy that might exist among specific polymorphic HERV-K in disease [45]. Our results provide a reference of global population diversity in HERV-K proviruses at all currently known polymorphic loci in the human genome and demonstrate that there are notable differences in the prevalence of HERV-Ks in different global populations and in the total number of HERV-Ks currently known to be polymorphic within a person's genome.

A model to estimate polymorphic HERV-K from whole genome sequence data
The goal of this research was to develop a computationally efficient and easy to use tool that could accurately report the status of all reported insertionally polymorphic HERV-Ks with coding potential (provirus) from whole genome sequence (WGS) data. We use the KGP database, which represents individuals in five super-populations and 26 populations, to establish the diversity in global populations at each known polymorphic HERV-K proviral locus and the total number of these polymorphic HERV-K in individual genomes to provide a foundation to study the role of HERV-K in human disease. Our reference set consists of all HERV-K sequences that are available in public databases and that can be unambiguously assigned a location in hg19. Sequences of HERV-K that are not present in hg19 but that were generated by PCR primers to the host flanking regions are included in the reference HERV-K set. From these HERV-K reference sequences, we generate a set of k-mers (see S2 Fig for optimizing k) that are unique to all HERV-Ks at each locus. The analysis of subject data starts with a data mining step that recovers all whole genome sequence reads that map to identified HERV-K elements in hg19. The rationale here is that polymorphic HERV-K that are not present in hg19 are greater than 80% homologous to those in the human reference genome and will map on existing elements. The recovered reads from a query WGS data set are then reduced to k-mers and mapped, requiring 100% match, to the reference set of k-mers (T), which represents all unique sites for HERV-K at each locus. The output is a ratio (n/T) of subject k-mers (n) that are 100% match to the reference k-mers (T) (see Methods for full details; the value of T for each HERV-K is in S1 Dataset:virus).
Our preliminary analysis of the KGP data demonstrated that our k-mer-based approach is sensitive to sequence depth; some HERV-K loci are represented by an almost continuous range of n/T values from 0-1 (S1 Fig), making presence/absence classification difficult. However, the majority of the KGP data is approximately 6x depth and thus to make use of this important resource, we developed a mixture model to statistically assign the n/T values from genomes to a cluster considering the sequence depth. K was optimized to 50 because this value improved our model computational efficiency and output ( Fig 1B, S1 Text, S2 Fig). The affect of sequence depth on n/T can be seen by comparing the sequence data of 28 individuals in the KGP data that have both low and high sequence depth data (Fig 1 shows a subset of eight individuals for clarity). If read depth is greater than 20, there is less dispersion of n/T values, most likely because more reads from the query WGS data are recovered from the mapped intervals. The states, 'provirus', 'solo LTR', and 'absent' are preliminarily assigned to each cluster based on the high depth data (data in Fig 1B used for description below). Individuals with n/T = 1 have the reference allele (represented by the yellow cluster of low depth data) and n/T = 0 (red cluster) indicates that the HERV-K is absent (no k-mers to unique sites in the HERV-K at this locus were recovered from mapped sequence reads). The k-mers derived from persons with low (green) and intermediate (blue) n/T values were mapped to the HERV-K reference for this locus to determine whether they localized only in the LTR (assign 'solo LTR' to green cluster) or in the coding region (assign 'provirus' to blue cluster) (S3 Fig).

Prevalence of polymorphic HERV-K in each KGP super-population
The WGS data of each individual in the KGP dataset were evaluated using our optimized analysis workflow. HERV-Ks on chrY were not considered. Twenty sites, omitting one at chr1:73594980 [see Methods] that have been reported to be polymorphic for presence/absence [10,11,34,46] were identified as polymorphic for a HERV-K provirus by our analysis (S1 Dataset:virus). Polymorphic HERV-Ks greater than 6 kbp in length cluster together in a phylogenetic analysis indicating that they are closely related (S4 Fig). The prevalence (proportion of individuals in a given population with a provirus present at a given locus) of the 20 polymorphic HERV-K proviruses varied from 0.9% to 99.5% when averaged across the entire KGP dataset (Table 1). However, there were notable differences in prevalence at each HERV-K site among the five super-populations (AFR, EAS, AMR, EUR, SAS; see Methods for key to abbreviations). Of the 20, the prevalence of seven polymorphic HERV-Ks was greater than 90% and the difference between populations with the lowest and highest prevalence was less than 6.5% (Table 1). There was 100% occupancy for six of the seven high prevalence polymorphic HERV-Ks (98.8% for the seventh), indicating that the rate of conversion to solo LTR is low for viruses at these sites (see S1 Text for occupancy and S2 Dataset:KGP(absence, solo, presence) for model prediction of solo LTR prevalence). Two polymorphic HERV-Ks had an overall prevalence of less than 10% in any population (Table 1) and were found in individuals of AFR origin; we found no evidence of a solo LTR at these two sites. Nine of the remaining 11 HERV-Ks are of interest because the difference between super-populations with the highest and lowest prevalence is between 28 and 80 percentage points (Table 1). Of note, for the three HERV-Ks with the largest difference among super-populations, the prevalence is lowest in EAS populations.
Individuals from African populations differ significantly from the other four super-populations in the prevalence of ten of the polymorphic HERV-K, three of which occur in close proximity on chr19. (Table 1, S2 Dataset:compare_prevalence). EUR and AFR super-populations are significantly different in the prevalence at all but one of the 20 polymorphic HERV-K based on adjusted p-values (S2 Dataset:compare_prevalence).

The number of polymorphic HERV-Ks per individual
The HERV-K genome is close to 10 kbp. As there are 20 known HERV-K loci with the potential to encode a provirus that are polymorphic in human populations, we asked if there is a difference in the burden of these repetitive, and potentially functional, viral elements among individuals. This was indeed the case. Of the 20 polymorphic HERV-K proviruses assessed, the number per person's genome ranges from 7-18 (  B) The result of the mixture model on the same data with k optimized to 50. The model returns four clusters each indicated by a unique color and eight of the 28 individuals that have both low and high depth sequence data are shown (see S1 Dataset:KGP for identification). The n/T ratio is 1 for persons with high depth data [red numbers, #6 and 12] who have the reference allele, while the corresponding low depth data [black numbers, yellow cluster] from the same individuals have n/T ranging from 0.7 to 0.9. There is less of an effect of sequence depth for individuals who do not have the HERV-K (n/T = 0, red cluster, #23 and 28). However optimizing k improves separation of the solo LTR (green cluster; #4 and #16) from the blue cluster (#25 and #14), which represents a state where some unique k-mers in the set T are missing in the query data (this is likely an allele; see S3 Fig   data suggest that a comprehensive investigation of polymorphic HERV-Ks may be a more productive means to advance studies of their potential disease impact.

Co-occurrence of polymorphic HERV-Ks
Our data provide a comprehensive picture of sites occupied by HERV-K provirus in each genome. Although most previous studies investigating a role of HERV-K in human disease assessed the prevalence of the HERV-K at a given locus, it is possible that, for example, two HERV-Ks each at 40% prevalence in a population rarely co-occur in an individual genome. By providing the status of all known polymorphic HERV-K in the genome, our tools facilitate such assessment and can advance investigation of HERV-K and human disease. We assessed combinations of three, four and five polymorphic HERV-Ks in KGP data and found that there are many combinations of co-occurring viruses that are population-specific (S3 Dataset). To facilitate exploration of HERV-K combinations among KGP populations, we developed a D3.j visualization tool (see Methods) that allows a user to choose any combination of the 20 polymorphic HERV-K proviruses and display the co-occurrence prevalence among the 26 populations represented in the KGP data. As an example, we show a combination of four HERV-Ks to represent the variation that occurs in KGP individuals, which in this case ranges from 3% in EAS to 59% in EUR (Fig 3A). We also determine that the three polymorphic HERV-Ks found on chr19 co-occur only from three AFR populations and in less than 2% of individuals ( Fig 3B).

KGP super-populations are distinguished by HERV-K status
Because there are clearly population-specific differences in both individual HERV-K prevalence and in the prevalence of HERV-K co-occurrence, we explored whether the presence or absence of these 20 documented polymorphic HERV-Ks is sufficient to distinguish populations using Fisher's linear discriminant analysis (LDA) [47]. Based on the status 'provirus', 'solo LTR', or 'absence', there is little resolution of AFR, EUR, and EAS super-populations ( Fig  4A). However, there is sufficient signature to separate AFR, EUR, and EAS if we utilize the n/T ratio of the 20 polymorphic HERV-Ks (S5 Fig) and we further improve population separation if we use the n/T ratio for all 96 HERV-Ks ( Fig 4B). This indicates that we are losing information by reducing the data to three states and that fixed HERV-K also contain signal for population of origin. An n/T = 1 indicates that the query set contains all k-mers that map to the reference set T for a specific HERV-K. If there is a HERV-K allele that has not been reported in any database but that is common in a population, we expect n/T <1 because we require 100% match to reference set T and k-mers covering allelic sites will be excluded (see Fig 1B, blue cluster for an example). We assessed the density distributions of n/T plots for each of the 96 HERV-Ks for evidence of population-specific alleles (S1 Text, S7 Fig). Five HERV-Ks have some indication of population specific distributions (S1 Dataset:virus). The HERV-K at chr1:155596457-155605636, which we report as fixed, is notable because the reference allele (n/T = 1) is only found in AFR (Fig 5A, S7 Fig). Individuals from other populations have n/T near 0.5. We mapped k-mers from individuals with n/T near 0.5 to the reference HERV-K sequences and confirmed that there is a loss of k-mers at several sites covered by the unique reference k-mers for this virus (S8 Fig). There are also cases where the reference allele is found in all populations except AFR (Fig 5B and see S7 Fig for additional examples).

Discussion
Our research provides a tool to mine whole genome sequence data to collectively evaluate the status of HERV-K provirus at known polymorphic and fixed sites in the human genome. The tool incorporates a statistical clustering algorithm to accommodate low depth sequence data and a visualization tool to explore the co-occurrence of known polymorphic HERV-K in the global populations represented in the KGP data. There are numerous significant differences in the prevalence of individual and co-occurring known polymorphic HERV-K among the five KGP super-populations. It is notable that individuals from EAS carry a lower total burden of the 20 polymorphic HERV-K than other represented populations. These data provide a comprehensive framework of genomic diversity among 20 documented polymorphic HERV-K proviruses to advance studies on potential roles for HERV-K in human disease, which have been alluring yet difficult to establish [21,22,24].
Tools developed to interrogate ERV insertional polymorphism typically exploit the unique signature created by the host-virus junction [11,48,49]. These approaches indicate that a site is occupied by an ERV but not whether there is a provirus associated with the site, which is more difficult to accomplish with short read sequence data. Our analysis tool provides an efficient means to detect occupancy and provirus status in one step. We decrease computational time by analyzing only the set of reads that map to existing HERV-K loci in the reference genome. This approach is justified because the known polymorphic HERV-K that are missing from the human reference are closely related to those in the reference genome assembly (see S4 Fig) and hence reads derived from them map to a related HERV-K in the reference. We employ k-mer counting methods, which also increase computational efficiency. A reference set of k-mers that is unique to each HERV-K is generated for each location in the genome and the proportion of reads (n/T) from the query set that maps to the k-mer reference set is reported as a continuous variable; there is no threshold of read count or coverage imposed for classification. Instead we utilize a mixture model to statistically cluster values based on n/T and sequence depth and assign the same HERV-K status to all individuals in a cluster. Clusters representing n/T of 1 consist of individuals from whom all the unique k-mers identified in the HERV-K reference set were recovered from their mapped WGS data. We classify other clusters by determining if k-mers mapped on the reference allele are distributed at sites in the coding portion of the genome or only in the LTR; reads mapping only in the LTRs are classified as solo LTR. This approach demonstrated that the k-mers derived from some individuals only covered a subset of the unique sites and led to the interesting finding that several HERV-K loci could have population specific alleles.
Wildschutte et al [11] have conducted the most comprehensive study of HERV-K prevalence in the KGP data to date. The goal of that paper was to identify new polymorphic insertions, either provirus or solo LTR, based on detecting reads containing the host virus junction. However, they implemented an additional step to detect provirus and provide the prevalence of some polymorphic HERV-K provirus for comparison with our results (see S1 Dataset:virus for comparison of prevalence values reported in Wildschutte et al [11]). There are five HERV-K previously reported in Subramanian et al 2011 [10] that were not included in Wildschutte et al [11]; all are polymorphic in our analysis (range 43-99%, see Table 1 and S1 Dataset:virus-column N). Seven polymorphic HERV-K, which Wildschutte et al [11] indicate occur in greater than 98% of KGP individuals, are fixed in our study. Our estimated prevalence for 14 HERV-K differs from that reported in Wildschutte et al [11] by 5% or more. Of these 14, the prevalence estimates at chr1:155596457-155605636 are most divergent. Our data show this site is fixed for provirus and Wildschutte et al [11] report that only 14% of the KGP data, all from AFR, have a HERV-K provirus integration. Our plots for chr1:155596457-155605636 show that AFR individuals carry the reference allele at this site (n/T near 1, Fig 5A) and all other individuals have n/T near 0.5. The k-mers from individuals with low n/T values for chr1:155596457-155605636 map to only a subset of sites marked by unique k-mers in the coding region (S8 Fig), which is consistent with sequence polymorphism or a deletion at these positions. The reference set T is small for this HERV-K and therefore overall coverage of the genome is low. Because Wildschutte et al [11] used a minimum coverage threshold for their kmer mapping method, it is possible that alleles present in non-AFR populations do not meet their inclusion criteria. There is a similar signal for alleles, represented by lower n/T values, at the other 13 HERV-K sites although the differences between our prevalence estimates and those of Wildschutte et al [11] are small (S1 Dataset:virus). In most cases these putative alleles are found in all populations at different frequencies but in five there is some degree of population specificity (Fig 5, S7 Fig, S1 Dataset:virus). Our results indicate that there could be considerably more sequence variation in HERV-K among human populations than previously appreciated. These data also suggest that using a HERV-K consensus sequence to study pathogenic potential could miss important features of HERV-K proviral polymorphism, which can be characterized by both the site occupancy status (presence/absence) and, when present, by sequence differences among individuals.
HERV-Ks are the youngest family of endogenous retroviruses in humans and consequently they share considerable sequence identity. This has the effect of limiting the number of unique sites associated with some HERV-K, which decreases the size of the reference set T (S1 Dataset:virus). The set T is small for near identical HERV-K such as HERV-Ks involved in a duplication event. The HERV-Ks at chr1:13458305-13467826 and chr1:13678850-13688242 are identical and cannot be distinguished. We report n/T for only one of these HERV-K (see S1 Dataset:virus, column M). We treat the two HERV-K proviruses spanning chr7:4622057-4640031 as a single virus with n/T = 1 reflecting the tandem arrangement found in the hg19. In this case, n/T<1 can mean either that both proviruses are present but with substitutions at a unique k-mer site or that one provirus converted to a solo LTR. Thus although an n/T ratio of 0 or 1 reliably indicates absence and presence of reference HERV-Ks, respectively, when T is small, sequence polymorphism and a deletion event can be difficult to distinguish from a solo LTR. However, because our mixture model statistically clusters similar n/T values based on sequence depth, all individuals in a cluster have the same status (e.g allele or solo LTR) even if we do not know what that state is. The ability of our tools to resolve the status of closely related HERV-K provirus sequences will improve as more empirical sequence data becomes available.
Our approach provides researchers with a rapid means to determine if the prevalence, and overall burden of the 96 HERV-K proviruses evaluated differ between a patient data set and the population represented in KGP to which they trace ancestry. The visualization tool will facilitate investigation of combinations of HERV-Ks in certain clinical conditions. The potential that HERV-K has multiple allelic forms in different populations is worthy of further analysis because a sequence allele could also contribute to a disease condition.

HERV-K proviruses
The 96 HERV-K proviruses previously reported [10,11,34,46] were supplemented with HERV-K alleles present in the NCBI nt database (November 2016 release) (92 in hg19, and 4 from the NCBI nt database). We required that any allele of a HERV-K from the nt database have at least 2kb of hg19 reference-matching host flanking sequence to confirm genome location. In total, 234 alleles were collected at the 96 known HERV-K loci. The location information and virus features are summarized in S1 Dataset: virus.

Developing a k-mer based detection model
We identified the k-mers that correspond to unique sequence characterizing each HERV-K. Kmers are substrings (subsequences) of length k that exist in a string (DNA sequence). The length k is determined empirically (S1 Text). Each k-mer is labeled with the corresponding viruses in which it is observed.
Only those k-mers referring to a single virus locus, unique k-mers, are selected for the set T. Where multiple alleles of a HERV-K are available, k-mers unique to all alleles at that location populations. In this case, all populations except AFR have the reference allele and all super-populations have an alternative allele that is not present in our reference set.
https://doi.org/10.1371/journal.pcbi.1006564.g005 comprise T. Multiple 2bps different k-mers (such as SNPs) corresponding to the same location on the virus, are merged into a single entry for the purposes of computing T. We map unique k-mers back to the corresponding alleles to determine coverage of the HERV-K and whether k-mers are located in LTRs (S3 Fig; S1 Dataset: virus).

Analysis of 1000 genome project (KGP) data
To develop a method to recover sequences containing information on HERV-K we leverage the fact that HERV-Ks are closely related. Thus, most sequence reads obtained from an individual with a polymorphic HERV-K that is absent in the human reference, hg19, will map to the location of a closely related HERV-K that is present the human genome reference. (As we show in S4 Fig, the known polymorphic HERV-K proviruses are closely related.) A file with the coordinates for all reported HERV-K insertions is used to extract mapped reads from a genome sequence file (S1 Dataset:bed, which provides the coordinates for both hg19 and hg38). Note that the KGP data were mapped to GRCh37, which includes the decoy sequence hs37d5. This decoy contains the HERV-K at chr1:73594980_73595948, which is not present in hg19. Thus, we did not recover any reads for this HERV-K, which is polymorphic but reportedly at high prevalence in most populations [11].
Our computational framework to indicate the status of each known HERV-K provirus is based on the n/T ratio, which is the proportion of k-mers in the data mined from WGS of each individual that are identical to the reference set T for each HERV-K provirus. Sequence reads are extracted from a mapped file of whole human genome sequence data based on coordinates corresponding to each annotated HERV-K. The reads are k-merized and mapped to the set T, which represents all unique k-mers assigned to each HERV-K in the reference set. We use exact match to map the k-mer data set to the unique k-mer references. The n/T ratio is an indicator of the presence of each HERV-K; n/T = 1 indicates that the individual has the HERV-K in our reference dataset documented to be at that locus while n/T = 0 indicates that no k-mers unique to a HERV-K locus were recovered (see Fig 1 for more explanation). Using a hash table (S1 Text), it takes 15 minutes to generate the n/T matrix for 100 files. The source code for the entire process is at https://github.com/lwl1112/polymorphicHERV

Dirichlet process Gaussian mixture model (DPGMM)
We utilized a statistical model to account for the dependency of the number of k-mers obtained from a person's sequence data (denoted by n ik for the ith subject and kth HERV-K, with i = 1,. . .,I,k = 1,. . .,96) that maps to the reference set T for each HERV-K on sequencing depth. Thus for each HERV-K we could statistically cluster those n ik /T values for i = 1,. . .,I based on the sequence depth of the WGS data for each individual for subsequent biological classification (provirus, solo LTR, absence, see Fig 1). More specifically in our analysis, for each k HERV-K, k = 1,. . .,96, consider a sample of size I measurements x i (i = 1:I), where each x i is a vector of length 2 x i = (x i1 ,x i2 ) with x i1 being the n ik /T measurement and x i2 the log function of depth. Here, for notation simplification, we use x i instead of x ik . To perform clustering analysis, we utilize the mixture model approach, which is arguably the most widely used statistical method for clustering. Specifically, we follow the work proposed by Lin et al. [52] that employs a Gaussian Mixture Model (GMM) with density function given by . N(μ j ,S j ) is the Gaussian density for the jth component parameterized by the 2-dimensional mean vector μ j and 2x2 covariance matrix S j . π {1:M} are the mixture components prior probabilities summing to 1. To allow a flexible modeling approach, we employ the standard Bayesian (truncated) Dirichlet Process prior for the parameters θ = (π j , μ j , S j , j = 1:M) [53,54]. The idea is that some of the mixture probabilities (π j ) can be zero, hence the actual number of mixture components needed may be smaller than the upper bound M. This mechanism allows automatic determination of the number of mixture components needed by the data set at hand.

Co-occurrence of polymorphic HERV-K
We consider that both the individual prevalence of a HERV-K and the co-occurrence of multiple HERV-Ks could differ among populations. The time of a brute-force approach for finding all combinations C m of size m from p polymorphic HERV-K is P p m¼1 p m À � ¼ 2 p À 1 À � , which is not efficient and is redundant. We employed the Apriori algorithm [55], which is commonly used for finding frequent pattern sets; in our case indicating which of the known polymorphic HERV-K frequently appear together. It first generates combinations C m (initialized to 1). In the optimization, frequent combinations F m are returned from candidates C m when prevalence exceeds the minimum threshold of co-occurrence. F m are then self-joined to generate combinations C m+1 of size m +1 and out of which F m+1 satisfy the minimum co-occurrence. In each pass, candidate combinations are pruned so as to avoid generating all combinations, which reduces running time significantly.

Statistical analysis of HERV-K frequencies across populations
We made statistical comparisons across 5 super-populations for the following three problems. For each problem, there are 5 2 À � = 10 families of 1-to-1 comparisons conducted. The 'prop-test' function in R is used to test whether the proportions for two super-populations are the same. 3. The co-occurrence for combinations of polymorphic HERV-K.
Therefore, multiple hypotheses would be conducted on frequencies F across super-populations P 1. . .5 as follows: Null hypothesis, H 0 : Alternative hypothesis, H A : F P i 6 ¼ F P j , where i6 ¼j.
A separate P-value is computed for each test and the Benjamini-Hochberg procedure [56] is used to account for multiple comparisons.

Visualization in D3.js
We utilized D3.js (Data Driven Documents) [57], an open-source java script library to create an interactive visualization to display co-occurrence of polymorphic HERV-Ks in human populations. Our visualization system includes two modules, a welcome page and a result page. Input JSON data include locations of polymorphic HERV-K, population information, and the 0/1 (absence / presence) matrix. (See S1 Text). Source code is available at: https://github.com/ lwl1112/polymorphicHERV/tree/master/visualization and a searchable tool with the data reported here is at: http://pages.iu.edu/~wli6/visualization/ All k-mers derived from the data mining step from each individual are mapped to the reference set of unique kmers, T, requiring 100% identity, to generate the set 'n' The first row shows the coverage of the set T on the HERV-K. The following plots show the mapping of the k-mer set 'n' from 8 individuals for the HERV-K at chr12: 55727215. # 6, 12, 14, and 25 (see S1 Dataset: KGP, column D for identification information) are labeled as 'provirus'. Note the drop out of the peaks near 3500 and 5000bp for #14 and #25, which accounts for a decrease in n/T in these individuals. #4 and 16 have low n/T and k-mers map to the LTR region indicated above the diagram; these are labeled as 'solo LTR'. #23, and 28 are labeled as 'absent'. For individuals with states 'solo LTR' and 'absent', there are some peaks in the coding region. This is most likely the result of assigning unique k-mers to this HERV-K that are shared with those from a HERV-K that is absent from the reference HERV-K dataset. (TIF)

S4 Fig. Maximum likelihood phylogenetic tree of fixed and polymorphic HERV-K.
To improve the alignment, only > = 6,500 bp HERV-Ks were included except for the HERV-K at chr1:75,842,771, which has a long deletion but aligns well in other regions. Maximum likelihood tree was generated using PhyML [4] using GTR with a gamma distribution. Node support was calculated using the alpha likelihood ratio test. Nodes with less than 0.9 alpha likelihood ratio test support were collapsed and colored in grey. HERV-K taxa are named after their genomic location in hg19. Polymorphic HERV-Ks identified in this study are indicated in red text. The chr8:146086169 HERV-K was identified in one individual in Wildschutte et al [5] but not found in this analysis. There is improved resolution of EAS from EUR and AFR using n/T compared to reducing the data to the three states 'provirus', 'solo LTR', 'absent' (Fig 4) for these 20 HERV-Ks. However, there is still substantial overlap of EUR and AFR based on n/T of the 20 polymorphic HERV-K studied. We assessed the density plots of all 96 HERV-K to determine if any peaks were specific to one of the super-populations. Shown are examples of candidate alleles specific to a population. In others several or all populations have the alleles but the prevalence is skewed. For example, the candidate allele for chr3:112743479-112752282 (the peak near n/T~0.7) appears to be more common in SAS individuals (pink trace). Similarly, EAS individuals (green trace) have a lower prevalence of the chr12:58721242-58730698 reference allele (n/T peak near 1) than do EUR (blue trace). Population-specific variation in HERV-K sequence could lead to under-estimation of proviral prevalence with mapping methods that require a coverage threshold.