Популяционная генетика

Исследование полногеномной ассоциации
Ищем SNP…
связанный с фенотипом.
Цель:
Объяснять
Понимание
Механизмы
Терапия
Предсказывать
Вмешательство
Профилактика
Так что же такое GWAS?

Определение
Любая связь между двумя измеренными величинами, которая делает их статистически зависимыми.
Наследственность
Доля дисперсии, объясняемая генетикой
P = G + E + G * E
Наследственность> 0
Ассоциация

Определение
Любая связь между двумя измеренными величинами, которая делает их статистически зависимыми.
Наследственность
Доля дисперсии, объясняемая генетикой
P = G + E + G * E
Наследственность> 0
Ассоциация

Почему?
Окружающая среда, взаимодействие генов и окружающей среды
Сложные черты, небольшие эффекты, редкие варианты
Уровни экспрессии генов
Методология GWAS?

Кейс-контроль
Четко определенный «случай»
Известная наследственность
Вариации
Количественные фенотипические данные
Например: Рост, концентрация биомаркеров
Явные модели
Например. Доминантный или рецессивный
Дизайн GWAS

Процесс
Визуализация
Филогенетика
PCA
Коррекция данных
Геномный контроль
Регрессия по основным компонентам PCA
Стратификация населения II

Что же там, в наших данных? Применим PCA -
Метод главных компонент
Для начала рассмотрим простые примеры
Новая ось PC1 соответствует направлению наибольшего изменения в данных.

Пример: планеты Солнечной системы
Как вы думаете, есть какая корреляция между колонками?

Пример: планеты Солнечной системы
Корреляция между колонками?
Три измерения явно излишни

Что показывает PCA?
Сколько компонент оставить? Два подхода
Все со стандартными отклонениями больше 1
Совокупная пропорция дисперсии больше 0.9

Результаты
Расстояние
Плотность
Диаметр
Saturn
PC1
PC2

Эффект нормализации: даем равный шанс разным группам измерений
С нормализацией 2 компоненты
Без нормализации 1 главная компонента

Теперь попробуйте сами
# загрузите файл с данными планет
# planets.csv
Planets=read.csv('Planets.csv', row.names = 1); Planets
cor(Planets)
pcaPlanets=princomp(Planets, cor=T)
summary(pcaPlanets)
biplot(pcaPlanets)
loadings(pcaPlanets)

Повторим упражнение с данными о странах мира
Данные
Вопросы:
Нужно ли проводить нормализацию?
Сколько главных компонент нужно сохранить?

Результаты

Простая карта
install.packages(“maptools”)
library(maptools)
data(wrld_simpl)
myCountries = wrld_simpl@data$NAME %in% row.names(Countries)
plot(wrld_simpl, col = c(gray(.80), "red")[myCountries+1])

Не совсем простая карта
MAX_INCOME=max(Countries$average.income)
# добавили колонку
Countries$index=round(10*Countries$average.income/MAX_INCOME); Countries
COL=rainbow(10) # задали цвета
ALL_COUNTRIES=as.data.frame(cbind( as.character(wrld_simpl@data$NAME ), gray(.80)))
head(ALL_COUNTRIES) # создали новый объект
names(ALL_COUNTRIES)=c('Name', "Color")
row.names(ALL_COUNTRIES)=ALL_COUNTRIES$Name
for(i in 1:dim(Countries)[1]){
ALL_COUNTRIES[row.names(Countries)[i],'Color']=COL[Countries[i,'index']]
}
head(ALL_COUNTRIES)# цвет зависит от дохода
unique(ALL_COUNTRIES$Color)# сколько категорий?
plot(wrld_simpl, col = ALL_COUNTRIES$Color) # нарисовали карту
# добавили легенду
sort(unique(Countries$index))*MAX_INCOME/10, COL[sort(unique(Countries$index))]
legend('bottomleft',
legend = paste("$",sort(unique(Countries$index))*MAX_INCOME/10),
col= COL[sort(unique(Countries$index))], pch=15)

K-MEANS
install.packages(“cluster”)
library(cluster)
Countries_clusters=kmeans(Countries[,1:3], centers=4, nstart=25)
clusplot(Countries, Countries_clusters$cluster,labels=3, color=TRUE)
ccs=data.frame(sapply(Countries[,1:3], scale)) ##нормализация
rownames(ccs)=rownames(Countries)
Countries_clusters_scaled=kmeans(ccs, centers=4, nstart=25)
clusplot(ccs, Countries_clusters_scaled$cluster,labels=3, color=TRUE)

Возвращаемся к нуклеотидам
Скачайте данные в вашу рабочую директорию
sativas413.ped
sativas413.fam
sativas413.map
sativas413.pheno
sativas413.csv
Мастерская GWAS

Библиотеки
#install packages
install.packages(c("poolr","qqman","BGLR","rrBLUP","DT", "dplyr"))
install.packages(c("rnaturalearth",'rnaturalearthdata','rgeos','ggspatial'))
devtools::install_github("dkahle/ggmap", ref = "tidyup")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SNPRelate")

Библиотеки
library(rrBLUP)
library(BGLR)
library(DT)
library(SNPRelate)
library(dplyr)
library(qqman)
library(poolr)
library(OpenStreetMap)
library(rjson)
library(rgdal)
library(RgoogleMaps)
library(mapproj)
library(sf)
library(OpenStreetMap)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

Шаг 1. Подготовка данных SNP в R
rm (список = ls ())
setwd («# ваш рабочий каталог, содержащий файл»)
Geno<- read_ped ("sativas413.ped")
head(Geno)
# объединяет данные маркерного аллеля в формате 0, 1, 2 и 3; 2 представляет отсутствующие данные
p = Geno$p; p
n = Geno$n; n
Geno = Geno$x; Geno
# Информация о присоединении
FAM <- read.table("sativas413.fam"); head( ​​FAM )
# Информация о положении снипов на геноме
MAP <- read.table("sativas413.map"); head( MAP )
# Перекодировать данные в файл ped
Geno[Geno == 2] <- NA # Преобразование отсутствующих данных в NA
Geno[Geno == 0] <- 0 # Преобразование 0 данных в 0
Geno[Geno == 1] <- 1 # Преобразование 1 в 1
Geno[Geno == 3] <- 2 # Преобразование 3 в 2
# Преобразование данных маркера в матрицу, транспонирование и проверка размерa
Geno <- matrix (Geno, nrow = p, ncol = n, byrow = TRUE)
Geno<- t (Geno)
dim(Geno)
Мастерская GWAS
0 00 Homozygote "1"/"1"
1 01 Heterozygote
2 10 Missing genotype
3 11 Homozygote "2"/"2"

Шаг 2. Считайте данные фенотипа в R
## прочитать фенотип
ris.pheno <- read.table ("sativas413.pheno", header = TRUE, stringsAsFactors = FALSE, sep = "\t")
# Просмотр первых нескольких столбцов и строк данных
ris.pheno [1: 5, 1: 5]
dim( ris.pheno)
# сравнить с фенотипическим файлом
rownames (Geno) <- FAM$V2; head(Geno)
table(rownames (Geno) == ris.pheno$NSFTVID)
# Теперь давайте извлечем первый фенотип припишем его объекту y
y <- matrix (ris.pheno$Flowering.time.at.Arkansas) # использовать первый фенотип
rownames (y) <- ris.pheno$NSFTVID
index<-!is.na (y)
y <- y [index, 1, drop = FALSE] # 374
Geno <- Geno [index,] # 374 x 36901
table(rownames (Geno) == rownames (y))
Мастерская GWAS

Шаг 3. Фильтрация данных SNP
# Здесь мы будем использовать цикл for, чтобы определить недостающие и преобразовать их в средние значения по столбцам.
for (j in 1: ncol (Geno)) {
Geno [, j] <- ifelse (is.na (Geno [, j]), mean (Geno [, j], na.rm = TRUE), Geno [, j])
}
# Фильтр редких аллелей. Здесь мы будем убирать аллели с частотой <5% и сохранять объект как Geno1
# Затем мы сопоставим и отбросим маркеры из файла информации о карте и сохраним его как MAP1
p <- colSums (Geno) / (2 * nrow (Geno))
maf <- ifelse (p> 0.5, 1 - p, p)
maf.index <- который (maf <0,05)
Geno1 <- Geno [, -maf.index]
dim(Geno1)
dim(Geno)
# подмножество на основе сохраненных SNP
MAP1 <- MAP [-maf.index,]; dim(MAP1)
Мастерская GWAS

Шаг 4. Структура популяции
# Создать файл геноматрицы и назначить имена строк и столбцов из файлов fam и map
Geno1 <- as.matrix (Geno1)
sample<- row.names (Geno1)
length(sample)
colnames(Geno1) <- MAP1 $ V2
snp.id <- colnames (Geno1)
length(snp.id)
snpgdsCreateGeno ("44k.gds", genmat = Geno1, sample.id = sample, snp.id = snp.id,
snp.chromosome = MAP1 $ V1, snp.position = MAP1 $ V4, snpfirstdim = FALSE)
# Теперь откройте файл 44k.gds
geno_44k <- snpgdsOpen ("44k.gds")
snpgdsSummary ("44k.gds")
Мастерская GWAS

Шаг 5. PCA
#PCA анализ
pca <- snpgdsPCA (geno_44k, snp.id = colnames (Geno1))
# график результатов PCA
pca <- data.frame (sample.id = row.names (Geno1), EV1 = pca $ eigenvect [, 1], EV2 = pca $ eigenvect [, 2], EV3 = pca $ eigenvect [, 3], EV4 = pca $ eigenvect [, 4], stringsAsFactors = FALSE)
plot(pca $ EV2, pca $ EV1, xlab = "собственный вектор 3", ylab = "собственный вектор 4")
# добавить информацию о населении на участок
pca_1 <- read.csv ("sativas413.csv", header = TRUE, stringsAsFactors = FALSE)
pca_2 <- pca_1 [match (pca $ sample.id, pca_1 $ NSFTV.ID),]
# Извлечь информацию о населении и добавить выходной файл PCA
pca_population <- cbind (pca_2 $ Sub.population, pca)
colnames (pca_population) [1] <- "Population"
# Постройте и добавьте названия населения
plot(pca_population $ EV1, pca_population $ EV2, xlab = "PC1", ylab = "PC2", col = c (1: 6) [factor (pca_population $ Population)])
legend(x = "right", legend= levels(factor (pca_population $Population)), col = c (1: 6), pch = 1, cex = 0.6)
Мастерская GWAS

Шаг 5а. Визуализация PCA
# Извлечь информацию о местоположении и добавить выходной файл PCA
pca_country <- as.data.frame(cbind(as.numeric(pca$EV1),as.numeric(pca$EV2),as.numeric(pca$EV3),
as.numeric(pca$EV4),as.numeric(pca_2$Latitude),
as.numeric(pca_2$Longitude)));head(pca_country)
names(pca_country)=c('PC1','PC2','PC3','PC4','LAT','LONG')
names(pca_country)
pca_country=na.omit(pca_country)
# корреляция с географией
cor (pca_country [, 1: 6], use= 'pairwise.complete.obs')
#map
map()
# определить цвета для карты
NCOLS = 1+max (unique (round (50 * (pca_country $ PC1 - min(pca_country $ PC1)))))
NCOLS
COL = rainbow(NCOLS)
pca_country$ind1=round (50 * (pca_country $ PC1 -min((pca_country $ PC1)))) +1
points(pca_country$LONG,pca_country$LAT, col=COLS[pca_country$ind1], pch=16)
Мастерская GWAS

Анализ PCA: цвет по населению и по географии
Мастерская GWAS

Аллели в отдельных локусах зависимы друг от друга
Проблема? Да и нет
Слишком много LD - проблема ? шум >> сигнал
Некоторые (предсказуемые) LD могут быть полезны ? позволяет использовать «маркерные» SNP
Нарушение равновесия по сцеплению (LD)

Стандартный GWAS
Одномерные методы
Использования взаимодействий
Многовариантные методы
Методы штрафной регрессии (LASSO)
Факториальные методы (ФС на основе DAPC)
Методы тестирования ассоциации

Вариации
Тестирование
Точный критерий Фишера, критерий тенденции Кохрана-Армитиджа, критерий хи-квадрат, дисперсионный анализ
Золотой стандарт - точный тест Фишера
Исправление
Бонферрони
Золотой стандарт - FDR
Одномерные методы
Статистика индивидуальных тестов
Поправка на множественное тестирование

Сильные стороны
Недостатки
Простота
Вычислительно быстро
Консервативный
Легко интерпретировать
Многомерная система, одномерная структура
Размер эффекта отдельных SNP может быть слишком маленьким.
Предельные эффекты отдельных SNP ≠ комбинированные эффекты
Одномерный GWAS- сильные и слабые стороны

А как насчет взаимодействия?
Тестирование на ассоциацию

Многовариантные методы
Тестирование на ассоциацию
Штрафная регрессия LASSO
Регрессия хребта (ридж Регрессия)
Нейронные сети
Штрафная регрессия
Байесовские подходы
Факторные методы
Сопоставление байесовской ассоциации эпистаза
Логические деревья
Модифицированное программирование экспрессии генов логической регрессии
Генетическое программирование для ассоциативных исследований
Выбор логической функции
Монте-Карло
Логическая регрессия
Логическая регрессия
Контролируемый PCA
Sparse-PCA
ФС на основе DAPC (snpzip)
Байесовское разбиение
Эластичная сетка
Байесовская логистическая регрессия с выбором переменной случайного поиска
Отношение шансов-
на основе MDR
Метод многофакторного уменьшения размерности
Нейронные сети, оптимизированные для генетического программирования
Параметрический
убывающий метод
Ограниченный
метод разделения
Комбинаторный метод разбиения
Случайные леса
Установить связь
подход
Непараметрические методы

GWAS
Ридж регрессия (Метод регуляризации Тихонова) с использованием пакета rrBLUP
Ридж регрессия BLUP (Meuwissen et al. 2001)
https://rdrr.io/cran/rrBLUP/man/GWAS.html
Модель маркерных эффектов
Y - вектор данных
µ - общее среднее
1 - вектор единиц
Xi - матрица расчета
g - генетический эффект i-го маркера
е - ошибка
у = µ1 + ∑Xi gi + e
Ридж регрессия - это способ создания модели, когда количество переменных-предикторов в наборе превышает количество наблюдений (ситуация с недостаточным объемом данных).

Код R: установите и загрузите необходимые библиотеки
install.packages (c («poolr», «qqman», «BGLR», «rrBLUP», «DT», «SNPRelate», «dplyr», «RgoogleMaps», «ggmap», «mapproj», «sf», "OpenStreetMap", "инструменты разработчика"))
install.packages (c ("rnaturalearth", 'rnaturalearthdata', 'rgeos', 'ggspatial'))
library(rrBLUP)
library(BGLR)
library(DT)
library(SNPRelate)
library(rgeos)
library(ggspatial)
library(qqman)
library(пуллер)
library(OpenStreetMap)
library(rjson)
library(rgdal)
library(RgoogleMaps)
library(mapproj)
library(sf)
library(dplyr)
Библиотека (инструменты разработчика)
devtools :: install_github ("dkahle / ggmap", ref = "tidyup")

GWAS
# создаем файл geno для анализа GWAS пакета rrBLUP
geno_final <- data.frame (marker = MAP1 [, 2], chrom = MAP1 [, 1], pos = MAP1 [, 4], t (Geno1 - 1), check.names = FALSE)
dim(Geno1)
# создаем фенофайл
pheno_final <- data.frame (NSFTV_ID = rownames (y), y = y)
# Запустить анализ GWAS
myGWAS <- GWAS (pheno_final, geno_final, min.MAF = 0,05, P3D = TRUE, plot = TRUE)
Мастерская GWAS

Мастерская GWAS

Использовать поиск по SNP
https://snp-seek.irri.org/
Мастерская GWAS

Есть много новых инструментов
Однако, если вы освоите старый добрый rrBLUP, это поможет понять другие.
Теперь упражнения
Посмотрите на доступные фенотипы, выберите один и повторите упражнения.
"Flowering.time.at.Aberdeen" "Flowering.time.at.Faridpur" "Flowering.time.at.Aberdeen" "FT.ratio.of.Arkansas.Aberdeen" "FT.ratio.of.Faridpur.Aberdeen" " Culm.habit "" Лист. Опушение Flag.leaf.length "" Flag.leaf.width "" Awn.presence "
«Число метелок на растение» «Высота растения» «Длина метелки» «Первичное число ветвей» «Число семян на одну ветвь» «Цветки на одну ветвь» «Метелка. Плодородие» Seed.length "" Seed.width "" Seed.volume "" Seed.surface.area "" Brown.rice.seed.length "" Brown.rice.seed.width "" Brown.rice.surface.area "" Brown .rice.volume "" Seed.length.width.ratio "" Brown.rice.length.width.ratio "" Seed.color "" Pericarp.color "" Straighthead.suseptability "" Blast.resistance "
"Amylose.content" "Alkali.spreading.value" "Protein.content"
"Year07Flowering.time.at.Arkansas" "Year06Flowering.time.at.Arkansas"
Мастерская GWAS

Подход к прогнозированию географической структуры населения (GPS)
Сделать вывод о происхождении человека из полногеномной коллекции маркеров, информативных о происхождении.
Мастерская GWAS

Чтобы сделать вывод о структуре популяции на основе данных генотипа, необходимо сначала уменьшить размерность набора данных из-за тысяч SNP, которые он включает.
От SNP к добавке
Тысячи SNP
ДОБАВКА
Пропорции примесей в географически смежных популяциях, таких как итальянцы и греки, а также в популяциях с похожей историей, таких как британцы и немцы, одинаковы.
Мастерская GWAS

Мастерская GWAS

Прогноз био-происхождения
Зная связь между географическими и генетическими расстояниями, можно ли определить географическое происхождение человека с известным генотипом?
Мы решили попробовать простой подход
Мастерская GWAS

Неизвестные образцы
Мастерская GWAS

GPS точно назначен
~ 100% всех людей в свои континентальные регионы
80% всех людей в страну происхождения
60% всех людей в своем внутреннем регионе страны
Мастерская GWAS

Моделирующая добавка

1001 Genomes - Каталог генетической изменчивости Arabidopsis thaliana.
Образцы из 22 стран
Массив SNP генотипа 250К
Координаты точки сбора
Мы взяли 450 инбредных штаммов A. thaliana из 22 стран, секвенированные и проанализированные Институтом биологии развития Max-Plank и компанией Monsanto из проекта 1001 Genomes (http://1001genomes.org/projects/MPICWang2013/)
Антон Елисеев / Кристина Кривоносова
Мастерская GWAS

Структура населения
Мэтью В. Хортон и др. Полногеномные паттерны генетической изменчивости во всем мире образцов A. thaliana из панели RegMap. Nat. Genet. 44 (2)
Мастерская GWAS

Фильтрация SNP
Мы взяли файлы вариантов, доступные на http://1001genomes.org/data/MPI/MPICWang2013/releases/current/.
отфильтрованы инделы и неаутосомные варианты. Обнаружено около 8 млн вариантов.
Мы выбрали 80 000 SNP из исходного набора данных и создали файл PLINK.
Используя PLINK, мы удалили SNP с частотой минорного аллеля <0,05.
и далее удалил LD с r2> 0,2, используя окно размером 50
Количество оставшихся SNP - 40 518.
Мастерская GWAS

ДОБАВКА
Мы выполнили несколько анализов ADMIXTURE с числом предковых популяций K от 3 до 20, а затем для последующих анализов было выбрано K = 14, так как это давало наименьшее медианное расстояние при проверке исключения по одному.
Затем мы разделили некоторые из исходных популяций на 41 генетически однородную субпопуляцию с помощью K-средних.
Мастерская GWAS

Примесь
Мастерская GWAS

Проверка GPS по одному разу
Процент популяций, которые точно нанесены на карту 60%
Среднее расстояние до правильного населения 200 км
Мастерская GWAS

Исходная гипотеза
Географическое положение Arabidopsis должно быть связано с холодоустойчивостью. Более холодный климат в Северном полушарии наблюдается в более высоких широтах. Растения могут адаптироваться к более холодному климату, регулируя время цветения. Следовательно, мы ожидаем увидеть связь с генами, контролирующими время и широту цветения.
Мастерская GWAS

GPS-анализ O. sativa
Мастерская GWAS

Точность для риса
Среднее расстояние: 4043 км
Мастерская GWAS

Прогнозирование GPS после SNP и географической фильтрации
Среднее расстояние: 1141 км
Мастерская GWAS

Почему не работает с рисом?
Рис не может выбрать свою вторую половинку, а арабидопсис может.
Мастерская GWAS

Наборы данных: WorldClime и SoilDB
Компоненты примеси больше зависят от климата, чем от параметров почвы.
Климат и почва не являются независимыми (p-значение теста Мантеля 0,03)
Затем спрогнозируйте класс компонента примеси по климату и почве.
Точность 0,85
Мастерская GWAS

Часть 1 заключение
Мы умеем моделировать дикие виды - не так уж и много с одомашненными
Климат и почва влияют на генетику (да!)
Мастерская GWAS

Семейство бобовых (бобовых) растений
Возможность установления клубенькового симбиоза
Высокая синтения к бобовым культурам
Небольшой диплоидный геном (500 МБ) секвенирован Tang et al., 2014 (BMC Genomics)
Многочисленные геномные ресурсы
RIL популяции
Коллекции мутантов
Большие базы данных EST и RNASEQ
Ресурсы HAPMAP Branca et al., 2011 (PNAS)
288 геномов уже секвенированы
Инструменты биотехнологии для валидации
Большое биоразнообразие Gentzbittel et al., 2015 (Front. Pl. Sci)
Medicago truncatula
M. truncatula - модельное бобовое растение.
Мастерская GWAS

Проект NSF HapMap
Текущее состояние: 288 последовательных образцов.
Итого: 16.516.721 SNP
30 строк при> 20X
Остальные -> 5X
Бранка и др., 2011 г. (PNAS)
Мастерская GWAS

12.11.2020
Зерновые бобовые
Кормовые бобовые
Люцерна
Клевер
Соя
Нута
Фасоль
Высокая сельскохозяйственная ценность
Горох
Арахис Чечевица
Зерновые Бобовые очень богаты белком (горох: в 3 раза больше, чем злаки)
Фото: Б.Жулье.
Бобовые - это основные сельскохозяйственные виды, используемые в питании человека и животных.
Мастерская GWAS

Данные HapMap для M. truncatula
Текущее состояние: 288 последовательных образцов.
# SNP
Chr1: 1,508,346
Chr2: 1.964.419
Chr3: 2.472.365
Chr4: 2.225.537
Chr5: 3.251.290
Chr6: 1.172.980
Chr7: 2.262.094
Chr8: 1.659.690
Итого: 16.516.721
~ 40% SNP в кодирующих регионах
18 сестринских таксонов
Вот данные, доступные для 288 образцов Medicago, которые уже были секвенированы, что привело к общему количеству более 16,5 миллионов SNP. Среди них 18 образцов, которые сильно отличались от других, были переклассифицированы в некоторые сестринские таксоны, несмотря на полное ботаническое сходство.
Мастерская GWAS

M. truncatula спонтанно встречается по всему Средиземноморскому бассейну.
Доступно несколько коллекций: DZ, TN, FR, AU, US,…. на тысячи образцов
Фото: А. Абдельгерфи, INA, Алжир, Алжир
Фото: J-M Prospéri, INRA Montpellier, Франция.
Фото: J-M Prospéri, INRA Montpellier, Франция.
Фото: Э. Ауани, CBBC, Тунис
Фото: Л. Генцбиттель, ENSAT, Тулуза, Франция.
Фото: Л. Генцбиттель, ENSAT, Тулуза, Франция.
Мастерская GWAS

Некоторые наследственные геномы могут быть адаптированы к местным биоклиматическим переменным.
Мы искали предполагаемую связь с 19 биоклиматическими переменными, как определено WorldClim (http://www.worldclim.org).
Геном «Южного побережья Туниса»
положительно коррелирует со средней годовой T ° C (BIO1) и средней T ° C более холодного квартала (BIO11)
и отрицательно коррелировал с годовым количеством осадков (BIO12).
Геном «испанского побережья»
отрицательно коррелирует с сезонностью T ° C (BIO4),
и годовой диапазон T ° C (BIO7).
«Греческий» геном
положительно коррелирует с осадками самого влажного месяца (BIO13),
сезонность осадков (BIO15),
осадки самого влажного квартала (BIO16)
и осадки самого холодного квартала (BIO19).
Испанский прибрежный
Греческий
BIO1: r = 0,26 ***
BIO11: r = -0,24 **
BIO12: r = - 0,28 ***
BIO4: r = -0,42 ***
BIO7: r = -0,31 ***
BIO13: r = 0,39 ***
BIO15: r = 0,45 ***
BIO16: r = 0,33 ***
BIO19: r = 0,29 ***
Геном «Северного побережья Туниса» не имеет существенной корреляции с какими-либо биоклиматическими переменными.
Мастерская GWAS

PCA определяет две субпопуляции Medicago среди 262 образцов.
40 000 случайно выбранных SNP
В соответствии с результатами Bonhomme et al. 2014 г.
Мастерская GWAS

Уточненная популяционная структура видов Medicago с использованием алгоритмов Admixture и GPS
Medicago, вероятно, будет иметь 8 наследственных геномов.
Gentzbittel et al., Genome biol. , 2019
Мастерская GWAS

Уточненная популяционная структура видов Medicago с использованием алгоритмов Admixture и GPS
Gentzbittel et al., Genome biol. , на пересмотре
Мастерская GWAS

Предсказано местонахождение 17 неизвестных образцов, включая эталонный геном Jemalong-A17.
Точность процедуры «Leave-one Out» <300 километров.
Линейная связь между географическими и генетическими расстояниями
при d <950 км
(R = 0,8, значение p = 10-4, критерий Мантеля)
Gentzbittel et al., Genome biol. , на пересмотре
Мастерская GWAS

Уточненная популяционная структура видов Medicago с использованием алгоритмов Admixture и GPS
Мастерская GWAS

Geographic distribution of the putative 8 ancestral genomes of M. truncatula
Spanish coastal
Algiers
North Tunisian coastal
Atlas
French
Greek
Spanish-Moroccan Inland
South Tunisian coastal
5 out of 8 populations are present in Mahgreb
Algiers (K1)
Spanish coastal (K2)
North Tunisian coastal (K3)
Atlas (K4)
South Tunisian coastal (K5)
French (K6)
Greek (K7)
Spanish-Moroccan Inland (K8)
Maghreb as a putative center of diversity for
M. truncatula species
French
Atlas
Gentzbittel et al., Genome biol. , 2019
GWAS workshop

Last glaciations & main glacial refugia may explain M. truncatula population structure
-20.000 BP
GWAS workshop

Admixture components are correlated with current conditions
Gentzbittel et al., Genome biol. , 2019
GWAS workshop

Verticillium sp. are soil-borne fungal pathogens for Legumes
GWAS workshop

Ancestral genomes may present different levels of resistance to Verticillium wilt
V. alfalfae symptoms on susceptible M.truncatula (Ben et al. 2013)
GWAS workshop

Ancestral genomes present different levels of resistance to Verticillium wilt
Linear Model (242 Mt HAPMAP accessions)
MSSCorLS ~ K2 + K5 + K7 + K8
r²=0.31, P-value < 2.2 x 10-16
Gentzbittel et al., Genome biol. , 2019
GWAS workshop

Can phenotype be predicted using WhoGem?
Use plant model (Medicago truncatula) to test applicability of admixture framework to predict phenotypes.
GPS was used to predict provenance of 59 unknown M. truncatula accessions, and MSS was inferred based on the similarity of admixture profiles to the training set accessions using WhoGem model.
M. truncatula was infected with V. alfalfae and wilting symptoms recorded.
MSS (Maximum symptom score) ranges from 0 (healthy) to 4 (dead), 2 threshold.
V. alfalfae symptoms on susceptible M.truncatula (Ben et al. 2013)
We thank Jean-Marie Prospéri (INRA Montpellier, France) for providing seeds of the accessions
GWAS workshop

Experimental validation of predicted resistance levels
21/29 samples that are predicted to be resistant were resistant
24/30 samples predicted to be susceptible were susceptible
Chi-square P-value = 1.5.10-4
Gentzbittel et al., Genome biol. , under revision
GWAS workshop

WhoGem models can be significant predictors of quantitative functional traits in plants
Comparisons of accuracy of prediction for several quantitative traits of adaptive relevance in M. truncatula using 5 Genomic Selection algorithms and our WhoGEM method.
Accuracy is computed using 50 runs of five-fold cross-validation for all methods.
Gentzbittel et al., Genome biol. , 2019
GWAS workshop

What is missing?
Distance between two points should not be just geometric distance. Add:
Mode of propagation (e.g. Medicago seeds stick to goats)
Soil
Climate
Geography
GWAS workshop

GWAS workshop

GWAS workshop

GWAS workshop

This methodology may serve as a basis for analyses in other plant species and for other functional quantitative
traits of interest.
Part 2 conclusions
Large proportion of phenotypic variation between individuals may be best explained by population admixture.
Variation in genome admixture proportion explains most of phenotypic variation for quantitative functional traits.
We experimentally confirm the prediction of differences in quantitative disease resistance levels in the wild model legume Medicago truncatula.
Admixture components were found to be significantly related to climate and geography, also positive selection at the species level might not explain current adaptation.
Phenotypes can be predicted using genome-wide patterns of admixture, when incorporating covariates such as individuals' provenance.
This insight contributes to the understanding of adaptation, and can accelerate plant and animal breeding, and biomedical research programs.
GWAS workshop

Practice
Admixture.csv
Pheno.csv
GWAS workshop

ADM=read.csv("Admixture.csv",row.names = 1)
# Instruction to name bars:
barNaming <- function(vec) {
retVec <- vec
for (k in 2:length(vec)) {
if (vec[k - 1] == vec[k])
retVec[k] <- ""
}
return(retVec)
}
# Size of the plot - these are margins, also related to the end appearance of the plot
par(mar = c(1, 1, 1, 1))
x=barplot(t(as.matrix(ADM[, 7:14])),
col = rainbow(8),
horiz = FALSE,
las = 2,
cex.axis = 0.1,
cex.names = 0.65,
)
text(x,0.5,barNaming(ADM$Country),srt=45, col='white')
GWAS workshop

Clustered vectors
GWAS workshop

Merge datasets
ADM=read.csv('Admixture.csv',row.names = 1); TEST=subset(ADM, is.na(ADM$Lat))
REF=subset(ADM, !is.na(ADM$Lat)); Phen=read.csv('pheno.csv',row.names = 1); head(Phen) #merge reference adm with pheno
AMDPH=na.omit(merge(REF, Phen, by='row.names'))
GWAS workshop

DATA=AMDPH[,c(4:6,8:15,17)];
cor(DATA)
GWAS workshop

##GPS
M=25;
for(i in 1:dim(TEST)[1]){
Y=TEST[i,7:14];Y
DIST=c()
for(j in 1:dim(AMDPH)[1]){
Z=AMDPH[j,8:15];Z
d=sum((Y-Z)^2);d DIST=c(DIST,d)
}
I=order(DIST)[1:M];I
Ph=AMDPH$MaxSymptomScore[I];Ph
Dist=DIST[I]+1.E-15;Dist
w=min(Dist)/Dist;w
predPh=sum(w*Ph)/sum(w);predPh
predLat=sum(w*AMDPH$Lat[I])/sum(w);predLat; predLon=sum(w*AMDPH$Long[I])/sum(w);predLon
output=paste(row.names(TEST[i,]),TEST$Pop[i],TEST$Country[i],predPh,predLat,predLon, sep=',')
write.table(output,"GPS_RES.csv",append = T, row.names = F, col.names = F, quote=T)
}
GWAS workshop

Results
GWAS workshop

#GLM
library(lmtest); library(dplyr); library(car); library(rcompanion)
#set up formulas
formula=formula(MaxSymptomScore~K1+K2+K3+K4+K5+K6+K7+K8)
formula0=MaxSymptomScore ~ 1
#initialize models
model.null = glm(formula0, data=AMDPH); summary(model.null)
model.full = glm(formula, data=AMDPH); summary(model.full )
#perform stepwise reduction
S=step(model.null, scope = list(upper=model.full), direction="both", test="Chisq", data=AMDPH)
summary(S)
#final model
model.final = glm(S$formula, data=AMDPH);model.final
summary(model.final)
AMDPH$predy = predict.glm(model.final, newdata = AMDPH, type="response" )

Laurent Gentzbittel, Ecolab, Toulouse, France/Skoltech, Moscow, Russia
Nevin Young's lab, Univ. of Minnesota, USA
Sergey Nuzdhin's lab, USC, USA
Chris Town's lab, JCVI, USA
Mounawer Badri & Naceur Djebali, CBBC, Tunisia
Jean-Marie Prospéri, INRA Mauguio, France
Nick Alexandrov & Dmitro Chebotarov, IRRI, Philippines
Collaborations and funding
Laurent Gentzbittel
GWAS workshop