Script.R

####-Importing Libraries and Set Up ----
# Avery Murphy
# October 25, 2022
# Assignment 2

# First the working directory is set using "session", "set working directory" and set it with the "to source file location" option 

# Loading Packages 
library(tidyverse)
library(ggplot2)
library(rentrez)
library(seqinr)
library(dplyr)
library(BiocManager)
library(Biostrings)
library(randomForest)
library(styler)

####Section1: Data Acquisition and Cleaning ----

#####_1.Data Acquisition ----
# Data was taken from the NCBI nucleotide database on October 22, 2024
# The data was extracted for each Genus separately so that they were sampled equally.
# The search was specified to only take COI gene sequence's. The search was done initially on each genus without a max set to see how many sample each genus returned.
Anop_data <- entrez_search(db = "nuccore", term = "(Anopheles[ORGN] AND COI[GENE])")
Aedes_data <- entrez_search(db = "nuccore", term = "(Aedes[ORGN] AND COI[GENE])")

# The max sample was set to 600 sequences for each genus so that there is room for filtering and enough data for the training set. 
Anop_data <- entrez_search(db = "nuccore", term = "(Anopheles[ORGN] AND COI[GENE])", retmax = 600, use_history=TRUE)
Aedes_data <- entrez_search(db = "nuccore", term = "(Aedes[ORGN] AND COI[GENE])", retmax = 600, use_history=TRUE)

#####_2.Organizing Data----

# Anopheles Data - Malaria Vector 
# The fasta sequences were extracted for each unique ID
Anop_fetch <- entrez_fetch(db = "nuccore", web_history = Anop_data$web_history, rettype = "fasta", retmax = 600)
# Then the data was written to file and put in the 'Assignment 2' folder inside the 'data' folder.
# write(Anop_fetch, "Anop_fetch.fasta", sep = "\n") 
# The fasta file was imported from the data file:
stringSet <- readDNAStringSet("../data/Anop_fetch.fasta")
# The fasta sequence were made into a DNA strings set
df_Anop <- data.frame(Identifier = names(stringSet), Nucelotide_Sequence = paste(stringSet))
# The species name was taken from each ID and put into a new column 
df_Anop$Species_Name <- word(df_Anop$Identifier, 2L, 3L)
df_Anop <- df_Anop[ , c("Identifier", "Species_Name", "Nucelotide_Sequence")]
view(df_Anop)
# Before the identifier names were cleaned, each was checked to make sure that the sequences were only from COI gene and the mitochondria 
# Then the identifier column was edited to only keep the unique identifier that will be needed for the random forest later in order to separate samples so they are not repeated between the training and validation set.
df_Anop$Identifier <- sub("^(\\S+).*", "\\1", df_Anop$Identifier)
# The names were edited to only keep the Genus names
df_Anop$Species_Name <- sub("^Anopheles.*", "Anopheles", df_Anop$Species_Name)
view(df_Anop)
# Then the column name was changed to Genus name instead of species name
names(df_Anop)[names(df_Anop) == "Species_Name"] <- "Genus_Name"
view(df_Anop)

# The same steps were now completed for the Aedes Data - Zika Vector 
Aedes_fetch <- entrez_fetch(db = "nuccore", web_history = Aedes_data$web_history, rettype = "fasta", retmax = 600)
# write(Aedes_fetch, "Aedes_fetch.fasta", sep = "\n") 
stringSet <- readDNAStringSet("../data/Aedes_fetch.fasta")
df_Aedes <- data.frame(Identifier = names(stringSet), Nucelotide_Sequence = paste(stringSet))
df_Aedes$Species_Name <- word(df_Aedes$Identifier, 2L, 3L)
df_Aedes <- df_Aedes[ , c("Identifier", "Species_Name", "Nucelotide_Sequence")]
view(df_Aedes)
df_Aedes$Identifier <- sub("^(\\S+).*", "\\1", df_Aedes$Identifier)
df_Aedes$Species_Name <- sub("^Aedes.*", "Aedes", df_Aedes$Species_Name)
view(df_Aedes)
names(df_Aedes)[names(df_Aedes) == "Species_Name"] <- "Genus_Name"
view(df_Aedes)

# Then the 2 genus data frames were merged vertically so all the rows for each variable were put together. 
full_df <- bind_rows(df_Anop, df_Aedes)

#####_3.Distribution and Filtering of Missing Sequence Data----

# First the sequence's will be looked at for missing data. A new variable was made to look at the amount of missing nucleotide's for each sequence 
full_df <- full_df %>% mutate(N_count = str_count(full_df$Nucelotide_Sequence, 'N'))
summary(full_df$N_count)
# There is not very much missing data in these sequences and the max is only 5.
view(full_df)
# When looking through the sequences, the N's are located within the sequence and not on the ends. Only one sequence has the max value and the N's are located close together within the sequence. 
full_df %>% filter(N_count >= 1)
# Only 36 sequences having missing nucleotide's and it is a very minimal amount. 

# Then the amount of gaps in the sequences were looked at:
full_df <- full_df %>% mutate(Missing_data = str_count(full_df$Nucelotide_Sequence, '-'))
view(full_df)
# There is no gaps in any of the sequences so no filtering step for this will be required.

# Now a column with each sequences length will added, so that the percentage of N's out of total base pairs for each sequence can be calculated. 
full_df <- full_df %>% mutate(Length = nchar(full_df$Nucelotide_Sequence))
full_df <- full_df %>% mutate(N_count_percentage = (N_count/Length) *100)
view(full_df)
# All the sequences with N's have a total of less than 1% N's. The small percentage should not impact the accuracy of sequence features, with that only the sequence with the largest percentage was removed because the N's were close together showing a poorly sequenced area, so the row 834 was removed. 
full_df <- full_df [ -(834), ]
dim(full_df)
# Check to see an observation was removed, there is now 1,999 rows. 

#####_4.Distribution and Filtering of Sequence Length----

#Next the sequence length is looked at after filtering for missing data. 
summary(full_df$Length)
# The summary shows a very large range in sequence length. 
view(full_df)
# Some of the lengths were over 15,000 bp for the Anopheles’s samples, at closer look by searching the ID's back in the NCBI nucleotide browser, these samples contained the whole mitochondrial genome and not just the COI gene. These samples can be removed by the row number.  The rows with over 15,000 bp were identified and removed. 
which(full_df$Length > 15000)
full_df <- full_df [ -c(445,446,447,448,449,450,451,452,453,454,455), ]
summary(full_df$Length)
# The summary is re-run and the max is now 700 bp

# Now the sequence lengths can be looked at for each Genus In Figure 1: 
boxplot(Length ~ Genus_Name, data= full_df, main = "Boxplot of Sequnce Length", ylab = "Base Pair Lengths", xlab = 'Genus Name', col = c( "lightblue", "lightgreen"))

# Figure 1: Examining the box plots, the 2 genus have over lapping lengths, both with some lengths much smaller and under 500 bp. Aedes genus has some longer sequence that are considered outliers compared to the other sample sequence lengths. 650 bp of the COI gene is typically used to distinguish organisms by the COI gene (Herbet et al., 2003).  The sequences with lengths between the q1 (572 bp) and the q3 (658 bp) will be kept for machine learning. This will remove sequence's with abnormally short and long lengths for both genus. 

full_df <- full_df %>% filter(Length >= 572 & Length <= 658)
# Confirm that the filtering worked through the summary and histogram 
summary(full_df$Length)
ggplot(full_df, aes(x=Length)) + geom_histogram()

table(full_df$Genus_Name)
# Now the data has been filtered for lengths and missing data, there are 427 Aedes and 436 Aanopheles samples that will be used to develop a classifier to separate these 2 genus based on the COI sequence.

####Section2: Machine Learning ----

#####_1.Calculating Sequence Properties----
# To separate the genus based on sequences, the properties of the sequences will need to first be determined. The K-mer frequency will be calculated.

# First a k-mer with length 3 will be calculated.  This calculates the frequency of 3 nucleotide’s appearing in a row within a sequences. Starting with a 3K-mer over a 2K-mer offers more specificity and more sequence features to train the classifier with. If the model can separate the sequences based on a 3 length K-Mer than a 2 K-Mer could be explored. First the nucleotide's need to be converted into biostrings so that not only character functions can be used. 
full_df <- as.data.frame(full_df)
full_df$Nucelotide_Sequence <- DNAStringSet(full_df$Nucelotide_Sequence)
class(full_df$Nucelotide_Sequence)

# Calculating the 3 K-Mer Frequency, as.prob = TRUE was used so it accounts for the variability in sequence length as they range from 200 bp differences.
full_df <- cbind(full_df, as.data.frame(trinucleotideFrequency(full_df$Nucelotide_Sequence, as.prob = TRUE)))
# Columns 8 to 71 are the 3K-Mer data 

# Then a new dataframe was made to look at the average of each 3 Length-Kmer across all sequence's in each genus. 
KMer3 <- full_df[, c(2, 8:71)]
genus_means <- KMer3 %>%
  group_by(Genus_Name) %>%
  summarise(across(everything(), mean))
view(genus_means)

# Looking through the average's some are similar but many are quite different such as AGA and ATG for example.  It is predicted at this point the model should be able to classify a genus based on the 3 K-Mer frequencies. 

#####_2.Developing Validation and Training Data----

# The nucleotides were converted back into character data for future use 
full_df$Nucelotide_Sequence <- as.character(full_df$Nucelotide_Sequence)

# First a validation data set must be made, this data will be remain hidden during the training process so that the algorithm cannot see the sequences and they will be used to validate the correct classification.  It is important that a random sample of the same number of sequences are taken from each genus. Around 25% of sequences for each genus will be set aside for the validation set. The data will grouped by genus and 105 random samples from each will be taken.  The seed will be set so that this can be reproduced and the same samples will always be taken. 
set.seed(217)
df_Validation <- full_df %>%
  group_by(Genus_Name) %>%
  sample_n(105)
# Checking there is 105 for each Genus
sort(table(df_Validation$Genus_Name))

# Now the training data for the random forest will be made, for this it is crucial that the data that was already taken for the validation set is not included in the training. The data is filtered with "!" so that it knows not to choose ones already in the validation data. This is why the unique NCBI ID was kept for each sequence so they could be separated. The training data consisted of 322 samples from each Genus (Took 427 the minimum - 105 from validation set to get 322)  
set.seed(13)
df_Training <- full_df %>%
  filter(!Identifier %in% df_Validation$Identifier) %>%
  group_by(Genus_Name) %>%
  sample_n(322)
sort(table(df_Training$Genus_Name))
# Confirmed that each genus had 322 samples

#####_2.Random Forest----

# Now the classifier will be built off the training data to classify genus based on the 3-Kmer Frequency.  The number of tress to grow was set to 100 initially to make sure every input row gets predicted and reduce the OOB error.
genus_classifier_1 <- randomForest::randomForest(x = df_Training[, 8:71], y = as.factor(df_Training$Genus_Name), ntree = 100, importance = TRUE)
genus_classifier_1
genus_classifier_1$confusion
# The classifier was not able to perfectly separate the samples into the correct genus, with 3 samples incorrectly classified. There was one false positive for Anopheles and 2 for Aedes. An error rate of 0.47%, meaning that 0.4 % of the time the model miss classified a sample. This is a very low error rate and the model classified quite well, but I decided to re-run with more decision trees, as this can increase randomization and reduces varießnce in the models prediction to see if this leads to perfect classification (Liaw & Wiener, 2002). 

genus_classifier_2 <- randomForest::randomForest(x = df_Training[, 8:71], y = as.factor(df_Training$Genus_Name), ntree = 150, importance = TRUE) 
genus_classifier_2
# Now the overall error rate dropped to 0.16, meaning the model incorrectly classified a sample 0.16% of the time which is lower then when only 100 trees were used. There was only 1 sample from the Anopheles that was classified incorrectly. This model will be used to test the validation data. 

genus_classifier_2$confusion
# The rows show the true classes and the columns shows prediction, in this training each sample was assigned to the correct genus. There was confusion and error for 1 sample of the Anopheles genus that was miss classified. 

genus_classifier_2$oob.times
summary(genus_classifier_2$oob.times)
# This shows the amount of times each sample was left out when building a decision tree. This increases the generality of the model.  The sample can only be predicted when it is not included. The average time each element could be predicted is 55 times. 

genus_classifier_2$votes
# Not every row has agreement for only one genus meaning that not all trees had 100% agreement when trying to classify each observation. Despite this the classifier was able to separate the training data samples correctly in each genus close to perfect.  

#Now the classifier will be tested on the validation set aside:
predictValidation <- predict(genus_classifier_2, df_Validation[, c(2, 8:71)])
predictValidation
table(observed = df_Validation$Genus_Name, predicted = predictValidation)

# The table will be made into a confusion matrix to display the results in Figure 2: 
confusion_matrix <- table(observed = df_Validation$Genus_Name, predicted = predictValidation)
confusion_df <- as.data.frame(confusion_matrix)
ggplot(data = confusion_df, aes(x = predicted, y = observed)) +
  geom_tile(aes(fill = Freq), color = "white") +
  scale_fill_gradient(low = "lightblue", high = "blue") +
  geom_text(aes(label = Freq), color = "white") +
  labs(title = "Confusion Matrix - 3Kmer", x = "Predicted", y = "Observed") +
  theme_minimal() + theme(plot.title = element_text(hjust = 0.5))

# Figure 2: The classifier was able to perfectly classify the samples into each genus, the first 105 as Aedes and the last 105 as Anopheles with 0 miss classified. Using 3-Kmer sequence properties of the COI gene worked to distinguish samples of 2 the Culicade genus. Now the classification will be done using a Kmer with length 2 instead to see if the model can still classify based on a less specific property. 

####Section3: Machine Learning with a Smaller K-mer ----

#####_1.Recalculating Sequence Properties----

full_df <- as.data.frame(full_df)
full_df$Nucelotide_Sequence <- DNAStringSet(full_df$Nucelotide_Sequence)
class(full_df$Nucelotide_Sequence)
# Calculating the 2 K-Mer frequency for each sequence:
full_df <- cbind(full_df, as.data.frame(dinucleotideFrequency(full_df$Nucelotide_Sequence, as.prob = TRUE)))
# Columns 72-87 have the 2 K-mer frequency. 

# Now the averages of each 2 K-Mer frequency for the sequences will be looked at:
KMer2 <- full_df[, c(2, 72:87)]
genus_means_2 <- KMer2 %>%
  group_by(Genus_Name) %>%
  summarise(across(everything(), mean))
view(genus_means_2)

# Looking through the averages of the 2 length K-Mers frequency, they appear more similar between genus's than what was calculated in the 3-Kmer.  Some were quite different like TC and AA.  It is predicted that the classification should still work with this sequence feature.

#####_2.Developing Validation and Training Data----

full_df$Nucelotide_Sequence <- as.character(full_df$Nucelotide_Sequence)

# Validation Data
set.seed(217)
df_Validation_2 <- full_df %>%
  group_by(Genus_Name) %>%
  sample_n(105)
sort(table(df_Validation_2$Genus_Name))
# Check each has 105 samples 

# Training Data 
set.seed(13)
df_Training_2 <- full_df %>%
  filter(!Identifier %in% df_Validation$Identifier) %>%
  group_by(Genus_Name) %>%
  sample_n(322)
sort(table(df_Training_2$Genus_Name))
# Check each has 322 samples

# Setting up the random forest and this time starting with 150 decision trees
genus_classifier_3 <- randomForest::randomForest(x = df_Training_2[, 72:87], y = as.factor(df_Training_2$Genus_Name), ntree = 150, importance = TRUE)
genus_classifier_3
genus_classifier_3$confusion
# The 2-Kmer had an over all error of 0.7%, a higher error rate than what was seen when a 3 K-Mer was used to train the model. It incorrectly identified 5 samples altogether with, 4 Anopheles genus were classified as Aedes and 1 Aedes was classified as Anophles.

genus_classifier_3$oob.times
summary(genus_classifier_3$oob.times)
#Overall the median time each sample could be predicted was the same as in the 3 K-Mmer classification. 

genus_classifier_3$votes
#Not always agreement as seen with the 3 K-Mer classification. 

# Testing out the classification on the validation data:
predictValidation_2 <- predict(genus_classifier_3, df_Validation_2[, c(2, 72:87)])
predictValidation_2
table(observed = df_Validation_2$Genus_Name, predicted = predictValidation_2)

# The table was again made into a confusion matrix to display the results in Figure 3:
confusion_matrix_2 <- table(observed = df_Validation_2$Genus_Name, predicted = predictValidation_2)
confusion_df_2 <- as.data.frame(confusion_matrix_2)
ggplot(data = confusion_df_2, aes(x = predicted, y = observed)) +
  geom_tile(aes(fill = Freq), color = "white") +
  scale_fill_gradient(low = "lightblue", high = "blue") +
  geom_text(aes(label = Freq), color = "white") +
  labs(title = "Confusion Matrix - 2Kmer", x = "Predicted", y = "Observed") +
  theme_minimal() + theme(plot.title = element_text(hjust = 0.5))

# Figure 3: The classifier was not able to assign all samples to the correct genus, with 3 samples miss classified.  The model incorrectly classified Aedes as Anopheles twice and incorrectly classified Anopheles as Aedes one time. The deep blue square forming the diagonal represent perfect classification.