GBLUP naive
BS
2025-04-23
GBLUP (naive) Анализ племенной ценности на основе линейной несмещенной модели предсказания с использованием генотипа (наивный пример)
Для начала создадим искусственные данные для демонстрации анализа. Используем следующие параметры:
3 параметра среды (например):
- Температура (°C)
- Количество осадков (мм)
- Уровень азота в почве (ppm)
3 фенотипических признака (например):
- Урожайность (кг/га)
- Высота растения (см)
- Содержание белка в зерне (%)
10 полиморфизмов SNP (биаллельные маркеры: 0, 1, 2). В данном примере, мы предположим, что 0, 1, 2 кодируют количество альтернативных аллелей у особи:
- 0: Гомозигота по референсному аллелю (например, AA).
- 1: Гетерозигота (например, Aa).
- 2: Гомозигота по альтернативному аллелю (например, aa).
set.seed(42) # Для воспроизводимости
# Параметры среды (для 50 генотипов)
n <- 50
environment <- data.frame(
Temperature = rnorm(n, mean = 20, sd = 5),
Precipitation = rnorm(n, mean = 100, sd = 30),
Soil_Nitrogen = rnorm(n, mean = 50, sd = 10)
)
head(environment)## Temperature Precipitation Soil_Nitrogen
## 1 26.85479 109.65776 62.00965
## 2 17.17651 76.48483 60.44751
## 3 21.81564 147.27183 39.96791
## 4 23.16431 119.28698 68.48482
## 5 22.02134 102.69282 43.33227
## 6 19.46938 108.29652 51.05514# Генотипы (10 SNP, 0/1/2)
snp <- matrix(sample(0:2, size = 10*n, replace = TRUE), nrow = n, ncol = 10)
colnames(snp) <- paste0("SNP", 1:10)
head(snp)## SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10
## [1,] 0 2 2 1 2 2 1 2 1 1
## [2,] 1 0 0 0 2 0 1 2 1 0
## [3,] 0 2 2 2 2 2 0 1 1 1
## [4,] 0 2 1 0 1 2 1 1 2 2
## [5,] 1 0 2 1 1 2 2 2 0 0
## [6,] 2 1 0 1 0 0 0 0 2 2# Фенотипы (с добавлением шума)
genetic_effect <- rowSums(snp) * 0.5 # Пример генетического эффекта
phenotype <- data.frame(
Yield = 50 + 0.8*environment$Temperature - 0.2*environment$Precipitation +
0.5*environment$Soil_Nitrogen + genetic_effect + rnorm(n, sd=2),
Height = 100 + 0.5*environment$Temperature + genetic_effect*0.3 + rnorm(n, sd=3),
Protein = 10 + 0.1*environment$Soil_Nitrogen + genetic_effect*0.2 + rnorm(n, sd=1)
)
head(phenotype)## Yield Height Protein
## 1 86.76015 116.4766 18.12844
## 2 81.58387 107.0911 16.04479
## 3 65.63853 115.8771 14.85348
## 4 81.96611 106.5557 16.81530
## 5 73.80500 110.1956 13.07007
## 6 74.71550 112.4060 16.76385Итак, мы сгенерировали 50 образцов (например, растений но могли бы быть и животные):
- 50 измерений среды: Каждая строка в environment описывает условия для соответствующего образца (температура, осадки и т.д.).
- 50 измерений фенотипов: Каждая строка в phenotype соответствует одному образцу (напр., урожайность растения №1, №2, …, №50).
- 50x10 SNP: Каждая строка в snp — это генотип одного образца по 10 маркерам (SNP1–SNP10).
Важно: Все данные выровнены по образцам (например, строка 5 в phenotype, environment и snp относится к одному и тому же растению).
Также нужно пояснить, что в коде выше:
genetic_effect <- rowSums(snp) * 0.5 — суммарный
эффект всех SNP для каждого образца (предполагается, что каждый SNP
вносит вклад в признак). Этот эффект добавлен в фенотипы (например,
Yield = ... + genetic_effect + шум).
Следующий шаг - рассчет матрицы геномных отношений между особями на основе SNP-маркеров. Эта матрица используется в GBLUP для учета генетической схожести между особями.
G-матрица отражает степень генетического родства между особями и вычисляется по формуле (VanRaden, 2008):
\[ G = \frac{(M - P)(M - P)^T}{2 \sum p_j(1 - p_j)} \] где:
M — матрица генотипов (0, 1, 2), Это исходная матрица генотипов размером n×m, где: n — число особей (например, 50 растений), m — число SNP (например, 10 маркеров). Каждая ячейка содержит значение 0, 1 или 2 — количество альтернативных аллелей у особи в данном SNP.
P — матрица ожидаемых генотипов 2(pj), где (pj) — частота альтернативного аллеля для SNP (j), (pj(1 - pj)) — дисперсия аллельных частот.
Это матрица, где каждая ячейка содержит ожидаемое количество альтернативных аллелей для SNP при условии случайного скрещивания (равновесие Харди-Вайнберга). Рассчитывается как 2pⱼ, где pⱼ — частота альтернативного аллеля для SNP j в популяции.
Как считать pⱼ:
Для SNP j: \[ p_j = \frac{\text{сумма альтернативных аллелей у всех особей}}{2 \times \text{число особей}} \]
Показатель дисперсия аллельных частот: pⱼ(1 - pⱼ)
Это генетическая дисперсия для SNP j, которая показывает, насколько разнообразны аллели в популяции. Максимальна при (pj = 0.5) (максимальное разнообразие), Равна нулю, если (pj = 0) или (pj = 1) (все особи одинаковы по этому SNP).
Например, если (pj = 0.5), то дисперсия: \[ (0.5 \times (1 - 0.5) = 0.25). \] Далее вычисляем частоты и матрицы:
# Частота альтернативного аллеля для каждого SNP
p <- colMeans(snp) / 2 # т.к. 0,1,2 кодируют количество аллелей
p## SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10
## 0.48 0.47 0.53 0.47 0.53 0.52 0.55 0.49 0.54 0.49# Создаем матрицу P (ожидаемые генотипы)
P <- matrix(2 * p, nrow = n, ncol = 10, byrow = TRUE)
head(P)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98
## [2,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98
## [3,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98
## [4,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98
## [5,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98
## [6,] 0.96 0.94 1.06 0.94 1.06 1.04 1.1 0.98 1.08 0.98# Центрированная матрица: W = M - P
# Центрирование матрицы M - P убирает смещение, вызванное разной частотой аллелей в популяции.
W <- snp - P
head(W)## SNP1 SNP2 SNP3 SNP4 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10
## [1,] -0.96 1.06 0.94 0.06 0.94 0.96 -0.1 1.02 -0.08 0.02
## [2,] 0.04 -0.94 -1.06 -0.94 0.94 -1.04 -0.1 1.02 -0.08 -0.98
## [3,] -0.96 1.06 0.94 1.06 0.94 0.96 -1.1 0.02 -0.08 0.02
## [4,] -0.96 1.06 -0.06 -0.94 -0.06 0.96 -0.1 0.02 0.92 1.02
## [5,] 0.04 -0.94 0.94 0.06 -0.06 0.96 0.9 1.02 -1.08 -0.98
## [6,] 1.04 0.06 -1.06 0.06 -1.06 -1.04 -1.1 -0.98 0.92 1.02Пример расчета (для 3 растений и 1 SNP)
Дано: \[ (M = [0, 1, 2]^T) \] (pj= 0.5), \[ (P = [1, 1, 1]^T) \] Вычисляем \((M - P): [ W = M - P = [-1, 0, 1]^T ]\)
Числитель \((W \cdot W^T): [ W \cdot W^T = \begin{bmatrix} 1 & 0 & -1 \ 0 & 0 & 0 \ -1 & 0 & 1 \end{bmatrix} ]\)
Знаменатель \((2 \sum p_j(1 - p_j))\):
Для 1 SNP: \((2 \times 0.5 \times 0.5 = 0.5)\).
Итоговая G-матрица: \([ G = \frac{W \cdot W^T}{0.5} = \begin{bmatrix} 2 & 0 & -2 \ 0 & 0 & 0 \ -2 & 0 & 2 \end{bmatrix} ]\)
G-матрица позволяет:
- Учесть генетическое родство между особями
- Предсказать селекционную ценность (breeding value) на основе SNP
- Устранить влияние популяционной структуры
# Пример данных: 3 растения, 1 SNP
M <- matrix(c(0, 1, 2), nrow = 3, ncol = 1) # Матрица генотипов (0,1,2)
colnames(M) <- "SNP1"
rownames(M) <- c("Plant1", "Plant2", "Plant3")
# 1. Рассчет частоты альтернативного аллеля (p_j)
p <- colMeans(M) / 2 # p = (сумма альтернативных аллелей) / (2 * n)
cat("Частота альтернативного аллеля (p_j):", p, "\n")## Частота альтернативного аллеля (p_j): 0.5# 2. Матрица ожидаемых генотипов (P)
P <- matrix(2 * p, nrow = 3, ncol = 1, byrow = TRUE)
cat("\nМатрица P (ожидаемые генотипы):\n")##
## Матрица P (ожидаемые генотипы):print(P)## [,1]
## [1,] 1
## [2,] 1
## [3,] 1# 3. Центрированная матрица (W = M - P)
W <- M - P
cat("\nЦентрированная матрица W = M - P:\n")##
## Центрированная матрица W = M - P:print(W)## SNP1
## Plant1 -1
## Plant2 0
## Plant3 1# 4. Числитель G-матрицы (W %*% t(W))
numerator <- W %*% t(W)
cat("\nЧислитель (W * W^T):\n")##
## Числитель (W * W^T):print(numerator)## Plant1 Plant2 Plant3
## Plant1 1 0 -1
## Plant2 0 0 0
## Plant3 -1 0 1# 5. Знаменатель (2 * сумма p_j(1-p_j))
denominator <- 2 * sum(p * (1 - p))
cat("\nЗнаменатель:", denominator, "\n")##
## Знаменатель: 0.5# 6. Итоговая G-матрица
G <- numerator / denominator
cat("\nМатрица G:\n")##
## Матрица G:print(G)## Plant1 Plant2 Plant3
## Plant1 2 0 -2
## Plant2 0 0 0
## Plant3 -2 0 2heatmap(G, symm = TRUE, main = "Матрица геномного родства (G)")
G-матрица показывает:
Plant1 и Plant3 имеют максимальное “отрицательное родство” (-2), так как их генотипы противоположны (0 vs 2).
Plant2 не имеет родства ни с кем, так как его генотип (1) равен ожидаемому (1) в популяции.
Интерпретация: Значения матрицы отражают генетическую схожесть между растениями. Например, G[1,1] = 2 означает, что Plant1 имеет вдвое большую генетическую вариацию относительно популяции.
Какие ввыводы можно сделать на этом микропримере: 1. Интерпретация диагональных элементов (родство особи с самой собой): Значения на диагонали (например, G[1,1] = 2, G[3,3] = 2) отражают генетическую вариацию особи относительно популяции. Чем больше значение, тем сильнее геном данной особи отклоняется от среднего по популяции. В примере Plant1 и Plant3 имеют максимальную вариацию, так как их генотипы (0 и 2) максимально отличаются от ожидаемого значения (1).
Интерпретация недиагональных элементов (родство между разными особями): G[1,3] = -2 указывает на отрицательное генетическое родство между Plant1 и Plant3. Это означает, что их генотипы противоположны (0 vs 2), и они менее схожи, чем в среднем по популяции. Нулевые значения (например, G[1,2] = 0) говорят об отсутствии родства между особями в рамках учтенных SNP.
Случай Plant2: Все элементы, связанные с Plant2, равны нулю, потому что его генотип (1) совпадает с ожидаемым значением (1) в популяции. Это показывает, как центрирование данных (вычитание матрицы P) устраняет влияние среднего по популяции.
Важность множества SNP в реальных анализах: В примере использован 1 SNP, что привело к экстремальным значениям в G-матрице. В реальных исследованиях используют тысячи SNP, и матрица получается более “гладкой”, а отрицательные значения встречаются реже. Чем больше SNP, тем точнее оценивается генетическое родство.
Применение G-матрицы в GBLUP: G-матрица позволяет учесть генетическую схожесть между особями при прогнозировании их селекционной ценности (например, урожайности). Отрицательные значения могут указывать на то, что особи менее подходят для скрещивания, но их интерпретация требует осторожности (особенно при малом числе SNP).
Нормировка матрицы: Деление на сумму дисперсий аллелей \(2 \sum p_j(1-p_j)\) делает G-матрицу масштабируемой, что позволяет сравнивать результаты разных исследований.
Рекомендации:
- Главная цель - выбор пар с низким генетическим родством для того, чтобы снизить инбридинг и увеличить генетическое разнообразие.
Для достижения цели можно рекомендовать следующее:
- Скрещивать особи с отрицательными или близкими к нулю значениями в G-матрице.
Например, в предложеном микропримере:
- Plant1 и Plant3 имеют G[1,3] = -2 — их скрещивание может дать гетерозиготное потомство с повышенной устойчивостью.
- Plant2 можно скрещивать с любым, так как его родство с другими нейтрально (G[2,1] = 0, G[2,3] = 0).
Продолжим с первоначальным примером:
Итак, наша цель оценить вклад генетических (SNP) и средовых факторов в фенотипы, а также предсказать селекционную ценность (breeding value) особей.
Модель GBLUP имеет вид: \[ [ \mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{Z}\mathbf{g} + \mathbf{e}, ] \] где:
- (\(\mathbf{y}\)) — вектор фенотипов,
- (\(\mathbf{X}\)) — матрица фиксированных эффектов (среда),
- (\(\mathbf{b}\)) — коэффициенты средовых эффектов,
- (\(\mathbf{Z}\)) — матрица случайных эффектов (обычно единичная матрица),
- (\(\mathbf{g}\)) — геномные эффекты (нормальное распределение: \(\mathbf{g} \sim N(0, \mathbf{G}\sigma^2_g)\)),
- (\(\mathbf{e}\)) — ошибки (нормальное распределение: \(\mathbf{e} \sim N(0, \mathbf{I}\sigma^2_e)\)).
# установим пакет sommer
#install.packages("sommer")
library(sommer)## Loading required package: Matrix## Loading required package: MASS## Loading required package: crayonЧуть-чуть исправим данные:
# 1. Генерация данных (с сохранением порядка ID)
set.seed(42)
n <- 50
ids <- paste0("Plant", 1:n) # Создаем вектор ID заранее
# Параметры среды
environment <- data.frame(
ID = factor(ids, levels = ids), # Фиксируем порядок уровней
Temperature = rnorm(n, 20, 5),
Precipitation = rnorm(n, 100, 30),
Soil_Nitrogen = rnorm(n, 50, 10)
)
# Генотипы (сохраняем порядок ID)
snp <- matrix(
sample(0:2, 10*n, TRUE),
nrow = n, ncol = 10,
dimnames = list(ids, paste0("SNP", 1:10)) # Используем заранее созданные ID
)
# Фенотипы (сохраняем порядок ID)
genetic_effect <- rowSums(snp) * 0.5
phenotype <- data.frame(
ID = factor(ids, levels = ids), # Фиксируем порядок
Yield = 50 + 0.8*environment$Temperature - 0.2*environment$Precipitation +
0.5*environment$Soil_Nitrogen + genetic_effect + rnorm(n, sd=2),
Height = 100 + 0.5*environment$Temperature + genetic_effect*0.3 + rnorm(n, sd=3),
Protein = 10 + 0.1*environment$Soil_Nitrogen + genetic_effect*0.2 + rnorm(n, sd=1)
)
# 2. Объединение данных БЕЗ изменения порядка
data <- merge(phenotype, environment, by = "ID", sort = FALSE) # Ключевое исправление!
# 3. Расчет G-матрицы после объединения
p <- colSums(snp)/(2*n)
W <- snp - matrix(2*p, n, 10, byrow = TRUE)
G <- tcrossprod(W)/(2*sum(p*(1-p)))
# Убедимся, что порядок ID совпадает
rownames(G) <- colnames(G) <- data$ID
# Проверка
identical(rownames(G), as.character(data$ID)) # Теперь TRUE## [1] TRUE# 4. Построение модели
library(sommer)
model <- mmer(
Yield ~ Temperature + Precipitation + Soil_Nitrogen,
random = ~ vs(ID, Gu = G),
rcov = ~ units,
data = data
)## This function has been deprecated. Please start using 'mmes' and its auxiliary functions (e.g., 'vsm', 'usm', 'dsm', 'ism', etc.). This function will be no longer maintained.
## This function has been deprecated. Please start using 'mmes' and its auxiliary functions (e.g., 'vsm', 'usm', 'dsm', 'ism', etc.). This function will be no longer maintained.## iteration LogLik wall cpu(sec) restrained
## 1 45.9008 0:6:4 0 0
## 2 46.2765 0:6:4 0 0
## 3 46.3641 0:6:4 0 0
## 4 46.3709 0:6:4 0 0
## 5 46.371 0:6:4 0 0# Результаты
summary(model)## This function has been deprecated. Please start using 'mmes' and its auxiliary functions (e.g., 'vsm', 'usm', 'dsm', 'ism', etc.). This function will be no longer maintained.## $groups
## Yield
## u:ID 50
##
## $varcomp
## VarComp VarCompSE Zratio Constraint
## u:ID.Yield-Yield 1.633952 1.0649190 1.534344 Positive
## units.Yield-Yield 3.855194 0.9059327 4.255497 Positive
##
## $betas
## Trait Effect Estimate Std.Error t.value
## 1 Yield (Intercept) 52.6212745 2.35075188 22.38487
## 2 Yield Temperature 0.8596639 0.05448487 15.77803
## 3 Yield Precipitation -0.2029211 0.01117551 -18.15766
## 4 Yield Soil_Nitrogen 0.5182167 0.03262717 15.88298
##
## $method
## [1] "NR"
##
## $logo
## logLik AIC BIC Method Converge
## Value 46.37096 -84.74191 -77.09382 NR TRUE
##
## attr(,"class")
## [1] "summary.mmer" "list"В выводе исеется предупреждение об устаревшей функции
Сообщение This function has been deprecated… указывает, что использована старая версия пакета sommer, где функция mmer() заменена на mmes(). Несмотря на это, модель отработала корректно и далле мы рассмотрим результат моделирования.
Моделирование показывает влияние различных эффектов:
Фиксированные эффекты (Fixed effects - $betas):
Показывают влияние параметров среды на урожайность (Yield):
- Intercept (52.62) — базовый уровень урожайности при нулевых значениях всех факторов среды (гипотетическая точка отсчета).
- Temperature (0.86) — увеличение температуры на 1°C повышает урожайность на 0.86 кг/га (при прочих равных).
- Precipitation (-0.20) — увеличение осадков на 1 мм снижает урожайность на 0.20 кг/га.
- Soil_Nitrogen (0.52) — увеличение азота в почве на 1 ppm повышает урожайность на 0.52 кг/га.
Стандартные ошибки (Std.Error) указывают на точность оценок. Например, для температуры оценка (0.86 Std.Error=0.05) означает, что истинный эффект, вероятно, лежит в диапазоне 0.81–0.94
Поскольку все значения |t.value| > 2, можно быть уверенным, что все факторы среды значимы .
Случайные эффекты (Random effects - $varcomp):
Отражают вклад генетических и неучтенных факторов в изменчивость урожайности:
ID (VarComp = 1.63) — генетическая дисперсия, объясняемая SNP (через G-матрицу). Чем выше значение, тем сильнее генетические различия между растениями влияют на урожайность. Residual (VarComp = 3.86) — остаточная дисперсия, которую нельзя объяснить ни генами, ни учтенными факторами среды (шум, погрешности измерений, неучтенные факторы).
Соотношение генетической и остаточной дисперсии позволяет оценить наследуемость признака: \[ h^2 = \frac{\text{Генетическая дисперсия}}{\text{Генетическая дисперсия} + \text{Остаточная дисперсия}} = \frac{1.634}{1.634 + 3.855} \approx 0.3 \] Это означает, что 30% изменчивости урожайности обусловлено генетическими факторами (SNP), остальное — средой и ошибкой.
Визуализируем селекционную ценность
# Генетические эффекты (breeding values)
breeding_values <- model$U$`u:ID`$Yield
# Преобразуем в вектор и создаем график
breeding_values <- unlist(breeding_values) # На случай, если это список
plot(breeding_values,
type = "h",
col = "blue",
xlab = "Особи (ID)",
ylab = "Селекционная ценность",
main = "Генетический потенциал урожайности")
То есть, для каждого растения мы теперь можем определить генетический потенциал урожайности, но что делать дальше?
Переходим к ключевым этапам: отбору селекционных пар и интерпретации генотипов. Далее разберем оба аспекта на примере наших данных.
Подбор селекционных пар
Для этой задачи следует использовать селекционную ценность (Breeding Values, BV)
Из модели GBLUP мы получили BV для каждой особи (например, model\(U\)u:ID$Yield). Это
генетический потенциал особи по признаку (урожайности).
Далее можно отобрать топ-10% особей, т.е. растения с наибольшими BV для скрещивания.
top_10 <- head(sort(breeding_values, decreasing = TRUE), n = round(0.1 * n))
cat("Лучшие особи для скрещивания:", names(top_10))## Лучшие особи для скрещивания: Plant25 Plant44 Plant7 Plant21 Plant48При этом нужно избегать инбридинга, т.е. выбирать пары с низким генетическим родством (на основе G-матрицы).
Вот пример для двух особей:
G["Plant25", "Plant44"] ## [1] 0.2116967# Значение близкое к 0 — хороший вариантG["Plant1", "Plant50"] ## [1] -0.4385662# Отрицательное значение — максимальное генетическое различиеЭтап оптимизации признаков
Далее, если нужно улучшить несколько признаков (например, Yield и Protein), можно использовать индекс селекции:
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Рассчитайте селекционные ценности для Protein (если еще не сделано)
protein_breeding_values <- data$Protein # или из модели GBLUP
# Масштабируем и преобразуем в вектор
yield_scaled <- as.vector(scale(breeding_values))
protein_scaled <- as.vector(scale(protein_breeding_values))
# Интегральный индекс
total_index <- 0.7*yield_scaled + 0.3*protein_scaled
# Ранжирование
ranked_data <- data.frame(
ID = data$ID,
Yield = breeding_values,
Protein = protein_breeding_values,
Index = total_index
) %>% arrange(desc(Index))
head(ranked_data)## ID Yield Protein Index
## Plant25 Plant25 2.469108 17.19685 1.607753
## Plant34 Plant34 1.629577 17.95521 1.310312
## Plant7 Plant7 1.804959 17.48540 1.308912
## Plant44 Plant44 1.975294 16.46500 1.192444
## Plant1 Plant1 1.029471 18.12844 1.022433
## Plant48 Plant48 1.789229 15.93651 0.984388%>% (pipe) — оператор из пакета dplyr, который передает результат предыдущей операции в следующую, формирую pipeline arrange(desc(Index)) сортирует данные по убыванию колонки Index.
Теперь на основе этих данных мы можем отобрать наиболее подходящих для селекции кандидатов. К примеру, мы можем отобрать топ 10% кандидатов для скрещивания:
top_selection <- ranked_data[1:round(0.1 * nrow(ranked_data)), ]
print(top_selection)## ID Yield Protein Index
## Plant25 Plant25 2.469108 17.19685 1.607753
## Plant34 Plant34 1.629577 17.95521 1.310312
## Plant7 Plant7 1.804959 17.48540 1.308912
## Plant44 Plant44 1.975294 16.46500 1.192444
## Plant1 Plant1 1.029471 18.12844 1.022433А чтобы избежать инбридинга, необходимо убедиться, что выбранные растения не слишком генетически близки. Для этого нужно использовать G-матрицу:
G[c("Plant25","Plant34", "Plant7", "Plant44", "Plant1"),c("Plant25","Plant34", "Plant7", "Plant44", "Plant1")]## Plant25 Plant34 Plant7 Plant44 Plant1
## Plant25 1.3556778 0.942238992 0.548870068 0.211696705 0.55689801
## Plant34 0.9422390 1.130895516 0.536828162 -0.001043632 0.34415767
## Plant7 0.5488701 0.536828162 1.146951391 0.006984305 0.75358247
## Plant44 0.2116967 -0.001043632 0.006984305 1.074699956 0.01501224
## Plant1 0.5568980 0.344157669 0.753582467 0.015012243 1.16300727Интерпретация: Значение -0.00104 указывает на низкое родство — это хорошая пара для скрещивания. Пороговое значение: для исключения рекомендовано избегать пар с родством > 0.5.
Проинтерпретируем данные:
print(top_selection)## ID Yield Protein Index
## Plant25 Plant25 2.469108 17.19685 1.607753
## Plant34 Plant34 1.629577 17.95521 1.310312
## Plant7 Plant7 1.804959 17.48540 1.308912
## Plant44 Plant44 1.975294 16.46500 1.192444
## Plant1 Plant1 1.029471 18.12844 1.022433Возьмем первую строку:
- Yield (2.47): Генетический потенциал урожайности на 2.47 кг/га выше среднего по популяции.
- Protein (17.20%): Содержание белка на 17.20% выше среднего.
- Index (1.60): Максимальный в выборке — растение идеально подходит для селекции.
В завершении, визуализируем наши результаты:
library(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:crayon':
##
## %+%# График селекционного индекса
ggplot(ranked_data, aes(x = Yield, y = Protein, color = Index)) +
geom_point(size = 3) +
geom_point(data = top_selection, color = "red", size = 8) + # Топ-кандидаты
scale_color_gradient(low = "blue", high = "red") +
labs(
title = "Распределение селекционного индекса",
x = "Селекционная ценность (Yield)",
y = "Селекционная ценность (Protein)"
)