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)))
## 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)))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.07058881Start 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,] -85.93596
##
## $stats
## [,1]
## [1,] 1.50768
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.3529408
##
## $ves
## [,1]
## [1,] 1.602349e-05
##
## $beta
## [,1]
## [1,] 5.298648e-15
## [2,] 2.132361e-01tr1.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,] 9.210247e-15
## [2,] -2.516358e-01tes.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,] 7.257034e-15
## [2,] 2.547853e-01knob180.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.350835e-15
## [2,] -4.254099e-02Include 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.1350998
##
## $REMLs
## [,1]
## [1,] -85.30971
##
## $stats
## [,1]
## [1,] 1.507682
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.3529407
##
## $ves
## [,1]
## [1,] 1.602348e-05
##
## $beta
## [,1]
## [1,] -3.290052e-05
## [2,] 2.132364e-01F.tr1.alt <- emma.REML.t_beta ( tr1 , alt , X0 = matrix (F, 93 , 1 ), K = A ) #no
F.tr1.alt## $ps
## [,1]
## [1,] 0.2034592
##
## $REMLs
## [,1]
## [1,] -125.4279
##
## $stats
## [,1]
## [1,] -1.28097
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.4764013
##
## $ves
## [,1]
## [1,] 0.3940113
##
## $beta
## [,1]
## [1,] -1.0069027
## [2,] -0.2327138F.tes.alt <- emma.REML.t_beta ( tes , alt , X0 = matrix (F, 93 , 1 ), K = A ) #no
F.tes.alt## $ps
## [,1]
## [1,] 0.1525412
##
## $REMLs
## [,1]
## [1,] -105.8068
##
## $stats
## [,1]
## [1,] 1.44269
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.5266867
##
## $ves
## [,1]
## [1,] 0.02785989
##
## $beta
## [,1]
## [1,] 0.0348864
## [2,] 0.2512215F.knob180.alt <- emma.REML.t_beta ( knob180 , alt , X0 = matrix (F, 93 , 1 ),K = A ) #no
F.knob180.alt## $ps
## [,1]
## [1,] 0.8504986
##
## $REMLs
## [,1]
## [1,] -127.5842
##
## $stats
## [,1]
## [1,] -0.1890193
##
## $dfs
## [,1]
## [1,] 91
##
## $vgs
## [,1]
## [1,] 0.8937359
##
## $ves
## [,1]
## [1,] 4.057555e-05
##
## $beta
## [,1]
## [1,] 4.283125e-05
## [2,] -4.254135e-02Alternate hypothesis. Drop lowest three elevation that show highgest inbreeding and introgression.
hidt<-subset(dt,pheno$Altitude>1900)
hiA<-A.mat(hidt)
#filter for hi altitude
hipheno<-filter(pheno,pheno$Altitude>1900)
## GENOME bp, standardized
higs <- t(as.matrix (scale(hipheno$GS_bp/1E9)))
## knob bp, standardized
hiknob180 <- t ( as.matrix(scale(hipheno$X180knobbp/1E9)))
# tr1 bp, standardized
hitr1 <- t ( as.matrix(scale(hipheno$TR1bp/1E9)))
#TEs
hites <- t ( as.matrix(scale(hipheno$TotalTebp/1E9)))
#altitude
hialt <- t ( as.matrix(scale(hipheno$Altitude)))
gs.hialt <- emma.REML.t_beta ( higs , hialt , K = hiA ) #no
gs.hialt## $ps
## [,1]
## [1,] 1
##
## $REMLs
## [,1]
## [1,] NA
##
## $stats
## [,1]
## [1,] NA
##
## $dfs
## [,1]
## [1,] NA
##
## $vgs
## [,1]
## [1,] NA
##
## $ves
## [,1]
## [1,] NA
##
## $beta
## function (a, b)
## .Internal(beta(a, b))
## <bytecode: 0x7faf388a5750>
## <environment: namespace:base>tr1.hialt <- emma.REML.t_beta ( hitr1 , hialt , K = hiA ) #no
tr1.hialt## $ps
## [,1]
## [1,] 1
##
## $REMLs
## [,1]
## [1,] NA
##
## $stats
## [,1]
## [1,] NA
##
## $dfs
## [,1]
## [1,] NA
##
## $vgs
## [,1]
## [1,] NA
##
## $ves
## [,1]
## [1,] NA
##
## $beta
## function (a, b)
## .Internal(beta(a, b))
## <bytecode: 0x7faf388a5750>
## <environment: namespace:base>tes.hialt <- emma.REML.t_beta ( hites , hialt , K = hiA ) #no
tes.hialt## $ps
## [,1]
## [1,] 1
##
## $REMLs
## [,1]
## [1,] NA
##
## $stats
## [,1]
## [1,] NA
##
## $dfs
## [,1]
## [1,] NA
##
## $vgs
## [,1]
## [1,] NA
##
## $ves
## [,1]
## [1,] NA
##
## $beta
## function (a, b)
## .Internal(beta(a, b))
## <bytecode: 0x7faf388a5750>
## <environment: namespace:base>knob180.alt <- emma.REML.t_beta ( hiknob180 , hialt , K = hiA ) #no
knob180.alt## $ps
## [,1]
## [1,] 1
##
## $REMLs
## [,1]
## [1,] NA
##
## $stats
## [,1]
## [1,] NA
##
## $dfs
## [,1]
## [1,] NA
##
## $vgs
## [,1]
## [1,] NA
##
## $ves
## [,1]
## [1,] NA
##
## $beta
## function (a, b)
## .Internal(beta(a, b))
## <bytecode: 0x7faf388a5750>
## <environment: namespace:base>