################################################### ### chunk number 1: ################################################### options(width=90) ################################################### ### chunk number 2: eval=FALSE ################################################### ## install.packages( c("xtable","combinat","gdata","gplots","mvtnorm"), dep=TRUE) ################################################### ### chunk number 3: eval=FALSE ################################################### ## repos <- c("http://www.warnes.net/bioc2007/", "http://cran.fhcrc.org") ## install.packages(c("GeneticsBase", ## #"GeneticsQC", ## "GeneticsDesign", "fbat"), ## repos=repos, type="source", dep=TRUE) ################################################### ### chunk number 4: ################################################### library(GeneticsBase) library(GeneticsDesign) library(fbat) ################################################### ### chunk number 5: ################################################### hm.a<-readGenes(gfile="hmCEU_YRI_chr22_ALLfbat.ped", gformat="fbat") print(hm.a) ################################################### ### chunk number 6: ################################################### hm.f<-readGenes(gfile="hmCEU_YRI_chr22_Foundersfbat.ped", gformat="fbat") ################################################### ### chunk number 7: ################################################### hm.a2<-hm.a[1:100,] ################################################### ### chunk number 8: ################################################### mG<-missGFreq(hm.a2, founderOnly=TRUE, quiet=FALSE) head(mG$nMissSubjects) ################################################### ### chunk number 9: ################################################### head(mG$nMissMarkers) ################################################### ### chunk number 10: ################################################### par(mfrow=c(2,1)) hist(mG$nMissSubjects[,1], main="", xlab="number of missing genotypes") plot(1:nrow(mG$nMissSubjects),mG$nMissSubjects[,1], xlab="subject ID", ylab="number of missing genotypes", type="h") title("Counts of missing genotypes") par(mfrow=c(1,1)) ################################################### ### chunk number 11: ################################################### ################ cM<-checkMarkers(hm.a2) head(cM) ################################################### ### chunk number 12: ################################################### cMend<-checkMendelian(hm.a2, quiet = FALSE) ################################################### ### chunk number 13: ################################################### head(cMend$nMerrMarker) ################################################### ### chunk number 14: ################################################### head(cMend$nMerrFamily) ################################################### ### chunk number 15: ################################################### par(mfrow=c(2,1)) hist(cMend$nMerrMarker, main="", xlab="number of Mendelian errors") plot(1:length(cMend$nMerrFamily),cMend$nMerrFamily, xlab="family ID", ylab="number of Mendian Errors") par(mfrow=c(1,1)) ################################################### ### chunk number 16: ################################################### t1<-alleleSummary(hm.a2[1:10]) t1 ################################################### ### chunk number 17: ################################################### t2<-genotypeSummary(hm.a2[1:10], founderOnly=TRUE) t2 ################################################### ### chunk number 18: ################################################### hwe<-HWE(hm.a2[1:10]) hwe ################################################### ### chunk number 19: ################################################### ld.small<-LD(hm.a2[1:20]) plot(ld.small) ################################################### ### chunk number 20: ################################################### ld <- LD(hm.a2) LDView(ld@"X^2") ################################################### ### chunk number 21: ################################################### res.A<-Armitage(hm.a2, method="A") head(res.A) ################################################### ### chunk number 22: ################################################### res.R<-Armitage(hm.a2, method="R") head(res.R) ################################################### ### chunk number 23: ################################################### res.D<-Armitage(hm.a2, method="D") head(res.D) ################################################### ### chunk number 24: ################################################### sampleInfo(hm.f)$race <- sampleInfo(hm.f)$affected raceval <- sampleInfo(hm.f)$race-1 sampleInfo(hm.f)$Norm0.0 <- rnorm(nobs(hm.f), mean=0.0*raceval) sampleInfo(hm.f)$Norm0.5 <- rnorm(nobs(hm.f), mean=0.5*raceval) sampleInfo(hm.f)$Norm1.0 <- rnorm(nobs(hm.f), mean=1.0*raceval) sampleInfo(hm.f)$Norm1.5 <- rnorm(nobs(hm.f), mean=1.5*sampleInfo(hm.f)$race) doSample <- function(raceval, mult) { prob <- c( 0.33-raceval*mult, 0.33+(raceval*mult)/2, 0.33+(raceval*mult)/2) factor(sample( x=c("Red","Green","Blue"), size=1, p=prob, rep=T)) } sampleInfo(hm.f)$Cat0.0 <- sapply( raceval, doSample, mult=0.0 ) sampleInfo(hm.f)$Cat0.1 <- sapply( raceval, doSample, mult=0.1 ) sampleInfo(hm.f)$Cat0.2 <- sapply( raceval, doSample, mult=0.2 ) sampleInfo(hm.f)$Cat0.3 <- sapply( raceval, doSample, mult=0.3 ) ################################################### ### chunk number 25: ################################################### model <- function( markerName ) { # extract requested genetic marker genotype <- genotypes(hm.f,marker=markerName) # get data frame to use for fitting the model mframe <- model.frame(race ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + genotype, data=sampleInfo(hm.f) ) # To test significance of a term, best method is to do anova of # the full model against a submodel omitting the particular term. # This avoids issues with changes in names of factor levels, # presence or absence of covaraiates, etc. result <- try( { fit.with <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5 + as.factor(genotype), data=mframe, family="binomial") fit.without <- glm( race==1 ~ sex + Norm0.0 + Norm0.5 + Norm1.0 + Norm1.5, data=mframe, family="binomial") anova(fit.with, fit.without, test="Chisq")$"P(>|Chi|)"[2] } ) if(class(result)=="try-error") return(NA) # or return(result) to see the error messages else result # full result. Usually we want to specify exactly which # parameters and stats get returned so the format is consistent # across all markers. } ################################################### ### chunk number 26: ################################################### t1 <- unix.time( fits <- sapply( markerNames(hm.f)[1:50], model ) )[3] t1 ################################################### ### chunk number 27: ################################################### t2 <- unix.time( fits <- sapply( markerNames(hm.f)[1:100], model ) )[3] t2 ################################################### ### chunk number 28: ################################################### t1 + (t2 - t1)/50 * (nmarker(hm.f)-50) / 60 ################################################### ### chunk number 29: ################################################### fits.sorted <- sort(fits) plot(-log(fits), type="h", col="blue") labels <- c(0.05, 0.01, 0.001, 0.0001, 0.00001, 0.00001) abline(h=-log(labels), lty=2, col="red") mtext(text=as.character(labels), side=4, at=-log(labels), col="red", las=1) title("Per-marker statistical significance") ################################################### ### chunk number 30: ################################################### f<-fbat(hm.a2) ################################################### ### chunk number 31: ################################################### summaryPvalue(f) ################################################### ### chunk number 32: ################################################### viewstat(f, "22_14884399_rs4911642") ################################################### ### chunk number 33: ################################################### gregorius(freq=0.15, N=20) ################################################### ### chunk number 34: ################################################### gregorius(freq=0.15, missprob=1-0.95) ################################################### ### chunk number 35: ################################################### GeneticPower.Quantitative.Numeric( N=50, freq=0.1, minh="recessive", alpha=0.05 ) GeneticPower.Quantitative.Factor( N=50, freq=0.1, minh="recessive", alpha=0.05 ) ################################################### ### chunk number 36: ################################################### power.range <- function( N, ... ) { sapply( N, function(n) GeneticPower.Quantitative.Numeric( N=n, ... ) ) } power.range( N=c(25,50,100,200,500), freq=0.1, minh="recessive", alpha=0.05 ) ################################################### ### chunk number 37: ################################################### fun <- function(p) power.range(freq=p, N=seq(100,1000,by=100), alpha=5e-2, minh='recessive') m <- sapply( X=seq(0.1,0.9, by=0.1), fun ) colnames(m) <- seq(0.1,0.9, by=0.1) rownames(m) <- seq(100,1000,by=100) print(round(m,2)) ################################################### ### chunk number 38: ################################################### res1<-GPC(pA=0.05, pD=0.1, RRAa=1.414, RRAA=2, r2=0.9, pB=0.06, nCase=500, ratio=1, alpha=0.05, quiet=FALSE) ################################################### ### chunk number 39: ################################################### res2<-GPC.default(pA=0.05, pD=0.1, RRAa=1.414, RRAA=2, Dprime=0.9, pB=0.06, nCase=500, ratio=1, alpha=0.05, quiet=FALSE)