Hi all,

I used Arlequin software to detect loci under selection in two populations. Actually, it worked well and it detected these loci. Now I want to plot these loci in a plot by using R. They provided R function, lociSelection.r, that could be used to generate the required plot (Page number 113 in Arlequin manual; http://cmpg.unibe.ch/software/arlequin35/man/Arlequin35.pdf). I tried to use this function but it did not work well with my data.

I do not know what are outfile, lociSelection, fStat, timeAttr and .png in the function down. I have two Arelquin output files, fdist2_ObsOut.txt and fdist2_simOut.txt.

I do not how can I adjust this function to fit my data. Any help well be appreciated.

Here is the function

################################################################################
# lociSelection function
#
# Author: Heidi Lischer
# Date: 8.2009
################################################################################

lociSelection <- function(pValList, FStatCiList, timeAttr, outfile, fStat){
   if(fStat == "FST"){
    fst <- 3
    fstP <- 4
  } else {
    fst <- 5
    fstP <- 6
  }
   MaxHet <- max(pValList[,2])
  MaxF <- max(pValList[,fst])
   MinHet <- 0
  MinF <- 0
   if(min(pValList[,2]) < 0){
    MinHet <- min(pValList[,2])
  }
  if(min(pValList[,fst]) < 0){
    MinF <- min(pValList[,fst])
  }            
   FST_h0.05 <- list()
  FST_0.05 <- list()
  FST_0.01 <- list()
  lociNames_0.01 <- list()
  a <- 1
  b <- 1
  c <- 1
  for(i in 1:nrow(pValList)){
    if(pValList[i, fstP] > 0.05){
      FST_h0.05[[a]] <- pValList[i, c(2, fst)]
      a <- a + 1   
    } else {
      if(pValList[i, fstP] > 0.01){
        FST_0.05[[b]] <- pValList[i, c(2, fst)]
        b <- b + 1
      } else{
        FST_0.01[[c]] <- pValList[i, c(2, fst)]
        lociNames_0.01[[c]] <- pValList[i, 1]
        c <- c + 1
      }
    }       
  }
  FST_h0.05 <- t(as.matrix(as.data.frame(FST_h0.05)))
  FST_0.05 <- t(as.matrix(as.data.frame(FST_0.05)))
  FST_0.01 <- t(as.matrix(as.data.frame(FST_0.01)))
  lociNames_0.01 <- as.vector(as.character(lociNames_0.01))

    # graphic---------------------------------------------------------------------
  outfileGraphic <- paste(**outfile, "lociSelection_", fStat," ", timeAttr,".png**",
                          sep="")
   #save graphic
  png(outfileGraphic.png, width = 1300, height = 1300, res=144)

  plot(FST_h0.05, cex=0.8, xlab="Heterozygosity", ylab="", main="",
       xlim=c(MinHet, MaxHet), ylim=c(MinF, MaxF), cex.axis=0.9)

  if(fStat == "FST"){
    mtext(text=expression(bold(Detection~of~loci~under~selection~from~genome~scans~based~on~F[ST])),
          side=3, line=1.5, cex=1.2)
    mtext(text=expression(F[ST]), side=2, line=2.5, cex=1.2)
  } else {
    mtext(text=expression(bold(Detection~of~loci~under~selection~from~genome~scans~based~on~F[CT])),
          side=3, line=1.5, cex=1.2)
    mtext(text=expression(F[CT]), side=2, line=2.5, cex=1.2)
  }
   if(length(FST_0.05) != 0){
    points(FST_0.05, pch=16, cex=0.8, col="blue")
  }
  if(length(FST_0.01) != 0){
    points(FST_0.01, pch=16, cex=0.8, col="red3")
  }
   lines(FStatCiList[,1], FStatCiList[,5])
   lines(FStatCiList[,1], FStatCiList[,4], lty=5)
  lines(FStatCiList[,1], FStatCiList[,6], lty=5)
   lines(FStatCiList[,1], FStatCiList[,3], lty=2, col="blue")
  lines(FStatCiList[,1], FStatCiList[,7], lty=2, col="blue")
   lines(FStatCiList[,1], FStatCiList[,2], lty=3, col="red3")
  lines(FStatCiList[,1], FStatCiList[,8], lty=3, col="red3")
   legend("topleft",
         c("", "", "", "50% quantile", "10% quantile", "5% quantile", "1% quantile"),
         lty=c(1, 1, 1, 1, 5, 2, 3),
         col=c("white", "white", "white", "black", "black", "blue", "red3"),
         cex=0.8,  bty="n")

  legend("topleft", legend=c("p > 5%", "p <= 5%", "p <= 1%"),
         pch=c(1, 16, 16), col=c("black", "blue", "red3"), cex=0.8, bty="n")

  dev.off()
    # graphic---------------------------------------------------------------------
  outfileGraphic <- paste(outfile, "lociSelection_Names_", fStat," ", timeAttr,".png",
                          sep="")
    #save graphic
  png(outfileGraphic, width = 1300, height = 1300, res=144)