| LD.decay | R Documentation |
This function calculates the LD decay based on a marker matrix and a map with distances between markers in cM or base pairs.
LD.decay(markers,map,silent=FALSE,unlinked=FALSE,gamma=0.95)
markers |
a numeric matrix of markers (columns) by individuals (rows) in -1, 0, 1 format. |
map |
a data frame with 3 columns "Locus" (name of markers), "LG" (linkage group or chromosome), and "Position" (in cM or base pairs). |
silent |
a TRUE/FALSE value statement indicating if the program should or should not display the progress bar. silent=TRUE means that will not be displayed. |
unlinked |
a TRUE/FALSE value statement indicating if the program should or should not calculate the alpha(see next argument) percentile of interchromosomal LD. |
gamma |
a percentile value for LD breakage to be used in the calculation of interchromosomal LD extent. |
a list with 3 elements; "by.LG", "all.LG", "LDM". The first element (by.LG) is a list with as many elements as chromosomes where each contains a matrix with 3 columns, the distance, the r2 value, and the p-value associated to the chi-square test for disequilibrium. The second element (all.LG) has a big matrix with distance, r2 values and p-values, for each point from all chromosomes in a single data.frame. The third element (LDM) is the matrix of linkage disequilibrium between pairs of markers.
If unlinked is selected the program should return the gamma percentile interchromosomal LD (r2) for each chromosome and average.
Laido, Giovanni, et al. Linkage disequilibrium and genome-wide association mapping in tetraploid wheat (Triticum turgidum L.). PloS one 9.4 (2014): e95211.
nInds=300
nMarks=500
#random population of nInds lines with nMarks markers
M <- matrix(rep(0,nInds*nMarks),nInds,nMarks)
for (i in 1:nInds) {
M[i,] <- ifelse(runif(nMarks)<0.5,-1,1)
}
rownames(M) <- paste0("ind",1:nInds)
colnames(M) <- paste0("mark",1:nMarks)
mapM = data.frame(Locus=colnames(M),Position=1:nMarks,LG=1)
res <- LD.decay(M, mapM)
#### subset only markers with significant LD
res$all.LG <- res$all.LG[which(res$all.LG$p < .001),]
#### plot the LD decay
with(res$all.LG, plot(r2~d,col=transp("cadetblue"),
xlim=c(0,55), ylim=c(0,1),
pch=20,cex=0.5,yaxt="n",
xaxt="n",ylab=expression(r^2),
xlab="Distance in cM")
)
axis(1, at=seq(0,55,5), labels=seq(0,55,5))
axis(2,at=seq(0,1,.1), labels=seq(0,1,.1), las=1)
#### if you want to add the loess regression lines
#### this could take a long time!!!
# mod <- loess(r2 ~ d, data=res$all.LG)
# par(new=TRUE)
# lilo <- predict(mod,data.frame(d=1:55))
# plot(lilo, bty="n", xaxt="n", yaxt="n", col="green",
# type="l", ylim=c(0,1),ylab="",xlab="",lwd=2)
# res3 <- LD.decay(markers=CPgeno, map=mapCP,
# unlinked = TRUE,gamma = .95)
# abline(h=res3$all.LG, col="red")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.