## ---- echo=FALSE, message=FALSE------------------------------------------ library(GENESIS) library(GWASTools) # file path to GDS file gdsfile <- system.file("extdata", "HapMap_ASW_MXL_geno.gds", package="GENESIS") # read in GDS data HapMap_geno <- GdsGenotypeReader(filename = gdsfile) # create a GenotypeData class object HapMap_genoData <- GenotypeData(HapMap_geno) # load saved matrix of KING-robust estimates data("HapMap_ASW_MXL_KINGmat") # run PC-AiR mypcair <- pcair(genoData = HapMap_genoData, kinMat = HapMap_ASW_MXL_KINGmat, divMat = HapMap_ASW_MXL_KINGmat) # run PC-Relate mypcrel <- pcrelate(genoData = HapMap_genoData, pcMat = mypcair$vectors[,1], training.set = mypcair$unrels) # generate a phenotype set.seed(4) pheno <- 0.2*mypcair$vectors[,1] + rnorm(mypcair$nsamp, mean = 0, sd = 1) ## ------------------------------------------------------------------------ # mypcair contains PCs from a previous PC-AiR analysis # mypcrel contains Kinship Estimates from a previous PC-Relate analysis # pheno is a vector of Phenotype values # make a data.frame mydat <- data.frame(scanID = mypcrel$sample.id, pc1 = mypcair$vectors[,1], pheno = pheno) head(mydat) # make ScanAnnotationDataFrame scanAnnot <- ScanAnnotationDataFrame(mydat) scanAnnot ## ---- eval=FALSE--------------------------------------------------------- # geno <- MatrixGenotypeReader(genotype = genotype, snpID = snpID, chromosome = chromosome, # position = position, scanID = scanID) # genoData <- GenotypeData(geno) ## ---- eval=FALSE--------------------------------------------------------- # geno <- GdsGenotypeReader(filename = "genotype.gds") # genoData <- GenotypeData(geno) ## ---- eval=FALSE--------------------------------------------------------- # snpgdsBED2GDS(bed.fn = "genotype.bed", bim.fn = "genotype.bim", fam.fn = "genotype.fam", # out.gdsfn = "genotype.gds") ## ---- eval=FALSE--------------------------------------------------------- # # read in GDS data # gdsfile <- system.file("extdata", "HapMap_ASW_MXL_geno.gds", package="GENESIS") # HapMap_geno <- GdsGenotypeReader(filename = gdsfile) ## ------------------------------------------------------------------------ # create a GenotypeData class object with paired ScanAnnotationDataFrame HapMap_genoData <- GenotypeData(HapMap_geno, scanAnnot = scanAnnot) HapMap_genoData ## ------------------------------------------------------------------------ myGRM <- pcrelateMakeGRM(mypcrel) myGRM[1:5,1:5] ## ------------------------------------------------------------------------ # fit the null mixed model nullmod <- fitNullMM(scanData = scanAnnot, outcome = "pheno", covars = "pc1", covMatList = myGRM, family = gaussian) ## ---- eval=FALSE--------------------------------------------------------- # nullmod <- fitNullMM(scanData = scanAnnot, outcome = "pheno", covars = c("pc1","pc2","sex","age"), # covMatList = myGRM, family = gaussian) ## ---- eval=FALSE--------------------------------------------------------- # nullmod <- fitNullMM(scanData = scanAnnot, outcome = "pheno", covars = "pc1" # covMatList = list("GRM" = myGRM, "House" = H), family = gaussian) ## ---- eval=FALSE--------------------------------------------------------- # nullmod <- fitNullMM(scanData = scanAnnot, outcome = "pheno", covars = "pc1", covMatList = myGRM, # family = gaussian, group.var = "race") ## ---- eval=FALSE--------------------------------------------------------- # nullmod <- fitNullMM(scanData = scanAnnot, outcome = "pheno", covars = "pc1", covMatList = myGRM, # family = binomial) ## ------------------------------------------------------------------------ assoc <- assocTestMM(genoData = HapMap_genoData, nullMMobj = nullmod, test = "Wald") ## ---- eval = FALSE------------------------------------------------------- # # mysnps is a vector of snpID values for the SNPs we want to test # assoc <- assocTestMM(genoData = HapMap_genoData, nullMMobj = nullmod, test = "Wald", # snp.include = mysnps) ## ---- eval = FALSE------------------------------------------------------- # assoc <- assocTestMM(genoData = HapMap_genoData, nullMMobj = nullmod, test = "Wald", # chromosome = 22) ## ------------------------------------------------------------------------ head(assoc) ## ------------------------------------------------------------------------ varCompCI(nullMMobj = nullmod, prop = TRUE) ## ---- echo=FALSE--------------------------------------------------------- close(HapMap_genoData)