Version:	2.3.1
Date:	2026-04-13
Depends:	R (≥ 4.0.0)
Imports:	mvtnorm, parallel, survival, tidyr, plyr, ggplot2, poisbinom, rms, methods
Suggests:	testthat, knitr, rmarkdown, biomaRt, VariantAnnotation, GenomicRanges, IRanges, S4Vectors, org.Hs.eg.db, TxDb.Hsapiens.UCSC.hg19.knownGene
Title:	SNPs-Based Whole Genome Association Studies
Description:	Functions to perform most of the common analysis in genome association studies are implemented. These analyses include descriptive statistics and exploratory analysis of missing values, calculation of Hardy-Weinberg equilibrium, analysis of association based on generalized linear models (either for quantitative or binary traits), and analysis of multiple SNPs (haplotype and epistasis analysis). Permutation test and related tests (sum statistic and truncated product) are also implemented. Max-statistic and genetic risk-allele score exact distributions are also possible to be estimated. The methods are described in Gonzalez JR et al., 2007 <doi:10.1093/bioinformatics/btm025>. This version includes internal copies of functions from the archived 'haplo.stats' package to maintain functionality.
URL:	https://github.com/isglobal-brge/SNPassoc
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
VignetteBuilder:	knitr
RoxygenNote:	7.2.2
NeedsCompilation:	yes
Packaged:	2026-04-13 08:38:52 UTC; dpelegri
Author:	Victor Moreno [aut], Juan R Gonzalez [aut], Dolors Pelegri [aut, cre]
Maintainer:	Dolors Pelegri <dolors.pelegri@isglobal.org>
Repository:	CRAN
Date/Publication:	2026-04-13 10:10:09 UTC

Bonferroni correction of p values

Description

This function shows the SNPs that are statistically significant after correcting for the number of tests performed (Bonferroni correction) for an object of class "WGassociation"

Usage

 Bonferroni.sig(x, model = "codominant", alpha = 0.05, 
      include.all.SNPs=FALSE)

Arguments

Details

After deciding the genetic model, the function shows the SNPs that are statistically significant at alpha level corrected by the number of performed tests.

Value

A data frame with the SNPs and the p values for those SNPs that are statistically significant after Bonferroni correction

See Also

WGassociation

Examples

data(SNPs)
datSNP<-setupSNP(SNPs,6:40,sep="")
ans<-WGassociation(protein~1,data=datSNP,model="all")
Bonferroni.sig(ans, model="codominant", alpha=0.05, include.all.SNPs=FALSE)

Population substructure

Description

This function estimates an inflation (or deflation) factor, lambda, as indicated in the paper by Devlin et al. (2001) and corrects the p-values using this factor.

Usage

GenomicControl(x, snp.sel)

Arguments

Details

This method is only valid for 2x2 tables. This means that the object of class 'WGassociation' might not have fitted the codominant model.

See reference for further details.

Value

The same object of class 'WGassociation' where the p-values have been corrected for genomic control.

References

B Devlin, K Roeder, and S.A. Bacanu. Unbiased Methods for Population Based Association Studies. Genetic Epidemiology (2001) 21:273-84

See Also

qqpval, WGassociation

Examples

data(SNPs) 
datSNP<-setupSNP(SNPs,6:40,sep="")
res<-WGassociation(casco,datSNP,model=c("do","re","log-add"))

# Genomic Control 
resCorrected<-GenomicControl(res)

Compute Generalized Inverse of Input Matrix

Description

Singular value decomposition (svd) is used to compute a generalized inverse of input matrix.

Usage

Ginv(x, eps = 1e-06)

Arguments

Details

The svd function uses the LAPACK standard library to compute the singular values of the input matrix. The rank of the matrix is determined by the number of singular values that are at least as large as max(svd)*eps.

Value

A list with the following components:

References

Press WH, Teukolsky SA, Vetterling WT, Flannery BP. (1992). Numerical recipes in C. The art of scientific computing. 2nd ed. Cambridge University Press, Cambridge. page 61.

Anderson, E., et al. (1994). LAPACK User's Guide, 2nd edition, SIAM, Philadelphia.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

svd

Examples

x <- matrix(c(1, 2, 1, 2, 3, 2), ncol = 3)
save <- Ginv(x)
ginv.x <- save$Ginv
rank.x <- save$rank

SNPs from HapMap project

Description

Information about 9307 SNPs from the HapMap project belonging to 22 chromosomes. Information about two different population is available: European population (CEU) and Yoruba (YRI). The genomic information (names of SNPs, chromosomes and genetic position) is also available in a data frame called 'HapMap.SNPs.pos'.

Usage

data(HapMap)

Format

A data frame with 120 observations on the 9808 variables (SNPs) and one variable called 'group' indicating the population.

Source

HapMap project (http://www.hapmap.org)

Examples

data(HapMap)

SNPs from HapMap project

Description

Information about 9307 SNPs from the HapMap project belonging to 22 chromosomes. Information about two different population is available: European population (CEU) and Yoruba (YRI). The genomic information (names of SNPs, chromosomes and genetic position) is also available in a data frame called 'HapMap.SNPs.pos'.

Usage

data(HapMap.SNPs.pos)

Format

A data frame with 120 observations on the 9808 variables (SNPs) and one variable called 'group' indicating the population.

Source

HapMap project (http://www.hapmap.org)

Examples

    data(HapMap.SNPs.pos)

max-statistic for a 2x3 table

Description

Compute pairwise linkage disequilibrium between genetic markers

Usage

LD(g1, ...)

## S3 method for class 'snp'
LD(g1, g2, ...)

## S3 method for class 'setupSNP'
LD(g1, SNPs, ...)


LDplot(x, digits = 3, marker, distance, which = c("D", "D'",
    "r", "X^2", "P-value", "n", " "), ...)

LDtable(x, colorcut = c(0, 0.01, 0.025, 0.05, 0.1, 1), 
       colors = heat.colors(length(colorcut)),
       textcol = "black", digits = 3, show.all = FALSE, 
       which = c("D", "D'", "r", "X^2", "P-value", "n"), 
       colorize = "P-value", cex, ...)

Arguments

Value

None

Author(s)

functions adapted from LD, LDtable and LDplot in package genetics by Gregory Warnes et al. (warnes@bst.rochester.edu)

References

genetics R package by Gregory Warnes et al. (warnes@bst.rochester.edu)

See Also

setupSNP snp

Internal SNPstat functions

Description

Internal SNPassoc functions

Usage

association.fit(var, dep, adj, quantitative, type, level, 
       nIndiv, genotypingRate = 0, ...) 
extractPval(x)
extractPval.i(i,x,pos,models)
SNPHWE(x)
GenotypeRate(x)
haplo.inter.fit(geno, var2, dep, adj = NULL, fam, 
          haplo.freq.min, ...)
crea.lab(x,pos.ini,cex,dist)
orderChromosome(x)
togeno(f,sep=sep,lab=lab)
expandsetupSNP(o)
pvalTest(dataX,Y,quantitative,type,genotypingRate)
modelTest(X,Y,quantitative,type,genotypingRate)
assoc(y,x,test="lrt",quantitative)
trim(s)
interleave(..., append.source=TRUE, sep=": ", drop=FALSE)
## Default S3 method:
codominant(o)
## Default S3 method:
dominant(o)
## Default S3 method:
recessive(o)
## Default S3 method:
overdominant(o)
## Default S3 method:
additive(o)

Details

These are not to be called by the user

Value

No return value, internal calls

SNPs in a case-control study

Description

SNPs data.frame contains selected SNPs and other clinical covariates for cases and controls in a case-control study

SNPs.info.pos data.frame contains the names of the SNPs included in the data set 'SNPs' including their chromosome and their genomic position

Usage

data(SNPs)

Format

'SNPs' data.frame contains the following columns:

Source

Data obtained from our department. The reference and details will be supplied after being published.

SNPs in a case-control study

Description

SNPs data.frame contains selected SNPs and other clinical covariates for cases and controls in a case-control study

SNPs.info.pos data.frame contains the names of the SNPs included in the data set 'SNPs' including their chromosome and their genomic position

Usage

data(SNPs.info.pos)

Format

'SNPs.info.pos' data.frame contains the following columns: A data frame with 35 observations on the following 3 variables.

snp

name of SNP

chr

name of chromosome

pos

genomic position

Source

Data obtained from our department. The reference and details will be supplied after being published.

Descriptive sample size and percentage

Description

This function computes sample size and percentage for each category of a categorical trait (e.g. case-control status) for each genotype (or combination of genotypes).

Usage

Table.N.Per(var, dep, subset = !is.na(var))

Arguments

Value

See Also

Table.mean.se

Examples

data(SNPs)
#sample size and percentage of cases and controls for each genotype 
Table.N.Per(SNPs$snp10001,SNPs$casco)

# The same table for a subset (males)
Table.N.Per(SNPs$snp10001,SNPs$casco,SNPs$sex=="Male")

# The same table assuming a dominant model
Table.N.Per(dominant(snp(SNPs$snp10001,sep="")),SNPs$casco,SNPs$sex=="Male")

Descriptive sample size, mean, and standard error

Description

This function computes sample size, mean and standard error of a quantitative trait for each genotype (or combination of genotypes)

Usage

Table.mean.se(var, dep, subset = !is.na(var))

Arguments

Value

See Also

Table.N.Per

Examples

data(SNPs)
# sample size, mean age and standard error for each genotype
Table.mean.se(SNPs$snp10001,SNPs$protein)

# The same table for a subset (males)
Table.mean.se(SNPs$snp10001,SNPs$protein,SNPs$sex=="Male")

# The same table assuming a dominant model
Table.mean.se(dominant(snp(SNPs$snp10001,sep="")),SNPs$protein,SNPs$sex=="Male")

Whole genome association analysis

Description

This function carries out a whole genome association analysis between the SNPs and a dependent variable (phenotype) under five different genetic models (inheritance patterns): codominant, dominant, recessive, overdominant and log-additive. The phenotype may be quantitative or categorical. In the second case (e.g. case-control studies) this variable must be of class 'factor' with two levels.

Usage

   WGassociation(formula, data, model = c("all"), 
                 quantitative = is.quantitative(formula, data),
                 genotypingRate = 80, level = 0.95, ...)

Arguments

Details

This function assesses the association between the response variable included in the left side in the 'formula' and the SNPs included in the 'data' argument adjusted by those variables included in the right side of the 'formula'. Different genetic models may be analyzed using 'model' argument.

Value

An object of class 'WGassociation'.

'summary' returns a summary table by groups defined in info (genes/chromosomes).

'WGstats' returns a detailed output, similar to the produced by association.

'pvalues' and 'print' return a table of p-values for each genetic model for each SNP. The first column indicates whether a problem with genotyping is present.

'plot' produces a plot of p values in the -log scale. See plot.WGassociation for further details.

'labels' returns the names of the SNPs analyzed.

The functions 'codominat', 'dominant', 'recessive', 'overdominant' and 'additive' are used to obtain the p values under these genetic models.

See examples for further illustration about all previous issues.

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

getSignificantSNPs association WGstats setupSNP plot.WGassociation

Examples

data(SNPs)
datSNP<-setupSNP(SNPs,6:40,sep="")
ansAll<-WGassociation(protein~1,data=datSNP,model="all")

# In that case the formula is not required. You can also write:
# ansAll<-WGassociation(protein,data=datSNP,model="all")


#only codominant and log-additive
ansCoAd<-WGassociation(protein~1,data=datSNP,model=c("co","log-add"))

#for printing p values
print(ansAll)
print(ansCoAd)

#for obtaining a matrix with the p palues
pvalAll<-pvalues(ansAll)
pvalCoAd<-pvalues(ansCoAd)

# when all models are fitted and we are interested in obtaining 
# p values for different genetic models

# codominant model
pvalCod<-codominant(ansAll)

# recessive model
pvalRec<-recessive(ansAll)

# and the same for additive, dominant or overdominant


#summary
summary(ansAll)

#for a detailed report
WGstats(ansAll)

#for plotting the p values
plot(ansAll)

Association analysis between a single SNP and a given phenotype

Description

This function carries out an association analysis between a single SNP and a dependent variable (phenotype) under five different genetic models (inheritance patterns): codominant, dominant, recessive, overdominant and log-additive. The phenotype may be quantitative or categorical. In the second case (e.g. case-control studies) this variable must be of class 'factor' with two levels.

Usage

association(formula, data, model=c("all"), model.interaction=
       c("codominant"), subset, name.snp = NULL, quantitative = 
       is.quantitative(formula,data), genotypingRate= 0, 
       level = 0.95, ...)

Arguments

Details

This function should be called by the user when we are interested in analyzing an unique SNP. It is recommended to use WGassociation function when more than one SNP is studied.

Value

For each genetic model (codominant, dominant, recessive, overdominant, and log-additive) the function gives a matrix with sample size and percentages for each genotype, the Odds Ratio and its 95% confidence interval (taking the most frequent homozygous genotype as the reference), the p-value corresponding to the likelihood ratio test obtained from a comparison with the null model, and the Akaike Information Criterion (AIC) of each genetic model. In the case of analyzing a quantitative trait, the function returns a matrix with sample size, mean and standard errors for each genotype, mean difference and its 95% confidence interval with respect to the most frequent homozygous genotype, the p-value obtained from an overall gene effect and the Akaike Information Criterion (AIC) of each genetic model.

When an interaction term (a categorical covariate with an SNP) is included in the model, three different tables are given. The first one correponds to the full interaction matrix where the ORs (or mean differences if a quantitative trait is analyzed) are expressed with respect to the non variant genotype and the first category of the covariate. The other two tables show the ORs and their 95% confidence intervals for both marginal models. P values for interaction and trend are also showed in the output.

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

Iniesta R, Guino E, Moreno V. Statistical analysis of genetic polymorphisms in epidemiological studies. Gac Sanit. 2005;19(4):333-41.

Elston RC. Introduction and overview. Statistical methods in genetic epidemiology. Stat Methods Med Res. 2000;9:527-41.

See Also

WGassociation

Examples

data(SNPs)

# first, we create an object of class 'setupSNP'
datSNP<-setupSNP(SNPs,6:40,sep="")

# case-control study, crude analysis
association(casco~snp10001, data=datSNP)

# case-control study, adjusted by sex and arterial blood pressure
association(casco~sex+snp10001+blood.pre, data=datSNP)


# quantitative trait, crude analysis
association(log(protein)~snp10001,data=datSNP)
# quantitative trait, adjusted by sex
association(log(protein)~snp10001+sex,data=datSNP)


#
# Interaction analysis
#

# Interaction SNP and factor
association(log(protein)~snp10001*sex+blood.pre, data=datSNP, 
          model="codominant")

# Interaction SNP and SNP (codominant and codominant)
association(log(protein)~snp10001*factor(snp10002)+blood.pre, 
          data=datSNP, model="codominant")

# Interaction SNP and SNP (dominant and recessive)
association(log(protein)~snp10001*factor(recessive(snp100019))+blood.pre, 
          data=datSNP, model="dominant")

SNP data on asthma case-control study

Description

data.frame with 51 SNPs and 6 epidemiological variables: country, gender, age, bmi, smoke ans case/control status.

Usage

data("asthma")

Format

The asthma data frame has 1578 rows (individuals) and 57 columns (variables) of data from a genetic study on asthma.

Value

An data.frame object.

Examples

data(asthma)
dim(asthma)
str(asthma)

Get gene symbol from a list of SNPs

Description

Get gene symbol from a list of SNPs

Usage

    getGeneSymbol(
      x,
      snpCol = 1,
      chrCol = 2,
      posCol = 3,
      db = TxDb.Hsapiens.UCSC.hg19.knownGene
    )

Arguments

Value

a data.frame having initial information and gene symbol

Examples

    
        snps = c('rs58108140','rs189107123','rs180734498','rs144762171')
        chr = c('chr1','chr1','chr1','chr1')
        pos = c(10583, 10611, 13302, 13327)
        
        x <- data.frame(snps, chr, pos )
        if (requireNamespace("VariantAnnotation", quietly = TRUE)) {
            getGeneSymbol(x)
        }
    

Get Latex output

Description

Create Latex output from association analyses

Usage

getNiceTable( x )

Arguments

Value

The R output of specific association analyses exported into LaTeX

Extract significant SNPs from an object of class 'WGassociation'

Description

Extract significant SNPs from an object of class 'WGassociation' when genomic information is available

Usage

getSignificantSNPs(x, chromosome, model, sig = 1e-15)

Arguments

Value

A list with the following components:

...

See Also

WGassociation

Examples

data(resHapMap)
# resHapMap contains the results for a log-additive genetic model

# to get the significant SNPs for chromosome 12
getSignificantSNPs(resHapMap,chromosome=12)
# to get the significant SNPs for chromosome 5
getSignificantSNPs(resHapMap,5)
# to get the significant SNPs for chromosome X at level 1e-8
getSignificantSNPs(resHapMap,5,sig=1e-8)

EM Computation of Haplotype Probabilities, with Progressive Insertion of Loci

Description

For genetic marker phenotypes measured on unrelated subjects, with linkage phase unknown, compute maximum likelihood estimates of haplotype probabilities. Because linkage phase is unknown, there may be more than one pair of haplotypes that are consistent with the oberved marker phenotypes, so posterior probabilities of pairs of haplotypes for each subject are also computed. Unlike the usual EM which attempts to enumerate all possible pairs of haplotypes before iterating over the EM steps, this "progressive insertion" algorithm progressively inserts batches of loci into haplotypes of growing lengths, runs the EM steps, trims off pairs of haplotypes per subject when the posterior probability of the pair is below a specified threshold, and then continues these insertion, EM, and trimming steps until all loci are inserted into the haplotype. The user can choose the batch size. If the batch size is chosen to be all loci, and the threshold for trimming is set to 0, then this algorithm reduces to the usual EM algorithm.

Usage

haplo.em(geno, locus.label=NA, miss.val=c(0, NA), weight, control= 
           haplo.em.control())

Arguments

Details

The basis of this progressive insertion algorithm is from the sofware snphap by David Clayton. Although some of the features and control parameters of this S-PLUS version are modeled after snphap, there are substantial differences, such as extension to allow for more than two alleles per locus, and some other nuances on how the alogrithm is implemented.

Value

list with components:

Note

Sorted order of haplotypes with character alleles is system-dependent, and can be controlled via the LC_COLLATE locale environment variable. Different locale settings can cause results to be non-reproducible even when controlling the random seed. This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

setupGeno, haplo.em.control

Examples

data(hla.demo)
attach(hla.demo)
geno <- hla.demo[,c(17,18,21:24)]
label <-c("DQB","DRB","B")
keep <- !apply(is.na(geno) | geno==0, 1, any)

save.em.keep <- haplo.em(geno=geno[keep,], locus.label=label)

# warning: output will not exactly match

print(save.em.keep)

Create the Control Parameters for the EM Computation of Haplotype Probabilities, with Progressive Insertion of Loci

Description

Create a list of parameters that control the EM algorithm for estimating haplotype frequencies, based on progressive insertion of loci. Non-default parameters for the EM algorithm can be set as parameters passed to haplo.em.control.

Usage

haplo.em.control(loci.insert.order=NULL, insert.batch.size = 6,
                             min.posterior = 1e-09, tol = 1e-05,
                             max.iter=5000, random.start=0, n.try = 10,
                             iseed=NULL, max.haps.limit=2e6, verbose=0)

Arguments

Details

The default is to use n.try = 10. If this takes too much time, it may be worthwhile to decrease n.try. Other tips for computing haplotype frequencies for a large number of loci, particularly if some have many alleles, is to decrease the batch size (insert.batch.size), increase the memory (max.haps.limit), and increase the probability of trimming off rare haplotypes at each insertion step (min.posterior).

Value

A list of the parameters passed to the function.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.em, haplo.score

Examples

# This is how it is used within haplo.score
#    > score.gauss <- haplo.score(resp, geno, trait.type="gaussian", 
#    >           em.control=haplo.em.control(insert.batch.size = 2, n.try=1))

GLM Regression of Trait on Ambiguous Haplotypes

Description

Perform glm regression of a trait on haplotype effects, allowing for ambiguous haplotypes. This method performs an iterative two-step EM, with the posterior probabilities of pairs of haplotypes per subject used as weights to update the regression coefficients, and the regression coefficients used to update the posterior probabilities.

Usage

haplo.glm(formula=formula(data), family=gaussian, data=parent.frame(),
           weights, na.action="na.geno.keep", start=NULL, 
           locus.label=NA, control=haplo.glm.control(), 
           method="glm.fit", model=TRUE, x=FALSE, y=TRUE, 
           contrasts=NULL, ...)

Arguments

Details

To properly prepare the data frame, the genotype matrix must be processed by setupGeno, and then included in the data frame with the response and other variables.

For binomial family, the initialization of values gives warnings if non-integer number of successes, which is a concern in these models because of the weights of posterior probability of each haplotype pair per subject. We supress the warnings by defining a haplo.binomial family, which we use if family=binomial is used.

Value

An object of class "haplo.glm" is returned. The output object from haplo.glm has all the components of a glm object, with a few more. It is important to note that some of the returned components correpond to the "expanded" version of the data. This means that each observation is expanded into the number of terms in the observation's posterior distribution of haplotype pairs, given the marker data. For example, when fitting the response y on haplotype effects, the value of y[i], for the ith observation, is replicated m[i] times, where m[i] is the number of pairs of haplotypes consistent with the observed marker data. The returned components that are expanded are indicated below by [expanded] in the definition of the component.

These expanded components may need to be collapsed, depending on the objective of the user. For example, when considering the influence of an observation, it may make sense to examine the expanded residuals for a single observation, perhaps plotted against the haplotypes for that observation. In contrast, it would not be sensible to plot all residuals against non-genetic covariates, without first collapsing the expanded residuals for each observation. To collapse, one can use the average residual per observation, weighted according to the posterior probabilities. The appropriate weight can be computed as wt = weight.expanded * haplo.post.info[[post]]. Then, the weighted average can be calculated as with(fit, tapply(residuals * wt, haplo.post.info[["indx"]], sum).

References

Lake S, Lyon H, Silverman E, Weiss S, Laird N, Schaid D (2002) Estimation and tests of haplotype-environment interaction when linkage phase is ambiguous. Human Heredity 55:56-65.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.glm.control, haplo.em, haplo.model.frame

Examples

cat(" FOR REGULAR USAGE, DO NOT DISCARD GENOTYPES WITH MISSING VALUES\n")
cat(" WE ONLY SUBSET BY keep HERE SO THE EXAMPLES RUN FASTER\n")

 data(hla.demo)
 geno <- as.matrix(hla.demo[,c(17,18,21:24)])
 keep <- !apply(is.na(geno) | geno==0, 1, any) # SKIP THESE THREE LINES
 hla.demo <- hla.demo[keep,]                   # IN AN ANALYSIS
 geno <- geno[keep,]                           # 
 attach(hla.demo)
 label <-c("DQB","DRB","B")
 y <- hla.demo$resp
 y.bin <- 1*(hla.demo$resp.cat=="low")

# set up a genotype array as a model.matrix for inserting into data frame
# Note that hla.demo is a data.frame, and we need to subset to columns
# of interest. Also also need to convert to a matrix object, so that
# setupGeno can code alleles and convert geno to 'model.matrix' class.

 geno <- setupGeno(geno, miss.val=c(0,NA))

  # geno now has an attribute 'unique.alleles' which must be passed to
  # haplo.glm as allele.lev=attributes(geno)$unique.alleles, see below

 my.data <- data.frame(geno=geno, age=hla.demo$age, male=hla.demo$male,
                      y=y, y.bin=y.bin)

 fit.gaus <- haplo.glm(y ~ male + geno, family = gaussian,  na.action=
               "na.geno.keep",allele.lev=attributes(geno)$unique.alleles, 
               data=my.data, locus.label=label,
               control = haplo.glm.control(haplo.freq.min=0.02))
 fit.gaus
 

Create list of control parameters for haplo.glm

Description

Create a list of control pararameters for haplo.glm. If no parameters are passed to this function, then all default values are used.

Usage

haplo.glm.control(haplo.effect="add", haplo.base=NULL,
                  haplo.min.count=NA, haplo.freq.min=.01,
                  sum.rare.min=0.001, haplo.min.info=0.001, 
                  keep.rare.haplo=TRUE,
                  eps.svd=sqrt(.Machine$double.eps),
                  glm.c=glm.control(maxit=500), 
                  em.c=haplo.em.control())

Arguments

Value

the list of above components

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.glm, haplo.em.control, glm.control

Examples

# NOT RUN
# using the data set up in the example for haplo.glm,
# the control function is used in haplo.glm as follows
#  > fit <- haplo.glm(y ~ male + geno, family = gaussian,  
#  >          na.action="na.geno.keep",
#  >          data=my.data, locus.label=locus.label,
#  >          control = haplo.glm.control(haplo.min.count=5,
#  >          em.c=haplo.em.control(n.try=1)))

Haplotype interaction with a covariate

Description

This function computes the ORs (or mean differences if a quantitative trait is analyzed) and their 95% confidence intervals corresponding to an interaction between the haplotypes and a categorical covariate

Usage

haplo.interaction(formula, data, SNPs.sel, quantitative = 
  is.quantitative(formula, data), haplo.freq.min = 0.05, ...)

Arguments

Details

The function estimates the haplotypes for the SNPs indicated in the 'SNPs.sel' argument. Then, usign 'haplo.glm' function (from 'haplo.stats' library) estimates the interaction between these haplotypes and the covariate indicated in the formula by means of 'interaction' function.

Value

Three different tables are given. The first one corresponds to the full interaction matrix where the ORs (or mean differences if a quantitative trait is analyzed) are expressed with respect to the most frequent haplotype and the first category of the covariate. The other two tables show the ORs (or mean differences if a quantitative trait is analyzed) and their 95% confidence intervals for both marginal models. P values for interaction are also showed in the output.

Examples

# not Run
library(SNPassoc)

data(SNPs)
datSNP<-setupSNP(SNPs,6:40,sep="")
res <- haplo.interaction(log(protein)~int(sex), data=datSNP,
                    SNPs.sel=c("snp100019","snp10001","snp100029"))
res

Sets up a model frame for haplo.glm

Description

For internal use within the haplo.stats library

Usage

haplo.model.frame(m, locus.label=NA, control=haplo.glm.control())

Arguments

Details

See haplo.glm description in help file and user manual

Value

A model frame with haplotypes modeled as effects

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

Score Statistics for Association of Traits with Haplotypes

Description

Compute score statistics to evaluate the association of a trait with haplotypes, when linkage phase is unknown and diploid marker phenotypes are observed among unrelated subjects. For now, only autosomal loci are considered.

Usage

haplo.score(y, geno, trait.type="gaussian", offset = NA, x.adj = NA,
            min.count=5, skip.haplo=min.count/(2*nrow(geno)),
            locus.label=NA, miss.val=c(0,NA), haplo.effect="additive",
            eps.svd=1e-5, simulate=FALSE, sim.control=score.sim.control(),
            em.control=haplo.em.control())

Arguments

Details

Compute the maximum likelihood estimates of the haplotype frequencies and the posterior probabilities of the pairs of haplotypes for each subject using an EM algorithm. The algorithm begins with haplotypes from a subset of the loci and progressively discards those with low frequency before inserting more loci. The process is repeated until haplotypes for all loci are established. The posterior probabilities are used to compute the score statistics for the association of (ambiguous) haplotypes with traits. The glm function is used to compute residuals of the regression of the trait on the non-genetic covariates.

Value

List with the following components:

References

Schaid DJ, Rowland CM, Tines DE, Jacobson RM, Poland GA. "Score tests for association of traits with haplotypes when linkage phase is ambiguous." Amer J Hum Genet. 70 (2002): 425-434.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.em, haplo.em.control, score.sim.control

Examples

  # establish all hla.demo data, 
  # remove genotypes with missing alleles just so haplo.score runs faster 
  # with missing values included, this example takes 2-4 minutes
  # FOR REGULAR USAGE, DO NOT DISCARD GENOTYPES WITH MISSING VALUES

  data(hla.demo)
  geno <- as.matrix(hla.demo[,c(17,18,21:24)])
  keep <- !apply(is.na(geno) | geno==0, 1, any)
  hla.demo <- hla.demo[keep,]
  geno <- geno[keep,]
  attach(hla.demo)
  label <- c("DQB","DRB","B")
 
# For quantitative, normally distributed trait:

  score.gaus <- haplo.score(resp, geno, locus.label=label, 
                            trait.type = "gaussian")
  print(score.gaus)

# For ordinal trait:
  y.ord <- as.numeric(resp.cat)
  score.ord <- haplo.score(y.ord, geno, locus.label=label,
                           trait.type="ordinal")
  print(score.ord)

# For a  binary trait and simulations,
# limit simulations to 500 in score.sim.control, default is 20000
  y.bin <-ifelse(y.ord==1,1,0)
  score.bin.sim <- haplo.score(y.bin, geno, trait.type = "binomial",
                     locus.label=label, simulate=TRUE,
                     sim.control=score.sim.control(min.sim=200,max.sim=500))

  print(score.bin.sim)

# For a binary trait, adjusted for sex and age:

  x <- cbind(male, age)
  score.bin.adj <- haplo.score(y.bin, geno, trait.type = "binomial", 
                               locus.label=label, x.adj=x)
  print(score.bin.adj)

Score Statistics for Association of Traits with Haplotypes

Description

Used to identify sub-haplotypes from a group of loci. Run haplo.score on all contiguous subsets of size n.slide from the loci in a genotype matrix (geno). From each call to haplo.score, report the global score statistic p-value. Can also report global and maximum score statistics simulated p-values.

Usage

haplo.score.slide(y, geno, trait.type="gaussian", n.slide=2,
                  offset = NA, x.adj = NA, min.count=5,
                  skip.haplo=min.count/(2*nrow(geno)),
                  locus.label=NA, miss.val=c(0,NA),
                  haplo.effect="additive", eps.svd=1e-5,
                  simulate=FALSE, sim.control=score.sim.control(),
                  em.control=haplo.em.control())

Arguments

Details

Haplo.score.slide is useful for a series of loci where little is known of the association between a trait and haplotypes. Using a range of n.slide values, the region with the strongest association will consistently have low p-values for locus subsets containing the associated haplotypes. The global p-value measures significance of the entire set of haplotypes for the locus subset. Simulated maximum score statistic p-values indicate when one or a few haplotypes are associated with the trait.

Value

List with the following components:

References

Schaid DJ, Rowland CM, Tines DE, Jacobson RM, Poland GA. "Score tests for association of traits with haplotypes when linkage phase is ambiguous." Amer J Hum Genet. 70 (2002): 425-434.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.score, score.sim.control

Examples

  data(hla.demo)

# Continuous trait slide by 2 loci on all 11 loci, uncomment to run it.
# Takes > 20 minutes to run
#  geno.11 <- hla.demo[,-c(1:4)]
#  label.11 <- c("DPB","DPA","DMA","DMB","TAP1","TAP2","DQB","DQA","DRB","B","A")
#  slide.gaus <- haplo.score.slide(hla.demo$resp, geno.11, trait.type = "gaussian",
#                                  locus.label=label.11, n.slide=2)

#  print(slide.gaus)
#  plot(slide.gaus)

# Run shortened example on 9 loci 
# For an ordinal trait, slide by 3 loci, and simulate p-values:
#  geno.9 <- hla.demo[,-c(1:6,15,16)]
#  label.9 <- c("DPA","DMA","DMB","TAP1","DQB","DQA","DRB","B","A")

#  y.ord <- as.numeric(hla.demo$resp.cat)

# data is set up, to run, run these lines of code on the data that was
# set up in this example. It takes > 15 minutes to run
#  slide.ord.sim <-  haplo.score.slide(y.ord, geno.9, trait.type = "ordinal",
#                      n.slide=3, locus.label=label.9, simulate=TRUE,
#                      sim.control=score.sim.control(min.sim=200, max.sim=500))

  # note, results will vary due to simulations
#  print(slide.ord.sim)
#  plot(slide.ord.sim)
#  plot(slide.ord.sim, pval="global.sim")
#  plot(slide.ord.sim, pval="max.sim")

HLA Loci and Serologic Response to Measles Vaccination

Description

A data frame with genotypes at eleven HLA-region loci genotyped for 220 subjects, phase not known. Contains measles vaccination response with covariate data.

Usage

data(hla.demo)

Format

A data frame with 220 observations on the following 26 variables.

resp

numeric, Quantitative response to Measles Vaccination

resp.cat

Category of vaccination response, a factor with levels high low normal

male

numeric, indicator of gener, 1=male, 0=female

age

numeric, subject's age

DPB.a1

first allele of genotype

DPB.a2

second allele of genotype

DPA.a1

first allele of genotype

DPA.a2

second allele of genotype

DMA.a1

first allele of genotype

DMA.a2

second allele of genotype

DMB.a1

first allele of genotype

DMB.a2

second allele of genotype

TAP1.a1

first allele of genotype

TAP1.a2

second allele of genotype

TAP2.a1

first allele of genotype

TAP2.a2

second allele of genotype

DQB.a1

first allele of genotype

DQB.a2

second allele of genotype

DQA.a1

first allele of genotype

DQA.a2

second allele of genotype

DRB.a1

first allele of genotype

DRB.a2

second allele of genotype

B.a1

first allele of genotype

B.a2

second allele of genotype

A.a1

first allele of genotype

A.a2

second allele of genotype

Source

Data set kindly provided by Gregory A. Poland, M.D. and the Mayo Clinic Vaccine Research Group for illustration only, and my not be used for publication.

References

Schaid DJ, Rowland CM, Tines DE, Jacobson RM, Poland GA. "Score tests for association of traits with haplotypes when linkage phase is ambiguous." Amer J Hum Genet. 70 (2002): 425-434.

Examples

data(hla.demo)

Collapsing (or recoding) genotypes into different categories (generally two) depending on a given genetic mode of inheritance

Description

codominant function recodifies a variable having genotypes depending on the allelic frequency in descending order.

dominant, recessive, and overdominant functions collapse the three categories of a given SNP into two categories as follows: Let 'AA', 'Aa', and 'aa' be the three genotypes. After determining the most frequent allele (let's suppose that 'A' is the major allele) the functions return a vector with to categories as follows. dominant: 'AA' and 'Aa-aa'; recessive: 'AA-Aa' and 'aa'; overdominant: 'AA-aa' vs 'Aa'.

additive function creates a numerical variable, 1, 2, 3 corresponding to the three genotypes sorted out by descending allelic frequency (this model is referred as log-additive).

Usage

codominant(o)
dominant(o)
recessive(o)
overdominant(o)
additive(o)

Arguments

Value

A snp object collapsing genotypes into different categories depending on a given genetic mode of inheritance

Examples

data(SNPs)
dominant(snp(SNPs$snp10001,sep=""))
overdominant(snp(SNPs$snp10001,sep=""))

Identify interaction term

Description

This is a special function used for 'haplo.interaction' function. It identifies the variable that will interact with the haplotype estimates. Using int() in a formula implies that the interaction term between this variable and haplotypes is included in 'haplo.glm' function.

Usage

int(x)

Arguments

Value

x

See Also

haplo.interaction

Examples

library(SNPassoc)

data(SNPs)
datSNP<-setupSNP(SNPs, 6:40, sep = "")
mod <- haplo.interaction(casco~int(sex)+blood.pre, data = datSNP, 
                         SNPs.sel = c("snp10001","snp10004","snp10005"))

Two-dimensional SNP analysis for association studies

Description

Perform a two-dimensional SNP analysis (interaction) for association studies with possible allowance for covariate

Usage

interactionPval(formula, data, quantitative = 
          is.quantitative(formula, data), model = "codominant")

Arguments

Details

The 'interactionPval' function calculates, for each pair of SNPs (i,j), the likelihood underling the null model L0, the likelihood under each of the single-SNP, L(i) and L(j), the likelihood under an additive SNP model La(i,j), and the likelihood under a full SNP model (including SNP-SNP interaction), Lf(i,j).

The upper triangle in matrix from this function contains the p values for the interaction (epistasis) log-likelihood ratio test, LRT, LRTij = -2 (log Lf(i,j) - log La(i,j))

The diagonal contains the p values from LRT for the crude effect of each SNP, LRTii = -2 (log L(i) - log L0)

The lower triangle contains the p values from LRT comparing the two-SNP additive likelihood to the best of the single-SNP models, LRTji = -2 (log La(i,j) - log max(L(i),L(j)))

In all cases the models including the SNPs are adjusted by the covariates indicated in the 'formula' argument. This method is used either for quantitative traits and dicotomous variables (case-control studies).

Value

The 'interactionPval' function returns a matrix of class 'SNPinteraction' containing the p values corresponding to the different likelihood ratio tests above describe.

Methods defined for 'SNPinteraction' objects are provided for print and plot. The plot method uses 'image' to plot a grid of p values. The upper triangle contains the interaction (epistasis) p values from LRT. The content in the lower triangle is the p values from the LRT comparing the additive model with the best single model. The diagonal contains the main effects pvalues from LRT. The 'plot.SNPinteraction' function also allows the user to plot the SNPs sorted by genomic position and with the information about chromosomes as in the 'plotMissing' function.

Note

two-dimensional SNP analysis on a dense grid can take a great deal of computer time and memory.

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

setupSNP

Examples

    data(SNPs)
    datSNP <- setupSNP(SNPs, 6:40, sep="")[1:10,]
    
    ansCod <- interactionPval(log(protein)~sex,datSNP)
    print(ansCod)
    plot(ansCod)

Print ORs and 95% confidence intervals for an object of class 'haplo.glm'

Description

Print ORs and confidence intervals for an object of class 'haplo.glm'

Usage

intervals(o, level=.95, ...)

Arguments

Value

intervals object with ORs and 95% confidence intervals for an object of class 'haplo.glm'

Examples

library(SNPassoc)

data(asthma, package = "SNPassoc")

asthma.s <- setupSNP(data=asthma, colSNPs=7:ncol(asthma), sep="")
trait <- asthma.s$casecontrol
snpsH <- c("rs714588", "rs1023555",  "rs898070")
genoH <- make.geno(asthma.s, snpsH)

mod <- haplo.glm( trait ~ genoH,
                  family="binomial", 
                  locus.label=snpsH,
                  allele.lev=attributes(genoH)$unique.alleles,
                  control = haplo.glm.control(haplo.freq.min=0.05))   
intervals(mod)
summary(mod)

Check whether a SNP is Monomorphic

Description

This function verifies when a SNP is Monomorphic

Usage

is.Monomorphic(x) 

Arguments

Value

A logical value TRUE if the SNP is Monomorphic, otherwise a FALSE

Examples

data(SNPs)
is.Monomorphic(SNPs$snp10001)
is.Monomorphic(SNPs$snp100020)
apply(SNPs[,20:30],2,is.Monomorphic)

Louis Information for haplo.glm

Description

For internal use within the haplo.stats library's haplo.glm function

Usage

louis.info(fit, epsilon=1e-8)

Arguments

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

Create a group of locus objects from some SNPs, assign to 'model.matrix' class.

Description

This function prepares the CRITICAL element corresponding to matrix of genotypes necessary to be included in 'haplo.glm' function.

Usage

make.geno(data, SNPs.sel)

Arguments

Value

the same as 'setupGeno' function, from 'haplo.stats' library, returns

See Also

snp

Examples

## Not run:
data(SNPs)
# first, we create an object of class 'setupSNP'
datSNP<-setupSNP(SNPs,6:40,sep="")
geno<-make.geno(datSNP,c("snp10001","snp10002","snp10003"))
## End(Not run)

max-statistic for a 2x3 table

Description

Computes the asymptotic p-value for max-statistic for a 2x3 table

Usage

maxstat(x, ...)

## Default S3 method:
maxstat(x, y, ...)

## S3 method for class 'table'
maxstat(x, ...)

## S3 method for class 'setupSNP'
maxstat(x, y, colSNPs=attr(x,"colSNPs"), ...)

## S3 method for class 'matrix'
maxstat(x, ...)

Arguments

Value

A matrix with the chi-square statistic for dominant, recessive, log-additive and max-statistic and its asymptotic p-value.

References

Gonzalez JR, Carrasco JL, Dudbridge F, Armengol L, Estivill X, Moreno V. Maximizing association statistics over genetic models (2007). Submitted

Sladek R, Rocheleau G, Rung J et al. A genome-wide association study identifies novel risk loci for type 2 diabetes (2007). Nature 445, 881-885

See Also

setupSNP

Examples

# example from Sladek et al. (2007) for the SNP rs1111875 
 tt<-matrix(c(77,298,310,122,316,231),nrow=2,ncol=3,byrow=TRUE)
 maxstat(tt)

 data(SNPs)
 maxstat(SNPs$casco,SNPs$snp10001) 
 myDat<-setupSNP(SNPs,6:40,sep="")
 maxstat(myDat,casco)

 

Remove rows with NA in covariates, but keep genotypes with NAs

Description

Removes rows with NA in response or covariates, but keeps subjects with NAs in their genotypes if not missing all alleles.

Usage

na.geno.keep(m)

Arguments

Value

a model matrix with rows removed if exclusion criteria requires it

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

Extract odds ratios, 95% CI and pvalues

Description

Extract odds ratios, 95

Usage

 odds(x, model=c("log-additive", "dominant", "recessive", "overdominant", "codominant"), 
      sorted=c("no","p-value","or"))

Arguments

Value

A matrix with OR 95% CI (lower, upper) and P value for the selected model. For codominant model, the OR and 95%CI are given for heterozygous and homozigous.

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

Examples

 data(SNPs)
 datSNP<-setupSNP(SNPs,6:40,sep="")
 ans<-WGassociation(casco~1,data=datSNP,model="all")
 odds(ans)

Permutation test analysis

Description

This function extract the p values for permutation approach performed using scanWGassociation function

Usage

   permTest(x, method="minimum", K)

Arguments

Details

This function extract the p values from an object of class 'WGassociation'. This object migth be obtained using the funcion called 'scanWGassociation' indicating the number of permutations in the argument 'nperm'.

Value

An object of class 'permTest'.

'print' returns a summary indicating the number of SNPs analyzed, the number of valid SNPs (those non-Monomorphic and that pass the calling rate), the p value after Bonferroni correction, and the p values based on permutation approach. One of them is based on considering the empirical percentil for the minimum p values, and the another one on assuming that the minimum p values follow a beta distribution.

'plot' produces a plot of the empirical distribution for the minimum p values (histogram) and the expected distribution assuming a beta distribution. The corrected p value is also showed in the plot.

See examples for further illustration about all previous issues.

References

Dudbridge F, Gusnanto A and Koeleman BPC. Detecting multiple associations in genome-wide studies. Human Genomics, 2006;2:310-317.

Dudbridge F and Koeleman BPC. Efficient computation of significance levels for multiple associations in large studies of correlated data, including genomewide association studies. Am J Hum Genet, 2004;75:424-435.

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

scanWGassociation

Examples

library(SNPassoc)

data(asthma, package = "SNPassoc")
asthma.s <- setupSNP(data=asthma, colSNPs=7:ncol(asthma), sep="")

ans <- WGassociation(casecontrol, data=asthma.s)

Function to plot -log p values from an object of class 'WGassociation'

Description

Function to plot -log p values from an object of class 'WGassociation'

Usage

## S3 method for class 'WGassociation'
plot(x, ...)

Arguments

Details

A panel with different plots (one for each mode of inheritance) are plotted. Each of them represents the -log(p value) for each SNP. Two horizontal lines are also plotted. One one them indicates the nominal statistical significance level whereas the other one indicates the statistical significance level after Bonferroni correction.

Value

No return value, just the plot

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

association setupSNP WGassociation

Examples

library(SNPassoc)

data(asthma, package = "SNPassoc")
asthma.s <- setupSNP(data=asthma, colSNPs=7:ncol(asthma), sep="")

ans <- WGassociation(casecontrol, data=asthma.s)

plot(ans)

Plot Haplotype Frequencies versus Haplotype Score Statistics

Description

Method function to plot a class of type haplo.score

Usage

## S3 method for class 'haplo.score'
plot(x, ...)

Arguments

Details

This is a plot method function used to plot haplotype frequencies on the x-axis and haplotype-specific scores on the y-axis. Because haplo.score is a class, the generic plot function can be used, which in turn calls this plot.haplo.score function.

Value

Nothing is returned.

References

Schaid DJ, Rowland CM, Tines DE, Jacobson RM, Poland GA. "Score tests for association of traits with haplotypes when linkage phase is ambiguous." Amer J Hum Genet. 70 (2002): 425-434.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.score

Examples

  data(hla.demo)
  geno <- as.matrix(hla.demo[,c(17,18,21:24)])
  keep <- !apply(is.na(geno) | geno==0, 1, any)
  hla.demo <- hla.demo[keep,]
  geno <- geno[keep,]
  attach(hla.demo)
  label <- c("DQB","DRB","B")
 
# For quantitative, normally distributed trait:

  score.gaus <- haplo.score(resp, geno, locus.label=label, 
                            trait.type = "gaussian")

  ## try: locator.haplo(1)

Plot of missing genotypes

Description

Plot a grid showing which genotypes are missing

Usage

plotMissing(x, print.labels.SNPs = TRUE, 
        main = "Genotype missing data", ...) 

Arguments

Details

This function uses 'image' function to plot a grid with black pixels where the genotypes are missing.

Value

No return value, just the plot

See Also

setupSNP

Examples

 data(SNPs)
 data(SNPs.info.pos) 
 ans<-setupSNP(SNPs,colSNPs=6:40,sep="")
 plotMissing(ans)
 
 # The same plot with the SNPs sorted by genomic position and 
 # showing the information about chromosomes

 ans<-setupSNP(SNPs,colSNPs=6:40,sort=TRUE,SNPs.info.pos,sep="") 
 plotMissing(ans)

Print a nice banner

Description

Print a centered banner that carries to multiple lines

Usage

printBanner(str, banner.width=options()$width, char.perline=.75*banner.width, border="=")

Arguments

Details

This function prints a nice banner in both R and S-PLUS

Value

nothing is returned

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

options

Examples

printBanner("This is a pretty banner", banner.width=40, char.perline=30)

# the output looks like this:
# ========================================
#         This is a pretty banner
# ========================================

Functions for inspecting population substructure

Description

This function plots ranked observed p values against the corresponding expected p values in -log scale.

Usage

qqpval(p, pch=16, col=4, ...)

Arguments

Value

No return value, just the plot

See Also

GenomicControl, WGassociation

Examples

data(SNPs)
datSNP<-setupSNP(SNPs,6:40,sep="")
res<-WGassociation(casco,datSNP,model=c("do","re","log-add"))

# observed vs expected p values for recessive model
qqpval(recessive(res))

Get related samples

Description

Get related samples

Usage

related(x)

Arguments

Value

A matrix with related individuals.

Examples

  library(SNPassoc)
  data(SNPs)

SNPs from HapMap project

Description

Information about 9307 SNPs from the HapMap project belonging to 22 chromosomes. Information about two different population is available: European population (CEU) and Yoruba (YRI). The genomic information (names of SNPs, chromosomes and genetic position) is also available in a data frame called 'HapMap.SNPs.pos'.

Usage

data(resHapMap)

Format

A data frame with 120 observations on the 9808 variables (SNPs) and one variable called 'group' indicating the population.

Source

HapMap project (http://www.hapmap.org)

Examples

    data(resHapMap)

Whole genome association analysis

Description

This function is obsolete due to some problems with gfotran compiler. Use 'WGassociation' function instead or send an e-mail to the maintainer for receiving a version including this function

Create the list of control parameters for simulations in haplo.score

Description

In the call to haplo.score, the sim.control parameter is a list of parameters that control the simulations. This list is created by this function, score.sim.control, making it easy to change the default values.

Usage

score.sim.control(p.threshold=0.25, min.sim=1000, max.sim=20000.,verbose=FALSE)

Arguments

Details

In simulations for haplo.score, employ the simulation p-value precision criteria of Besag and Clifford (1991). The criteria ensures both the global and the maximum score statistic simulated p-values be precise for small p-values. First, perform min.sim simulations to guarantee sufficient precision for the score statistics on individual haplotypes. Then continue simulations as needed until simulated p-values for both the global and max score statistics meet precision requirements set by p.threshold.

Value

A list of the control parameters:

References

Besag, J and Clifford, P. "Sequential Monte Carlo p-values." Biometrika. 78, no. 2 (1991): 301-304.

Note

This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.score

Examples

# it would be used in haplo.score as appears below
#
# score.sim.500 <- haplo.score(y, geno, trait.type="gaussian", simulate=T, 
#                sim.control=score.sim.control(min.sim=500, max.sim=2000)

Create a group of locus objects from a genotype matrix, assign to 'model.matrix' class.

Description

The function makes each pair of columns a locus object, which recodes alleles to numeric and saves the original alleles as an attribute of the model.matrix.

Usage

setupGeno(geno, miss.val=c(0,NA), locus.label=NULL)

Arguments

Details

This function contains the essential parts of the loci function, which is no longer within haplo.stats

Value

A 'model.matrix' object with the alleles recoded to numeric values, and the original values are stored in the 'unique.alleles' attribute. The ith item of the unique.alleles list is a vector of unique alleles for the ith locus.

Note

A matrix that contains all elements of mode character will be sorted in alphabetic order. This order may differ across platforms according to your setting of LC_COLLATE. See the note in haplo.em about how this sort order affects results. This function is a modified version of the original from the haplo.stats package (Copyright 2003 Mayo Foundation for Medical Education and Research). It has been included in SNPassoc to ensure continued functionality following the archival of haplo.stats from CRAN.

Author(s)

Original code by Daniel J. Schaid. Adapted for SNPassoc by Dolors Pelegrí-Sisó.

See Also

haplo.glm, haplo.em

Examples

# Create some loci to work with
a1 <- 1:6
a2 <- 7:12

b1 <- c("A","A","B","C","E","D")
b2 <-c("A","A","C","E","F","G")

c1 <- c("101","10","115","132","21","112")
c2 <- c("100","101","0","100","21","110")

myGeno <- data.frame(a1,a2,b1,b2,c1,c2)
myGeno <- setupGeno(myGeno)
myGeno

attributes(myGeno)$unique.alleles

Convert columns in a dataframe to class 'snp'

Description

setupSNP Convert columns in a dataframe to class 'snp'

summary.setupSNP gives a summary for an object of class 'setupSNP' including allele names, major allele frequencie, an exact thest of Hardy-Weinberg equilibrium and percentage of missing genotypes

Usage

setupSNP(data, colSNPs, sort = FALSE, info, sep = "/", ...)

Arguments

Value

a dataframe of class 'setupSNP' containing converted SNP variables. All other variables will be unchanged.

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

snp

Examples

 data(SNPs)
 myDat<-setupSNP(SNPs,6:40,sep="")


#sorted SNPs and having genomic information
 data(SNPs.info.pos)
 myDat.o<-setupSNP(SNPs,6:40,sep="",sort=TRUE, info=SNPs.info.pos)

# summary
 summary(myDat.o)

# plot one SNP
  plot(myDat,which=2)

SNP object

Description

snp creates an snp object

is returns TRUE if x is of class 'snp'

as attempts to coerce its argument into an object of class 'snp'

reorder change the reference genotype

summary gives a summary for an object of class 'snp' including genotype and allele frequencies and an exact thest of Hardy-Weinberg equilibrium

plot gives a summary for an object of class 'snp' including genotype and allele frequencies and an exact thest of Hardy-Weinberg equilibrium in a plot. Barplot or pie are allowed

[.snp is a copy of [.factor modified to preserve all attributes

Usage

  snp(x, sep = "/", name.genotypes, reorder="common", 
    remove.spaces = TRUE, allow.partial.missing = FALSE) 

  is.snp(x)
 
  as.snp(x, ...)

  ## S3 method for class 'snp'
additive(o)

Arguments

Details

SNP objects hold information on which gene or marker alleles were observed for different individuals. For each individual, two alleles are recorded.

The snp class considers the stored alleles to be unordered , i.e., "C/T" is equivalent to "T/C". It assumes that the order of the alleles is not important.

When snp is called, x is a character vector, and it is assumed that each element encodes both alleles. In this case, if sep is a character string, x is assumed to be coded as "Allele1<sep>Allele2". If sep is a numeric value, it is assumed that character locations 1:sep contain allele 1 and that remaining locations contain allele 2.

additive.snp recodes the SNPs for being analyzed as a linear covariate (codes 0,1,2)

Value

The snp class extends "factor" where the levels is a character vector of possible genotype values stored coded by paste( allele1, "", allele2, sep="/")

References

JR Gonzalez, L Armengol, X Sole, E Guino, JM Mercader, X Estivill, V Moreno. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007;23(5):654-5.

See Also

association

Examples

# some examples of snp data in different formats

dat1  <- c("21", "21", "11", "22", "21",
                    "22", "22", "11", "11", NA)
ans1  <- snp(dat1,sep="")
ans1

dat2 <- c("A/A","A/G","G/G","A/G","G/G",
                    "A/A","A/A","G/G",NA)
ans2  <- snp(dat2,sep="/")
ans2

dat3 <- c("C-C","C-T","C-C","T-T","C-C",
                    "C-C","C-C","C-C","T-T",NA)
ans3 <- snp(dat3,sep="-")
ans3


dat4 <- c("het","het","het","hom1","hom2",
                    "het","het","hom1","hom1",NA)
ans4 <- snp(dat4,name.genotypes=c("hom1","het","hom2"))
ans4


# summary 
summary(ans3)

# plots

plot(ans3)
plot(ans3,type=pie)
plot(ans3,type=pie,label="SNP 10045")

Sort a vector of SNPs by genomic position

Description

This function sorts a vector with the position of SNPs in a data frame using another data frame which contains information about SNPs, their chromosome, and their genomic position

Usage

sortSNPs(data, colSNPs, info)

Arguments

Details

First of all, the function obtains a vector with the SNPs sorted using the data frame with the genomic positions (see 'dataSNPs.pos' argument). Then, the columns which indicate where the information about the genotypes is in our data frame, are sorted using this vector.

This information is useful when WGassociation function is called since it allow the user to summaryze the results with the SNPs sorted by genomic position

Value

a vector indicating the colums where the SNPs are recorded in our data frame ('data' argument), sorted using the genomic positions listed in 'dataSNPs.pos' argument)

Examples

#
# data(SNPs)
# data(SNPs.info.pos)
# colSNPs.order<-sortSNPs(SNPs,c(6:40),SNPs.info.pos)
#

Test for Hardy-Weinberg Equilibrium

Description

Test the null hypothesis that Hardy-Weinberg equilibrium holds in cases, controls and both populations.

print print the information. Number of digits, and significance level can be changed

Usage

tableHWE(x, strata, ...)

Arguments

Details

This function calculates the HWE test for those variables of class 'snp' in the object x of class 'setupSNP'

Value

A matrix with p values for Hardy-Weinberg Equilibrium

Author(s)

This function is based on an R function which computes an exact SNP test of Hardy-Weinberg Equilibrium written by Wigginton JE, Cutler DJ and Abecasis GR available at http://csg.sph.umich.edu/abecasis/Exact/r_instruct.html

References

Wigginton JE, Cutler DJ and Abecasis GR (2005). A note on exact tests of Hardy-Weinberg equilibrium. Am J Hum Genet 76:887-93

See Also

setupSNP

Examples

data(SNPs)
ans<-setupSNP(SNPs,6:40,sep="")
res<-tableHWE(ans)
print(res)
#change the significance level showed in the flag column
print(res,sig=0.001)


#stratified analysis
res<-tableHWE(ans,ans$sex)
print(res)