Post date: Dec 17, 2014 12:3:59 AM
x<-matrix(scan("betasStripe.txt",n=340687347*2,sep=" "),nrow=340687347,ncol=2,byrow=T)
nl<-3743817
nc<-91
beta<-matrix(x[,2],nrow=nl,ncol=91,byrow=F)
pip<-matrix(x[,1],nrow=nl,ncol=91,byrow=F)
b1<-apply(beta[,1:45],1,mean)
b2<-apply(beta[,46:91],1,mean)
#data: b1 and b2
#t = 1818.858, df = 3743815, p-value < 2.2e-16
#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# 0.6843833 0.6854588
#sample estimates:
# cor
#0.6849215
pip1<-apply(pip[,1:45],1,mean)
pip2<-apply(pip[,46:91],1,mean)
#Pearson's product-moment correlation
#data: pip1 and pip2
#t = 3340.014, df = 3743815, p-value < 2.2e-16
#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# 0.8650365 0.8655455
#sample estimates:
# cor
#0.8652912
bm1<-b1*pip1
bm2<-b2*pip2
cor.test(bm1,bm2)
# Pearson's product-moment correlation
#
#data: bm1 and bm2
#t = 5717.588, df = 3743815, p-value < 2.2e-16
#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# 0.9471263 0.9473345
#sample estimates:
# cor
#0.9472305
## finals
b<-(b1+b2)/2
pip<-(pip1+pip2)/2
bm<-(bm1+bm2)/2
## pip > 0.1
lg8_ord15_scaf270-440236, C, A
pip = 0.109, b = -0.877, bm = -0.0936
lg8_ord37_scaf262-658900, T, C
pip = 0.221, b = -1.692, bm = -0.381
grep -n lg8_ord15_scaf270-440236 maf0.05LocList*txt
#maf0.05LocListA.txt:3421209:lg8_ord15_scaf270-440236, C, A,
#maf0.05LocListComb.txt:3226824:lg8_ord15_scaf270-440236, C, A
#maf0.05LocListC.txt:3418331:lg8_ord15_scaf270-440236, C, A,
grep -n lg8_ord37_scaf262-658900 maf0.05LocList*txt
#maf0.05LocListA.txt:3464055:lg8_ord37_scaf262-658900, T, C,
#maf0.05LocListComb.txt:3267745:lg8_ord37_scaf262-658900, T, C
#maf0.05LocListC.txt:3461399:lg8_ord37_scaf262-658900, T, C,
grep lg8_ord15_scaf270-440236 maf05timemaA*/output/pimass_timemaArep*.mcmc.txt | cut -f 4-5 | perl -p -i -e 's/^ //;s/\s+/ /' > betaAlg8_ord15_scaf270-440236.txt
bAs270_440236<-read.table("../output_pimassMaf5/betaAlg8_ord15_scaf270-440236.txt",header=F)
mean(bAs270_440236[,2])
[1] -0.0305688
grep lg8_ord15_scaf270-440236 maf05timemaC*/output/pimass_timemaArep*.mcmc.txt | cut -f 4-5 | perl -p -i -e 's/^ //;s/\s+/ /' > betaClg8_ord15_scaf270-440236.txt
bCs270_440236<-read.table("../output_pimassMaf5/betaClg8_ord15_scaf270-440236.txt",header=F)
mean(bCs270_440236[,2])
[1] -0.0200815
grep lg8_ord37_scaf262-658900 maf05timemaA*/output/pimass_timemaArep*.mcmc.txt | cut -f 4-5 | perl -p -i -e 's/^ //;s/\s+/ /' > betaAlg8_ord37_scaf262-658900.txt
bAs262_658900<-read.table("../output_pimassMaf5/betaAlg8_ord37_scaf262-658900.txt",header=F)
mean(bAs262_658900[,2])
#[1] -0.0073228
grep lg8_ord37_scaf262-658900 maf05timemaC*/output/pimass_timemaCrep*.mcmc.txt | cut -f 4-5 | perl -p -i -e 's/^ //;s/\s+/ /' > betaClg8_ord37_scaf262-658900.txt
bCs262_658900<-read.table("../output_pimassMaf5/betaClg8_ord37_scaf262-658900.txt",header=F)
mean(bCs262_658900[,2])
[1] [1] -0.016062