Mexicana_Emma
Main hypothesis: selection for genome size along an altitudinal gradient across the entire range.
library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.3.2## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr## Warning: package 'ggplot2' was built under R version 3.3.2## Warning: package 'tidyr' was built under R version 3.3.2## Conflicts with tidy packages ----------------------------------------------## filter(): dplyr, stats
## lag(): dplyr, statslibrary (emma )
library(ggExtra)## Warning: package 'ggExtra' was built under R version 3.3.2library(rrBLUP)
# return beta function for emma
source("~/Projects/results/bilinski_gsize/SelectionTests/jri_emma_returnbeta.txt")Phenotypes
#phenotype data in matching order
pheno <- read.csv("~/Projects/results/bilinski_gsize/Mexicana/Master_mexnucleo_pruned.csv")
## GENOME bp, standardized
gs <- t(as.matrix (scale(pheno$GS_bp/1E9,scale=FALSE,center=FALSE)))
## knob bp, standardized
knob180 <- t ( as.matrix(scale(pheno$X180knobbp/1E9)))
# tr1 bp, standardized
tr1 <- t ( as.matrix(scale(pheno$TR1bp/1E9)))
#TEs
tes <- t ( as.matrix(scale(pheno$TotalTebp/1E9)))
#altitude
alt <- t ( as.matrix(scale(pheno$Altitude,scale=FALSE,center=TRUE)+0.5))Genotypes
geno <- read.csv("~/Projects/results/bilinski_gsize/SNP_data/Mexicana_conversion/GBS_final.txt")
geno$X <- NULL
#kinship matrix
dt <- t(geno)
A <- A.mat(dt)Calculate inbreeding
#Usese #F=\sum_{i=1}^{loci}{\frac{\frac{homo_true-p_i^2}{p_i}}{nloci}}
altgeno=geno+1
homo_1<-altgeno==0
homo_2<-altgeno==2
p=rowSums(altgeno,na.rm=T)/(2*(ncol(altgeno)-rowSums(is.na(altgeno))))
F=colSums(homo_1-(1-p)^2/(1-p)+(homo_2-p^2)/p,na.rm=T)/nrow(altgeno)
F=as.matrix(F,ncol=1)
#check inbreeding coeffs vs GS
Ftest<-data.frame(as.vector(gs),F)
Fplot<-ggplot(Ftest)+geom_point(aes(y=as.vector.gs.,x=F))+ylab("standardized gs")
ggExtra::ggMarginal(Fplot, type = "histogram")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#check vs altitude
Atest<-data.frame(as.vector(alt),F)
Aplot<-ggplot(Atest)+geom_point(aes(y=as.vector.alt.,x=F))+ylab("standardized altitude")
ggExtra::ggMarginal(Aplot, type = "histogram")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#get mean
Fprime=F; Fprime[Fprime<0]=0
mean(Fprime)## [1] 0.07058881gtest<-data.frame(as.vector(alt),as.vector(gs),F)
bob<-group_by(gtest,as.vector.alt.) %>% summarize(mean(as.vector.gs.),mean(F))
head(bob)## # A tibble: 6 × 3
## as.vector.alt. `mean(as.vector.gs.)` `mean(F)`
## <dbl> <dbl> <dbl>
## 1 -749.51613 5.347704 0.029486637
## 2 -675.51613 5.452893 0.117640732
## 3 -459.51613 6.046213 0.130199964
## 4 -87.51613 6.304432 0.018017327
## 5 12.48387 6.276587 -0.003365026
## 6 67.48387 6.319102 0.013608780colnames(bob)=c("popalt","popgs","popF")
cor(bob$popF,bob$popalt)## [1] -0.04392102Start Testing
Gs and components vs altitude
gs.alt <- emma.REML.t_beta ( gs , alt , K = A ) #no
gs.alt## $ps
## [,1]
## [1,] 0.1351003
##
## $REMLs
## [,1]
## [1,] 16.78149
##
## $stats
## [,1]
## [1,] 1.50768
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.03692075
##
## $ves
## [,1]
## [1,] 1.6762e-06
##
## $beta
## [,1]
## [1,] 6.0412176011
## [2,] 0.0001764108tr1.alt <- emma.REML.t_beta ( tr1 , alt , K = A ) #no
tr1.alt## $ps
## [,1]
## [1,] 0.2140261
##
## $REMLs
## [,1]
## [1,] -126.5884
##
## $stats
## [,1]
## [1,] -1.251317
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.6501334
##
## $ves
## [,1]
## [1,] 0.224181
##
## $beta
## [,1]
## [1,] 0.000321827
## [2,] -0.000643654tes.alt <- emma.REML.t_beta ( tes , alt , K = A ) #no
tes.alt## $ps
## [,1]
## [1,] 0.1407168
##
## $REMLs
## [,1]
## [1,] -106.3851
##
## $stats
## [,1]
## [1,] 1.486071
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.503268
##
## $ves
## [,1]
## [1,] 0.05239913
##
## $beta
## [,1]
## [1,] -0.0003258550
## [2,] 0.0006517101knob180.alt <- emma.REML.t_beta ( knob180 , alt , K = A ) #no
knob180.alt## $ps
## [,1]
## [1,] 0.8504999
##
## $REMLs
## [,1]
## [1,] -128.2105
##
## $stats
## [,1]
## [1,] -0.1890177
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.8937361
##
## $ves
## [,1]
## [1,] 4.057556e-05
##
## $beta
## [,1]
## [1,] 5.440736e-05
## [2,] -1.088147e-04Include per individual inbreeding as a covariate, still nothing associated with altitude.
F.gs.alt <- emma.REML.t_beta ( gs , alt , X0 = matrix (F, 93 , 1 ), K = A )
F.gs.alt #still no## $ps
## [,1]
## [1,] 0.5593477
##
## $REMLs
## [,1]
## [1,] -278.941
##
## $stats
## [,1]
## [1,] 0.5859679
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.001221887
##
## $ves
## [,1]
## [1,] 26.91386
##
## $beta
## [,1]
## [1,] 2.629761e+01
## [2,] 8.114309e-04F.tr1.alt <- emma.REML.t_beta ( tr1 , alt , X0 = matrix (F, 93 , 1 ), K = A ) #no
F.tr1.alt## $ps
## [,1]
## [1,] 0.2046628
##
## $REMLs
## [,1]
## [1,] -125.4324
##
## $stats
## [,1]
## [1,] -1.277536
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.4757504
##
## $ves
## [,1]
## [1,] 0.3948001
##
## $beta
## [,1]
## [1,] -1.0057379815
## [2,] -0.0005934622F.tes.alt <- emma.REML.t_beta ( tes , alt , X0 = matrix (F, 93 , 1 ), K = A ) #no
F.tes.alt## $ps
## [,1]
## [1,] 0.1529214
##
## $REMLs
## [,1]
## [1,] -105.8091
##
## $stats
## [,1]
## [1,] 1.441339
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.5264829
##
## $ves
## [,1]
## [1,] 0.02810037
##
## $beta
## [,1]
## [1,] 0.0306105010
## [2,] 0.0006420067F.knob180.alt <- emma.REML.t_beta ( knob180 , alt , X0 = matrix (F, 93 , 1 ),K = A ) #no
F.knob180.alt## $ps
## [,1]
## [1,] 0.8490635
##
## $REMLs
## [,1]
## [1,] -127.5837
##
## $stats
## [,1]
## [1,] -0.1908561
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.8937291
##
## $ves
## [,1]
## [1,] 4.057524e-05
##
## $beta
## [,1]
## [1,] 0.0008878460
## [2,] -0.0001098901Alternate hypothesis. Drop lowest three elevation that show highgest inbreeding and introgression.
higeno <- read.csv("~/Downloads/GBS_alt_threshold_final.txt")
higeno$X <- NULL
#kinship matrix
hidt <- t(higeno)
hiA <- A.mat(hidt)
#hi altitude only
hipheno <- read.csv("~/Downloads/Pheno_alt_threshold_ordered.csv")
## GENOME bp, standardized
higs <- t(as.matrix (scale(hipheno$GS_bp/1E9)))
## knob bp, standardized
hiknob180 <- t ( as.matrix(scale(hipheno$X180knobbp/1E9,scale=FALSE)))
# tr1 bp, standardized
hitr1 <- t ( as.matrix(scale(hipheno$TR1bp/1E9,scale=FALSE,center=TRUE)))
#TEs
hites <- t ( as.matrix(scale(hipheno$TotalTebp/1E9,scale=FALSE)))
#altitude
hialt <- t ( as.matrix(scale(hipheno$Altitude,scale=FALSE)+0.5))
gs.hialt <- emma.REML.t_beta ( higs , hialt , K = hiA ) #no
gs.hialt$ps## [,1]
## [1,] 0.0006754742#try with pop as covariate
pop=matrix(c(rep(1,9),rep(2,9),rep(3,8),rep(4,9),rep(5,9),rep(6,9),rep(7,9),rep(8,8)),ncol=1)
gs.pophialt <- emma.REML.t_beta ( higs , rbind(hialt,t(pop)), K = hiA ) #no
gs.pophialt$ps## [,1]
## [1,] 0.0006754742
## [2,] 1.0000000000tr1.hialt <- emma.REML.t_beta ( hitr1 , hialt , X0=t(higs), K = hiA ) #no
tr1.hialt$ps## [,1]
## [1,] 0.07002969tes.hialt <- emma.REML.t_beta ( hites , hialt , X0=t(higs), K = hiA ) #no
tes.hialt$ps## [,1]
## [1,] 0.1810733knob180.alt <- emma.REML.t_beta ( hiknob180 , hialt , K = hiA ) #no
knob180.alt$ps## [,1]
## [1,] 0.7337353