SPADEmissing.R

#!/usr/bin/env Rscript

args = commandArgs(trailingOnly=TRUE)

if (length(args) == 0) {
  stop("At least one argument should be provided (input file address missing).", call.=FALSE)
} else if (length(args) == 1) {
  args[2] = "out.txt"
} else if (length(args) > 2) {
  stop("No more than two arguments should be provided.", call.=FALSE)
}

if (!file.exists(args[1])) {
  cat(paste0("Error: file not found at ", args[1], "\n"))
  quit(save = "no", status = 0)
}

d=read.table(args[1],header=T,check.names=F)

cF=d$dF_A1A1+d$dF_A1A2+d$dF_A2A2
cM=d$dM_A1A1+d$dM_A1A2+d$dM_A2A2
hF=d$hF_A1A1+d$hF_A1A2+d$hF_A2A2
hM=d$hM_A1A1+d$hM_A1A2+d$hM_A2A2

d$OR_Missing_Case_FemaleVsMale=NA
d$log10P_Missing_Case_FemaleVsMale=NA
d$OR_Missing_Control_FemaleVsMale=NA
d$log10P_Missing_Control_FemaleVsMale=NA

for(i in 1:nrow(d)){
  a1=d$dF_MISSING_CT[i]; b1=cF[i]; c1=d$dM_MISSING_CT[i]; d1=cM[i]
  m1=matrix(c(a1,b1,c1,d1),2,2)
  p1=fisher.test(m1)$p.value
  o1=if(a1==0|c1==0){(a1+0.5)/(b1+0.5)/((c1+0.5)/(d1+0.5))}else{(a1/b1)/(c1/d1)}
  d$OR_Missing_Case_FemaleVsMale[i]=o1
  d$log10P_Missing_Case_FemaleVsMale[i]=-log10(p1)
  
  a2=d$hF_MISSING_CT[i]; b2=hF[i]; c2=d$hM_MISSING_CT[i]; d2=hM[i]
  m2=matrix(c(a2,b2,c2,d2),2,2)
  p2=fisher.test(m2)$p.value
  o2=if(a2==0|c2==0){(a2+0.5)/(b2+0.5)/((c2+0.5)/(d2+0.5))}else{(a2/b2)/(c2/d2)}
  d$OR_Missing_Control_FemaleVsMale[i]=o2
  d$log10P_Missing_Control_FemaleVsMale[i]=-log10(p2)
}

write.table(d,args[2],row.names=F,sep="\t",quote=F)

cat(paste0("Done! Output written to ", args[2], "\n"))