(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



Detectability of runs of homozygosity is influenced by analysis parameters and population-specific demographic history [1]

['Gabriel A. A. Silva', 'College Of Forestry', 'Wildlife', 'Environment', 'Auburn University', 'Auburn', 'Alabama', 'United States Of America', 'Avril M. Harder', 'Kenneth B. Kirksey']

Date: 2024-12

For each simulated population, all analyses were based on 50 randomly sampled individuals. For the declining population, extremely low heterozygosity in two individuals precluded estimation of HMM transition probabilities via Viterbi training, and all downstream analyses for this demographic scenario were based on the remaining 48 individuals. Mean heterozygosity (± SD) varied across scenarios: large population mean = 3.17 x10 -5 ± 1.26 x 10 −5 ; bottlenecked population mean = 3.05 x 10 −5 ± 5.81 x 10 −6 ; small population mean = 2.87 x 10 −5 ± 1.22 x 10 −5 ; declining population mean = 2.44 x 10 −5 ± 8.81 x 10 −6 ( S17 Fig ). Mean F ROH generally increased with decreasing mean heterozygosity: large population mean = 0.21; bottlenecked population mean = 0.48; small population mean = 0.62; declining population mean = 0.65 ( Fig 1B , 1D , 1F , and 1H ). For all demographic scenarios, length-specific F ROH ( Methods ) mostly decreased with increasing ROH length with some notable exceptions in the bottlenecked and declining population scenarios ( Fig 1C , 1E , 1G , and 1I ). For the bottlenecked and declining population scenarios, a few individuals had high proportions of their genomes located in very long ROHs ( Fig 1G and 1I ), consistent with very recent inbreeding events in each population. Following downsampling and SNP filtering, the final mean coverage across all demographic scenarios was 4.87, 9.77, 14.69, and 29.20, and 44.99 for the 5X, 10X, 15X, 30X, and 50X VCF files, respectively.

ROH calling accuracy across methods.

Before comparing results from the three main methods—BCFtools Genotypes, BCFtools Likelihoods, and PLINK—we examined the effect of using either default or Viterbi-trained HMM transition probabilities on the results produced by each BCFtools method. Within each method, we compared overall F ROH values obtained when using default HMM transition probabilities (default values:--hw-to-az = 6.7 x 10−8,--az-to-hw = 5 x 10−9) against those calculated when using probabilities estimated via Viterbi training (mean values across all groups:--hw-to-az = 9.2 x 10−7,--az-to-hw = 3.8 x 10−6; S3 Table). (In practice, mean probabilities were calculated and applied for each demographic scenario and BCFtools method combination, but note the small magnitude of variation in Viterbi-trained probabilities across all groups relative to the difference between these probabilities and the program defaults.) The slope and intercept parameters never significantly differed between the two sets of results across all demographic scenarios and coverage levels, but for most comparisons, the slope and intercept parameters of the model built using Viterbi-trained results were closer to 1 and 0, respectively, than for the default values model (S4 Table and S1–S4 Figs). Using Viterbi-trained transition probabilities slightly increased false positive rates but also decreased false negative rates (S5–S8 Figs). Because of these slight, if not significant, improvements in overall F ROH estimation accuracy and false negative rates, we chose to proceed with Viterbi-trained HMM transition probabilities for both BCFtools methods and present only these results below (we also describe one additional benefit of Viterbi-trained probabilities below).

We used our simulated data set and linear models to determine whether each approach tends to over- or underestimate true F ROH . All three methods underestimated F ROH , with all model intercepts across coverage levels negative and different from zero, except for PLINK when calling ROHs in the large population at 5X coverage (i.e., no other 95% CIs for intercepts included zero; Fig 2 and S5 Table). Within each demographic scenario, the three methods produced similar results with respect to overall F ROH accuracy. For example, within all four scenarios, intercepts did not differ among methods or coverage levels, with two exceptions at 5X coverage when using PLINK: intercepts differed between this group and others within the large, declining, and the majority of coverage levels for small populations. Perhaps the most visually striking points of comparison across demographic scenarios are the slope estimates specific to each. Within each method, demographic scenarios were consistently ordered by slope estimates, with the declining population having the largest slope estimate (for BCFtools, all slopes > 1 with no 95% CIs including one), followed by the small, bottlenecked, and large populations (for the large population, all BCFtools slopes < 1, all PLINK slope 95% CIs include one; Fig 2 and S5 Table). Thus, accuracy of overall F ROH estimates calculated for the declining population varied the most with respect to true F ROH (i.e., of all demographic scenarios, slope parameter estimates for this population most greatly differed from one), followed by the large population. With respect to true F ROH , accuracy went up with increasing true F ROH for the declining population and went down with increasing true F ROH for the large population.

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 2. The relationship between true and inferred F ROH values depends on inference method and population demographic scenario. Each regression line represents linear model results for a single level of coverage with the shaded areas representing 95% confidence intervals. Each point represents data for a single simulated individual. Panels display outcomes using BCFtools in Genotypes mode (A-D), BCFtools in Likelihoods mode (E-H) and, PLINK (I-L), as well as by population scenarios including large (A, E, I), small (B, F, J), bottlenecked (C, G, K), and declining (D, H, L) populations. Dashed line is 1:1 line and x- and y-axes are consistent within each demographic scenario. Note the differing slopes across demographic scenarios (e.g., among panels A-D) and differing overall accuracies across methods (e.g., differing distances between regression lines and 1:1 line among panels D, H, and L). Another version of this figure with consistent axis limits across panels and colorized by sequencing depth is available in S19 Fig. https://doi.org/10.1371/journal.pcbi.1012566.g002

We calculated false negative (i.e., failing to call a ROH present in an individual) and false positive (i.e., calling a ROH that was not present in an individual) rates to further assess each method’s accuracy. With respect to false positive rates, PLINK performed poorly relative to the other two methods, with false positive rates consistently higher than BCFtools Genotypes and BCFtools Likelihoods (Fig 3A). Across all demographic scenarios and coverage levels, median false positive rates ranged from 6.14 x 10−7 to 0.031 for BCFtools Genotypes, from 4.54 x 10−8 to 0.007 for BCFtools Likelihoods, and from 0.008 to 0.112 for PLINK. For all three methods, increasing coverage to 10X corresponded to decreasing false positive rates, but rates did not continue to improve at coverages above 10X (i.e., 50% quantiles around median values overlap at 10-50X within each method and demographic scenario; results for all coverage levels provided in S9 and S10 Figs). Variation in false positive rates among samples at each coverage level was smallest for BCFtools Likelihoods, followed by BCFtools Genotypes, with PLINK showing the greatest variation across samples (summary statistics provided in S6 Table).

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 3. PLINK outperforms BCFtools with respect to false negative rates, but underperforms with respect to false positive rates. A) False positive (i.e., incorrectly calling a base position as being located in a ROH) and B) false negative (i.e., failing to identify a base position as being located in a ROH) rates across demographic scenarios and methods. Horizontal lines indicate median values and shaded boxes are 50% quantiles. Note the difference in scale of y-axis between panels A and B. Both BCFtools approaches outperform PLINK with respect to false positive rates but the reverse is true for false negative rates. Increasing coverage corresponds to decreasing false positive rates and to increasing false negative rates. Values displayed for 5X and 50X coverages; data for all coverage levels presented in S9 and S10 Figs, as well as a scatter plot in S18 Fig. https://doi.org/10.1371/journal.pcbi.1012566.g003

Patterns in false negative rates were generally in the opposite direction and magnitude to those we observed with false positives: both BCFtools methods performed poorly relative to PLINK, with BCFtools Genotypes producing slightly lower rates (median range across all demographic scenarios and coverage levels = 0.36–0.74) than BCFtools Likelihoods (median range = 0.49–0.78; Fig 3B). PLINK exhibited lower false negative rates than the BCFtools approaches (median range = 0.32–0.53) and less variation among samples at each coverage level. All three methods produced false negative rates that increased slightly at 10X coverage relative to 5X, with rates leveling off at 10X and higher coverage levels. Examples of false negative and false positive scenarios can be seen in S10 Fig, which illustrates a full chromosome of true and called ROHs for one exemplar individual from each demographic scenario.

We also examined how true and called values of F ROH varied for ROHs of different lengths. For the simulated data, all three methods almost always underestimated the proportion of the genome located in short ROHs across demographic scenarios, with the 95% CI less than zero for all tests other than PLINK at 5X coverage (Fig 4; all levels of coverage presented in S12–S15 Figs). For longer ROHs, results were variable and dependent on demographic scenario. For example, all three methods appeared to approach accuracy for longer ROHs in the large population, but this is likely due to the relative paucity of longer ROHs in this demographic scenario (Fig 1). For the bottlenecked and declining populations (wherein some individuals had very long ROHs), both PLINK and BCFtools Genotypes tended to overestimate F ROH for very long ROHs. PLINK generally produced the highest overestimates of F ROH and the most variation across samples of the three approaches, followed by BCFtools Genotypes. Finally, one coverage-related trend emerged across ROH length categories, methods, and demographic scenarios, with F ROH estimates calculated at 5X coverage often exceeding estimates calculated at higher coverage levels. Across all length bins and demographic scenarios, individual estimates of length-specific F ROH calculated at 5X were greater than those calculated at 10X for 37%, 53%, and 53% of BCFtools Likelihoods, BCFtools Genotypes, and PLINK estimates, respectively.

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 4. Increasing true ROH length corresponds to increasing detection. Called F ROH −True F ROH displayed by length bin (short, intermediate, long, very long) and demographic scenario (A: large population; B: small population; C: bottlenecked population; D: declining population) at 15X (results for all coverage levels presented in S9–S10 Figs). For BCFtools Genotypes and PLINK, F ROH for short ROHs is consistently underestimated whereas F ROH for very long ROHs is overestimated when these ROHs are present. BCFtools Likelihoods does not overestimate ROHs in any length bin. https://doi.org/10.1371/journal.pcbi.1012566.g004

To further investigate how called ROHs correspond to true ROHs, we identified regions of overlap between true and called ROHs within each individual and at each coverage level using a unique identifier for each true and called ROH. We found no instances of true ROHs being split into multiple called ROHs, but multiple true ROHs were often lumped together into a single called ROH. This pattern held true for all three methods at most coverage levels (Figs 5 and S16). For BCFtools Genotypes and PLINK, increasing coverage did not appear to ameliorate this problem (i.e., the mean number of true ROHs lumped into a single called ROH changed very little with increasing coverage; Fig 5C and 5D). However, for BCFtools Likelihoods, the number of true ROHs contained in a single called ROH decreased with increasing coverage, reaching a 1:1 ratio at 30X. Across methods, the mean number of true ROHs combined into a single called ROH increased with increasing ROH length, except for BCFtools Likelihoods at coverage levels ≥ 30X (S16 Fig). This lumping can be seen in Fig 5B.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1012566

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.