make phenotype residuals from covariates

make_residuals.R

library(data.table)

arguments <- commandArgs(T)

phenfile <- arguments[1]
covfile <- arguments[2]
outfile <- arguments[3]
lmodel <- arguments[4]

# expect lmodel to be "logistic" or "linear"

stopifnot(lmodel %in% c("logistic", "linear"))

phen <- fread(phenfile)
names(phen)[1:2] <- c("FID", "IID")

cov <- fread(covfile)
names(cov)[1:2] <- c("FID", "IID")

phen <- subset(phen, FID %in% cov$FID)
cov <- subset(cov, FID %in% phen$FID)

index <- match(cov$FID, phen$FID)
cov <- cov[index, ]

stopifnot(all(phen$FID == cov$FID))

dat <- data.frame(pxxxxxx=phen[,3], cov[,-c(1:2)])
dat$pxxxxxx[1] <- NA

f <- paste0("pxxxxxx ~ ", paste(names(cov)[-c(1:2)], collapse=" + "))

if(lmodel == "linear")
{
	dat$resid <- residuals(lm(f, data=dat, na.action="na.exclude"))
} else {
	dat$resid <- residuals(glm(f, data=dat, na.action="na.exclude", family="binomial"))
}

write.table(data.frame(FID=phen$FID, IID=phen$IID, y=dat$resid), file=outfile, row=FALSE, col=FALSE, qu=FALSE)