R script for computation of G values from multiple simulation files

R script for computation of G values from multiple simulation files

This script embeds the functions fromArpToSFS() and G2DcalcMultipop() (see the corresponding files) and allows one compute G values from multiple .arp files, for individual contigs and from one pair of populations. All scripts are provided in their current state of development with no warranty that they will work in all cases, and only for reference purposes. They can be used and modified by citing the Author (see script header). Example files can be obtained from the Author upon request.

#this script was written by Ivan Scotti, INRA, URFM
#Site Agroparc, Domaine Saint-Paul
#84914 AVIGNON Cedex 9
#you are free to use and modify it, provided you quote the author
#for questions, please contact Ivan Scotti at
#ivan.scotti[at]paca.inra.fr
#15 July 2015
#
#embedding the fromArpToSFS() and G2DcalcMultiPop() functions to intake
#multiple .arp files
#useful when multiple simulations are produced from the same parameter set
#recommended by Niealsen eet al. Genome res. 19, 838 (2009)
#example arguments
#root="crops113_maxLmodif_1_" #1.arp
#NbSims=100
#nChroms=104
#sampleSizes=5
multiG2DcalcMultiPop<-function(root,nPops,NbSims,nChroms,sampleSizes,Pop1,Pop2)
{
  #vector to stock all values
  multiGval<-numeric()
  #getting .arp file names
  fileNames<-paste(root,1:NbSim,".arp",sep="")
  for (i in fileNames)
  {
    #embedding the two functions
    fromArpToSFS(infile=i,nChrom=nChroms,nPop=nPops,sampleSize=sampleSizes)
    #NB G2Dcalc() uses defaults because it
    #inherits values from fromArpToSFS()
    multiGval<-c(multiGval,G2DcalcMultiPop(,firstPop=Pop1,secondPop=Pop2,sampleSize=sampleSizes))
  }
  #output sent to connection
  multiGval  
}
# #usage:
# wapaSimGvaluesPops1and2<-multiG2DcalcMultiPop(root="waSNP10_maxLmodif_1_",
#                                               nPops=4,
#                                               NbSims=100,
#                                               nChroms=2883,
#                                               sampleSizes=40,
#                                               Pop1=1,
#                                               Pop2=2)
# hist(wapaSimGvaluesPops1and2)

Date de modification : 22 juin 2023 | Date de création : 16 juillet 2015 | Rédaction : Ivan Scotti