fbatc

Description

Family based test for a group of markers conditional on another group of markers (typically conditional on a single marker). To start the graphical interface, provide no options, i.e. type fbatc() and press return.

Usage

fbatc( ped=NULL, phe=NULL, data=mergePhePed( ped, phe),
       trait="AffectionStatus", traitType="auto",
       markerAnalyze=NULL, markerCondition=NULL,
       offset=NULL,
       tempPrefix="temp",
       MAXITER=1000, TOL=sqrt(.Machine$double.eps),
       verbose=FALSE )

Arguments

Details

Implements the test as described in Hoffmann et. al. (please see References).

The results returned are a data.frame object. The column ‘pvalue’ and ‘rank’ are the pvalue and rank of the empirical covariance matrix of the model-based test (dichotomous or normal). The column ‘pvalueR’ and ‘rankR’ are the pvalue and rank of the robust test. The model-based test has considerable more power over the robust test, but must assume a disease model. Please see Hoffmann et. al. for more details.

This requires that FBAT be installed. If it is not, then the routine will attempt to automatically install it when given permission to do so by the user.

References

Hoffmann, Thomas J. and Laird, Nan M. Parsing the Effects of Individual SNPs in Candidate Genes in Families. Submitted.

Examples

## Not run: 
  set.seed(13)

  ## We simulate NO LD HERE, and a completely random trait!
  ## Data here is only to show how to use the function

  ###################
  ## IGNORE START: ##
  ###################

  ##  You can safely ignore how the data is generated,
  ##  and just see how to use it afterward.
  NUM_FAMILIES <- 500
  AFREQ <- c(0.2,0.2)
  ped <- as.ped( data.frame( pid = kronecker(1:NUM_FAMILIES,c(1,1,1)),
                            id = kronecker( rep(1,NUM_FAMILIES), c(1,2,3) ),
                            idfath = kronecker( rep(1,NUM_FAMILIES), c(0,0,1) ),
                            idmoth = kronecker( rep(1,NUM_FAMILIES), c(0,0,2) ),
                            sex = rep(0,NUM_FAMILIES*3),
                            AffectionStatus = kronecker( rep(1,NUM_FAMILIES), c(0,0,2) ),
                            m0.a = rep(0,NUM_FAMILIES*3),      ## missing for now
                            m0.b = rep(0,NUM_FAMILIES*3),
                            m1.a = rep(0,NUM_FAMILIES*3),
                            m1.b = rep(0,NUM_FAMILIES*3) ) )
  CUR_FAMILY <- 1
  while( CUR_FAMILY<=NUM_FAMILIES ) {
    ## Indexing: start=father, (start+1)=mother, (start+2)=child
    start <- CUR_FAMILY*3-2

    ## Draw the parental genotypes from the population
    ped$m0.a[start:(start+1)] <- rbinom( 1, 1, AFREQ[1] ) + 1
    ped$m0.b[start:(start+1)] <- rbinom( 1, 1, AFREQ[1] ) + 1
    ped$m1.a[start:(start+1)] <- rbinom( 1, 1, AFREQ[2] ) + 1
    ped$m1.b[start:(start+1)] <- rbinom( 1, 1, AFREQ[2] ) + 1

    ## Draw the children's genotype from the parents
    ma <- rbinom( 1, 1, 0.5 )
    mb <- rbinom( 1, 1, 0.5 )
    if( rbinom( 1, 1, 0.5 ) == 0 ) {
      ped$m0.a[start+2] <- ped$m0.a[start]
      ped$m1.a[start+2] <- ped$m1.a[start]
    }else{
      ped$m0.a[start+2] <- ped$m0.b[start]
      ped$m1.a[start+2] <- ped$m1.b[start]
    }
    if( rbinom( 1, 1, 0.5 ) == 0 ) {
      ped$m0.b[start+2] <- ped$m0.a[start+1]
      ped$m1.b[start+2] <- ped$m1.a[start+1]
    }else{
      ped$m0.b[start+2] <- ped$m0.b[start+1]
      ped$m1.b[start+2] <- ped$m1.b[start+1]
    }

    CUR_FAMILY <- CUR_FAMILY + 1
  }

  ## Create a completely random phenotype as well
  phe <- as.phe( data.frame( pid=ped$pid, id=ped$id, qtl=rnorm(nrow(ped)) ) )

  ################
  ## IGNORE END ##
  ################

  ## Look at the first part of the pedigree
  print( head( ped ) )
  ## Look at the first part of the phenotype
  print( head( phe ) )

  ## Binary trait
  ## -- fbatc default trait is AffectionStatus
  ## -- fbatc default trait type is 'auto'
  ## - Test marker m1 conditional on m0
  print(  fbatc( ped=ped, markerAnalyze="m1", markerCondition="m0" )  )
  ## - Do the test the other way around, m0 conditional on m1
  print(  fbatc( ped=ped, markerAnalyze="m0", markerCondition="m1" )  )
  ## - Otherwise, we could have done this in one step;
  ##    if markerCondition is not specified, each member
  ##     of markerAnalyze is used.
  print(  fbatc( ped=ped, markerAnalyze=c("m0","m1") )  )

  ## QTL
  print(  fbatc( ped=ped, phe=phe, trait="qtl", markerAnalyze="m1", markerCondition="m0" )  )
  print(  fbatc( ped=ped, phe=phe, trait="qtl", markerAnalyze="m0", markerCondition="m1" )  )

  ## Additionally, we could write out the data that we
  ##  generated to disk so that we can then use it.
  write.ped( "simulated", ped ) ## to simulated.ped
  write.phe( "simulated", phe ) ## to simulated.phe

## End(Not run)

FBAT-C Stepwise Strategy

Description

Apply the FBAT-C test in a stepwise fashion using fbatcStrategyStep (which does forward selection with fbatcStrategyStepUp, followed by backwards selection with fbatcStrategyStepDown) and get the results ready for publication with fbatcStrategyStepLatex.

Usage

fbatcStrategyStepUp(ped, phe, markers=pedMarkerNames(ped), trait="trait",
  traitType="auto", alphaMMarker=0.05, alphaStep=alphaMMarker, sortByCorrelation=TRUE,
  tempPrefix="temp_", sim=FALSE, debug=FALSE )
fbatcStrategyStepDown(ped, phe,  markers=pedMarkerNames(ped),
  markersChosen=pedMarkerNames(ped), markersChosenR=markersChosen, trait="trait",
  traitType="auto", alphaMMarker=0.05, alphaStep=alphaMMarker, sortByCorrelation=TRUE,
  tempPrefix="temp_", sim=FALSE, debug=FALSE )
fbatcStrategyStep(ped, phe, markers=pedMarkerNames(ped), trait="trait",
  traitType="auto", alphaMMarker=0.05, alphaStep=alphaMMarker, sortByCorrelation=TRUE,
  tempPrefix="temp_", sim=FALSE, debug=FALSE )
fbatcStrategyStepLatex(res, digits=4, ffile="", preamble=FALSE, build=preamble, pdf="")
## S3 method for class 'fbatcSStep'
print(x,...)

Arguments

Details

fbatcStrategy returns a list with the following components.

mmarkerPvalue: p-value of the multi-marker test on those markers (Rakovski et. al 2008).

correlation: correlation matrix of the markers

univariate: univariate results

step: (model-based test) list of components pvalue (ith pvalue of the conditional test of markersAnalyze[i] on all markersCondition), numInf (number of informative families in the ith test), markersAnalyze, and markersCondition

markersChosen: (model-based test) results from applying step-up strategy

stepR, markersChosenR: (model-free test) results similar to step and markersChosen.

fbatge

Description

Family based test for gene-environment interaction utilizing arbitrary family structures and multiple affected offspring. This method is recommended over the fbati routine in most scenarios.

If no arguments are passed, then a friendly graphical interface is presented to the user.

fbatge [GxE test], fbatj (see fbatj help) [G,GxE test], fbatme (see fbatme help) [G test] generally have more options than fbatgeAll. fbatgeAll runs all three tests, and gives results of all of them, and so uses only the options that are common to all three functions.

Usage

fbatge( ped=NULL, phe=NULL,
        env=NULL, cov=NULL,
        trait="AffectionStatus", geno=NULL,
        strategy="hybrid", model="additive" )
fbatgeAll( ped=NULL, phe=NULL, env=NULL, trait="AffectionStatus" )

Arguments

Details

Implements the test as described in Hoffmann et. al. (please see References).

NOTE: The allele frequency is simply based on the allele frequency in all genotyped individuals, and is not the best choice.

References

Hoffmann, Thomas J., and Laird, Nan M. Combining Multiple Disease Models for a More Powerful Gene-Environment Interaction Test in Nuclear Families.

Examples

  ## Not run: 
    example( fbati ) ## See fbati, creates a dataset for us in 'phe' and 'ped'
    print( fbatge( ped=ped, phe=phe, env="env" ) )
    ## The results are very close to the FBAT-I function, which
    ##  we would expect for trios.
  
## End(Not run)

fbati

Description

Family based test for gene-environment interaction for bi-allelic snps, command/line or GUI (provide no options to start the graphical interface, i.e. just type fbati() and press return).

Usage

fbati( ped=NULL, phe=NULL,
       data=mergePhePed(ped,phe),
       marker=NULL, ## pairs...
       env,
       method="fbati",
       model="additive",
       iter=10000,
       seed=7,
       maxSib=3,
       fixNames=TRUE,
       debug=FALSE )

Arguments

Details

Returns a data.frame object with the results. The columns entitled GX...X indicate the number of informative families in each strata for the given marker. If these columns do not show up, it indicates there was only one type of strata.

The parents need not be in the dataset if they have completely missing genotypes (they will be inserted), but the snps must currently be bi-allelic (or you will get error messages).

fread.ped and fread.phe are suggested to enforce loading the whole dataset.

NOTE 1: The fbati test was developed for families with at least one affected, so if there is more than one affected individual per family, only a random affected one will be used, and a random unaffected to reduce strata, unless strataFix is disabled. This is done on a per marker basis, thus the seed is set before every marker to obtain reproducible results.

NOTE 2: The data is converted into nuclear families. This is done by a call to ‘nuclifyMerged’ to the passed in dataset to enforce this consistency.

References

Hoffmann, Thomas J., Lange, Christoph, Vansteelandt, Stijn, and Laird, Nan M. Gene-Environment Interaction Test for Dichotomous Traits in Trios and Sibships. Submitted.

S. L. Lake and N. M. Laird. Tests of gene-environment interaction for case-parent triads with general environmental exposures. Ann Hum Genet, 68(Pt 1):55-64, Jan 2004.

Examples

## Not run: 
## Data is simulated according to the formula in the
##  paper (you can see it from the code).

## Set the seed so you get the same results
set.seed(13)

## Constants (you can vary these)
NUM_FAMILIES <- 500
AFREQ <- 0.1  ## Allele frequency
BG <- -0.25   ## main effect of gene
BE <- 0       ## main effect of environment
BGE <- 0.75   ## main gene-environment effect
ENV <- 0.2    ## environmental exposure frequency

## (but don't modify this one)
MAX_PROB <- exp( BG*2 + BE*1 + BGE*2*1 )

#####################################
## Create a random dataset (trios) ##
#####################################

## -- genotypes are set to missing for now,
##    everyone will be affected
ped <- as.ped( data.frame( pid = kronecker(1:NUM_FAMILIES,c(1,1,1)),
                           id = kronecker( rep(1,NUM_FAMILIES), c(1,2,3) ),
                           idfath = kronecker( rep(1,NUM_FAMILIES), c(0,0,1) ),
                           idmoth = kronecker( rep(1,NUM_FAMILIES), c(0,0,2) ),
                           sex = rep(0,NUM_FAMILIES*3),
                           AffectionStatus = kronecker( rep(1,NUM_FAMILIES), c(0,0,2) ),
                           m0.a = rep(0,NUM_FAMILIES*3),      ## missing for now
                           m0.b = rep(0,NUM_FAMILIES*3) ) )   ## missing for now
## -- envioronment not generated yet
phe <- as.phe( data.frame( pid = ped$pid,
                           id = ped$id,
                           env = rep(NA,NUM_FAMILIES*3) ) )   ## missing for now

## 50/50 chance of each parents alleles
mendelTransmission <- function( a, b ) {
  r <- rbinom( length(a), 1, 0.75 )
  return( a*r + b*(1-r) )
}

## Not the most efficient code, but it gets it done;
##  takes < 5 sec on pentium M 1.8Ghz
CUR_FAMILY <- 1
while( CUR_FAMILY<=NUM_FAMILIES ) {
  ## Indexing: start=father, (start+1)=mother, (start+2)=child
  start <- CUR_FAMILY*3-2

  ## Draw the parental genotypes from the population
  ped$m0.a[start:(start+1)] <- rbinom( 1, 1, AFREQ ) + 1
  ped$m0.b[start:(start+1)] <- rbinom( 1, 1, AFREQ ) + 1

  ## Draw the children's genotype from the parents
  ped$m0.a[start+2] <- mendelTransmission( ped$m0.a[start], ped$m0.b[start] )
  ped$m0.b[start+2] <- mendelTransmission( ped$m0.a[start+1], ped$m0.b[start+1] )

  ## Generate the environment
  phe$env[start+2] <- rbinom( 1, 1, ENV )

  ## and the affection status
  Xg <- as.integer(ped$m0.a[start+2]==2) + as.integer(ped$m0.b[start+2]==2)
  if( rbinom( 1, 1, exp( BG*Xg + BE*phe$env[start+2] + BGE*Xg*phe$env[start+2] ) / MAX_PROB ) == 1 )
    CUR_FAMILY <- CUR_FAMILY + 1
  ## otherwise it wasn't a valid drawn individual
}


##############
## Analysis ##
##############

## Print the first 4 families
print( head( ped, n=12 ) )
print( head( phe, n=12 ) )

## NOTE: We could have just put all of this info into a single dataframe otherwise,
##  that would look like just the results of this
data <- mergePhePed( ped, phe )
print( head( data, n=12 ) )

## And run the analysis on all the markers
fbati( ped=ped, phe=phe, env="env" )

## Or do it via the merged data.frame object
##  7 and 8 correspond to the marker columns
fbati( data=data, env="env", marker=c(7,8) )

## You may also want to up the number of Monte-Carlo
##  iterations from the default

## And we could also run a joint test instead
##  (see fbatj)
## fbatj temporarily removed from namespace
#fbatj( ped=ped, phe=phe, env="env" )
#fbatj( data=data, env="env", marker=c(7,8) )

## This won't be run, but we could do this with the gui.
## It requires the data to be written to disk, so we do so:
write.ped( ped, "simulated" )
write.phe( phe, "simulated" )
## Then start the GUI -- specify the options as before,
##  but for the first two, navigate to the 'simulated.ped' and 'simulated.phe' files.
fbati()

## End(Not run)

FBAT Main effects Test

Description

Family based test for the main genetic effect, using the variance based on Mendelian transmissions. The null hypothesis is that there is no linkage and no association.

Usage

fbatme( ped=NULL, phe=NULL,
        data=mergePhePed(ped,phe),
        marker=NULL,
        trait="AffectionStatus",
        model="additive",
        fixNames=TRUE,
        verbose = FALSE )

Arguments

Details

Returns a data.frame object with the results. Uses the variance based on Mendelian transmissions.

NOTE: The allele frequency is simply based on the allele frequency in all genotyped individuals, and is not the best choice.

Launchpad

Description

Provides a GUI launchpad for routines in the fbati (i.e. this package) and pbatR (a dependency of this package) R packages.

Usage

launchpad()

See Also

fbati

Nuclify and Merge

Description

mergePhePed merges a phenotype and pedigree object into a single data.frame object.

nuclifyMerged chops a merged object into nuclear families of a dataset, generally a necessary preprocessing option for tests.

nuclify chops instead a ‘ped’ and ‘phe’ object separately.

Usage

mergePhePed(ped, phe)
nuclifyMerged(data, OUT_MULT=2)
nuclify(ped, phe)

Arguments

Details

mergePhePed and nuclifyMerged both return data.frame objects. nuclify returns a list that contains the ‘phe’ object and the ‘ped’ object with those respective names (see pbatR documentation, both objects extend data.frame objects, and can be used for the most part as if data.frame objects). When the data is nuclified, the parents of the nuclified families parents are lost.

NOTE: nuclifyMerged will modify the pedigree id (pid) to be [100*(previous pid) + (nuclear family index)]. This should make it easy to observe the results of this call to your dataset.

Examples

## Create some pedigree structure
##
##   100 --- 101
##        |
##       201---202
##           |
##     -------------
##     |   |   |   |
##    301 302 303 304
ped <- as.ped( data.frame( pid    = rep(1,8),
                           id     = c(100,101, 201,202, 301,302,303,304),
                           idfath = c(0,0,     100,0,   201,201,201,201),
                           idmoth = c(0,0,     101,0,   202,202,202,202),
                           sex    = c(1,2,     1,2,     2,2,2,2),
                           AffectionStatus = rep(0,8),
                           m0.a = rep(2,8),
                           m0.b = rep(2,8) ) )

## Which should chop up into
##   100 --- 101    201---202
##        |             |
##       201      -------------
##                |   |   |   |
##               301 302 303 304
nuclifyMerged( ped )

## NOTE: We could have merged the ped with a phe object,
##        via the 'mergePhePed' routine before running.

Strata Reduction

Description

Reduces the number of strata in nuclear pedigrees for testing (for use with FBAT-I). For nuclear families with both parents, a random affected child is drawn. For nuclear families with at least one parent missing, a random affected and another random sib is used (parents ignored).

Usage

strataReduce( data, envCol, m0, m1=m0+1, maxSib=3 )

Arguments

Examples

## Function creates a family with the specified markers and statuses
## NOTE: affection is false/true, whereas it is coded 1/2 in the ped file
createFam <- function( pa=c(0,0), pb=c(0,0),
                       ca, cb,
                       caffected=rep(TRUE,length(ca)),
                       env=1:length(ca) ) {
  ## pid, id, idfath, idmoth, sex, affection, m0a, m0b
  numC <- length(ca)
  return( data.frame( pid=rep(1,2+numC),
                      id=1:(2+numC),
                      idfath=c(0,0,rep(1,numC)),
                      idmoth=c(0,0,rep(2,numC)),
                      sex=c(2,1,rep(0,numC)),
                      affection=c(0,0,as.integer(caffected)+1),
                      m0.a=c(pa,ca), m0.b=c(pb,cb),
                      env=c(NA,NA,env) ) )
}
## Function tests/exemplifies the strataReduce(...) routine
srFam <- function( ... ) {
  data <- createFam( ... )
  data2 <- strataReduce( data=data, envCol=9, m0=7, maxSib=2 )
  cat( "Original data:\n" )
  print( data )
  cat( "Reduced stratification data:\n" )
  print( data2 )
}

## Basic sib test
srFam( ca=c(1,1,2), cb=c(1,2,2) )

## Basic trio test
srFam( ca=c(1,1,2), cb=c(1,2,2), pa=c(1,1), pb=c(2,2) )

## a fairly comprehensive test here
## The affected should always be one of the first three,
##  the unaffected could be one the first eight
for( i in 1:10 )
  srFam( ca=c(1:8,0,0), cb=c(1:8,0,0),
         pa=c(1,1),
         caffected=c(rep(TRUE,6),rep(FALSE,4)),
         env=c(1:3,rep(NA,7)) )

## Now just to make sure, a full pedigree, rather than just one family
data <- createFam( ca=1:2, cb=1:2 )
for( i in 2:10 )
  data <- rbind( data, createFam( ca=1:2, cb=1:2 ) )
  cat( "Original data (full pedigree):\n" )
  print( data )
  cat( "Reduced stratification data (full pedigree), maxSib=3\n" )
  print( strataReduce( data=data, envCol=9, m0=7 ) )