rm(list=ls()) setwd("C:/Users/shuzhang/Desktop/Research/Zsl/PS/spatial/paper/Scientific_Report/program/") source("subfunction.R") #source("C:\\Users\\Shuanglin Zhang\\Desktop\\Research\\Zsl\\PS\\spatial\\kpop\\simulation.R") qq=scan("freq.txt") rep1=5 sample=1000 L=10000 Fst=0.01 mu=0.2 mm=4 pvper=matrix(0,L*rep1,mm) power=matrix(0,9,mm+3) nn0=0 for(k in c(2,10,20)) { if(k==2) npop=1 else npop=2 if(k==20) mu=5 else if (k==2) mu=1 else mu=2 name=paste("pc2a_f01_",k,"pop.txt",sep="") pc0=read.table(file=name,header=FALSE, skip=1) pc0=as.matrix(pc0[,-1]) for(kk in c(10)) { nn0=nn0+1 LL=L*kk for(j in 1:rep1) { print(paste("k=",k,"kk=",kk,"rep1=",j)) xx=gen_kpop_r(sample, 2*LL, Fst, k, qq, mu, npop) y=xx[[1]] y=quan_to_qual(y,0.16) xx=xx[[2]] pp=colMeans(xx)/2 xx=xx[,which(pp>=0.002)[1:LL]] m1=(j-1)*L+1 m2=j*L pvper[m1:m2,1]=uncorrect_region1(y,xx,kk) pvper[m1:m2,2]=GC(pvper[m1:m2,1]) pvper[m1:m2,3]=pc_linear_region1(y, xx, pc0[,1:10], kk) aa=sample(c(1:LL),L,replace=FALSE) npc=5 h1=new_choose_h(y,xx[,aa], pc0, 10) pvper[m1:m2,4]=pc_nonp_region1(y, xx, pc0, 10, h1, kk) } for(jj in 1:3) { nn=jj+(nn0-1)*3 if(jj==1) alpha0=0.05 else if(jj==2) alpha0=0.01 else alpha0=0.001 power[nn,1]=k power[nn,2]=kk power[nn,3]=alpha0 for(ii in 1:mm) { power[nn,(ii+3)]=sum(pvper[,ii]