Deal with cases where there less than two pops for Fst

tools/Fst_gds/Fst_gds.script.R

@@ -2,6 +2,7 @@ library(SNPRelate);
 library(SeqArray);
 library(ape)
 library(ggtree)
+library(ggplot2)
 
 popmap_file = snakemake@input[["popmap_file"]]
 gdsFile = snakemake@input[["gds"]]
@@ -34,49 +35,57 @@ if (popmap_file != '')
     if (nbpops > 1)
     {
 
-    fstMatrix     = matrix(0, nbpops, nbpops)  
-    meanfstMatrix = matrix(0, nbpops, nbpops)    
-    
-    for (i in seq(1,nbpops-1) )
-    {
-        # Two populations: 
-        for (j in seq(i+1, nbpops) )
+        fstMatrix     = matrix(0, nbpops, nbpops)  
+        meanfstMatrix = matrix(0, nbpops, nbpops)    
+        
+        for (i in seq(1,nbpops-1) )
         {
-            selected <- tab[tab$popcode %in% c(Allpops[i], Allpops[j]),]
-            v <- snpgdsFst(gds, sample.id=selected$sample.id, population=as.factor(selected$popcode), method="W&C84", autosome.only =F)
-            fstMatrix[i,j] = v$Fst
-            fstMatrix[j,i] = v$Fst
-            meanfstMatrix[i,j] = v$MeanFst
-            meanfstMatrix[j,i] = v$MeanFst
+            # Two populations: 
+            for (j in seq(i+1, nbpops) )
+            {
+                selected <- tab[tab$popcode %in% c(Allpops[i], Allpops[j]),]
+                v <- snpgdsFst(gds, sample.id=selected$sample.id, population=as.factor(selected$popcode), method="W&C84", autosome.only =F)
+                fstMatrix[i,j] = v$Fst
+                fstMatrix[j,i] = v$Fst
+                meanfstMatrix[i,j] = v$MeanFst
+                meanfstMatrix[j,i] = v$MeanFst
+            }
         }
-    }
-    colnames(meanfstMatrix) = Allpops;    rownames(meanfstMatrix) = Allpops
-    colnames(fstMatrix) = Allpops;    rownames(fstMatrix) = Allpops
+        colnames(meanfstMatrix) = Allpops;    rownames(meanfstMatrix) = Allpops
+        colnames(fstMatrix) = Allpops;    rownames(fstMatrix) = Allpops
+
+        write.table(meanfstMatrix, file = meanfstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
+        write.table(fstMatrix, file = fstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
+                
+        png(file=Fst_heatmap_png, width=1024, height=1024); 
+        heatmap(fstMatrix) ;
+                
+        dev.off();
 
-    write.table(meanfstMatrix, file = meanfstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
-    write.table(fstMatrix, file = fstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
-            
+        png(file=NJ_png, width=1024, height=1024);
+        forme   = "rectangular" 
+        title=paste0("NJ plot on pairewise Fst.")
+        
+        distance = fstMatrix
+        tree = nj(distance)
+        p <- ggtree(tree,layout=forme) + geom_text(aes(label=label), size=3, color="purple", hjust=-0.3) + theme_tree2("gray86") +ggtitle(title)
+        plot(p)
+        dev.off()
+    }   
+} else {
+    # need popmap with at least 2 pops
+    write.table(0, file = meanfstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
+    write.table(0, file = fstFile, append = FALSE, quote = F, sep = '\t', row.names = T, col.names = NA);
     png(file=Fst_heatmap_png, width=1024, height=1024); 
-    heatmap(fstMatrix) ;
-             
-    dev.off();
+        p = ggplot() + ggtitle("You have to provide a popmap with samples from at least 2 populations") + theme(panel.border = element_rect(colour = "black", fill=NA, size=2))
+        print(p)
+    dev.off()
 
-    png(file=NJ_png, width=1024, height=1024);
-    forme   = "rectangular" 
-	title=paste0("NJ plot on pairewise Fst.")
-	
-    distance = fstMatrix
-	tree = nj(distance)
-	p <- ggtree(tree,layout=forme) + geom_text(aes(label=label), size=3, color="purple", hjust=-0.3) + theme_tree2("gray86") +ggtitle(title)
-	plot(p)
+    png(file=NJ_png, width=1024, height=1024); 
+        p = ggplot() + ggtitle("You have to provide a popmap with samples from at least 2 populations") + theme(panel.border = element_rect(colour = "black", fill=NA, size=2))
+        print(p)
     dev.off()
-    }
-    else {
-        # need more than 1 pop
-    } 
-}
-else{
-    # need popmap
+
 }
 seqClose(gds);