Nothing
R/plot.pedigree.R
In pbatR: Pedigree/Family-Based Genetic Association Tests Analysis and Power
Defines functions
plotPed
plot.pedigree
Documented in
plotPed
## Taken from the kinship package (plot.pedigree.R, pedigree.R, kinship.R), licensed under (GPL >= 2), and modified as below to 1) display more information, 2) sensible margins so you can see the output, 3) handle AffectionStatus without crashing, and just generally play nicer with our data.
## FROM pedigree.R in kinship package
## From plot.pedigree.R in kinship
# $Id: plot.pedigree.s,v 1.10 2003/10/24 14:13:34 Atkinson Exp $
plot.pedigree <- function(x, id = x$id, sex = x$sex, status = x$status,
affected = x$affected,
cex = 1, col = rep(1, length(x$id)),
symbolsize = 1, branch = 0.6,
packed = T, align = packed, width = 8,
density=c(-1, 50,70,90), mar=c(4.1, 1, 4.1, 1),
angle=c(90,70,50,0), keep.par=F,
extraInfo = NULL, ## Tom added
...)
{
#library(kinship2)
#col <- rep("white",length(x$id))
#col[affected==1] <- "black"
col[is.na(affected)] <- "red" ## missing!!
maxlev <- max(x$depth) + 1
n <- length(x$depth)
## DATA CHECKS
sex <- as.numeric(sex)
sex[is.na(sex)] <- 3
if(any(sex < 1 | sex > 4))
stop("Invalid sex code")
if(length(sex) != n)
stop("Wrong length for sex")
if(is.null(status))
status <- rep(0, n)
else {
if(!all(status == 0 | status == 1 | status == 2))
stop("Invalid status code")
if(length(status) != n)
stop("Wrong length for status")
}
if(!is.null(id)) {
if(length(id) != n)
stop("Wrong length for id")
}
## Define plotting symbols for sex by affected
## in S, 0=square, 1= circle, 5=diamond, 2=triangle
## 15, 16, 18 and 17 are filled versions of them
symbol <- c(0, 1, 5, 2)[sex]
## CREATE SYMBOLS/PLOTS
points.sym <- function(x, y, symbol, status, affected, size, col, aspect,
angle, density, adj)
{
circle <- function(cx, cy, r, code=0, col=1, angle, density, ...)
{
# cx, cy, coordinates for center; r is radius
z <- (0:360 * pi)/180
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
if(sum(code,na.rm=T)==0) polygon(x, y, border=T,
density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1],...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
z <- (0:180 * pi)/180
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
z <- (180:360 * pi)/180
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
z <- (0:90 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3],...)
z <- (180:270 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
z <- (270:360 * pi)/180
x <- c(cx, sin(z) * r + cx)
y <- c(cy, cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
z <- (0:90 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
z <- (90:180 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4],...)
z <- (180:270 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
z <- (270:360 * pi)/180
x <- c(cx, sin(z) * r + cx)
y <- c(cy, cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
square <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- cx + c(-r,-r,r,r)
y <- cy + (1/adj)*c(-r,r,r,-r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0, col=col,...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0,0,r,r)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle[3], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
diamond <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- cx + c(-r,0,r,0)
y <- cy + (1/adj)*c(0,r,0,-r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(r,-r,0)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle=angle[3], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(-r,0,0)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
triangle <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
r <- r*1.25
a <- 3*r/sqrt(3)
b <- r/2
x <- cx + c((-1/2)*a, (1/2)*a, 0)
y <- cy + (1/adj)*c(-b, -b, r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1],angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c((-1/2)*a, 0, 0)
y <- cy + (1/adj)*c(-b, -b, r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0, (1/2)*a, 0)
y <- cy + (1/adj)*c(-b, -b, r)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
midx <- (r*(.5)*a)/(b+r)
x <- cx + c((-1/2)*a,-midx, 0, 0)
y <- cy + (1/adj)*c(-b, 0, 0, -b)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-midx, 0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,midx)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle=angle[3],...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
midx <- (r*(.5)*a)/(b+r)
x <- cx + c((-1/2)*a, -midx , 0, 0)
y <- cy + (1/adj)*c(-b, 0, 0, -b)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-midx, 0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,midx)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,midx,(1/2)*a,0)
y <- cy + (1/adj)*c(0,0,-b,-b)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
# in S, 0=square, 1= circle, 5=diamond, 2=triangle
# 15, 16, 18 and 17 are filled versions of them
for(i in (1:length(x))[!is.na(match(symbol,c(0,15)))]){
square(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i], angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(1,16)))]){
circle(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(5,18)))]){
diamond(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(2,17)))]){
triangle(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
# Draw slash if status=1
who <- (status == 1)
if(any(who)) {
deltax <- size*cos(pi/4) + aspect
deltay <- (1/adj)*(size*sin(pi/4) + aspect)
segments(x[who] - deltax, y[who] - deltay,
x[who] + deltax, y[who] + deltay)
}
}
#
# Now, get the structure of the plot, and lay out the main region
#
plist <- kinship2::align.pedigree(x, packed = packed, width = width, align =
align)
who <- (plist$nid > 0)
xp <- plist$pos[who]
yp <- - (row(plist$nid))[who]
np <- plist$nid[who]
oldpar <- par(mar = mar, err=-1)
par(usr = c(range(xp) + c(-0.5, 0.5), - (maxlev + 0.5), -0.5))
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
symbolsize <- symbolsize * adj
radius <- 0.08 * symbolsize
delta <- .1
aspect <- 0.05 * symbolsize #used for status variable
# textoff <- - (1/adj)*(radius) - .5*cex*(par()$"1em"[2])
textheight <- cex*par()$cin[2] #/4 ## ? div by 7 ??
textoff <- - (1/adj)*(radius) - .5*cex*(par()$"cin"[2]/7) #- textheight*0.1
totalTextHeight <- textheight
#cat("totalTextHeight before:", totalTextHeight, "\n")
#print(extraInfo)
#print(dim(extraInfo))
if(!is.null(extraInfo))
totalTextHeight <- textheight * (1 + ncol(extraInfo))
#cat("totalTextHeight after:", totalTextHeight, "\n")
## Tom addition to fix the margins
#print( xp )
#print( yp )
#print(totalTextHeight)
xmin <- min(xp)
xmax <- max(xp)
ymin <- min(yp)
ymax <- max(yp)
xlen <- xmax - xmin
ylen <- ymax - ymin
fudge <- max(0.1 * c(xlen, ylen)) + totalTextHeight
xmin <- xmin - fudge
xmax <- xmax + fudge
ymin <- ymin - fudge
ymax <- ymax + fudge
#cat("%%", xmin, xmax, ymin, ymax, "\n")
#plot(xp,yp,axes=F,type='n',xlab='',ylab='',...)
plot(c(xmin,xmax), c(ymin,ymax), axes=F,type='n',xlab='',ylab='',...)
points.sym(xp, yp, symbol[np], status[np], affected[np, ,drop=F],
radius, col[np], aspect, angle=angle, density=density, adj=adj)
#print(id)
#print(extraInfo)
#stop()
for(p in 1:length(np)){
textOut <- c()
if(!is.null(id))
textOut <- c(textOut, paste("ID:", id[np[p]]))
if(!is.null(extraInfo))
for(i in 1:ncol(extraInfo))
textOut <- c(textOut, paste(colnames(extraInfo)[i], ": ", extraInfo[np[p],i], sep=""))
if(length(textOut) > 0)
text(xp[p], yp[p] + textoff, paste(textOut, collapse="\n"), cex=cex, pos=1)
}
if(!is.null(id)) {
#if(is.null(extraInfo))
# text(xp, yp + textoff, paste("ID:",id[np]), cex = cex, col = col[np])
#else
# text(xp, yp + textoff, paste(paste("ID:",id[np]), extraInfo[np,1:ncol(extraInfo)], collapse="\n"), cex=cex, col=col[np])
###text(xp, yp + textoff, id[np], cex = cex, col = col[np])
### Tom's additions, so we get some more information here...
#text(xp, yp + textoff, paste("ID:",id[np]), cex = cex, col = col[np])
#
##text(xp, yp + textoff - textheight, "Bite me", cex=cex )
if(!is.null(extraInfo)){
##text(xp, yp + textoff, paste(extraInfo[np,],collapse="\n"))
#for(i in 1:ncol(extraInfo))
# text(xp, yp + textoff - textheight*i, paste(colnames(extraInfo)[i], ": ", extraInfo[np,i], sep=""))
}
}
# Draw in the family connections
for(i in 1:maxlev) {
# between spouses
if(any(plist$spouse[i, ]>0)) {
temp <- (1:ncol(plist$spouse))[plist$spouse[i, ]>0]
segments(plist$pos[i, temp] + radius, rep( - i, length(temp)),
plist$pos[i, temp + 1] - radius, rep(-i, length(temp)))
# to be added -- double segment for plist$spouse==2
}
}
for(i in 2:maxlev) {
# Children to parents, and across sibs
zed <- unique(plist$fam[i, ])
zed <- zed[zed > 0]
for(fam in zed) {
xx <- plist$pos[i - 1, fam + 0:1]
parentx <- mean(xx)
# Draw uplines from each subject
if(!is.null(plist$twins)){
tw.left <- (plist$twins[i, ] > 0 & plist$twins[i, ] < 4) & plist$fam[i,]==fam
mz.left <- (plist$twins[i, ] == 1) & plist$fam[i,]==fam
un.left <- (plist$twins[i, ] == 3) & plist$fam[i,]==fam
## figure out twins IN FAMILY (find next person in family)
tw.right <- rep(F,length(tw.left))
famlst <- plist$fam[i,]
twloc <- (1:length(tw.right))[tw.left]
famloc <- (1:length(tw.right))[famlst==fam]
flag <- NULL
for(j in twloc) flag <- c(flag,min(famloc[j<famloc]))
tw.right[flag] <- T
mz.right <- rep(F,length(mz.left))
mzloc <- (1:length(mz.right))[mz.left]
flag <- NULL
for(j in mzloc) flag <- c(flag,min(famloc[j<famloc]))
mz.right[flag] <- T
## unknown zygosity of the twin (type=3)
un.right <- rep(F,length(un.left))
unloc <- (1:length(un.right))[un.left]
flag <- NULL
for(j in unloc) flag <- c(flag,min(famloc[j<famloc]))
un.right[flag] <- T
}
if(is.null(plist$twins)){
twn <- length(plist$fam[i,]==fam)
tw.left <- rep(F,twn)
tw.right <- rep(F,twn)
mz.left <- rep(F,twn)
mz.right <- rep(F,twn)
un.left <- rep(F,twn)
un.right <- rep(F,twn)
}
tw <- tw.left|tw.right
##mz <- mz.left|mz.right ## Codetools 02/20/2014
##un <- un.left|un.right ## Codetools 02/20/2014
who <- (plist$fam[i, ] == fam) & !tw
xx <- plist$pos[i, who]
yy <- rep( - i, length = sum(who))
ww <- plist$nid[i, who]
xx.l <- plist$pos[i, tw.left]
yy.l <- rep( - i, length = sum(tw.left))
ww.l <- plist$nid[i, tw.left]
xx.r <- plist$pos[i, tw.right]
yy.r <- rep( - i, length = sum(tw.right))
ww.r <- plist$nid[i, tw.right]
segments(xx, yy + (1/adj) * radius, xx, yy + 3 * delta, col
= col[ww])
#####################################################
# draw diag line for twins
## left side of twin
segments(xx.l, yy.l + (1/adj) * radius, .5*(xx.l+xx.r), yy.l + 3 * delta,
col = col[ww.l])
## right side of twin
segments(xx.r, yy.r + (1/adj) * radius, .5*(xx.l+xx.r), yy.r + 3 * delta,
col = col[ww.r])
## draw midpoint MZ twin line
xx.lm <- plist$pos[i, mz.left]
yy.lm <- rep( - i, length = sum(mz.left))
xx.rm <- plist$pos[i, mz.right]
yy.rm <- rep( - i, length = sum(mz.right))
segments(.5*(xx.lm+.5*(xx.rm+xx.lm)),
.5*(yy.lm+(1/adj)*radius + yy.lm+3*delta),
.5*(xx.rm+.5*(xx.rm+xx.lm)),
.5*(yy.rm+(1/adj)*radius + yy.lm+3*delta))
## add question mark for unknown zygosity
##xx.lu <- plist$pos[i, un.left] ## Codetools 02/20/2014
##yy.lu <- rep( - i, length = sum(un.left)) ## Codetools 02/20/2014
##xx.ru <- plist$pos[i, un.right] ## Codetools 02/20/2014
##yy.ru <- rep( - i, length = sum(un.right)) ## Codetools 02/20/2014
## text('?',.5*(xx.lu+.5*(xx.ru+xx.lu)),
## .5*(yy.lu+(1/adj)*radius + yy.lu+3*delta), cex=adj*.3)
#####################################################
who <- (plist$fam[i, ] == fam)
xx <- plist$pos[i, who]
yy <- rep( - i, length = sum(who))
ww <- plist$nid[i, who]
# add connector
xx2 <- plist$pos[i,]
xx2 <- xx2[who]
## check to see if twins on either end of family
if(sum(tw)>=2){
tw.who <- tw[who]
n <- length(tw.who)
n2 <- sum(tw.who)
flagL <- sum(tw.who[1:2]) == 2
flagR <- sum(tw.who[c(n,n-1)]) == 2
xx2.save <- xx2
if(flagL) xx2[tw.who][1:2] <- rep(mean(xx2.save[tw.who][1:2]),2)
if(flagR) xx2[tw.who][c(n2,n2-1)] <- rep(mean(xx2.save[tw.who][c(n2,n2-1)]),2)
}
segments(min(xx2), 3 * delta - i, max(xx2), 3 * delta - i)
# Draw line to parents
x1 <- mean(range(xx))
y1 <- 3 * delta - i
if(branch == 0)
segments(x1, y1, parentx, - (i - 1))
else {
y2 <- 1 - i
x2 <- parentx
ydelta <- ((y2 - y1) * branch)/2
segments(c(x1, x1, x2), c(y1, y1 + ydelta, y2 - ydelta),
c(x1, x2, x2), c(y1 + ydelta, y2 - ydelta, y2))
}
}
}
# Connect up remote spouses
arcconnect <- function(x, y)
{
# draw an approx arc whose midpoint is 1/2 unit higher than its chord,
# with x[1], y[1], and x[2], y[2] as the ends of the chord
xx <- seq(x[1], x[2], length = 15)
yy <- seq(y[1], y[2], length = 15) + 0.5 - (seq(-7, 7))^2/98
lines(xx, yy, lty = 2)
}
xx <- table(np)
xx <- xx[xx > 1]
#those who appear multiple times
if(length(xx) > 0) {
multi <- as.numeric(names(xx))
for(i in 1:length(xx)) {
x2 <- xp[np == multi[i]]
y2 <- yp[np == multi[i]]
nn <- xx[i]
for(j in 1:(nn - 1))
arcconnect(x2[j + 0:1], y2[j + 0:1])
}
}
### print message if not plotting everyone
ckall <- x$id[is.na(match(x$id,x$id[plist$nid]))]
if(length(ckall>0)) cat('Did not plot the following people:',ckall,'\n')
if(!keep.par) par(oldpar)
## reorder xp and yp so that they are in the same order as plist
tmp <- plist$nid[plist$nid!=0]
xp2 <- xp[order(tmp)]
yp2 <- yp[order(tmp)]
invisible(list(plist=plist,object=x,x=xp2,y=yp2,textoff=textoff,
symbolsize=symbolsize))
}
## This was taken from my fbati R package...
mergePhePed <-
function (ped, phe)
{
if (is.null(phe) || is.null(ped)) {
return(NULL)
}
data <- merge(ped, phe, all.x = TRUE)
data <- data[order(data$pid, data$id), ]
return(data)
}
alignPhePed <-
function (ped, phe) ## modified from nuclify -- no longer nuclifies...
{
ped <- sort(ped)
phe <- sort(phe)
ndata <- mergePhePed(ped, phe)
nped <- ndata[, 1:ncol(ped)]
class(nped) <- c("ped", "data.frame")
nphe <- ndata[, c(1, 2, (ncol(ped) + 1):ncol(ndata))]
class(nphe) <- c("phe", "data.frame")
return(list(ped = nped, phe = nphe))
}
## Tom's code from pbatR:
## And then this is my own function to handle plotting a 'ped' class object
## geno \in "none", "percent", "each"
plotPed <- function( ped, sink=NULL, geno="percent", phe=NULL, pheCols=NULL ) {
##require( kinship2 ) ## replaced from 'library' for codetools...
## Bloody affectionStatus + kinship = hell on your pedigree
alev <- unique(ped$AffectionStatus)
if(any(is.na(alev)) || max(alev,na.rm=TRUE)==1){
## actually, this is what we want
}else{
## coded 0, 1, 2...
ped$AffectionStatus[ped$AffectionStatus == 0] <- NA
ped$AffectionStatus <- ped$AffectionStatus - 1
}
## is it symbolic? it can't be for these routines...
# if( is.sym(ped) )
# ped <- as.ped( ped, clearSym=TRUE )
## Handle genotype completions of individuals
extraInfo <- NULL
if(geno == "percent") {
extraInfo <- data.frame(G=paste(round(100 * rowSums(ped[,7:ncol(ped)]!=0) / (ncol(ped)-6), digits=2), "%", sep=""))
}
if(geno == "each") {
extraInfo <- data.frame(matrix("NA",nrow=nrow(ped),ncol=(ncol(ped)-6)/2),stringsAsFactors=FALSE)
colnames(extraInfo) <- pedMarkerNames(ped)
for(i in 1:ncol(extraInfo)){
a0 <- ped[[6+i*2-1]]
a1 <- ped[[6+i*2]]
#extraInfo[[i]] <- paste(min(a0,a1), "/", max(a0,a1),sep="")
extraInfo[[i]] <- paste(apply(ped[,6+i*2-c(0,1)],1,min), "/", apply(ped[,6+i*2-c(0,1)],1,max), sep="")
miss <- a0==0 | a1==0
if(any(miss))
extraInfo[[i]][miss] <- "?"
}
}
## New, put in some phenotype information?
pheInfo <- NULL
if(!is.null(phe)){
res <- alignPhePed(ped, phe)
ped <- res$ped
phe <- res$phe
}
if(!is.null(pheCols)) {
pheInfo <- phe[, match(pheCols, colnames(phe))]
if(is.null(extraInfo))
extraInfo <- pheInfo
else
extraInfo <- cbind(extraInfo, pheInfo)
}
#print("extraInfo")
#print(head(extraInfo))
## move it to their format
if( any( ped$sex==0 ) )
ped$sex[ped$sex==0] <- 3;
## Seriously, this program is just an absolute !!@#$% with affection status, I don't know what the is wrong with it's code, so we hack it below instead...
#ped$affection <- 0
#ped$affection[ped$AffectionStatus==2] <- 1
#affection <- matrix[ped$AffectionStatus]
## If sink = filename, sink each plot to a file!
## See if we should sink it to file
if( !is.null(sink) ) {
pdf( sink );
}else{
#par(ask=TRUE);
}
for( pid in unique(ped$pid) ) {
## pull out the pedigree piece
subPed <- ped[ ped$pid==pid, ]
subExtraInfo <- NULL
#if(!is.null(extraInfo))
# subExtraInfo <- as.data.frame(extraInfo[ ped$pid==pid, ]) ## Bloody R and it's vector matrix idiocy
if(!is.null(extraInfo))
subExtraInfo <- subset(extraInfo, ped$pid==pid)
## fix it so it's their program happy
pedigr <- kinship2::pedigree(id = subPed$id, dadid=subPed$idfath, momid=subPed$idmoth, sex=subPed$sex) #, affected=subPed$affection)
pedigr$affected <- matrix(subPed$AffectionStatus) ## Workaround the bug in the pedigree routine -- hey, it works!!!
#print( pedigr )
#print( str(pedigr) )
SUCCESS <- FALSE; ## sometimes it fails...
try( {
plot.pedigree(pedigr, extraInfo=subExtraInfo);
title(paste("Pedigree", pid));
SUCCESS <- TRUE;
} );
if( !SUCCESS )
print( paste( "Plotting pedigree", pid, "failed." ) );
}
## Close off the filename if necessary
if( !is.null(sink) )
dev.off()
}
Try the pbatR package in your browser
Any scripts or data that you put into this service are public.
pbatR documentation built on May 29, 2024, 10:27 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotPed plot.pedigree
Documented in plotPed
## Taken from the kinship package (plot.pedigree.R, pedigree.R, kinship.R), licensed under (GPL >= 2), and modified as below to 1) display more information, 2) sensible margins so you can see the output, 3) handle AffectionStatus without crashing, and just generally play nicer with our data.
## FROM pedigree.R in kinship package
## From plot.pedigree.R in kinship
# $Id: plot.pedigree.s,v 1.10 2003/10/24 14:13:34 Atkinson Exp $
plot.pedigree <- function(x, id = x$id, sex = x$sex, status = x$status,
affected = x$affected,
cex = 1, col = rep(1, length(x$id)),
symbolsize = 1, branch = 0.6,
packed = T, align = packed, width = 8,
density=c(-1, 50,70,90), mar=c(4.1, 1, 4.1, 1),
angle=c(90,70,50,0), keep.par=F,
extraInfo = NULL, ## Tom added
...)
{
#library(kinship2)
#col <- rep("white",length(x$id))
#col[affected==1] <- "black"
col[is.na(affected)] <- "red" ## missing!!
maxlev <- max(x$depth) + 1
n <- length(x$depth)
## DATA CHECKS
sex <- as.numeric(sex)
sex[is.na(sex)] <- 3
if(any(sex < 1 | sex > 4))
stop("Invalid sex code")
if(length(sex) != n)
stop("Wrong length for sex")
if(is.null(status))
status <- rep(0, n)
else {
if(!all(status == 0 | status == 1 | status == 2))
stop("Invalid status code")
if(length(status) != n)
stop("Wrong length for status")
}
if(!is.null(id)) {
if(length(id) != n)
stop("Wrong length for id")
}
## Define plotting symbols for sex by affected
## in S, 0=square, 1= circle, 5=diamond, 2=triangle
## 15, 16, 18 and 17 are filled versions of them
symbol <- c(0, 1, 5, 2)[sex]
## CREATE SYMBOLS/PLOTS
points.sym <- function(x, y, symbol, status, affected, size, col, aspect,
angle, density, adj)
{
circle <- function(cx, cy, r, code=0, col=1, angle, density, ...)
{
# cx, cy, coordinates for center; r is radius
z <- (0:360 * pi)/180
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
if(sum(code,na.rm=T)==0) polygon(x, y, border=T,
density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1],...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
z <- (0:180 * pi)/180
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
z <- (180:360 * pi)/180
x <- sin(z) * r + cx
y <- cos(z) * r * 1/adj + cy
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
z <- (0:90 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3],...)
z <- (180:270 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
z <- (270:360 * pi)/180
x <- c(cx, sin(z) * r + cx)
y <- c(cy, cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
z <- (0:90 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
z <- (90:180 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4],...)
z <- (180:270 * pi)/180
x <- c(cx,sin(z) * r + cx)
y <- c(cy,cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
z <- (270:360 * pi)/180
x <- c(cx, sin(z) * r + cx)
y <- c(cy, cos(z) * r * 1/adj + cy)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
square <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- cx + c(-r,-r,r,r)
y <- cy + (1/adj)*c(-r,r,r,-r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0, col=col,...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0,0,r,r)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle[3], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,-r,0,0)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(0,r,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,0,r,r)
y <- cy + (1/adj)*c(-r,0,0,-r)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
diamond <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
x <- cx + c(-r,0,r,0)
y <- cy + (1/adj)*c(0,r,0,-r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1], angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(r,-r,0)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle=angle[3], ...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,0,-r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-r,0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,0,r)
y <- cy + (1/adj)*c(-r,0,0)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
triangle <- function(cx, cy, r, code=0, col=1, angle, density, ...) {
# cx, cy, coordinates for center; r is radius
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
r <- r*1.25
a <- 3*r/sqrt(3)
b <- r/2
x <- cx + c((-1/2)*a, (1/2)*a, 0)
y <- cy + (1/adj)*c(-b, -b, r)
if(sum(code,na.rm=T)==0) polygon(x,y,border=T,density=0,col=col, ...)
else {
if(length(code)==1) polygon(x,y,border=T,col=col,
density=density[1],angle=angle[1], ...)
if(length(code)==2) {
polygon(x,y,border=T,density=0, ...)
x <- cx + c((-1/2)*a, 0, 0)
y <- cy + (1/adj)*c(-b, -b, r)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(0, (1/2)*a, 0)
y <- cy + (1/adj)*c(-b, -b, r)
polygon(x,y,border=T,col=col,density=density[2]*code[2],
angle=angle[2], ...)
}
if(length(code)==3) {
polygon(x,y,border=T,density=0, ...)
midx <- (r*(.5)*a)/(b+r)
x <- cx + c((-1/2)*a,-midx, 0, 0)
y <- cy + (1/adj)*c(-b, 0, 0, -b)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-midx, 0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,midx)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=density[3]*code[3],
angle=angle[3],...)
}
if(length(code)==4) {
polygon(x,y,border=T,density=0, ...)
midx <- (r*(.5)*a)/(b+r)
x <- cx + c((-1/2)*a, -midx , 0, 0)
y <- cy + (1/adj)*c(-b, 0, 0, -b)
polygon(x,y,border=T,col=col,density=code[1]*density[1],
angle=angle[1], ...)
x <- cx + c(-midx, 0,0)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[2]*density[2],
angle=angle[2], ...)
x <- cx + c(0,0,midx)
y <- cy + (1/adj)*c(0,r,0)
polygon(x,y,border=T,col=col,density=code[3]*density[3],
angle=angle[3], ...)
x <- cx + c(0,midx,(1/2)*a,0)
y <- cy + (1/adj)*c(0,0,-b,-b)
polygon(x,y,border=T,col=col,density=code[4]*density[4],
angle=angle[4], ...)
}
if(length(code)>4) stop('Can only plot up to 4 levels of codes')
}
invisible()
}
# in S, 0=square, 1= circle, 5=diamond, 2=triangle
# 15, 16, 18 and 17 are filled versions of them
for(i in (1:length(x))[!is.na(match(symbol,c(0,15)))]){
square(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i], angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(1,16)))]){
circle(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(5,18)))]){
diamond(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
for(i in (1:length(x))[!is.na(match(symbol,c(2,17)))]){
triangle(x[i],y[i],size, code=(affected[i,,drop=F]),
col=col[i],angle=angle, density=density)
}
# Draw slash if status=1
who <- (status == 1)
if(any(who)) {
deltax <- size*cos(pi/4) + aspect
deltay <- (1/adj)*(size*sin(pi/4) + aspect)
segments(x[who] - deltax, y[who] - deltay,
x[who] + deltax, y[who] + deltay)
}
}
#
# Now, get the structure of the plot, and lay out the main region
#
plist <- kinship2::align.pedigree(x, packed = packed, width = width, align =
align)
who <- (plist$nid > 0)
xp <- plist$pos[who]
yp <- - (row(plist$nid))[who]
np <- plist$nid[who]
oldpar <- par(mar = mar, err=-1)
par(usr = c(range(xp) + c(-0.5, 0.5), - (maxlev + 0.5), -0.5))
pin <- par()$pin
usr <- par()$usr
adj <- (pin[2]/pin[1])/((usr[4] - usr[3])/(usr[2] - usr[1]))
symbolsize <- symbolsize * adj
radius <- 0.08 * symbolsize
delta <- .1
aspect <- 0.05 * symbolsize #used for status variable
# textoff <- - (1/adj)*(radius) - .5*cex*(par()$"1em"[2])
textheight <- cex*par()$cin[2] #/4 ## ? div by 7 ??
textoff <- - (1/adj)*(radius) - .5*cex*(par()$"cin"[2]/7) #- textheight*0.1
totalTextHeight <- textheight
#cat("totalTextHeight before:", totalTextHeight, "\n")
#print(extraInfo)
#print(dim(extraInfo))
if(!is.null(extraInfo))
totalTextHeight <- textheight * (1 + ncol(extraInfo))
#cat("totalTextHeight after:", totalTextHeight, "\n")
## Tom addition to fix the margins
#print( xp )
#print( yp )
#print(totalTextHeight)
xmin <- min(xp)
xmax <- max(xp)
ymin <- min(yp)
ymax <- max(yp)
xlen <- xmax - xmin
ylen <- ymax - ymin
fudge <- max(0.1 * c(xlen, ylen)) + totalTextHeight
xmin <- xmin - fudge
xmax <- xmax + fudge
ymin <- ymin - fudge
ymax <- ymax + fudge
#cat("%%", xmin, xmax, ymin, ymax, "\n")
#plot(xp,yp,axes=F,type='n',xlab='',ylab='',...)
plot(c(xmin,xmax), c(ymin,ymax), axes=F,type='n',xlab='',ylab='',...)
points.sym(xp, yp, symbol[np], status[np], affected[np, ,drop=F],
radius, col[np], aspect, angle=angle, density=density, adj=adj)
#print(id)
#print(extraInfo)
#stop()
for(p in 1:length(np)){
textOut <- c()
if(!is.null(id))
textOut <- c(textOut, paste("ID:", id[np[p]]))
if(!is.null(extraInfo))
for(i in 1:ncol(extraInfo))
textOut <- c(textOut, paste(colnames(extraInfo)[i], ": ", extraInfo[np[p],i], sep=""))
if(length(textOut) > 0)
text(xp[p], yp[p] + textoff, paste(textOut, collapse="\n"), cex=cex, pos=1)
}
if(!is.null(id)) {
#if(is.null(extraInfo))
# text(xp, yp + textoff, paste("ID:",id[np]), cex = cex, col = col[np])
#else
# text(xp, yp + textoff, paste(paste("ID:",id[np]), extraInfo[np,1:ncol(extraInfo)], collapse="\n"), cex=cex, col=col[np])
###text(xp, yp + textoff, id[np], cex = cex, col = col[np])
### Tom's additions, so we get some more information here...
#text(xp, yp + textoff, paste("ID:",id[np]), cex = cex, col = col[np])
#
##text(xp, yp + textoff - textheight, "Bite me", cex=cex )
if(!is.null(extraInfo)){
##text(xp, yp + textoff, paste(extraInfo[np,],collapse="\n"))
#for(i in 1:ncol(extraInfo))
# text(xp, yp + textoff - textheight*i, paste(colnames(extraInfo)[i], ": ", extraInfo[np,i], sep=""))
}
}
# Draw in the family connections
for(i in 1:maxlev) {
# between spouses
if(any(plist$spouse[i, ]>0)) {
temp <- (1:ncol(plist$spouse))[plist$spouse[i, ]>0]
segments(plist$pos[i, temp] + radius, rep( - i, length(temp)),
plist$pos[i, temp + 1] - radius, rep(-i, length(temp)))
# to be added -- double segment for plist$spouse==2
}
}
for(i in 2:maxlev) {
# Children to parents, and across sibs
zed <- unique(plist$fam[i, ])
zed <- zed[zed > 0]
for(fam in zed) {
xx <- plist$pos[i - 1, fam + 0:1]
parentx <- mean(xx)
# Draw uplines from each subject
if(!is.null(plist$twins)){
tw.left <- (plist$twins[i, ] > 0 & plist$twins[i, ] < 4) & plist$fam[i,]==fam
mz.left <- (plist$twins[i, ] == 1) & plist$fam[i,]==fam
un.left <- (plist$twins[i, ] == 3) & plist$fam[i,]==fam
## figure out twins IN FAMILY (find next person in family)
tw.right <- rep(F,length(tw.left))
famlst <- plist$fam[i,]
twloc <- (1:length(tw.right))[tw.left]
famloc <- (1:length(tw.right))[famlst==fam]
flag <- NULL
for(j in twloc) flag <- c(flag,min(famloc[j<famloc]))
tw.right[flag] <- T
mz.right <- rep(F,length(mz.left))
mzloc <- (1:length(mz.right))[mz.left]
flag <- NULL
for(j in mzloc) flag <- c(flag,min(famloc[j<famloc]))
mz.right[flag] <- T
## unknown zygosity of the twin (type=3)
un.right <- rep(F,length(un.left))
unloc <- (1:length(un.right))[un.left]
flag <- NULL
for(j in unloc) flag <- c(flag,min(famloc[j<famloc]))
un.right[flag] <- T
}
if(is.null(plist$twins)){
twn <- length(plist$fam[i,]==fam)
tw.left <- rep(F,twn)
tw.right <- rep(F,twn)
mz.left <- rep(F,twn)
mz.right <- rep(F,twn)
un.left <- rep(F,twn)
un.right <- rep(F,twn)
}
tw <- tw.left|tw.right
##mz <- mz.left|mz.right ## Codetools 02/20/2014
##un <- un.left|un.right ## Codetools 02/20/2014
who <- (plist$fam[i, ] == fam) & !tw
xx <- plist$pos[i, who]
yy <- rep( - i, length = sum(who))
ww <- plist$nid[i, who]
xx.l <- plist$pos[i, tw.left]
yy.l <- rep( - i, length = sum(tw.left))
ww.l <- plist$nid[i, tw.left]
xx.r <- plist$pos[i, tw.right]
yy.r <- rep( - i, length = sum(tw.right))
ww.r <- plist$nid[i, tw.right]
segments(xx, yy + (1/adj) * radius, xx, yy + 3 * delta, col
= col[ww])
#####################################################
# draw diag line for twins
## left side of twin
segments(xx.l, yy.l + (1/adj) * radius, .5*(xx.l+xx.r), yy.l + 3 * delta,
col = col[ww.l])
## right side of twin
segments(xx.r, yy.r + (1/adj) * radius, .5*(xx.l+xx.r), yy.r + 3 * delta,
col = col[ww.r])
## draw midpoint MZ twin line
xx.lm <- plist$pos[i, mz.left]
yy.lm <- rep( - i, length = sum(mz.left))
xx.rm <- plist$pos[i, mz.right]
yy.rm <- rep( - i, length = sum(mz.right))
segments(.5*(xx.lm+.5*(xx.rm+xx.lm)),
.5*(yy.lm+(1/adj)*radius + yy.lm+3*delta),
.5*(xx.rm+.5*(xx.rm+xx.lm)),
.5*(yy.rm+(1/adj)*radius + yy.lm+3*delta))
## add question mark for unknown zygosity
##xx.lu <- plist$pos[i, un.left] ## Codetools 02/20/2014
##yy.lu <- rep( - i, length = sum(un.left)) ## Codetools 02/20/2014
##xx.ru <- plist$pos[i, un.right] ## Codetools 02/20/2014
##yy.ru <- rep( - i, length = sum(un.right)) ## Codetools 02/20/2014
## text('?',.5*(xx.lu+.5*(xx.ru+xx.lu)),
## .5*(yy.lu+(1/adj)*radius + yy.lu+3*delta), cex=adj*.3)
#####################################################
who <- (plist$fam[i, ] == fam)
xx <- plist$pos[i, who]
yy <- rep( - i, length = sum(who))
ww <- plist$nid[i, who]
# add connector
xx2 <- plist$pos[i,]
xx2 <- xx2[who]
## check to see if twins on either end of family
if(sum(tw)>=2){
tw.who <- tw[who]
n <- length(tw.who)
n2 <- sum(tw.who)
flagL <- sum(tw.who[1:2]) == 2
flagR <- sum(tw.who[c(n,n-1)]) == 2
xx2.save <- xx2
if(flagL) xx2[tw.who][1:2] <- rep(mean(xx2.save[tw.who][1:2]),2)
if(flagR) xx2[tw.who][c(n2,n2-1)] <- rep(mean(xx2.save[tw.who][c(n2,n2-1)]),2)
}
segments(min(xx2), 3 * delta - i, max(xx2), 3 * delta - i)
# Draw line to parents
x1 <- mean(range(xx))
y1 <- 3 * delta - i
if(branch == 0)
segments(x1, y1, parentx, - (i - 1))
else {
y2 <- 1 - i
x2 <- parentx
ydelta <- ((y2 - y1) * branch)/2
segments(c(x1, x1, x2), c(y1, y1 + ydelta, y2 - ydelta),
c(x1, x2, x2), c(y1 + ydelta, y2 - ydelta, y2))
}
}
}
# Connect up remote spouses
arcconnect <- function(x, y)
{
# draw an approx arc whose midpoint is 1/2 unit higher than its chord,
# with x[1], y[1], and x[2], y[2] as the ends of the chord
xx <- seq(x[1], x[2], length = 15)
yy <- seq(y[1], y[2], length = 15) + 0.5 - (seq(-7, 7))^2/98
lines(xx, yy, lty = 2)
}
xx <- table(np)
xx <- xx[xx > 1]
#those who appear multiple times
if(length(xx) > 0) {
multi <- as.numeric(names(xx))
for(i in 1:length(xx)) {
x2 <- xp[np == multi[i]]
y2 <- yp[np == multi[i]]
nn <- xx[i]
for(j in 1:(nn - 1))
arcconnect(x2[j + 0:1], y2[j + 0:1])
}
}
### print message if not plotting everyone
ckall <- x$id[is.na(match(x$id,x$id[plist$nid]))]
if(length(ckall>0)) cat('Did not plot the following people:',ckall,'\n')
if(!keep.par) par(oldpar)
## reorder xp and yp so that they are in the same order as plist
tmp <- plist$nid[plist$nid!=0]
xp2 <- xp[order(tmp)]
yp2 <- yp[order(tmp)]
invisible(list(plist=plist,object=x,x=xp2,y=yp2,textoff=textoff,
symbolsize=symbolsize))
}
## This was taken from my fbati R package...
mergePhePed <-
function (ped, phe)
{
if (is.null(phe) || is.null(ped)) {
return(NULL)
}
data <- merge(ped, phe, all.x = TRUE)
data <- data[order(data$pid, data$id), ]
return(data)
}
alignPhePed <-
function (ped, phe) ## modified from nuclify -- no longer nuclifies...
{
ped <- sort(ped)
phe <- sort(phe)
ndata <- mergePhePed(ped, phe)
nped <- ndata[, 1:ncol(ped)]
class(nped) <- c("ped", "data.frame")
nphe <- ndata[, c(1, 2, (ncol(ped) + 1):ncol(ndata))]
class(nphe) <- c("phe", "data.frame")
return(list(ped = nped, phe = nphe))
}
## Tom's code from pbatR:
## And then this is my own function to handle plotting a 'ped' class object
## geno \in "none", "percent", "each"
plotPed <- function( ped, sink=NULL, geno="percent", phe=NULL, pheCols=NULL ) {
##require( kinship2 ) ## replaced from 'library' for codetools...
## Bloody affectionStatus + kinship = hell on your pedigree
alev <- unique(ped$AffectionStatus)
if(any(is.na(alev)) || max(alev,na.rm=TRUE)==1){
## actually, this is what we want
}else{
## coded 0, 1, 2...
ped$AffectionStatus[ped$AffectionStatus == 0] <- NA
ped$AffectionStatus <- ped$AffectionStatus - 1
}
## is it symbolic? it can't be for these routines...
# if( is.sym(ped) )
# ped <- as.ped( ped, clearSym=TRUE )
## Handle genotype completions of individuals
extraInfo <- NULL
if(geno == "percent") {
extraInfo <- data.frame(G=paste(round(100 * rowSums(ped[,7:ncol(ped)]!=0) / (ncol(ped)-6), digits=2), "%", sep=""))
}
if(geno == "each") {
extraInfo <- data.frame(matrix("NA",nrow=nrow(ped),ncol=(ncol(ped)-6)/2),stringsAsFactors=FALSE)
colnames(extraInfo) <- pedMarkerNames(ped)
for(i in 1:ncol(extraInfo)){
a0 <- ped[[6+i*2-1]]
a1 <- ped[[6+i*2]]
#extraInfo[[i]] <- paste(min(a0,a1), "/", max(a0,a1),sep="")
extraInfo[[i]] <- paste(apply(ped[,6+i*2-c(0,1)],1,min), "/", apply(ped[,6+i*2-c(0,1)],1,max), sep="")
miss <- a0==0 | a1==0
if(any(miss))
extraInfo[[i]][miss] <- "?"
}
}
## New, put in some phenotype information?
pheInfo <- NULL
if(!is.null(phe)){
res <- alignPhePed(ped, phe)
ped <- res$ped
phe <- res$phe
}
if(!is.null(pheCols)) {
pheInfo <- phe[, match(pheCols, colnames(phe))]
if(is.null(extraInfo))
extraInfo <- pheInfo
else
extraInfo <- cbind(extraInfo, pheInfo)
}
#print("extraInfo")
#print(head(extraInfo))
## move it to their format
if( any( ped$sex==0 ) )
ped$sex[ped$sex==0] <- 3;
## Seriously, this program is just an absolute !!@#$% with affection status, I don't know what the is wrong with it's code, so we hack it below instead...
#ped$affection <- 0
#ped$affection[ped$AffectionStatus==2] <- 1
#affection <- matrix[ped$AffectionStatus]
## If sink = filename, sink each plot to a file!
## See if we should sink it to file
if( !is.null(sink) ) {
pdf( sink );
}else{
#par(ask=TRUE);
}
for( pid in unique(ped$pid) ) {
## pull out the pedigree piece
subPed <- ped[ ped$pid==pid, ]
subExtraInfo <- NULL
#if(!is.null(extraInfo))
# subExtraInfo <- as.data.frame(extraInfo[ ped$pid==pid, ]) ## Bloody R and it's vector matrix idiocy
if(!is.null(extraInfo))
subExtraInfo <- subset(extraInfo, ped$pid==pid)
## fix it so it's their program happy
pedigr <- kinship2::pedigree(id = subPed$id, dadid=subPed$idfath, momid=subPed$idmoth, sex=subPed$sex) #, affected=subPed$affection)
pedigr$affected <- matrix(subPed$AffectionStatus) ## Workaround the bug in the pedigree routine -- hey, it works!!!
#print( pedigr )
#print( str(pedigr) )
SUCCESS <- FALSE; ## sometimes it fails...
try( {
plot.pedigree(pedigr, extraInfo=subExtraInfo);
title(paste("Pedigree", pid));
SUCCESS <- TRUE;
} );
if( !SUCCESS )
print( paste( "Plotting pedigree", pid, "failed." ) );
}
## Close off the filename if necessary
if( !is.null(sink) )
dev.off()
}
Try the pbatR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.