.Rhistory

# assign potential group types
group_types <- c('groupA', 'groupB', 'groupC', 'groupD', 'groupE',
'groupF', 'groupG', 'groupH', 'groupI', 'groupJ', 'groupK',
'groupL', 'groupM', 'groupN', 'groupO')
# assign probabilites for each group type. For now give each group the same probability of occurring
prob_group_types <- c(0.0357, 0.0357, 0.0357, 0.0357, 0.0357, 0.0357, 0.0357, 0.0357, 0.0357, 0.0357,
0.0357, 0.0357, 0.0357, 0.0357, 0.0357)
# create data frame
comb_data <- data.frame(group_types, prob_group_types)
comb_data
?sample
dateString <- "2010/05/01"
date <- as.Date(dateString, format="%Y/%m/%d")
year <- format(date, "%Y")
month <- format(date, "%b")
year
month
path <- "/home/josephcrispell/Desktop/Research/EdgeArea_UK/Fastqs_06-19_n68/FASTQC/"
# Read in the FASTQC file summary table
summaryFile <- paste0(path, "Fastqc_summary_30-10-19.txt")
summary <- read.table(summaryFile, header=TRUE, sep="\t", stringsAsFactors=FALSE)
# Sort the file names
summary <- summary[order(summary$FileName), ]
# Look at GC distribution
hist(summary$GC, breaks=20, las=1)
summar[summar$GC == 68, ]
summary[summary$GC == 68, ]
summary[summary$GC == 68, c("FileName", "GC")]
summary[summary$GC != 65, c("FileName", "GC")]
boxplot(summary$LeftTrim, summary$RightTrim, names=c("LEFT", "RIGHT"), las=1, frame=FALSE,
main="Trimming suggestions", outcol=rgb(0,0,0, 0.1), pch=19)
# Check adapter flag
table(summary$AdapterContentFlag)
hist(summary$NumberReads, breaks=100, las=1)
summary[which(summary$NumberReads < 250000), c("FileName", "NumberReads")]
summary[summary$NumberReads > 5000000, c("FileName", "NumberReads")]
ui <- fluidPage(
# App title ----
titlePanel("Hello World!"),
# Sidebar layout with input and output definitions ----
sidebarLayout(
# Sidebar panel for inputs ----
sidebarPanel(
# Input: Slider for the number of bins ----
sliderInput(inputId = "bins",
label = "Number of bins:",
min = 5,
max = 50,
value = 30)
),
# Main panel for displaying outputs ----
mainPanel(
# Output: Histogram ----
plotOutput(outputId = "distPlot")
)
)
)
# Define server logic required to draw a histogram ----
server <- function(input, output) {
# Histogram of the Old Faithful Geyser Data ----
# with requested number of bins
# This expression that generates a histogram is wrapped in a call
# to renderPlot to indicate that:
#
# 1. It is "reactive" and therefore should be automatically
#    re-executed when inputs (input$bins) change
# 2. Its output type is a plot
output$distPlot <- renderPlot({
x    <- faithful$waiting
bins <- seq(min(x), max(x), length.out = input$bins + 1)
hist(x, breaks = bins, col = "#75AADB", border = "orange",
xlab = "Waiting time to next eruption (in mins)",
main = "Histogram of waiting times")
})
}
shinyApp(ui = ui, server = server)
install.packages("shiny")
# Load libraries
library(shiny)
# Construct the user interface object - controls layout and appearance
ui <- fluidPage(
# App title ----
titlePanel("Hello World!"),
# Sidebar layout with input and output definitions ----
sidebarLayout(
# Sidebar panel for inputs ----
sidebarPanel(
# Input: Slider for the number of bins ----
sliderInput(inputId = "bins",
label = "Number of bins:",
min = 5,
max = 50,
value = 30)
),
# Main panel for displaying outputs ----
mainPanel(
# Output: Histogram ----
plotOutput(outputId = "distPlot")
)
)
)
# Define server function - instructiosn for the computer to build the app
server <- function(input, output) {
# Histogram of the Old Faithful Geyser Data ----
# with requested number of bins
# This expression that generates a histogram is wrapped in a call
# to renderPlot to indicate that:
#
# 1. It is "reactive" and therefore should be automatically
#    re-executed when inputs (input$bins) change
# 2. Its output type is a plot
output$distPlot <- renderPlot({
x    <- faithful$waiting
bins <- seq(min(x), max(x), length.out = input$bins + 1)
hist(x, breaks = bins, col = "#75AADB", border = "orange",
xlab = "Waiting time to next eruption (in mins)",
main = "Histogram of waiting times")
})
}
shinyApp(ui = ui, server = server)
runApp('Desktop/RShinyApp.R')
?fileInput
?read.dna
?ape::read.dna
runApp('Desktop/RShinyApp.R')
?fileInput
runApp('Desktop/RShinyApp.R')
runApp('Desktop/RShinyApp.R')
library(basicPlotteR)
runApp('Desktop/RShinyApp.R')
plotFASTA
librarY(basicPlotteR)
remove.packages("basicPlotteR")
devtools::install_github("JosephCrispell/basicPlotteR")
runApp('Desktop/RShinyApp.R')
plotSamplingLocations <- function(xCoords, yCoords, scaleSize=10000, xPad=0.75, yPad=0.1, ...){
# Plot the sampling locations
plot(x=xCoords, y=yCoords, ...)
# Note the size of the axes
axisLimits <- par()$usr
xAxisSize <- axisLimits[2] - axisLimits[1]
yAxisSize <- axisLimits[4] - axisLimits[3]
# Note the start and end of scale box
xStart <- axisLimits[1] + (xAxisSize*xPad)
xEnd <- xStart + scaleSize
yStart <- axisLimits[3] + (yAxisSize*yPad)
yEnd <- yStart + scaleSize
# Add a scale box
rect(xleft=xStart, ybottom=yStart, xright=xEnd, ytop=yEnd)
# Add label for scale box
label <- bquote(paste(.(scaleSize/1000), "km")^2)
text(x=xStart+(0.5*scaleSize), y=yStart, labels=label, pos=1)
}
getTipInfo <- function(tipLabels, sampleInfo, propNs){
# Set the row names of the sample info table
rownames(sampleInfo) <- sampleInfo$Sample.Ref
# Create the tipInfo table
tipInfo <- sampleInfo[tipLabels, c("Sample.Ref", "Mapx", "Mapy", "Eartag", "CPHH", "Breakdown.ID", "County", "Year", "TYPING.RESULT")]
rownames(tipInfo) <- tipLabels
# Add the proportion Ns to tipInfo table
rownames(propNs) <- parseTipLabels(propNs$ID)
tipInfo$ProportionNs <- propNs[tipLabels, "ProportionNs"]
return(tipInfo)
}
parseTipLabels <- function(tipLabels){
# Check if already parsed
if(grepl(tipLabels[1], pattern="_") == FALSE){
return(tipLabels)
}
# Initialise a vector to store the parsed labels
parsedLabels <- c()
# Examine eahc tip label
for(index in seq_along(tipLabels)){
# Remove the ">" from the start of the tip label
label <- tipLabels[index]
if(grepl(tipLabels[1], pattern="^>")){
label <- substr(tipLabels[index], 2, nchar(tipLabels[index]))
}
# Split the label and select first part
label <- strsplit(label, split="_")[[1]][1]
# If it doesn't start with "AF" change "-" to "/"
if(grepl(label, pattern="^AF") == FALSE){
label <- gsub(pattern="-", replacement="/", label)
}
# Store the parsed label
parsedLabels[index] <- label
}
return(parsedLabels)
}
calculatePropNsInSequence <- function(sequenceIndex, sequences){
# Count the number of Ns in the sequence
nucleotideCounts <- table(sequences[sequenceIndex, ])
names(nucleotideCounts) <- toupper(names(nucleotideCounts))
# Note the number of Ns
numberMissing <- 0
if("N" %in% names(nucleotideCounts)){
numberMissing <- nucleotideCounts[["N"]]
}
return(numberMissing / ncol(sequences))
}
getNSitesInFASTA <- function(fastaFile){
# Open a connection to a file to read (open="r")
connection <- file(fastaFile, open="r")
# Get first line of file
firstLine <- readLines(connection, n=1)
# Close file connection
close(connection)
# Get the number of sites used in the FASTA file from first line
nSites <- as.numeric(strsplit(firstLine, " ")[[1]][2])
return(nSites)
}
runRAXML <- function(fastaFile, date, path, nBootstraps=100, nThreads=10, outgroup=NULL, model="GTRCAT"){
# Note the directory for the RAxML output files
directory <- paste(path, "RAxML_", date, sep="")
# Check if analyses already run
alreadyRun <- dir.exists(directory)
# If not already run, create output directory for RAxML
if(alreadyRun == FALSE){
suppressWarnings(dir.create(directory))
}
# Set the Working directory - this will be where the output files are dumped
setwd(directory)
# Build analysis name
analysisName <- paste("RaxML-R_", date, sep="")
# Check if already Run and just want to retrieve tree
if(alreadyRun == FALSE){
# Build the command
seeds <- sample(1:100000000, size=2, replace=FALSE) # For parsimony tree and boostrapping
if(is.null(outgroup)){
command <- paste("raxmlHPC",
" -f a", # Algorithm: Rapid boostrap inference
" -N ", nBootstraps,
" -T ", nThreads,
" -m ", model, " -V", # -V means no rate heterogenity
" -p ", seeds[1], " -x ", seeds[2], # Parsimony and boostrapping seeds
" -n ", analysisName,
" -s ", fastaFile, sep="")
}else{
command <- paste("raxmlHPC",
" -f a", # Algorithm: Rapid boostrap inference
" -N ", nBootstraps,
" -T ", nThreads,
" -m ", model, " -V", # -V means no rate heterogenity
" -p ", seeds[1], " -x ", seeds[2], # Parsimony and boostrapping seeds
" -n ", analysisName,
" -s ", fastaFile,
" -o ", outgroup, sep="")
}
system(command, intern=TRUE)
}
# Get the tree and read it in
treeBS <- getTreeFileWithSupportValues(analysisName)
return(treeBS)
}
getTreeFileWithSupportValues <- function(analysisName){
# Get files in current working directory
files <- list.files()
# Select the tree file with BS support values
treeBSFile <- files[grepl(files, pattern=paste("RAxML_bipartitions[.]", analysisName, sep="")) == TRUE]
# Open the file
treeBS <- read.tree(treeBSFile)
return(treeBS)
}
# Load libraries
library(ape) # Reading and phylogeny
library(basicPlotteR) # Because I made it! :-)
# Set the path
path <- file.path("~", "Desktop", "Research", "EdgeArea_UK")
# Get the current date
date <- format(Sys.Date(), "%d-%m-%y")
#### Read in sequencing data ####
# Read in the sample metadata
sampleInfoFile <- file.path(path, "SampleInformation_28-05-19.csv")
sampleInfo <- read.table(sampleInfoFile, header=TRUE, sep=",", stringsAsFactors=FALSE)
sampleInfo <- sampleInfo[sampleInfo$Sample.Ref != "", ]
# Read in the FASTA file
fastaFile <- file.path(path, "vcfFiles", "sequences_Prox-10_31-10-2019.fasta")
sequences <- read.dna(fastaFile, format="fasta", as.character=TRUE)
nSites <- ncol(sequences)
# Plot the FASTA file
fastaPlotFile <- paste0(substr(fastaFile, 0, nchar(fastaFile) - 6), ".pdf")
plotFASTA(sequences, pdfFileName=fastaPlotFile, pdfHeight=14, pdfWidth=21, lineForSequenceNames=-3)
# Calculate sequencing quality in the FASTA file
# Note average proportion Ns may be higher because of the extreme coverage in one genome - reaching areas with no coverage in other genomes
propNs <- data.frame("ID"=rownames(sequences), "ProportionNs"=sapply(1:nrow(sequences), calculatePropNsInSequence, sequences),
stringsAsFactors=FALSE)
rownames(propNs) <- parseTipLabels(propNs$ID)
uninstall.packages("basicPlotteR")
remove.packages("basicPlotteR")
devtools::install_github("JosephCrispell/basicPlotteR")
# Load libraries
library(ape) # Reading and phylogeny
library(geiger) # For the tips function in used in defineBranchColoursOfClade
library(basicPlotteR) # Because I made it! :-)
# Set the path
path <- file.path("~", "Desktop", "Research", "EdgeArea_UK")
# Get the current date
date <- format(Sys.Date(), "%d-%m-%y")
# Read in the sample metadata
sampleInfoFile <- file.path(path, "SampleInformation_28-05-19.csv")
sampleInfo <- read.table(sampleInfoFile, header=TRUE, sep=",", stringsAsFactors=FALSE)
sampleInfo <- sampleInfo[sampleInfo$Sample.Ref != "", ]
# Read in the FASTA file
fastaFile <- file.path(path, "vcfFiles", "sequences_Prox-10_31-10-2019.fasta")
sequences <- read.dna(fastaFile, format="fasta", as.character=TRUE)
nSites <- ncol(sequences)
# Plot the FASTA file
fastaPlotFile <- paste0(substr(fastaFile, 0, nchar(fastaFile) - 6), ".pdf")
plotFASTA(sequences, pdfFileName=fastaPlotFile, pdfHeight=14, pdfWidth=21, lineForSequenceNames=-3)
# Calculate sequencing quality in the FASTA file
# Note average proportion Ns may be higher because of the extreme coverage in one genome - reaching areas with no coverage in other genomes
propNs <- data.frame("ID"=rownames(sequences), "ProportionNs"=sapply(1:nrow(sequences), calculatePropNsInSequence, sequences),
stringsAsFactors=FALSE)
rownames(propNs) <- parseTipLabels(propNs$ID)
defineBranchColoursOfClade <- function(tree, nodeDefiningClade,
colour, defaultColour){
branchColours <- rep(defaultColour, dim(tree$edge)[1])
clade <- tips(tree, node=nodeDefiningClade)
branchesInClades <- which.edge(tree, clade)
branchColours[branchesInClades] <- colour
return(branchColours)
}
plotSamplingLocations <- function(xCoords, yCoords, scaleSize=10000, xPad=0.75, yPad=0.1, ...){
# Plot the sampling locations
plot(x=xCoords, y=yCoords, ...)
# Note the size of the axes
axisLimits <- par()$usr
xAxisSize <- axisLimits[2] - axisLimits[1]
yAxisSize <- axisLimits[4] - axisLimits[3]
# Note the start and end of scale box
xStart <- axisLimits[1] + (xAxisSize*xPad)
xEnd <- xStart + scaleSize
yStart <- axisLimits[3] + (yAxisSize*yPad)
yEnd <- yStart + scaleSize
# Add a scale box
rect(xleft=xStart, ybottom=yStart, xright=xEnd, ytop=yEnd)
# Add label for scale box
label <- bquote(paste(.(scaleSize/1000), "km")^2)
text(x=xStart+(0.5*scaleSize), y=yStart, labels=label, pos=1)
}
getTipInfo <- function(tipLabels, sampleInfo, propNs){
# Set the row names of the sample info table
rownames(sampleInfo) <- sampleInfo$Sample.Ref
# Create the tipInfo table
tipInfo <- sampleInfo[tipLabels, c("Sample.Ref", "Mapx", "Mapy", "Eartag", "CPHH", "Breakdown.ID", "County", "Year", "TYPING.RESULT")]
rownames(tipInfo) <- tipLabels
# Add the proportion Ns to tipInfo table
rownames(propNs) <- parseTipLabels(propNs$ID)
tipInfo$ProportionNs <- propNs[tipLabels, "ProportionNs"]
return(tipInfo)
}
parseTipLabels <- function(tipLabels){
# Check if already parsed
if(grepl(tipLabels[1], pattern="_") == FALSE){
return(tipLabels)
}
# Initialise a vector to store the parsed labels
parsedLabels <- c()
# Examine eahc tip label
for(index in seq_along(tipLabels)){
# Remove the ">" from the start of the tip label
label <- tipLabels[index]
if(grepl(tipLabels[1], pattern="^>")){
label <- substr(tipLabels[index], 2, nchar(tipLabels[index]))
}
# Split the label and select first part
label <- strsplit(label, split="_")[[1]][1]
# If it doesn't start with "AF" change "-" to "/"
if(grepl(label, pattern="^AF") == FALSE){
label <- gsub(pattern="-", replacement="/", label)
}
# Store the parsed label
parsedLabels[index] <- label
}
return(parsedLabels)
}
calculatePropNsInSequence <- function(sequenceIndex, sequences){
# Count the number of Ns in the sequence
nucleotideCounts <- table(sequences[sequenceIndex, ])
names(nucleotideCounts) <- toupper(names(nucleotideCounts))
# Note the number of Ns
numberMissing <- 0
if("N" %in% names(nucleotideCounts)){
numberMissing <- nucleotideCounts[["N"]]
}
return(numberMissing / ncol(sequences))
}
getNSitesInFASTA <- function(fastaFile){
# Open a connection to a file to read (open="r")
connection <- file(fastaFile, open="r")
# Get first line of file
firstLine <- readLines(connection, n=1)
# Close file connection
close(connection)
# Get the number of sites used in the FASTA file from first line
nSites <- as.numeric(strsplit(firstLine, " ")[[1]][2])
return(nSites)
}
runRAXML <- function(fastaFile, date, path, nBootstraps=100, nThreads=10, outgroup=NULL, model="GTRCAT"){
# Note the directory for the RAxML output files
directory <- paste(path, "RAxML_", date, sep="")
# Check if analyses already run
alreadyRun <- dir.exists(directory)
# If not already run, create output directory for RAxML
if(alreadyRun == FALSE){
suppressWarnings(dir.create(directory))
}
# Set the Working directory - this will be where the output files are dumped
setwd(directory)
# Build analysis name
analysisName <- paste("RaxML-R_", date, sep="")
# Check if already Run and just want to retrieve tree
if(alreadyRun == FALSE){
# Build the command
seeds <- sample(1:100000000, size=2, replace=FALSE) # For parsimony tree and boostrapping
if(is.null(outgroup)){
command <- paste("raxmlHPC",
" -f a", # Algorithm: Rapid boostrap inference
" -N ", nBootstraps,
" -T ", nThreads,
" -m ", model, " -V", # -V means no rate heterogenity
" -p ", seeds[1], " -x ", seeds[2], # Parsimony and boostrapping seeds
" -n ", analysisName,
" -s ", fastaFile, sep="")
}else{
command <- paste("raxmlHPC",
" -f a", # Algorithm: Rapid boostrap inference
" -N ", nBootstraps,
" -T ", nThreads,
" -m ", model, " -V", # -V means no rate heterogenity
" -p ", seeds[1], " -x ", seeds[2], # Parsimony and boostrapping seeds
" -n ", analysisName,
" -s ", fastaFile,
" -o ", outgroup, sep="")
}
system(command, intern=TRUE)
}
# Get the tree and read it in
treeBS <- getTreeFileWithSupportValues(analysisName)
return(treeBS)
}
getTreeFileWithSupportValues <- function(analysisName){
# Get files in current working directory
files <- list.files()
# Select the tree file with BS support values
treeBSFile <- files[grepl(files, pattern=paste("RAxML_bipartitions[.]", analysisName, sep="")) == TRUE]
# Open the file
treeBS <- read.tree(treeBSFile)
return(treeBS)
}
# Read in the sample metadata
sampleInfoFile <- file.path(path, "SampleInformation_28-05-19.csv")
sampleInfo <- read.table(sampleInfoFile, header=TRUE, sep=",", stringsAsFactors=FALSE)
sampleInfo <- sampleInfo[sampleInfo$Sample.Ref != "", ]
# Read in the FASTA file
fastaFile <- file.path(path, "vcfFiles", "sequences_Prox-10_31-10-2019.fasta")
sequences <- read.dna(fastaFile, format="fasta", as.character=TRUE)
nSites <- ncol(sequences)
# Plot the FASTA file
fastaPlotFile <- paste0(substr(fastaFile, 0, nchar(fastaFile) - 6), ".pdf")
plotFASTA(sequences, pdfFileName=fastaPlotFile, pdfHeight=14, pdfWidth=21, lineForSequenceNames=-3)
# Calculate sequencing quality in the FASTA file
# Note average proportion Ns may be higher because of the extreme coverage in one genome - reaching areas with no coverage in other genomes
propNs <- data.frame("ID"=rownames(sequences), "ProportionNs"=sapply(1:nrow(sequences), calculatePropNsInSequence, sequences),
stringsAsFactors=FALSE)
rownames(propNs) <- parseTipLabels(propNs$ID)
# Build a phylogeny using RAxML
tree <- runRAXML(fastaFile, date="31-10-19", file.path(path, "vcfFiles", ""), outgroup="\\>Ref-1997")
# Remove Reference
tree <- drop.tip(tree, tree$tip.label[grepl(tree$tip.label, pattern=">Ref-1997")])
# Convert the branch lengths to SNPs
if(mean(tree$edge.length) < 1){
tree$edge.length <- tree$edge.length * nSites
}
# Parse the tip labels
tree$tip.label <- parseTipLabels(tree$tip.label)
# Extract the densely sample clade
node <- 121
clade <- extract.clade(tree, node=node)
# Get the tip sampling and quality information
tipInfo <- getTipInfo(tree$tip.label, sampleInfo, propNs)
# Count number of tips that we have no data for
warning(length(which(is.na(tipInfo$Sample.Ref))), " sequence labels not found in sample information!")
# Plot the spatial sampling
plotSamplingLocations(tipInfo$Mapx, tipInfo$Mapy, bty="n", xaxt="n", yaxt="n", xlab="", ylab="", pch=19,
asp=1, main="Sampling locations",
col=ifelse(tipInfo$Sample.Ref %in% clade$tip.label, rgb(1,0,0, 0.25), rgb(0,0,0, 0.25)))
legend("right", legend=c("all", "in clade"), text.col=c(rgb(0.1,0.1,0.1, 0.75), rgb(1,0,0, 0.75)), bty="n")
# Plot the temporal sampling
plotMultipleHistograms(distributions=list(tipInfo$Year, tipInfo[tipInfo$Sample.Ref %in% clade$tip.label, "Year"]),
colours=c(rgb(0.1,0.1,0.1, 0.5), rgb(1,0,0, 0.5)), nBins=25, las=1, xaxt="n", xlab="Year",
main="Sampling dates")
yearRange <- range(tipInfo$Year, na.rm=TRUE)
axis(side=1, at=seq(from=yearRange[1], to=yearRange[2]+2, by=2), xpd=TRUE)
legend("top", legend=c("all", "in clade"), text.col=c(rgb(0.1,0.1,0.1, 0.75), rgb(1,0,0, 0.75)), bty="n")
plot.phylo(tree, show.tip.label=FALSE, edge.width=2, edge.color=branchColours)
branchColours <- defineBranchColoursOfClade(tree, nodeDefiningClade=node, colour="red", defaultColour="black")
plot.phylo(tree, show.tip.label=FALSE, edge.width=2, edge.color=branchColours)
addSNPScale(position="bottom", size=50)
plot.phylo(tree, show.tip.label=FALSE, edge.width=2, edge.color=branchColours)
addSNPScale(position="bottom", size=50)
addSNPScale(position="bottom", size=50, lwd=2)
addSNPScale(position="bottom", size=50, lineWidth=2)
plot.phylo(tree, show.tip.label=FALSE, edge.width=2, edge.color=branchColours)
addSNPScale(position="bottom", size=50, lineWidth=4)
?addSNPScale
remove.packages("basicPlotteR")
library("devtools")
library("roxygen2")
packageDirectory <- "/home/josephcrispell/Desktop/Research/basicPlotteR/"
setwd(packageDirectory)
document()