totajuliusd / topr

topr is a collection of plotting functions for visualizing and exploring genetic association results. Association results from multiple phenotypes can be viewed simultaneously, over the entire genome (Manhattan plot) or in the more detailed regional view.
Other
45 stars 13 forks source link

Multiple plots in a for loop #34

Closed hlnicholls closed 8 months ago

hlnicholls commented 8 months ago

Hi, thank you for developing this package, it's been very helpful and intuitive.

I've got multiple GWAS' that I've tried to make manhattan plots for in a for loop. The for loop runs without error but doesn't actually produce/save any plots. If I take my code out of the for loop and run each phenotype individually that works fine.

Is there something wrong with my for loop for using topr? Here is my code:

phenotypes <- c("phenotype1", "phenotype2", "phenotype3")

for (phenotype in phenotypes) {
  res <- filter(res_all, Phenotype == phenotype)
  gwas_filtered <- res[res$P < 5e-8, ]
  print('get genes')
  gwas_annotated <- annotate_with_nearest_gene(gwas_filtered)
  res[res$P < 5e-8, "Gene_Symbol"] <- gwas_annotated$Gene_Symbol
  print('plotting...')
  png_filename <- paste0("/user/plots/", phenotype, "_annotated.png")
  png(png_filename, width = 3000, height = 1800, res = 300)
  topr::manhattan(res, annotate = 5e-9, title=phenotype)
  dev.off()
}

The above runs without error but doesn't save any plots. When I try taking the same code out of the for loop my plot saves to a png file as expected:

phenotype <- "phenotype1"
res <- filter(res_all, Phenotype == phenotype)
gwas_filtered <- res[res$P < 5e-8, ]
print('get genes')
gwas_annotated <- annotate_with_nearest_gene(gwas_filtered)
res[res$P < 5e-8, "Gene_Symbol"] <- gwas_annotated$Gene_Symbol
print('plotting...')
png_filename <- paste0("/user/plots/", phenotype, "_annotated.png")
png(png_filename, width = 3000, height = 1800, res = 300)
topr::manhattan(res, annotate = 5e-9, title=phenotype)
dev.off()
totajuliusd commented 8 months ago

Hi, try putting your manhattan call within the plot() function. This works with the test data:

phenotypes <- c("CD_UKBB", "CD_FINNGEN")
for (phenotype in phenotypes){
  png_filename <- paste0( phenotype, "_annotated.png")
  png(png_filename, width = 3000, height = 1800, res = 300)
  plot(topr::manhattan(get(phenotype), annotate = 5e-9, title=phenotype))
  dev.off()
}
hlnicholls commented 8 months ago

This has worked, thank you!