Changes to R script

Scripts/ROHcircle.R

@@ -46,10 +46,6 @@ margin <- 10000000
 # Function to derive global position of locations
 glob_pos <- function(locations, chromosomes, margin) {
   
-
-  
-  chromosomes[chromosomes == "23"] <- "19"
-  chromosomes[chromosomes == "X"] <- "19"
   
   locations_chromosomes <- data.frame(
     loc=locations,
@@ -81,53 +77,52 @@ bcftools <- bcftools[bcftools$Quality > 95,]
 
 
 
-for (ind in samples){
-  het <- read.table(paste0("het_", ind), header=T)
-  het <- het[!het$PI == "PI",]
-  het$PI <- as.numeric(het$PI)
+
+het <- read.table(paste0("het_", ind), header=T)
+het <- het[!het$PI == "PI",]
+het$PI <- as.numeric(het$PI)
   
-  bcftools.ind <- bcftools[bcftools$Sample == ind,]
+bcftools.ind <- bcftools[bcftools$Sample == ind,]
 
-  plink.ind <- plink[plink$IID == ind,]
-  #plink.ind <- plink[plink$KB > 500,]
+plink.ind <- plink[plink$IID == ind,]
+#plink.ind <- plink[plink$KB > 500,]
 
-  plt <- ggplot() +
-    geom_segment(data=het, mapping=aes(x=glob_pos(BIN_START, CHROM), y=0, xend=glob_pos(BIN_START, CHROM), yend=PI)) +
-    geom_point(aes(x=end_position, y=0), size=5, alpha=0) +
-    coord_polar(start = 0)
+plt <- ggplot() +
+  geom_segment(data=het, mapping=aes(x=glob_pos(BIN_START, CHROM), y=0, xend=glob_pos(BIN_START, CHROM), yend=PI)) +
+  geom_point(aes(x=end_position, y=0), size=5, alpha=0) +
+  coord_polar(start = 0)
   
-  if (dim(plink.ind)[1] >= 1) {
-    plt <- plt + 
-      geom_segment(data=plink.ind, mapping=aes(x=glob_pos(POS1, CHR), y=-0.002, xend=glob_pos(POS2, CHR), yend=-0.002), color="blue", size=3, lineend="butt")
+if (dim(plink.ind)[1] >= 1) {
+  plt <- plt + 
+    geom_segment(data=plink.ind, mapping=aes(x=glob_pos(POS1, CHR), y=-0.002, xend=glob_pos(POS2, CHR), yend=-0.002), color="blue", size=3, lineend="butt")
       
-    }
-  if (dim(bcftools.ind)[1] >= 1) {
-    plt <- plt + 
-      geom_segment(data=bcftools.ind, mapping=aes(x=glob_pos(Start, Chromosome), y=-0.004, xend=glob_pos(End, Chromosome), yend=-0.004), color="red", size=3)
-    nrow(bcftools.ind)
-    sum(bcftools.ind$Length)
-    sum(bcftools.ind$Length) / nrow(bcftools.ind)
-
   }
-  
-  plt <- plt +
-    annotate("text", x=0, y=-0.1, vjust=-0.5, label=ind, size=20,) +
-    annotate("text", x=0, y=-0.1, label="plink", vjust=1, hjust=0, size=15, color="blue") +
-    annotate("text", x=0, y=-0.1, label="bcftools", vjust=1,hjust=1, size=15, color="red") +
-    #annotate("text", x=0, y=-0.1, label=sum(bcftools.ind$Length), vjust=3, hjust=0, size=10) +
-    #annotate("text", x=0, y=-0.1, label=sum(bcftools.ind$Length), vjust=3,hjust=1, size=10)
-  
-  png(file=paste0("het_roh_", ind, ".png"), units="cm", width=60, height=60, res=600) 
-  print(
-    plt + theme_nothing() + 
-      theme(plot.margin = unit(rep(0, 4), "null")) +
-      scale_x_continuous(expand=c(0,0)) +
-      scale_y_continuous(limits=c(-0.1, 0.048), expand=c(0,0)) +
-      labs(x = NULL, y = NULL)
+if (dim(bcftools.ind)[1] >= 1) {
+  plt <- plt + 
+    geom_segment(data=bcftools.ind, mapping=aes(x=glob_pos(Start, Chromosome), y=-0.004, xend=glob_pos(End, Chromosome), yend=-0.004), color="red", size=3)
+  nrow(bcftools.ind)
+  sum(bcftools.ind$Length)
+  sum(bcftools.ind$Length) / nrow(bcftools.ind)
 
-  )
-  dev.off()
 }
+  
+plt <- plt +
+  annotate("text", x=0, y=-0.1, vjust=-0.5, label=ind, size=20,) +
+  annotate("text", x=0, y=-0.1, label="plink", vjust=1, hjust=0, size=15, color="blue") +
+  annotate("text", x=0, y=-0.1, label="bcftools", vjust=1,hjust=1, size=15, color="red") +
+  #annotate("text", x=0, y=-0.1, label=sum(bcftools.ind$Length), vjust=3, hjust=0, size=10) +
+  #annotate("text", x=0, y=-0.1, label=sum(bcftools.ind$Length), vjust=3,hjust=1, size=10)
+
+png(file=paste0("het_roh_", ind, ".png"), units="cm", width=60, height=60, res=600) 
+print(
+  plt + theme_nothing() + 
+    theme(plot.margin = unit(rep(0, 4), "null")) +
+    scale_x_continuous(expand=c(0,0)) +
+    scale_y_continuous(limits=c(-0.1, 0.048), expand=c(0,0)) +
+    labs(x = NULL, y = NULL)
+  )
+dev.off()
+