Post date: Sep 21, 2018 10:5:2 PM
1. Combine AB/BA posterior probabilities then average across chains. For whichever ancestry most likely, assign that locus that ancestry. Format for R/qtl linkage map creation.
## L1 has 241 indvs, 34284 loci
ch1<-as.matrix(read.table("L1_struc_out_popinfo1.txt_ss", sep=' ',nrow=241*34284,
blank.lines.skip=TRUE,colClasses=c(rep("numeric",6),"NULL")))
ch2<-as.matrix(read.table("L1_struc_out_popinfo2.txt_ss", sep=' ',nrow=241*34284,
blank.lines.skip=TRUE,colClasses=c(rep("numeric",6),"NULL")))
ch3<-as.matrix(read.table("L1_struc_out_popinfo3.txt_ss", sep=' ',nrow=241*34284,
blank.lines.skip=TRUE,colClasses=c(rep("numeric",6),"NULL")))
ch4<-as.matrix(read.table("L1_struc_out_popinfo4.txt_ss", sep=' ',nrow=241*34284,
blank.lines.skip=TRUE,colClasses=c(rep("numeric",6),"NULL")))
ch5<-as.matrix(read.table("L1_struc_out_popinfo5.txt_ss", sep=' ',nrow=241*34284,
blank.lines.skip=TRUE,colClasses=c(rep("numeric",6),"NULL")))
ch1<-as.matrix(transform(ch1, het=ch1[,4]+ch1[,5])[-c(4:5)])
ch2<-as.matrix(transform(ch2, het=ch2[,4]+ch2[,5])[-c(4:5)])
ch3<-as.matrix(transform(ch3, het=ch3[,4]+ch3[,5])[-c(4:5)])
ch4<-as.matrix(transform(ch4, het=ch4[,4]+ch4[,5])[-c(4:5)])
ch5<-as.matrix(transform(ch5, het=ch5[,4]+ch5[,5])[-c(4:5)])
avg<-matrix(NA,nrow=8262444,ncol=5)
avg[,1:2]<-ch1[,1:2]
avg[,3]<-apply(cbind(ch1[,3],ch2[,3],ch3[,3],ch4[,3],ch5[,3]),1,mean)
avg[,4]<-apply(cbind(ch1[,5],ch2[,5],ch3[,5],ch4[,5],ch5[,5]),1,mean)
avg[,5]<-apply(cbind(ch1[,4],ch2[,4],ch3[,4],ch4[,4],ch5[,4]),1,mean)
out<-matrix(NA, ncol=34284+1, nrow=(241)+1)
ids<- read.table('L1_ids.txt', header=FALSE,)
loci<-read.table('locusids2.txt', header=FALSE)
out[2:242,1]<-as.character(ids[,1])
out[1,2:34285]<-as.character(loci[,1])
out[1,1]<-'ids'
max<-apply(avg[,3:5],1,which.max)
# colnames(avg)[3:5]<-c("L","H","M")
for(i in 1:241){
#for individual 1, get all loci, then move to ind 2
locus<-seq(0+i,8262444,241)
for(j in 1:34284){
if(max[locus[j]]==1){out[i+1,j+1]<-'A'}
else{if(max[locus[j]]==2){out[i+1,j+1]<-'H'}
else{if(max[locus[j]]==3){out[i+1,j+1]<-'B'}}}
}}
test<-rbind(out[1,],c("",rep(1,34284)),out[-1,])
write.table(test,'L1_map_in.csv',quote=FALSE,sep = ',', na='NA', row.names=FALSE, col.names=FALSE)
2. Create linkage map
library(qtl)
mapthis<-read.cross("csv",file="L1_map_in.csv",estimate.map=FALSE)
mapthis <- est.rf(mapthis)