## current file paths are different from the ones used in my own computer ## all scripts are kept in script folder ## final vcf and results are kept in folder data ## contact cjnankai@hotmail.com for further questions ################################################################################ ###### raw reads alignment # all runs have been performed in Uppmax super cluster in Uppsala University sbatch -o log.txt genome_mapping.sh fq_folder out_folder # fq_folder: folders for fastq, out_folder for output g.vcf folder # merge g.vcf files of one plate into one g.vcf file sbatch -o log.txt call-merge-gvcf.sh gvcf_in_folder gvcf_out_folder platexx 13093 # # call genotypes from g.vcf files sbatch -o log.txt call-genotype-gvcf.sh gvcf_folder output.vcf # apply recalibration sbatch -o log.txt ApplyRecalibration.sh output.vcf # final vcf file named "spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate.vcf.gz" ### filtering snp with mr >0.5, dp<2 ## LD remove SNPs of significant p<=0.05 and r2>=0.2 ## remove 4 Inds with genotype mr >0.9 ## only keep noncoding snps perl filterSNP-by-list.pl spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_snpList.txt \ /Users/chenjun/Desktop/exomeCapture/mapping/vcf/spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50.vcf.gz \ spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding.vcf # remove all individuals with more than 90% SNP failed R count_miss=function(x){ round(sum(grepl('\\.', x))/399801,2) } require(data.table) data<-fread('spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding.vcf') res.miss<-apply(data[,10:ncol(data)],2,count_miss) bad.ind<-which(res.miss>=0.9)+9 write.table(data[,-bad.ind],file='~/Desktop/exomeCapture/structure/input/spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90.vcf',quote=F,row.names=F,col.names=T,sep="\t") ###### spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90.vcf is the final VCF ## stitch chr ## in order to generate input file for Admixture, we joined contigs randomly by stretches of 'N' perl stitch-chr.pl spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90.vcf \ spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90_stitched.vcf.gz 1>noncoding_snp_stitched_map.txt # generate input file for admixture plink --vcf /Users/chenjun/Desktop/exomeCapture/structure/input/all_unlinked/spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_mr90_stitched.vcf.gz \ --recode --make-bed --out /Users/chenjun/Desktop/exomeCapture/EIG-6.1.4/input/all_unlinked/spruce_exomeCapture_allPlates_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_mr90_stitche # run admixture K=1 to 10 sbatch -o log.txt run_admixture.sh # run eigenSoft ./smartpca -p par.example #par.example input file for eigenSoft # run treemix ## bootstrap for 100 iterations with K=100 and get consensus tree and support score sbatch run_treemix_boot.sh spruce_exomeCapture_allPlates_FullDateSetRG_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90_stitched.treeMix.gz ### use the consensus tree topology to infer admixture events sbatch run_treemix.sh spruce_exomeCapture_allPlates_FullDateSetRG_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90_stitched.treeMix.gz ########################################## simulation ########################################## ## 1) subsampling allele frequency to create 3D observed SFS # we need at least 10 alleles from each pop and from outgroup perl convert-vcf2snp-spruce.pl spruce_exomeCapture_allPlatesNP92_RG0802_diploid_scaffolds_VarRecalibrate_AD2_MR50_unlinked_noncoding_mr90_stitched.vcf.gz allele_counts/ perl ../polarize-allele-freq.pl Omorika-allele-freq.txt SE9-allele-freq.txt RO-allele-freq.txt BY-allele-freq.txt ALP-allele-freq.txt Kirov-allele-freq.txt > combined-allele-freq.txt ### polarize allele freq source('~/Desktop/exomeCapture/simulations/functions.r') data<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/combined-allele-freq.txt",header=T) DAF<-apply(data[,-1],1,get_DAF) my_DAF<-data.frame(data$snp,t(DAF)) names(my_DAF)<-names(data) write.table(my_DAF,file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/subsampled-n10-combined-polarized-allele-freq.txt',quote=F,row.names=F,sep="\t") ######## get 3D SFS matrix # sfs1<-get_3D_SFS(my_DAF[,c(2,5,6)]) ## BY ROM SE !!! order is important for sfs ## not used for the publication sfs2<-get_3D_SFS(my_DAF[,c(3,5,6)]) # ALP, ROM, SE # write.table(t(sfs1),file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs3D/BYROMSE_DSFS.obs',sep="\t",quote=F,row.names=F,col.names=F) write.table(t(sfs2),file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs3D/ALPROMSE_DSFS.obs',sep="\t",quote=F,row.names=F,col.names=F) ###### convert 3D SFS to 2D # sfs1_2d<-get_2d_from_3d(sfs1) ## we remove the 1st and last cell in the sfs 3d vector because these two values cannot be estimated from 3pop simulation sfs2_2d<-get_2d_from_3d(sfs2) ## 2d joint SFS for pairwise 3 pop simulation #### only for 2 pop simulation # write.table(sfs1_2d[[1]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/BYROMSE_jointDAFpop1_0.obs',quote=F,row.names=T,col.names=T,sep="\t") # write.table(sfs1_2d[[2]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/BYROMSE_jointDAFpop2_0.obs',quote=F,row.names=T,col.names=T,sep="\t") # write.table(sfs1_2d[[3]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/BYROMSE_jointDAFpop2_1.obs',quote=F,row.names=T,col.names=T,sep="\t") write.table(sfs2_2d[[1]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/ALPROMSE_jointDAFpop1_0.obs',quote=F,row.names=T,col.names=T,sep="\t") write.table(sfs2_2d[[2]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/ALPROMSE_jointDAFpop2_0.obs',quote=F,row.names=T,col.names=T,sep="\t") write.table(sfs2_2d[[3]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/sfs2D/ALPROMSE_jointDAFpop2_1.obs',quote=F,row.names=T,col.names=T,sep="\t") ###### 2) run fastsimcoal2 SFS ##### we used 3D SFS to simulate and use 2D for plots. detail tpl and est files shown in the input folder fsc25221 -t bifurcate.tpl -e bifurcate.est -M0.05 -n100000 -N100000 -l10 -L20 -d -u obs1<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE_jointDAFpop1_0.obs",skip=1,row.names=1,header=T) ## alp- rom obs2<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE_jointDAFpop2_0.obs",skip=1,row.names=1,header=T) ## alp- SE obs3<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE_jointDAFpop2_1.obs",skip=1,row.names=1,header=T) ## rom- SE model1<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE/ALPROMSE_jointDAFpop1_0.txt",header=T,row.names=1) model2<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE/ALPROMSE_jointDAFpop2_0.txt",header=T,row.names=1) model3<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/ALPROMSE/ALPROMSE_jointDAFpop2_1.txt",header=T,row.names=1) # -u model<-as.matrix(read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/triadent/triadent/triadent_DSFS.txt",skip=2)) model_3d<-get_2d_from_3d(model) ## model[0]==0 and model[1331]=0 model1<-model_3d[[1]] model2<-model_3d[[2]] model3<-model_3d[[3]] ## for 2D joint SFS in 3pop simulation since no value estimated at [0,0] cell by fastsimcoal, we mask this cell as well ## we further mask the last cell [10,10] in comparison with outgroup res1<-eval_2d_sfs(model1,obs1,las=1,xlab='ALP',ylab='ROM',breaks=seq(-4.8,-0.6,length=20)) # we can actually fit [10,10] using SE as outgroup p1<-recordPlot() res2<-eval_2d_sfs(model2,obs2,las=1,xlab='ALP',ylab='SE',breaks=seq(-4.8,-0.6,length=20)) # unable to estimate fixation cell [10,10] without outgroup p2<-recordPlot() res3<-eval_2d_sfs(model3,obs3,las=1,xlab='ROM',ylab='SE',breaks=seq(-4.8,-0.6,length=20)) # unable to estimate fixation cell [10,10] without outgroup p3<-recordPlot() res1<-eval_2d_sfs(model1,obs1,las=1,xlab='ALP',ylab='ROM',breaks=seq(-4.8,-0.6,length=20),mask_mono=F) # for 3d matrix,we can fit the mono and fixation res2<-eval_2d_sfs(model2,obs2,las=1,xlab='ALP',ylab='SE',breaks=seq(-4.8,-0.6,length=20),mask_mono=F) # res3<-eval_2d_sfs(model3,obs3,las=1,xlab='ROM',ylab='SE',breaks=seq(-4.8,-0.6,length=20),mask_mono=F) # fst_by_alp<-apply(my_DAF[,2:3],1,calc_hdFst_separate) apply(fst_by_alp,1,sum) ### gain confident intervals by booststrap 100 runs fsc25221 -i bifurcate_maxL.par -n100 -j -d -s0 -x –I -u -q bst<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/bifurcate/test/bifurcate-3Dfinal/booststrap_100runs/3D-confident-interval/bifurcate_bst100_3dSFS_bestlhoods.txt",header=T) ## model selection mig<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/bifurcate/test/bifurcate-3Dfinal/model_comparison_100estimates/model_comparisons_100estimates_bifurcateMig.txt",header=T) nomig<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/bifurcate/test/bifurcate-3Dfinal/model_comparison_100estimates/model_comparisons_100estimates_bifurcateNoMig.txt",header=T) res<-matrix(numeric(100*2),nrow=100) for(i in 1:100){ aic_x=2*(6-mig$MaxEstLhood[i]) aic_y=2*(6-nomig$MaxEstLhood[i]) delta_aic_x=aic_x-min(aic_x,aic_y) delta_aic_y=aic_y-min(aic_x,aic_y) wi_x=exp(-0.5*delta_aic_x) wi_y=exp(-0.5*delta_aic_y) res[i,]=c(wi_x/(wi_x+wi_y),wi_y/(wi_x+wi_y)) } sum(res[,1]>res[,2]) ## 95% runs support migration != 0 [1] 95 est<-read.table("~/Desktop/exomeCapture/simulations/fsc_mac64/models/ALP-RO-SE/bifurcate/test/bifurcate-3Dfinal/model_comparison_100estimates/bifurcate_run9.bestlhoods",header=T) plot(c(10,2e7),c(1e4,1e7),log='xy',xlab='Time (year ago)',ylab=expression(N[e]),type='n',yaxt='n',xaxt='n') rect(1e7,1e3,2e7,2e7,col='gray',border='gray') axis(2,at=c(1e4,1e5,1e6,1e7)) axis(1,at=c(10,1e3,1e5,1e7),lab=c(0,'1e3', 1e5,1e7)) text(1e7,1e5,paste('Population','divergence',sep="\n")) plot_demo(bst,Ne=bst$N0,col='#FF999910',extend=T) # booststrap plot_demo(est,Ne=est$N1,lwd=3,col=4,lty=2) plot_demo(est,Ne=est$N2,lwd=3,col='green',lty=4) plot_demo(est,Ne=est$N0,lwd=3,col=2,lty=1,extend=T) legend('topleft', legend=c('Alp','Carp','Fenno'),col=c(2,4,'green'),lty=c(1,2,4)) #### calculate wi aic weight of evidence res<-matrix(numeric(100*2),nrow=100) for(i in 1:100){ aic_x=2*(6-mig$MaxEstLhood[i]) aic_y=2*(6-nomig$MaxEstLhood[i]) delta_aic_x=aic_x-min(aic_x,aic_y) delta_aic_y=aic_y-min(aic_x,aic_y) wi_x=exp(-0.5*delta_aic_x) wi_y=exp(-0.5*delta_aic_y) res[i,]=c(wi_x/(wi_x+wi_y),wi_y/(wi_x+wi_y)) } ######################### simulation for omorika, obovata and pabies perl count_alleles.pl allele_counts/Kirov/ > Kirov-allele-freq.txt perl count_alleles.pl allele_counts/CH-DE-DK/ > ALP-allele-freq.txt perl count_alleles.pl allele_counts/P.omorika/ > Omorika-allele-freq.txt perl count_alleles.pl allele_counts/RO/ > RO-allele-freq.txt perl count_alleles.pl allele_counts/SE-FI/ > SE9-allele-freq.txt perl count_alleles.pl allele_counts/P.obovata/ > Obovata-allele-freq.txt omk<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/Omorika-allele-freq.txt",header=T) obv<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/Obovata-allele-freq.txt",header=T) swe<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/SE9-allele-freq.txt",header=T) kirov<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/Kirov-allele-freq.txt",header=T) rom<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/RO-allele-freq.txt",header=T) alp<-read.table("~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/ALP-allele-freq.txt",header=T) omk$snp<-paste(omk$Contig, omk$Pos,sep=':') obv$snp<-paste(obv$Contig, obv$Pos,sep=':') swe$snp<-paste(swe$Contig, swe$Pos,sep=':') rom$snp<-paste(rom$Contig, rom$Pos,sep=':') kirov$snp<-paste(kirov$Contig, kirov$Pos,sep=':') alp$snp<-paste(alp$Contig, alp$Pos,sep=':') names(omk)[3:6]<-paste(names(omk)[3:6],1,sep='') names(obv)[3:6]<-paste(names(obv)[3:6],2,sep='') names(swe)[3:6]<-paste(names(swe)[3:6],3,sep='') names(kirov)[3:6]<-paste(names(kirov)[3:6],4,sep='') names(alp)[3:6]<-paste(names(alp)[3:6],5,sep='') names(rom)[3:6]<-paste(names(rom)[3:6],6,sep='') data<-merge(omk,obv[,3:7],by='snp') data<-merge(data,swe[,3:7],by='snp') data<-merge(data,kirov[,3:7],by='snp') data<-merge(data,alp[,3:7],by='snp') data<-merge(data,rom[,3:7],by='snp') data$MAI<-apply(data[,4:27],1,find_minor_allele,n.pop=6) data$omk<-apply(data[,c(4:7,28)],1,function(a){a[floor(a[5])]}) data$obv<-apply(data[,c(8:11,28)],1,function(a){a[floor(a[5])]}) data$swe<-apply(data[,c(12:15,28)],1,function(a){a[floor(a[5])]}) data$kirov<-apply(data[,c(16:19,28)],1,function(a){a[floor(a[5])]}) data$alp<-apply(data[,c(20:23,28)],1,function(a){a[floor(a[5])]}) data$rom<-apply(data[,c(24:27,28)],1,function(a){a[floor(a[5])]}) data$omk.sub<-apply(data[,c(4:7,29)],1,subsample_maf) # subsample MAF with size = 10 for all populations data$obv.sub<-apply(data[,c(8:11,30)],1,subsample_maf) data$swe.sub<-apply(data[,c(12:15,31)],1,subsample_maf) data$kirov.sub<-apply(data[,c(16:19,32)],1,subsample_maf) # subsample MAF with size = 10 for all populations data$alp.sub<-apply(data[,c(20:23,33)],1,subsample_maf) data$rom.sub<-apply(data[,c(24:27,34)],1,subsample_maf) write.table(data,file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/subsampled-n10-combined-MAF-omk-obv-swe-hyb-alp-rom.txt',quote=F,row.names=F,sep="\t",) sfs_2d<-get_jointMAF_SFS(data[,35:40],amb.index=data$MAI!=floor(data$MAI),size=c(10,10,10,10,10,10)) write.table(sfs_2d[[1]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop1_0.obs',quote=F,sep="\t") write.table(sfs_2d[[2]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop2_0.obs',quote=F,sep="\t") write.table(sfs_2d[[3]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop3_0.obs',quote=F,sep="\t") write.table(sfs_2d[[4]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop4_0.obs',quote=F,sep="\t") write.table(sfs_2d[[5]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop5_0.obs',quote=F,sep="\t") write.table(sfs_2d[[6]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop2_1.obs',quote=F,sep="\t") write.table(sfs_2d[[7]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop3_1.obs',quote=F,sep="\t") write.table(sfs_2d[[8]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop4_1.obs',quote=F,sep="\t") write.table(sfs_2d[[9]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop5_1.obs',quote=F,sep="\t") write.table(sfs_2d[[10]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop3_2.obs',quote=F,sep="\t") write.table(sfs_2d[[11]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop4_2.obs',quote=F,sep="\t") write.table(sfs_2d[[12]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop5_2.obs',quote=F,sep="\t") write.table(sfs_2d[[13]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop4_3.obs',quote=F,sep="\t") write.table(sfs_2d[[14]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop5_3.obs',quote=F,sep="\t") write.table(sfs_2d[[15]],file='~/Desktop/exomeCapture/simulations/input/obs_sfs/sfs/omk-obv-swe-hyb-alp-rom/all6_jointMAFpop5_4.obs',quote=F,sep="\t") ## run simulations fsc25221 -t omk-obv-swe.tpl -e omk-obv-swe.est -M0.001 -n100000 -N100000 -l10 -L20 -m -q fsc25221 -t all6.tpl -e all6.est -M0.001 -n100000 -N100000 -l10 -L20 -m -q # plot source('~/Desktop/exomeCapture/simulations/functions.r') obs.files<-list.files('~/Desktop/exomeCapture/simulations/fsc_mac64/models/omk-obv-swe',full=T,pattern='.obs$') model.files<-list.files('~/Desktop/exomeCapture/simulations/fsc_mac64/models/omk-obv-swe/all6/',pattern='.txt$',full=T) model<-list() obs<-list() for (i in 1:15){ obs[[i]]<-as.matrix(read.table(obs.files[i],header=T,skip=1)) model[[i]]<-as.matrix(read.table(model.files[i],header=T)) } p1<-recordPlot() p2<-recordPlot() p3<-recordPlot() res1<-eval_2d_sfs(model[[1]],obs[[1]],las=1,xlab='P.omorika',ylab='P.obovata',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res2<-eval_2d_sfs(model[[2]],obs[[2]],las=1,xlab='P.omorika',ylab='P.abes (FENNO)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res3<-eval_2d_sfs(model[[3]],obs[[3]],las=1,xlab='P.obovata',ylab='P.abes (FENNO)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res4<-eval_2d_sfs(model[[4]],obs[[4]],las=1,xlab='P.omorika',ylab='Hybrid',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res5<-eval_2d_sfs(model[[5]],obs[[5]],las=1,xlab='P.obovata',ylab='Hybrid',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res6<-eval_2d_sfs(model[[6]],obs[[6]],las=1,xlab='P.abies (FENNO)',ylab='Hybrid',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res7<-eval_2d_sfs(model[[7]],obs[[7]],las=1,xlab='P.omorika',ylab='P.abies (ALP)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res8<-eval_2d_sfs(model[[8]],obs[[8]],las=1,xlab='P.obovata',ylab='P.abies (ALP)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res9<-eval_2d_sfs(model[[9]],obs[[9]],las=1,xlab='P.abies (FENNO)',ylab='P.abies (ALP)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res10<-eval_2d_sfs(model[[10]],obs[[10]],las=1,xlab='Hybrid',ylab='P.abies (ALP)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res11<-eval_2d_sfs(model[[11]],obs[[11]],las=1,xlab='P.omorika',ylab='P.abies (ROM)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res12<-eval_2d_sfs(model[[12]],obs[[12]],las=1,xlab='P.obovata',ylab='P.abies (ROM)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res13<-eval_2d_sfs(model[[13]],obs[[13]],las=1,xlab='P.abies (FENNO)',ylab='P.abies (ROM)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res14<-eval_2d_sfs(model[[14]],obs[[14]],las=1,xlab='Hyrbid',ylab='P.abies (ROM)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T) res15<-eval_2d_sfs(model[[15]],obs[[15]],las=1,xlab='P.abies (ALP)',ylab='P.abies (ROM)',breaks=seq(-8.2,-0.45,length=20),mask_mono=T,max.freq=0.7,maf=T)