Background Color:
 
Background Pattern:
Reset
Thursday, November 23, 2017
menu click
Search

Tutorial

We will use some linkage data for Rheumatoid Arthritis provided during the 15th Genetic Analysis Workshop. The data were collected by the NARAC (North American Rheumatoid Arthritis Consortium) and consists mostly of Affected Sib pairs but there are a number of larger families with both affected and unaffected individuals. Several covariates were also gathered but we will concentrate on the 'ever smoked' status variable (HxCigE). The R code below assumes that all files used are present in the R working directory. All files required to run this analysis can be downloaded with the R code at the bottom of this page.

One initial step before actually running LinkGLM.R itself is to compute the covariance matrix for IBD sharing for each family in the dataset. Quinta Helmer wrote a C++ program calling upon Merlin which makes those computations by Monte Carlo simulations. For the present data, the command to be typed at the DOS prompt to call upon this program is "linkage_score_v0_3 --covpihat --reruns 100 --start 0 --stop 200 --grid 5 -p merlin.ped -d merlin.dat -m merlin.map". This will generate a '.var' file for each family in the ped file and these files will be used as input by LinkGLM.R. The computations can take a while so the actual output files for this call are provided in the downloadable .zip file at the bottom of this page.

We then compile the LinkGLM.R code and auxiliary functions,

source("LinkGLM.R")

we first analyze the data without any covariates

LinkGLM.func(
model="logistic",   # for binary traits
trait="RA",         # trait label as in .dat file
trait.col=6,        # trait column as when read by the 'read.table' R command
beta=-6,            # this should reflect the prevalence of RA in the population
rho = 0.8,          # sib-sib correlation of random effects
sigma2 = 2.5,       # variance of random effects
start.pos = 0,
stop.pos = 200,
grid = 5 )

Note: a method used to derive sensible values for beta, rho and sigma2 is described in Hum Hered. 64:5-15 (2007) and a similar analysis is described in fuller details in Lebrec et al. Adjusting for sex and anti-CCP levels in linkage analysis of rheumatoid arthritis (To appear in BMC Genetics).

We can read in the output, summarize the results using the res.func function and plot the LOD scores

res <- read.table("LinkGLM.sc", header=T)
res.plot <- by( res, INDICES=res$POSITION, FUN=res.func)
x <- as.numeric(dimnames(res.plot)$"res$POSITION")
y <- unlist(res.plot)
y <- y[seq(1,length(y), by=2)]
plot(x, y, type="l", xlab="Position (cM)", ylab="LOD", col=1, ylim=c(-1,3))
title("NARAC Data Chr.06")

We now perform the analysis using HxCigE as a covariate:

LinkGLM.func(
model="logistic",
trait="RA",
trait.col=6,
path.sc = "LinkGLM_HxCigE.sc",
path.wgt = "LinkGLM_HxCigE.wgt",
path.err = "LinkGLM_HxCigE.err",
path.out= "LinkGLM_HxCigE.out",
beta=c(-6.7, +1.1 ),
cov.lst = c("HxCigE"),
cov.cols=c( 43 ),
cov.typ=c("categorical"),
rho = 0.8,
sigma2 = 2.5,
start.pos = 0,
stop.pos = 200,
grid = 5 )

We read in the ouput and display the results on the same plot as the initial analysis without covariate:

res <- read.table("LinkGLM_HxCigE.sc", header=T)
res.plot <- by( res, INDICES=res$POSITION, FUN=res.func)
x <- as.numeric(dimnames(res.plot)$"res$POSITION")
y <- unlist(res.plot)
y <- y[seq(1,length(y), by=2)]
lines(x,y, lty=2)