---
title: "sar11blastp.Rmd"
output: html_document
date: "2023-04-27"
output: html_document
style: tidyverse

---

```{r setup, include = FALSE}
library(knitr)
library(tidyverse)
library(kableExtra)
library(DT)
library(reshape2)
library(Biostrings)
library(tm)
```

# Make a database

Comparing protein sequence 012276695protein.faa (for easier terms called 95) to 000012345.1protein.faa (for easier terms called 45). First creating a database for 95.  
```{r, engine='bash'}
/home/shared/ncbi-blast-2.11.0+/bin/makeblastdb \
-in ../data_raw/012276695protein.faa \
-dbtype prot \
-out ../blastdb/95-protein
```


Running BLAST. Comparing 95 to 45. 
```{r, engine='bash'}
/home/shared/ncbi-blast-2.11.0+/bin/blastp \
-query ../data_raw/000012345.1protein.faa \
-db ../blastdb/95-protein \
-out ../output/95-45_blastp.tab \
-evalue 1E-20 \
-num_threads 20 \
-max_target_seqs 1 \
-outfmt 6
```


Running BLAST. Comparing 95 to 85.
```{r, engine='bash'}
/home/shared/ncbi-blast-2.11.0+/bin/blastp \
-query ../data_raw/900177485.1protein.faa \
-db ../blastdb/95-protein \
-out ../output/95-85_blastp.tab \
-evalue 1E-20 \
-num_threads 20 \
-max_target_seqs 1 \
-outfmt 6
```


#Get the query sequence

```{r, engine='bash'}
head ../data_raw/012276695protein.faa
echo "How many sequences are there?"
grep -c ">" ../data_raw/012276695protein.faa
```

##Histogram for the BLAST file of 95-45

```{r, eval=TRUE}
# Read FASTA file

blast_file <- read.csv("../output/95-45_blastp.tab", sep = '\t', header = FALSE) 

#Extract columns desired

blast2.0_file <- blast_file[, c("V1", "V2", "V11", "V12")]

#Make into dataframe

blast2.0_file_df <- data_frame(blast2.0_file)
colnames(blast2.0_file_df) <- c("Query", "Database", "Eval", "Score")

# Plot histogram using ggplot2

ggplot(blast2.0_file_df, aes(x = Score)) +
  geom_histogram(binwidth = 50, fill = "blue", alpha = 0.5) +
  labs(title = "Histogram of Score", x = "Score", y = "Frequency")

```


#Joining blast table with annotation table

fasta to tab-delimited
```{r, engine='bash'}
perl -e '$count=0; $len=0; while(<>) {s/\r?\n//; s/\t/ /g; if (s/^>//) { if ($. != 1) {print "\n"} s/ |$/\t/; $count++; $_ .= "\t";} else {s/ //g; $len += length($_)} print $_;} print "\n"; warn "\nConverted $count FASTA records in $. lines to tabular format\nTotal sequence length: $len\n\n";' ../data_raw/012276695protein.faa > ../output/95protein.tab
```


Having the comparison of 95 to 45 into a CSV file.
```{r, eval=TRUE}
# Read out the blast file into csv
bltabl <- read.csv("../output/95-45_blastp.tab", sep = '\t', header = FALSE)


#Extract columns desired
 bltabl_2.0 <- bltabl[, c("V1", "V2", "V11", "V12")]

#Make into dataframe

bltabl2.0_df <- data_frame(bltabl_2.0)

colnames(bltabl2.0_df) <- c("Query", "Database", "Eval", "Score")
```

Reading the file from 95 protein to CSV file.
```{r, eval=TRUE}
spgo <- read.csv("../output/95protein.tab", sep = '\t', header = TRUE)
colnames(spgo) <- c("Database", "Protein", "Coding_Sequence")
```



```{r, eval=TRUE}
# Join the data frames
joined_df <- left_join(bltabl2.0_df, spgo, by = c("Database" = "Database"))

# Print the table with kable
kable(head(joined_df, n = 10), caption = "Joined data frame with kable formatting (first 10 rows)") 
  
```

```{r, eval=TRUE}
# Read joined blast table with annotation table
joined_df <- left_join(bltabl2.0_df, spgo, by = "Database")
#c("Database" = "Database"))

# Select the column of interest
column_name <- "Protein"  # Replace with the name of the column of interest
column_data <- joined_df[[column_name]]

# Count the occurrences of the strings in the column
string_counts <- table(column_data)

# Convert to a data frame, sort by count, and select the top 10
string_counts_df <- as.data.frame(string_counts)
colnames(string_counts_df) <- c("String", "Count")
string_counts_df <- string_counts_df[order(string_counts_df$Count, decreasing = TRUE), ]
top_10_strings <- head(string_counts_df, n = 10)

# Plot the top 10 most common strings using ggplot2
ggplot(top_10_strings, aes(x = reorder(String, -Count), y = Count, fill = String)) +
  geom_bar(stat = "identity", position = "dodge", color = "black") +
  labs(title = "Top 10 Protein hits",
       x = column_name,
       y = "Count") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip()

```

