Mercurial > repos > dereeper > roary_plots

comparison Roary/bin/roary-create_pan_genome_plots.R @ 0:c47a5f61bc9f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Fri, 14 May 2021 20:27:06 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c47a5f61bc9f
1 #!/usr/bin/env Rscript
2 # ABSTRACT: Create R plots
3 # PODNAME: create_plots.R
4 # Take the output files from the pan genome pipeline and create nice plots.
5 library(ggplot2)
6
7
8 mydata = read.table("number_of_new_genes.Rtab")
9 boxplot(mydata, data=mydata, main="Number of new genes",
10 xlab="No. of genomes", ylab="No. of genes",varwidth=TRUE, ylim=c(0,max(mydata)), outline=FALSE)
11
12 mydata = read.table("number_of_conserved_genes.Rtab")
13 boxplot(mydata, data=mydata, main="Number of conserved genes",
14 xlab="No. of genomes", ylab="No. of genes",varwidth=TRUE, ylim=c(0,max(mydata)), outline=FALSE)
15
16 mydata = read.table("number_of_genes_in_pan_genome.Rtab")
17 boxplot(mydata, data=mydata, main="No. of genes in the pan-genome",
18 xlab="No. of genomes", ylab="No. of genes",varwidth=TRUE, ylim=c(0,max(mydata)), outline=FALSE)
19
20 mydata = read.table("number_of_unique_genes.Rtab")
21 boxplot(mydata, data=mydata, main="Number of unique genes",
22 xlab="No. of genomes", ylab="No. of genes",varwidth=TRUE, ylim=c(0,max(mydata)), outline=FALSE)
23
24 mydata = read.table("blast_identity_frequency.Rtab")
25 plot(mydata,main="Number of blastp hits with different percentage identity", xlab="Blast percentage identity", ylab="No. blast results")
26
27
28 library(ggplot2)
29 conserved = colMeans(read.table("number_of_conserved_genes.Rtab"))
30 total = colMeans(read.table("number_of_genes_in_pan_genome.Rtab"))
31
32 genes = data.frame( genes_to_genomes = c(conserved,total),
33 genomes = c(c(1:length(conserved)),c(1:length(conserved))),
34 Key = c(rep("Conserved genes",length(conserved)), rep("Total genes",length(total))) )
35
36 ggplot(data = genes, aes(x = genomes, y = genes_to_genomes, group = Key, linetype=Key)) +geom_line()+
37 theme_classic() +
38 ylim(c(1,max(total)))+
39 xlim(c(1,length(total)))+
40 xlab("No. of genomes") +
41 ylab("No. of genes")+ theme_bw(base_size = 16) + theme(legend.justification=c(0,1),legend.position=c(0,1))+
42 ggsave(filename="conserved_vs_total_genes.png", scale=1)
43
44 ######################
45
46 unique_genes = colMeans(read.table("number_of_unique_genes.Rtab"))
47 new_genes = colMeans(read.table("number_of_new_genes.Rtab"))
48
49 genes = data.frame( genes_to_genomes = c(unique_genes,new_genes),
50 genomes = c(c(1:length(unique_genes)),c(1:length(unique_genes))),
51 Key = c(rep("Unique genes",length(unique_genes)), rep("New genes",length(new_genes))) )
52
53 ggplot(data = genes, aes(x = genomes, y = genes_to_genomes, group = Key, linetype=Key)) +geom_line()+
54 theme_classic() +
55 ylim(c(1,max(unique_genes)))+
56 xlim(c(1,length(unique_genes)))+
57 xlab("No. of genomes") +
58 ylab("No. of genes")+ theme_bw(base_size = 16) + theme(legend.justification=c(1,1),legend.position=c(1,1))+
59 ggsave(filename="unique_vs_new_genes.png", scale=1)