Pan-Genome R Notebook

library("ape")
library("vegan")

Loading required package: permute
Loading required package: lattice
This is vegan 2.4-6

library("philentropy")
library("ggplot2")
setwd("~/Desktop/Intersection_OMCL/")
pan_matrix <- read.table(file = "pangenome_matrix_t0.tab",stringsAsFactors = TRUE, header = TRUE, row.names = 1)
sp <- specaccum(pan_matrix,method="random",permutations=100)

Set of permutations < 'minperm'. Generating entire set.

plot(sp,ci.type="poly",col="blue",lwd=2,ci.lty=0,ci.col="lightblue",xlab="Genomes",ylab="Proteins",main="Gene accumulation plot")

# Rarefaction curve with box and whiskers for the resamplings
plot(sp,ci.type="poly",col="blue",lwd=2,ci.lty=0,ci.col="lightblue",xlab="Genomes",ylab="Proteins",main="Gene accumulation plot")
boxplot(sp,col="yellow",add=TRUE)

#to plot to file
jpeg('rarefaction_plot.jpg') 
plot(sp,ci.type="poly",col="blue",lwd=2,ci.lty=0,ci.col="lightblue",xlab="Genomes",ylab="Proteins",main="Gene accumulation plot")
boxplot(sp,col="yellow",add=TRUE)
#dev.off()

# Make a histogram based on gene presence absence
# Version one very awkward
n_row <- nrow(pan_matrix)   #count number of rows
n_col <- ncol(pan_matrix)   #Ditto columns
HG <- vector(mode="numeric", length = 0)
 
 for(col in 1:n_col){
     col_count = 0
     for(row in 1:n_row){
         if(pan_matrix[row,col] > 0){
             col_count = col_count + 1
         }
     }
     HG[col] <- col_count
 }
 t_HG <- t(HG) 
 hist(HG,col="red",xlab="Present in # Genomes",ylab="# HGs",main="Histogram of HG Frequenies")

NA

# Make a histogram based on gene presence absence
# Version two (after googling R and histogram)
qplot(HG, 
       geom="histogram",
      binwidth = 0.5,
      main = "Histogram of HG Frequenies",
      xlab="Present in # Genomes",
      ylab="# HGs",
      col=I("blue"), 
      fill=I("red"),
      alpha=I(.6)
       )  

To plot this to file

#to plot to file
jpeg('histogram_plot2.jpg') 
qplot(HG, 
       geom="histogram",
      binwidth = 0.5,
      main = "Histogram of HG Frequenies",
      xlab="Present in # Genomes",
      ylab="# HGs",
      col=I("blue"), 
      fill=I("red"),
      alpha=I(.6)
       )  
#dev.off()

#make a distance matrix tree based on gene presence absence
r_names <- row.names(pan_matrix)
matrix_dist <- distance(pan_matrix, method = "jaccard")
row.names(matrix_dist) <- r_names
tree_bionj <- bionj(matrix_dist)
plot(tree_bionj)

Using gene presence absence create a PCoA plot

pcoa <- wcmdscale(matrix_dist)
pcoa_names <- row.names(pcoa)
plot(pcoa,pch=1,cex=2, col=c("red"))
text(pcoa, pcoa_names,cex=0.9)

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.