R Under development (unstable) (2024-03-03 r86036 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > ## Expanded test script for pedgene package > devel=FALSE > if(devel) { + rfiles <- list.files(path="../R/", pattern="*.R$") + for(f in rfiles) source(paste("../R/",f,sep="")) + dfiles <- list.files(path="../data/", pattern="*.RData$") + for(d in dfiles) load(paste("../data/",d,sep="")) + library(kinship2) + library(survey) + } else { + + require(pedgene) + data(example.ped) + data(example.geno) + data(example.map) + data(example.relation) + } Loading required package: pedgene Loading required package: Matrix Loading required package: CompQuadForm Loading required package: survey Loading required package: grid Loading required package: survival Attaching package: 'survey' The following object is masked from 'package:graphics': dotchart Loading required package: kinship2 Loading required package: quadprog > > #require(survery) > > > ###################################### > ## From Dan Weeks, issues to check > ###################################### > ## 1) missid ="0" when the rest of the ids are character > ## 2) skip pedigree checking, checkpeds=TRUE/FALSE > ## 3) character alleles > ## 4) disconnected pedigrees > ## 5) "flipped" 0/2 geno counts > > ## base case, m-b weights > pg.out.m2 <- pedgene(ped=example.ped, geno=example.geno, map=example.map, male.dose=2, + weights.mb=TRUE,checkpeds=TRUE) > > # summary/print and plot methods for this object > print.pedgene(pg.out.m2,digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 80.1 0.389 -2.21 2 AX X 7 3 211.5 0.134 -3.02 pval.burden 1 0.02692 2 0.00251 > > ## with twins > example.ped[10,"sex"] <- 2 > #example.relation[2,4] <- 1 > #colnames(example.relation) <- c("ped","id1","id2", "code") > #data(example.relation) > #example.relation > pg.out.m2.twins <- pedgene(ped=example.ped, geno=example.geno, relation=example.relation, + map=example.map, male.dose=2, weights.mb=TRUE, checkpeds=TRUE) > > # summary/print and plot methods for this object > print.pedgene(pg.out.m2.twins,digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 80.1 0.389 -2.21 2 AX X 7 3 211.5 0.134 -3.02 pval.burden 1 0.02692 2 0.00251 > > ## base case, beta weights, no pedcheck > pg.beta.m2 <- pedgene(ped=example.ped, geno=example.geno, map=example.map, male.dose=2, verbose.return=TRUE) > names(pg.beta.m2) [1] "pgdf" "call" "save" > lapply(pg.beta.m2$save, dim) $geno [1] 18 20 $ped [1] 18 7 $map [1] 20 2 > > print(pg.beta.m2, digits=4) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 144.8 0.3498 -2.048 2 AX X 7 3 124.1 0.2038 -2.168 pval.burden 1 0.04061 2 0.03019 > > ## base case, mb weights, method=kounen, no pedcheck > pg.kounen.m2 <- pedgene(ped=example.ped, geno=example.geno, map=example.map, male.dose=2, weights.mb=TRUE,method="kounen") > print(pg.kounen.m2,digits=4) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 80.1 0.3894 -2.213 2 AX X 7 3 211.5 0.1342 -3.023 pval.burden 1 0.026925 2 0.002506 > > ## try making ped1 disconeeded by taking 2nd-generation parents away > ## this creates an error with latest kinship2, so removing > #pg.out.m2.rm34 <- pedgene(example.ped[-(3:4),], example.geno, example.map, male.dose=2, checkpeds=FALSE, weights.mb=TRUE) > #pg.out.m2.rm34 > > ## Test character ids, which is robust now because we're now making super-ids by > ## pasting ped-person together within the function > char.ped <- with(example.ped, data.frame(famid=as.character(famid), person=as.character(person), father=as.character(father), mother=as.character(mother), sex=sex, trait=trait)) > > > ## as long as subject and ped ids are character, not factor, this will work > ## pedgene makes sure to not treat character as factor > pg.out.m2.char <- pedgene(char.ped, example.geno, example.map, male.dose=2, checkpeds=FALSE) > pg.out.m2.char gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 144.8264 0.3498380 -2.047518 2 AX X 7 3 124.0772 0.2037578 -2.167566 pval.burden 1 0.04060727 2 0.03019169 > > ## show that it accepts 23 as X, but recodes 23 to X within the function > map23 <- example.map > map23$chrom[map23$chrom=="X"] <- 23 > pg.X23.m2 <- pedgene(ped=example.ped, geno=example.geno, map=map23, male.dose=2, + weights=NULL, checkpeds=TRUE) > > print(pg.X23.m2, digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 145 0.350 -2.05 2 AX X 7 3 124 0.204 -2.17 pval.burden 1 0.0406 2 0.0302 > > > ## geno row with all NA > geno.narow <- example.geno > geno.narow[4,3:ncol(example.geno)] <- NA > # to check if male dose>1 for males on X chrom -- works > #geno.narow[3,3:ncol(example.geno)] <- ifelse(geno.narow[3,2:ncol(example.geno)]==0,0,2) > pg.narow.m2 <- pedgene(ped=example.ped, geno=geno.narow, map=example.map, male.dose=2, + weights=NULL, weights.mb=TRUE,checkpeds=TRUE) > print(pg.narow.m2,digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 66.8 0.445 -2.08 2 AX X 7 3 169.6 0.197 -2.82 pval.burden 1 0.03733 2 0.00477 > > ## choose marker 4 as the 1-snp to represent the gene > ## single-snp genes reduce kernel test to burden, check that p-vals agree, stat.kern=stat.burd^2 > ## This also caused indexing problems in the past. > pg.g1.m2 <- pedgene(ped=example.ped, geno=example.geno[,c(1:2,6,13:22)], + map=example.map[c(1,11:20),], male.dose=2,weights.mb=TRUE) > pg.g1.m2 gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 1 0 1.694118 0.1930591 -1.301583 2 AX X 7 3 211.482122 0.1342214 -3.022656 pval.burden 1 0.193059061 2 0.002505672 > > # male dose=1 > pg.out.m1 <- pedgene(example.ped, example.geno, example.map, male.dose=1 ) > > print(pg.out.m1, digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 145 0.35 -2.05 2 AX X 7 3 31 0.33 -1.78 pval.burden 1 0.0406 2 0.0758 > > > ## test with no map arg (all variants in one gene columns 3:12) > pg.out.nomap <- pedgene(example.ped, example.geno[,1:12]) > pg.out.nomap gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 unknown unknown 7 3 144.8264 0.349838 -2.047518 pval.burden 1 0.04060727 > > ## test with extra subject in geno, make sure it is removed > example2.geno <- rbind(example.geno[1,],example.geno) > pg.out <- pedgene(ped=example.ped, geno=example2.geno, map=example.map, male.dose=2, + weights.mb=TRUE,checkpeds=TRUE, method=NA) > warnings() > example2.geno[1,1:2] <- c(10,5) > pg.out <- pedgene(ped=example.ped, geno=example2.geno, map=example.map, male.dose=2, + weights.mb=TRUE,checkpeds=TRUE) > warnings() > pg.out gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 80.10206 0.3894249 -2.212604 2 AX X 7 3 211.48212 0.1342214 -3.022656 pval.burden 1 0.026924950 2 0.002505672 > > ## Testing first gene with dose=2-dose > geno.recode <- cbind(example.geno[,1:2], 2-example.geno[,grep("AA", names(example.geno))]) > pg.recode.mb <- pedgene(example.ped, geno.recode, male.dose=2, weights.mb=TRUE) > ## note when map not given, assumes all 1 gene, and assigns "unknown" gene/chrom > pg.recode.mb gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 unknown unknown 7 3 80.10206 0.3894249 2.212604 pval.burden 1 0.02692495 > > pg.recode.beta <- pedgene(example.ped, geno.recode, male.dose=2) > ## note when map not given, assumes all 1 gene, and assigns "unknown" gene/chrom > pg.recode.beta gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 unknown unknown 7 3 2.472408e-49 0.9986723 0.001904113 pval.burden 1 0.9984807 > > > ## weights, Madsen-Browning > maf <- colMeans(example.geno[,-(1:2)]/2) > ## maf not correct for X matrix, b/c n-alleles for males is not 2 > ## so these results will be a little different for X-chrom > > pg.out.wts.m2 <- pedgene(example.ped, example.geno, map=example.map, + weights=1/sqrt((maf*(1-maf)))) > # note stat, pval for AX gene don't match pg.out.m2 > print(pg.out.wts.m2) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 80.10206 0.3894249 -2.212604 2 AX X 7 3 330.53325 0.1336853 -3.019160 pval.burden 1 0.026924950 2 0.002534764 > > > ## one column genotype > pg.out.1snp <- pedgene(example.ped, example.geno[,c(1,2,4),drop=FALSE], map=example.map[2,,drop=FALSE]) > pg.out.1snp gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 1 0 0 1 0 pval.burden 1 1 > > ## plot, consider using the unrelated kernel-clustering plot method to show > ## regions of clustering more than expected, > ## plot gene regions separately > > > ## Testing many genes at once: > > genobig <- example.geno > mapbig <- example.map > for(k in 2:10) { + genobig <- cbind(genobig, example.geno[,-(1:2)]) + mapbig <- rbind(mapbig, example.map) + mapbig$gene[((k-1)*20+1):(20*k)] <- paste(example.map$gene[1:20],k,sep="") + } > > ## Add two genes: one with 1 variant. Another with no markers with variance > genobig <- cbind(genobig, example.geno[,6], rep(1, nrow(example.geno)), rep(2, nrow(example.geno))) > mapbig <- rbind(mapbig, c(10, "onevar"), c(11,"novar"), c(11, "novar")) > > > pgbig.m2 <- pedgene(example.ped, genobig, mapbig, male.dose=2, acc.davies=1e-4) gene: ' novar ' has no markers after removing markers with all same genotype > pgbig.m1 <- pedgene(example.ped, genobig, mapbig, male.dose=1, acc.davies=1e-4) gene: ' novar ' has no markers after removing markers with all same genotype > > print(pgbig.m2, digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 144.83 0.350 -2.05 2 AX X 7 3 124.08 0.204 -2.17 3 AA2 1 7 3 144.83 0.350 -2.05 4 AX2 X 7 3 124.08 0.204 -2.17 5 AA3 1 7 3 144.83 0.350 -2.05 6 AX3 X 7 3 124.08 0.204 -2.17 7 AA4 1 7 3 144.83 0.350 -2.05 8 AX4 X 7 3 124.08 0.204 -2.17 9 AA5 1 7 3 144.83 0.350 -2.05 10 AX5 X 7 3 124.08 0.204 -2.17 11 AA6 1 7 3 144.83 0.350 -2.05 12 AX6 X 7 3 124.08 0.204 -2.17 13 AA7 1 7 3 144.83 0.350 -2.05 14 AX7 X 7 3 124.08 0.204 -2.17 15 AA8 1 7 3 144.83 0.350 -2.05 16 AX8 X 7 3 124.08 0.204 -2.17 17 AA9 1 7 3 144.83 0.350 -2.05 18 AX9 X 7 3 124.08 0.204 -2.17 19 AA10 1 7 3 144.83 0.350 -2.05 20 AX10 X 7 3 124.08 0.204 -2.17 21 onevar 10 1 0 1.69 0.193 -1.30 22 novar 11 0 2 NA NA NA pval.burden 1 0.0406 2 0.0302 3 0.0406 4 0.0302 5 0.0406 6 0.0302 7 0.0406 8 0.0302 9 0.0406 10 0.0302 11 0.0406 12 0.0302 13 0.0406 14 0.0302 15 0.0406 16 0.0302 17 0.0406 18 0.0302 19 0.0406 20 0.0302 21 0.1931 22 NA > print(pgbig.m1, digits=3) gene chrom n.variant n.noninform stat.kernel pval.kernel stat.burden 1 AA 1 7 3 144.83 0.350 -2.05 2 AX X 7 3 31.01 0.330 -1.78 3 AA2 1 7 3 144.83 0.350 -2.05 4 AX2 X 7 3 31.01 0.330 -1.78 5 AA3 1 7 3 144.83 0.350 -2.05 6 AX3 X 7 3 31.01 0.330 -1.78 7 AA4 1 7 3 144.83 0.350 -2.05 8 AX4 X 7 3 31.01 0.330 -1.78 9 AA5 1 7 3 144.83 0.350 -2.05 10 AX5 X 7 3 31.01 0.330 -1.78 11 AA6 1 7 3 144.83 0.350 -2.05 12 AX6 X 7 3 31.01 0.330 -1.78 13 AA7 1 7 3 144.83 0.350 -2.05 14 AX7 X 7 3 31.01 0.330 -1.78 15 AA8 1 7 3 144.83 0.350 -2.05 16 AX8 X 7 3 31.01 0.330 -1.78 17 AA9 1 7 3 144.83 0.350 -2.05 18 AX9 X 7 3 31.01 0.330 -1.78 19 AA10 1 7 3 144.83 0.350 -2.05 20 AX10 X 7 3 31.01 0.330 -1.78 21 onevar 10 1 0 1.69 0.193 -1.30 22 novar 11 0 2 NA NA NA pval.burden 1 0.0406 2 0.0758 3 0.0406 4 0.0758 5 0.0406 6 0.0758 7 0.0406 8 0.0758 9 0.0406 10 0.0758 11 0.0406 12 0.0758 13 0.0406 14 0.0758 15 0.0406 16 0.0758 17 0.0406 18 0.0758 19 0.0406 20 0.0758 21 0.1931 22 NA > > > proc.time() user system elapsed 1.64 0.23 1.87