load("/tmp/pairs_data.Rdata") #pairs
load(file="/framshare/fulldata.Rdata") #fulldata
paste(fulldata$idmoth,fulldata$idfath,sep="_")->fulldata$family_id
snplist<-read.csv("~/framshare/snplist.csv",as.is=TRUE)
loadnewblock<-(1:max(snplist$block))*5000


lm.fun<-function(temp.geno,pairs,fulldata) {
  #height_sd_12 ~ age_1+age_2+age_sd_12+sex_1+sex_2+genetic stuff+parental.geno
  dat<-merge(temp.geno,fulldata,by="id")
  dat[dat$generation==3,]->g3
  split(g3,g3$family_id)->fams
  sapply(fams,nrow)->index
  fams[index>1]->fams
  pair.fun<-function(x) {
    x[order(x$id),]->x
    x$snp->va
    x$id->id
    nrow(x) -> N
    tmp<-list()
    for (i in 1:(N-1)) {
      index2<-(i+1):N
      index1<-rep(i,length(index2))
      rowSums(cbind(va[index1],va[index2]))->s
      paste(id[index1],id[index2],sep="_")->id.tmp
      data.frame(pair.id=id.tmp,sum_snp=s)->tmp[[i]]
    }
    do.call("rbind",tmp)->df
  }
  lapply(fams,pair.fun)->tmp
  do.call("rbind",tmp)->tmp
  merge(pairs,tmp,by="pair.id")->pairs
  #
  subset(dat,select=c("id","snp"))->tmp
  names(tmp)[2]<-"snp.dad"
  merge(g3,tmp,by.x="idfath",by.y="id")->dat
  names(tmp)[2]<-"snp.mom"
  merge(dat,tmp,by.x="idmoth",by.y="id")->dat
  #
  subset(dat,select=c("snp.mom","snp.dad"))->tmp
  fun<-function(x) paste(sort(x),collapse=".")
  apply(tmp,1,fun)->tmp
  data.frame(family_id=dat$family_id,parent.gene=factor(tmp,levels=c("0.0","0.1","0.2","1.1","1.2","2.2")))->parent
  parent[!duplicated(parent),]->parent
  merge(pairs,parent,by="family_id")->pairs
  as.factor(pairs$sd_sex)->pairs$sd_sex
  lm(sd_height~age_1+age_2+sd_age+sd_sex+sum_snp+parent.gene,data=pairs)->mod
  summary(mod)$coef
}

#loadnewblock-length 53
library(parallel)
makeCluster(15,"SOCK")->cl
out<-list()
for (iii in 1:length(loadnewblock)) { #see here
  setwd("/framshare/snps/")
  geno<-read.csv(paste("block",iii,".csv",sep=""),as.is=TRUE)
  names(geno)[1]<-"id"
  geno<-geno[geno$id%in%fulldata$id,]
  names(geno)[-1]->nms
  outer.dat<-list()
  for (j in nms) { 
    temp.geno<-geno[,c("id",j)]
    names(temp.geno)[2]<-"snp"
    temp.geno->outer.dat[[j]]
  }
  print(round(iii/length(loadnewblock), digits=3))
  clusterApply(cl,outer.dat,lm.fun,pairs=pairs,fulldata=fulldata)->tmp
  nms->names(tmp)
  save(tmp,file=paste("~/tmp/sib_SD",iii,".Rdata",sep=""))
  tmp->out[[iii]]
}
stopCluster(cl)
do.call("c",out)->out
save(out,file="~/sib_SD.Rdata")