Commit 5f32ac2d authored by Helene Rimbert's avatar Helene Rimbert
Browse files

initial commit

parent 1fb01e23
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
############
##########
#GPS_opt
#############
##############
##First read in the arguments listed at the command line
args=(commandArgs(TRUE))
##args is now a list of character vectors
## First check to see if arguments are passed.
## Then cycle through each element of the list and evaluate the expressions.
if(length(args)==0){
print("No arguments supplied.")
##supply default values
#a = 1
#b = c(1,1,1)
run=1
scenario=1
}else{
for(i in 1:length(args)){
eval(parse(text=args[[i]]))
}
}
print(paste("run",run))
print(paste("scenario",scenario))
library(Rcpp)
Sys.setenv("PKG_CXXFLAGS"=paste("-std=c++11 -O3 -I", '/work/sbensadoun/comp_PS_GS/script/sor/stats/include' , sep=""))
sourceCpp("/work/sbensadoun/comp_PS_GS/script/sor/src/sor.cpp", verbose = T)
source('/work/sbensadoun/comp_PS_GS/script/sor4.R')
library(parallel)
#################
source("/work/sbensadoun/comp_PS_GS/script/function_opt.r")
library(rrBLUP)
library(data.table)
ncpu=1
outdir="/work/sbensadoun/comp_PS_GS/GPS/opt_cross"
geno.file="/work/sbensadoun/comp_PS_GS/donnees_entree/wp4_pheno_prune0.5_01_nohet.txt"
map.file="/work/sbensadoun/comp_PS_GS/donnees_entree/wp4_pheno_prune0.5_def.map"
design.file="/work/sbensadoun/comp_PS_GS/donnees_entree/design_190308.txt"
effect.file=paste("WP4_12K_100effect_run",run,".txt",sep="")
setwd(outdir)
design=read.table(design.file,header=T)
design$strategy="GPS"
design$crossing="opt"
#run=1
parameter=as.vector(as.matrix(design[scenario,]))
names(parameter)=colnames(design)
param.PS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"PS.NC","PS.N2","PS.N3","PS.N4",
"alpha4PS","alpha3PS","alphaPS",
"CTPS",
"PS.n2",
"PS.n3",
"PS.n4",
"years",
"CTPS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GPS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GPS.NC","GPS.N2","GPS.N3","GPS.N4",
"alpha4GPS","alpha3GPS","alphaGPS",
"CTGPS",
"GPS.n2",
"GPS.n3",
"GPS.n4",
"years",
"CTGPS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GPS.NC=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GPSn2.NC","GPSn2.N2","GPSn2.N3","GPSn2.N4",
"alpha4GPSn2","alpha3GPSn2","alphaGPSn2",
"CTGPSn2",
"GPSn2.n2",
"GPSn2.n3",
"GPSn2.n4",
"years",
"CTGPSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GS.NC","GS.N2","GS.N3","GS.N4",
"alpha4GS","alpha3GS","alphaGS",
"CTGS",
"GS.n2",
"GS.n3",
"GS.n4",
"years",
"CTGS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GS.NC=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GSn2.NC","GSn2.N2","GSn2.N3","GSn2.N4",
"alpha4GSn2","alpha3GSn2","alphaGSn2",
"CTGSn2",
"GSn2.n2",
"GSn2.n3",
"GSn2.n4",
"years",
"CTGSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
parameter=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","NC","K","CT","CTgs","lambda","alpha2",
"PS.NC","PS.N2","PS.N3","PS.N4",
"GPS.NC","GPS.N2","GPS.N3","GPS.N4",
"GPSn2.NC","GPSn2.N2","GPSn2.N3","GPSn2.N4",
"GS.NC","GS.N2","GS.N3","GS.N4",
"GSn2.NC","GSn2.N2","GSn2.N3","GSn2.N4",
####
"alpha4GPS","alpha3GPS","alphaGPS",
"alpha4GPSn2","alpha3GPSn2","alphaGPSn2",
"alpha4PS","alpha3PS","alphaPS",
"alpha4GS","alpha3GS","alphaGS",
"alpha4GSn2","alpha3GSn2","alphaGSn2",
"CTPS",
"CTGPS",
"CTGPSn2",
"CTGS",
"CTGSn2",
"PS.n2",
"GPS.n2",
"GPSn2.n2",
"GS.n2",
"GSn2.n2",
"PS.n3",
"GPS.n3",
"GPSn2.n3",
"GS.n3",
"GSn2.n3",
"PS.n4",
"GPS.n4",
"GPSn2.n4",
"GS.n4",
"GSn2.n4",
"years",
"years.gs",
"CTPS.year",
"CTGPS.year",
"CTGPSn2.year",
"CTGS.year",
"CTGSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param=param.GPS
names(param)=c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"NC","N2","N3","N4",
"alpha4","alpha3","alpha",
"CT",
"n2",
"n3",
"n4",
"years",
"CT.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)
NC=as.numeric(param["NC"])
NP=as.numeric(param["NP"])
K=as.numeric(param["K"])
N3=as.numeric(param["N3"])
N4=as.numeric(param["N4"])
bycross4=as.numeric(param["bycross4"])
bycrossp=as.numeric(param["bycrossp"])
size4=as.numeric(param["size4"])
sizep=as.numeric(param["sizep"])
crossing=param["crossing"]
strategy=param["strategy"]
h2=as.numeric(param["h2"])
header=paste(strategy,crossing,sep="_")
header2="rep1_WP4_12K_100QTL"
res=c("strategy","scenario","run","year","CG","K","bycross4","bycrossp","size4","sizep","h2","NC","N2","N3","N4","TBV10","TBV25","TBV50","TBV100","TBV200","rate10","rate25","rate50","rate100","rate200","rate","rateNP10","rateNP25","rateNP50","rateNP100","rateNP200","nbqtl10","nbqtl25","nbqtl50","nbqtl100","nbqtl200","nbal","accu","rmse","var.tbv")
write.table(t(data.frame(res)),paste(header2,"_",header,"_scenario",scenario,"_run",run,"_stat.txt",sep=""),quote=F,col.names=F,row.names=F,sep="\t")
####
y=as.numeric(param["years"])
years=seq(5,y,5)
maxy=max(years)
effects=as.vector(read.table(effect.file)[,1])
geno2=read.table(geno.file,check.names=F)
map2=read.table(map.file)
TBV=generate_BV(as.matrix(geno2),effects)
geno3=geno2[,effects==0]
map3=map2[effects==0,]
pheno=prod_noise(TBV,h2)
prog=geno2
if (bycross4==0)
{
n3s=balance(N3,NC)
}
names=1001:(1000+NP)
prog3=geno3
names(pheno)=rownames(geno2)
estimated_effects=estimate_effect(pheno,geno3)
map=as.matrix(map2)[,c(1,3)]
#map[,2]=round(map2[,4]/4000000,2)
map=apply(map,2,as.numeric)
map[,2]=round(map[,2],2)
year=0
nam=paste(header2,"_",header,"_scenario",scenario,"_run",run,"_year",year,"_",sep="")
write.table(data.frame(estimated_effects),paste(nam,"_estimated_effects.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
fwrite(data.frame(prog),paste(nam,"_progeny.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
fwrite(data.frame(pheno),paste(nam,"_pheno.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
nam2=paste(header2,"_",header,"_scenario",scenario,"_run",run,sep="")
fwrite(data.frame(prog),paste(nam2,"_pro.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
fwrite(data.frame(pheno),paste(nam2,"_pheno.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
prog=prog[rev(order(pheno)),]
pheno=pheno[rev(order(pheno)),]
prog3=prog[,effects==0]
for (i in 1:K){
for (j in 1:3){
print(i)
print(j)
#year=(12*(i-1))+(4*j)
year=(15*(i-1))+(5*j)
nam=paste(header2,"_",header,"_scenario",scenario,"_run",run,"_year",year,"_",sep="")
if (crossing=="random")
{
cross=cross_random2(Nind=NP,NC=NC)
cr=cbind(rownames(prog)[cross[,1]],rownames(prog)[cross[,2]])
cross=data.frame(cross)
}
if (crossing=="opt")
{
cr=cross_value_best_chr(estimated_effects,as.matrix(prog3),map3)
print("cr1 ok")
cr=as.matrix(cr[1:NC,1:2])
cr=apply(cr,2,as.numeric)
print("cr2 ok")
cross=data.frame(cr)
cr=cbind(rownames(prog)[as.numeric(as.vector(cr[,1]))],rownames(prog)[as.numeric(as.vector(cr[,2]))])
}
########
data=cbind.data.frame(rownames(prog),prog)
npro=n3s
nam=paste(header2,"_",header,"_scenario",scenario,"_run",run,"_year",year,"_",sep="")
filePrefix=nam
compressFile=TRUE
sor_parallel(map, cross, data, filePrefix=nam, compressFile=TRUE, ncpu=ncpu,npro=npro)
res=vector("list",nrow(cross))
nam=paste(header2,"_",header,"_scenario",scenario,"_run",run,"_year",year,"_",sep="")
nam2=paste(nam,as.vector(cross[,1]),"_x_",as.vector(cross[,2]),sep="")
for (l in 1:length(nam2)){
#print(l)
res[[l]]=sor_read_file(nam2[l])/2
}
file.remove(nam2)
##############
print("pro ok")
date()
result=cycle_GPS(header2=header2,header=header,year=year,cross=cr,pro.sp=res,parameters=param,scenario=scenario,run=run,estimated_effects=estimated_effects,effects=effects,names=names,n3s=n3s)
date()
print("cycle ok")
ped=result[[1]]
prog=result[[2]]
pheno4=result[[3]]
prog4=result[[4]]
nam=paste(header2,"_",header,"_scenario",scenario,"_run",run,"_year",year,"_",sep="")
colnames(prog)=as.vector(map2[,2])
fwrite(data.frame(prog),paste(nam,"_progeny.txt",sep=""),quote=F,row.names=T,col.names=T,sep="\t")
fwrite(data.frame(ped),paste(nam,"_pedigree.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
nam2=paste(header2,"_",header,"_scenario",scenario,"_run",run,sep="")
fwrite(data.frame(prog4),paste(nam2,"_pro.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t",append=T)
fwrite(data.frame(pheno4),paste(nam2,"_pheno.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t",append=T)
if (year!=maxy)
{
names=names+NP
prog.tmp=fread(paste(nam2,"_pro.txt",sep=""),data.table=F)
pheno.tmp=fread(paste(nam2,"_pheno.txt",sep=""),data.table=F)
estimated_effects=estimate_effect(as.vector(as.matrix(pheno.tmp)),as.matrix(prog.tmp[,effects==0]))
write.table(data.frame(estimated_effects),paste(nam,"_estimated_effects.txt",sep=""),quote=F,row.names=F,col.names=F,sep="\t")
prog3=prog[,effects==0]
print("read ok")
}
}
}
############
##########
#GPS_random
#############
##############
##First read in the arguments listed at the command line
args=(commandArgs(TRUE))
##args is now a list of character vectors
## First check to see if arguments are passed.
## Then cycle through each element of the list and evaluate the expressions.
if(length(args)==0){
print("No arguments supplied.")
##supply default values
#a = 1
#b = c(1,1,1)
run=1
scenario=1
}else{
for(i in 1:length(args)){
eval(parse(text=args[[i]]))
}
}
print(paste("run",run))
print(paste("scenario",scenario))
library(Rcpp)
Sys.setenv("PKG_CXXFLAGS"=paste("-std=c++11 -O3 -I", '/work/sbensadoun/comp_PS_GS/script/sor/stats/include' , sep=""))
sourceCpp("/work/sbensadoun/comp_PS_GS/script/sor/src/sor.cpp", verbose = T)
source('/work/sbensadoun/comp_PS_GS/script/sor4.R')
library(parallel)
#################
source("/work/sbensadoun/comp_PS_GS/script/function_opt.r")
library(rrBLUP)
library(data.table)
ncpu=1
outdir="/work/sbensadoun/comp_PS_GS/GPS/random"
geno.file="/work/sbensadoun/comp_PS_GS/donnees_entree/wp4_pheno_prune0.5_01_nohet.txt"
map.file="/work/sbensadoun/comp_PS_GS/donnees_entree/wp4_pheno_prune0.5_def.map"
design.file="/work/sbensadoun/comp_PS_GS/donnees_entree/design_190308.txt"
effect.file=paste("WP4_12K_100effect_run",run,".txt",sep="")
setwd(outdir)
design=read.table(design.file,header=T)
design$strategy="GPS"
design$crossing="random"
#run=1
parameter=as.vector(as.matrix(design[scenario,]))
names(parameter)=colnames(design)
param.PS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"PS.NC","PS.N2","PS.N3","PS.N4",
"alpha4PS","alpha3PS","alphaPS",
"CTPS",
"PS.n2",
"PS.n3",
"PS.n4",
"years",
"CTPS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GPS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GPS.NC","GPS.N2","GPS.N3","GPS.N4",
"alpha4GPS","alpha3GPS","alphaGPS",
"CTGPS",
"GPS.n2",
"GPS.n3",
"GPS.n4",
"years",
"CTGPS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GPS.NC=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GPSn2.NC","GPSn2.N2","GPSn2.N3","GPSn2.N4",
"alpha4GPSn2","alpha3GPSn2","alphaGPSn2",
"CTGPSn2",
"GPSn2.n2",
"GPSn2.n3",
"GPSn2.n4",
"years",
"CTGPSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GS=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GS.NC","GS.N2","GS.N3","GS.N4",
"alpha4GS","alpha3GS","alphaGS",
"CTGS",
"GS.n2",
"GS.n3",
"GS.n4",
"years",
"CTGS.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
param.GS.NC=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","K","CT","lambda","alpha2",
"GSn2.NC","GSn2.N2","GSn2.N3","GSn2.N4",
"alpha4GSn2","alpha3GSn2","alphaGSn2",
"CTGSn2",
"GSn2.n2",
"GSn2.n3",
"GSn2.n4",
"years",
"CTGSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"
)]
parameter=parameter[c("scenario","CC","CDH","CM2","CM3","CM4","CP3","CP4","CG","Nref","NP","NC","K","CT","CTgs","lambda","alpha2",
"PS.NC","PS.N2","PS.N3","PS.N4",
"GPS.NC","GPS.N2","GPS.N3","GPS.N4",
"GPSn2.NC","GPSn2.N2","GPSn2.N3","GPSn2.N4",
"GS.NC","GS.N2","GS.N3","GS.N4",
"GSn2.NC","GSn2.N2","GSn2.N3","GSn2.N4",
####
"alpha4GPS","alpha3GPS","alphaGPS",
"alpha4GPSn2","alpha3GPSn2","alphaGPSn2",
"alpha4PS","alpha3PS","alphaPS",
"alpha4GS","alpha3GS","alphaGS",
"alpha4GSn2","alpha3GSn2","alphaGSn2",
"CTPS",
"CTGPS",
"CTGPSn2",
"CTGS",
"CTGSn2",
"PS.n2",
"GPS.n2",
"GPSn2.n2",
"GS.n2",
"GSn2.n2",
"PS.n3",
"GPS.n3",
"GPSn2.n3",
"GS.n3",
"GSn2.n3",
"PS.n4",
"GPS.n4",
"GPSn2.n4",
"GS.n4",
"GSn2.n4",
"years",
"years.gs",
"CTPS.year",
"CTGPS.year",
"CTGPSn2.year",
"CTGS.year",
"CTGSn2.year",
"bycross4",
"bycrossp",
"size4",
"sizep",
"h2",
"strategy",
"crossing"