Hi all,
I want to associate mortality with ~100K SNPs, in 6,500 samples that are
divided up into 60 breeds.
So it's important to account for population stratification in my analysis.
I'm using egscore (the eigenstrat method) for the association (and I've
tried using the polygen and grammar packages too). The problem is that I
cannot get lambda (the inflation factor) to go below the accepted inflation
threshold of 1.1, it seems to converge around 1.5-1.6 when I include an
increasing number of PC axes into the analysis.
Then, I also tried to use the polygen and grammar packages, but the same
thing happens.
Here is the code that I am using:
library(GenABEL)
#Load the data as a gwaa.class
my.geno.data <- load.gwaa.data(phenofile =
"pheno.dat",genofile="youroutput.dat")
#Calculate the IBS matrix
kin<-ibs(my.geno.data, weight="freq", snpfreq=NULL)
diag(kin) <-hom(my.geno.data)$Var
# Estimation of polygenic model, This estimates the residuals when the
effect of covariates breeds are factored out.
polygen <-polygenic(mortality~breed,kin=kin,my.geno.data)
grammarobject <-qtscore(polygen$pgres,data=my.geno.data,clam=FALSE)
....so then you found that the lambda was still > 1.1
#So then I used egscore on the output from polygen
output <-egscore(polygen$pgresidualY,data=my.geno.data,kin=kin,naxes=X)
....where I iteratively included an increasing number of PC axes (naxes=X),
lambda still > 1.1, and doesn't change if I run egscore on the raw data
instead of the environmental residuals (again iterating through axes), and
it also doesn't change regardless of if I include breed as a co-variate, or
as a stratification variable.
output <-egscore(mortality~breed,data=my.geno.data,kin=kin,naxes=X)
Am I doing something wrong? What else can I try to properly account for
population stratification in an association between trait and SNPs?
Thanks
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
