Evolution of genomic sequence inhomogeneity at mid-range scales

Substitution and polymorphism inside MRI regions

¿From the whole-genome human-chimp-macaque alignment we extracted all the aligned sequences with inhomogeneous nucleotide compositions that satisfy the criteria for MRI (so-called X-rich MRI regions; see the Materials and Methods section) and also control regions with nucleotide compositions equal to the average values for the entire human genome. We used the default MRI region length of 100 nucleotides for all computations. Only SNPs located within these MRI and control regions were studied. We particularly focused on the single nucleotide substitutions that maintain or erode MRI features. For example, in GT-rich MRI regions we counted the total number of novel polymorphisms that erode the feature, i.e. G or T → C or A substitutions, denoted as N _GT→CA and also the total number of those that maintain the MRI features, i.e. C or A → G or T substitutions, denoted as N _CA→GT . In addition, the entire set of recent human substitutions; that is, those nucleotides that differed in human but were the same in chimp and macaque, were processed for the MRI and control regions and presented as "fixed substitutions". The substitution ratio, S _X (recall that: S _X = N _X→nonX /N _nonX→X ) for the numbers of substitutions that maintain and erode X-rich MRI features was calculated for each substitution subtype (rare, minor, medium, and major SNPs and 'fixed'--refer to the Methods section for a detailed explanation) and presented in Figures 1 and 2. With respect to X-rich or poor MRI regions, the X in Figure 1 represents a two base combination such as GC, AG, GT, etc. while in Figure 2X can be any single nucleotide, e.g., A, T, C, and G. If the S _X -ratio is equal to 1 the X-rich region does not tend towards a change in its X-base composition. When S _X > 1, the substitutions reduce the X-richness of the examined regions, whereas when S _X < 1, substitution rates elevate the X-richness of the regions. Figures 1 and 2 demonstrate clear linear trends for S _X -ratios with respect to increasing fixation of substitutions within human populations.

For the cases of GT-, AC-, AG-, and TC-rich MRI regions (Figure 1), all S-ratios for rare SNPs are close to 1.8 (showing erosion of the MRI features). For major SNPs and fixed mutations the S _GT and S _AC -ratios reach 1.0 (which means no change in the corresponding base composition) and S _GA and S _CT -ratios reach 1.2 respectively. As for the corresponding control GT-, AC-, AG-, and TC-average regions (all having 50% of corresponding base composition) these lines are flat with all S-ratios equal to 1. The latter result is highly expected because of the symmetry of (+) and (-) chromosomal strands for these particular base compositions. Figure 1 also demonstrates that in GC-rich MRI regions the S _GC -ratio change has the highest slope from 7.0 for rare SNPs to 1.6 for fixed substitutions. In AT-rich MRI regions (also referred to as nonGC-rich in the tables) the change of S _AT -ratio has the lowest slope starting from 1.7 (rare SNPs) and ending at 1.3 (fixed substitutions). The control regions with the average GC/AT compositions (40-42% GC and 58-60% AT) also demonstrate a clear change of S-ratios during substitution fixation. In the control GC-average regions, rare SNPs favor increasing AT-richness (S _GC -ratio of 1.3) whereas fixed mutations demonstrate the opposite effect (S _GC -ratio of 0.8).

The data for the S-ratios for single nucleotides (Figure 2) are very similar to the trends seen in GC- and AT-rich regions. As expected from (+/-) strand symmetry, S _G -ratios are equal to S _C -ratios and represent about a half of the GC trend. The minor differences between G- and C-rich regions are within the errors of measurement. In the same way the S _A -ratios are seen to be the same as the S _T -ratios and they comprise approximately half of the effect seen for AT-rich regions.

In order to estimate mutation rates for MRI regions versus their respective control regions, we counted the occurrence rates for rare SNPs. The frequency ratio of rare SNPs in MRI rich regions to those in the control regions was calculated. The smallest ratio observed was for A+C content (0.464). This means the frequency of rare SNPs within MRI AC-rich regions is approximately half that of control regions. The highest occurrence ratio for rare SNPs was observed in G- and C-rich MRI regions (1.16 and 1.17 respecitvely). Thus, the occurrence rates of rare SNPs is slightly lower in MRI regions than in the corresponding control regions with the exception of G- and C-rich MRI regions. The entire dataset for the SNPs occurrences in MRI and control regions is presented in Additional file 1. The prevalence of rare and minor SNPs over major SNPs was also observed, their proportion over every MRI and control regions being 5.79.

Insertions and deletions inside MRI regions

Tables 2 and 3 demonstrate that in the human genome there is a prevalence of deletions over insertions (i.e. negative values of ΔX and ΔnonX) for every type of nucleotide content studied and for every type of MRI and control region with the exception of GC-indels in GC-rich MRI regions. In the last case ΔGC is positive and equal to 1405 added nucleotides (over a total set of 1.8 million nucleotides). For all other cases of X except X = GC, short and medium indels cause gradual contraction of genomic regions in humans. This means that there is no nucleotide composition equilibrium to which the indels drive the genome in the indefinite future and, therefore, these equilibria have not been calculated. Table 2 shows that, for every X-rich region, indels result in the increasing the richness of corresponding MRI regions (positive net X% change for X-rich region and negative net X% change for nonX-rich region). In all X-control regions the net X% change is several times less than in the corresponding X-rich and nonX-rich regions.

Finally, we calculated the percentage of nucleotide composition changes in case of both substitutions and indels separately, that occurred in the human genome during last ten million years after the divergence of human and chimpanzee. These results are presented in Tables 2 and 3 and serve to measure the relative importance of substitutions versus indels to the nucleotide composition of MRI regions.

Genomic samples and computation of recent human mutations ("fixed substitutions")

Taking human-chimp (human build 36.1 and chimp build 2 version 1) and human-macaque (macaque build v1 edit4) whole-genome pairwise alignments from the UCSC Genome Browser [21]http://hgdownload.cse.ucsc.edu/downloads.html as input, we generated a Perl script for the identification of the common genomic regions for these three species. The process involved the usage of the human genomic sequence as the reference for the location with the chimp and macaque sequences being extracted only in areas where the sequences of all three species were represented. We then invoked the ClustalW (v1.83) program with default parameters to obtain a whole-genome human-chimpanzee-macaque triple alignment. The obtained alignment is available at our website http://bpg.utoledo.edu/human_chimp_macaque.html. This triple alignment was used to calculate the dataset of recent mutations in humans. We considered a recent substitution at a particular position (for example T → C at position 23456719 on chromosome 7) to be valid if the human genome has a C base while both chimp and macaque have a T base in the corresponding aligned positions. In addition, we required that the quality of the alignment in the vicinity of the mutation be reliable (more than 70% similarity between human and macaque in the 20 bp flanking region [-10, +10]). The frequency table of all inferred recent human mutations is presented in the Additional file 3. We analyzed these recent substitutions together with the SNP datasets and call the former mutations "fixed substitutions," assuming that the majority of them occurred less than 10 million years ago and were already fixed across all human populations. In the same manner we processed indels in the triple alignments and computed all unambiguous cases of human insertions and deletions with sizes from 1 to 49 nucleotides.

Processing of SNP data

Over 4.62 million human SNPs from all chromosomes were obtained (dbSNP build 128 [22], ftp://ftp.ncbi.nih.gov/snp/), filtered for completeness and correctness annotations (676499 records discarded total), and mapped onto the whole-genome human-chimpanzee alignment. SNP allele frequencies were averaged from the frequency data of all populations of that allele. However, only those SNPs that were successfully located within the alignment were processed further. For each SNP site we verified the existence of the particular polymorphic bases in the specified position of the human genome reference sequence and also in the corresponding aligned position on the chimp genomic sequence. If any of these two species had different bases than the SNP alleles, the SNP was discarded (20469 SNPs discarded total).

Otherwise, we defined the origin of the polymorphism based on the chimpanzee nucleotide. Consider the following example to illustrate this process: suppose one has an A/G polymorphism located at position 34567812 of chromosome 5 with an average A allele frequency of 0.6 and a G allele frequency of 0.4. Then at position 34567812 of chromosome 5 of the human genome reference sequence (Genbank build 36.1), we would first examine if the A or G allele is present at that position and discard the SNP if not. Next, using the flanking region of that SNP we could align the chimp genomic sequence. If the chimp nucleotide were T or C then the SNP would also be discarded because those alleles are not a part of the human haplotype at that position. However supposing that the chimp nucleotide were G, then the polymorphism would be declared as a G → A polymorphism with G being declared the ancestral allele that at some point in human evolution mutated into an A allele within some human population(s). From the frequency data we may finally characterize this example SNP more precisely as a 0.4G → 0.6A polymorphism.

Using this approach we successfully characterized 3.93 million human SNPs. This last group of SNPs was divided into four subgroups based on the abundance of the mutant allele in the given human populations:

I. rare polymorphisms with the frequency of the mutated allele being less than 3%;

II. minor polymorphisms with frequencies ranging from 3% to 20%;

III. medium polymorphisms with frequencies going from 20% to 80%; and

IV. major polymorphisms with the frequency being above 80%.

For our method, misidentification of the ancestral allele might arise when the site for the human SNP is also polymorphic in chimp populations (e.g. A/G polymorphism) or for the possible case that this site had a recent substitution in chimps (A → G) after their divergence from humans. Human and chimpanzee genomes only differ by 1.23% due to single nucleotide substitutions with 1.06% being due to fixed substitutions and the rest (0.17%) being due to polymorphisms in human and chimp [20]. Moreover, according to the Chimpanzee Sequencing and Analysis Consortium the average estimated error rate of human alleles being misidentified due to chimp polymorphisms is only ~1.6% across all typical SNPs, which is acceptably low. It is also observed, however, that in the mutational hotspot of the CpG dinucleotide, there is an increased error rate for ancestral misidentification. If the human alleles are C/T in the CpG and TpG context and the chimp allele is T (in the TpG context) then the estimated error rate is actually 9.8% [20]. Thus, in the context of studying our MRI regions, any substitution (especially in GC-rich MRI regions since they contain an overabundance of G and C) going from TpG → CpG could have the ancestral allele misidentified, which would mean that the substitution would actually be CpG → TpG, although in the case of GC-rich MRI regions where such dinucleotides are more likely, an error rate of 9.8% is not sufficient to change the trend or conclusion of our results.

X-rich MRI genomic regions and control regions with average base composition

Any base or combination of bases can be described by a parameter X. For example, X could be G-base; C+T-bases; or A+T+G bases, et cetera. It is also useful to refer to nonX base(s) as all bases not X. Thus, X + nonX must represent all four nucleotides A, G, T, and C. For the examples above, nonX are A+T+C-bases; G+A-bases, and C-base, respectively. MRI is characterized by a specific base composition within a region under analysis. We characterize X-rich MRI regions based on an overabundance of the X base(s) within a region of a certain length (the so-called window), where the percentage of X should be above a certain threshold (Bechtel et al 2008). We calculated MRI regions in the human genome for single nucleotides and various nucleotide combinations using a stretchy window of 100+ nucleotides with the following threshold parameters: for A or T the threshold was 49%; for G or C we used 40%; for G+C it was 70%; for A+T the threshold was at 80%; for G+T, C+A, G+A, and C+T were at 70%; nonA or nonT was 87%; and non G or non C the threshold was 93%. These thresholds were chosen experimentally in such a way that MRI regions should represent about 2% of the whole human genome. A stretchy window of N + nucleotides means that we scan genomic sequence with an N-size window to find a genomic MRI region that fits the threshold criterion, then we extend the window above the detected region by 10 nt steps until the criterion is no longer met. After registering the full MRI region we jump beyond the current MRI region and continue with the default N-size window. Using this approach we characterized all MRI regions in the triple human-chimp-macaque alignments using the human sequence for calculating nucleotide composition and MRI features. We also discarded those MRI regions in the alignments where the indel composition exceeded 50%. For the collection of control regions with average base compositions we used the same stretchy window approach with the nucleotide composition corresponding to the following average genomic frequencies: for A, T between 30 and 31% thresholds; for G, C between 20-21%; for G+C between 40-42%; A+T at 58-60%; G+T, C+A, G+A, or C+T were at 49-51%.

Note that control regions with genome-average X-composition also have genome-averaged nonX-composition. Therefore, their subsitution ratios are in inverse proportion to each other: S _X = 1/S _nonX . Due to this only one ratio for X and nonX pair is shown in Figures 1 and 2.

Calculation of the substitution ratios in MRI and control regions

Studying SNPs and fixed substitutions in X-rich MRI regions we measured the number of changes from X to nonX (denoted as N _X→nonX ) and also the number of changes from nonX to X (denoted as N _nonX→X ). The fluctuations in the observed distribution of N _X→nonX and N _nonX→X are well-known as Poisson noise. Thus, the standard deviation for the true values for N _X→nonX and N _nonX→X is calculated according to v the Poisson distribution, that is: σ _N = . For each X-rich MRI region we measured the substitution S _X -ratios with S _X = N _X→nonX /N _nonX→X shown in Figures 1 and 2. The propagation of uncertainty for a ratio f = A/B can be calculated using the formula (σ _f /f)² = (σ _A /A)² + (σ _B /B)² - 2(σ _A ·σ _B )/(A·B)·ρ _AB , where ρ _AB is the correlation coefficient for A and B variables. Because the observed frequency of having a SNP at a genomic site in humans is less than 1%, it is correct to assume that the correlation between N _X→nonX and N _nonX→X is negligible. Therefore, the standard deviation for the S _X ratio was calculated by the following formula:

Calculation of base composition equilibrium for the observed substitution rates

As described in the previous paragraphs, for studying SNPs and fixed substitutions in X-rich MRI regions of the human genome we measured the number of changes from X to nonX (denoted as N _X→nonX ) and also the number of changes from nonX to X (denoted as N _nonX→X ). These N _X→nonX and N _nonX→X helped us to estimate the frequencies of these two types of mutations per X or nonX site, named here as F _X→nonX and F _nonX→X , correspondingly. Suppose one has a sample of MRI regions with a total nucleotide sequence length of L and a composition of X with the region richness given as P _X being measured in numbers from 0 to 1. Then, the total number of X sites in this sample will be L·P _X , and the total number of nonX sites will be L·(1 - P _X ). During a certain time interval called ΔT there will be ΔN _X→nonX and ΔN _nonX→X substitutions. Therefore the frequency of substitutions per site is F _X→nonX = ΔN _X→nonX /(ΔT·L·P _X ) and F _nonX→X = ΔN _nonX→X /(ΔT·L·(1 - P _X )). It is impossible to measure directly these ΔN values for a specific time interval of ΔT. However, with a good approximation we can assume that the frequencies are proportional to the observed numbers N _X→nonX and N _nonX→X and can be represented by the simple formula: F _X→nonX = A·N _X→nonX /P _X and F _nonX→X = A·N _nonX→X /(1 - P _X ), where A is a scaling factor having the same value for F _X→nonX and F _nonX→X , since N _X→nonX and N _nonX→X are counted from the same sample. In a gedanken experiment, let's assume that the current F _X→nonX and F _nonX→X values will stay unchangeable forever for our MRI sample. Then, in time, mutations should alter the base composition of our sample until it reaches an equilibrium composition with a new percentage for X-bases denoted here as Q _X . This equilibrium composition Q _X can be computed using the observed parameters of P _X , N _X→nonX , and N _nonX→X . Indeed, under the equilibrium, the number of changes from X to nonX must be equal to the number of reverse changes from nonX to X, or:

(1)

We can compute these ΔN _X→nonX and ΔN _nonX→X values from frequencies in such a way:

also in a similar way

By putting these transformations into Equation 1 we get:

or

(2)

Finally, simple transformation of Equation 2 gives us the final Equation 3 for calculation of equilibrium percentage:

In the Results section, Formula 3 is used to compute the equilibrium percentage for X-bases in the studied MRI regions.