README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  message = FALSE,
  warning = FALSE
)

library(tidyverse)

ggplot2::theme_set(ggprism::theme_prism(palette = "viridis"))

```

# GLabR <a href='https://github.com/baynec2/GLabR'><img src='man/hex_GLabR.png' align="right" height="139" /></a>


<!-- badges: start -->
[![R-CMD-check](https://github.com/baynec2/GLabR/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/baynec2/GLabR/actions/workflows/R-CMD-check.yaml)
<!-- badges: end -->

The goal of GLabR is to provide a centralized locations to hold functions that are routinely useful in the Gonzalez Lab at UCSD

## Installation

You can install the development version of GLabR like so:

```{r,eval=FALSE}
devtools::install_github("baynec2/GlabR")
```

## Examples

### Data normalization (batch correction)

The first use case of GlabR is to normalize data that we export from proteome discoverer. Here, a text file containing PSMs is exported and subsequently needs to be processed. To accomplish this, there are a few different steps that need to be followed . These are described below. 

1. We need to combine PSMs from each of the different fractions (there are multiple fractions corresponding to a single sample). This is accomplished using the combine_psm_fractions() function. Essentially, this function filters out PSMs based on certain criteria and then sums the intensities for each protein that they map to.

2. Then we need to normalize our data to account for batch corrections. This is accomplished using the normalize_to_bridge() function. 

  * In some cases, there may not be a bridge channel included, and in these cases we will need to use the normalize_1plex() function instead

3. After this, we still might want to normalize the data further. To get data that is more normally distributed, we can use Leigh-ana's method of box cox normalization through the la_box_cox_norm() function. 


Here we can see an example of the overall workflow using a single plex example set

```{r example}
library(GLabR)

data = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_1plex()
 
head(data)                   
```


The final norm column has the completely normalized data. 



Note that we decided to set this up such that NAs are . Previous code converted NAs to 1s. This has the advantage that it is easier to deal with but is not entirely appropriate to do. As such, I have decided to make the user explicitly change NAs if they desire. 



Here we can see an example of the overall workflow using an example set with multiple plexes 

```{r}
data = read_delim("tests/testdata/normalize_to_bridge/PSM_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126)


head(data)
```



If you prefer data in a wide format ( as reported using the previous script), you can specify data_format = "wide" within the normalize to bridge function. This will give a column containing the final_norm values for each Sample/TMT combination. Note that this is also an option for the normalize_1plex function.  An example is shown below

```{r}
data = read_delim("tests/testdata/normalize_to_bridge/PSM_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126,data_format = "wide")


head(data)

```



If we wanted to then use the box cox norm method, we could use la_box_cox_norm() function. 

```{r}
data = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126) %>% 
  la_box_cox_norm()


head(data)
```





Note that when using this function, it is expecting your data to be in the long format. If not, it will not work. If you wish for the results reported from la_box_cox_nrom to be in the wide format, you can specify that like below.

```{r}
data = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126) %>% 
  la_box_cox_norm(data_format = "wide")


head(data)
```

Out of curiosity, what does the normalized data look like when comparing the two methods. 

```{r}

p1 = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126) %>% 
  la_box_cox_norm() %>% 
  ggplot(aes(box_cox_scaled_values))+
  geom_histogram()+
  ggtitle("box cox scaled values")



p2 = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_to_bridge(bridge_channel_plex = 126) %>% 
  ggplot(aes(final_norm))+
  geom_histogram()+
  ggtitle("batch corrected")

p3 = ggpubr::ggarrange(p1,p2)

p3
```
Here we can see the overall distribution of the data. Looks like the box cox does a pretty good job of normalizing the data! 

Note that when I was making this package, I discovered a problem in the script that was originally used to normalize the data. Essentially, the original script was only assigning the NAs from some columns as 1, ultimately causing results to deviate from what was intended (NA from all columns having 1)

Here I have kept the function that produces the same results as the output from the original script (nonnormalizedall.txt) as combine_psm_fractions_replica() for posterity. 


Here we can see the slight differences produced from these different functions. 
```{r}
#using function to replicate new process
new = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  mutate(method = "new")

#using function to replicate old process
old = read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions_replica() %>% 
  mutate(method = "old")

#Combining data
comb = bind_rows(new,old) %>% 
  pivot_wider(names_from = method, values_from = value)


#Plotting
p1 = comb %>% 
  ggplot(aes(new,old))+
  geom_point()+
  geom_smooth(method = "lm")

p1
```

As we can see above, while this data is correlated, it is not exactly the same. **Future work should only use the combine_psm_fractions() function.** 

### Determining Differential Protein Abundance

Usually when analyzing proteomics data, we will be interested in assessing the differences between two or more groups. 


There are a number of ways to do this, but the most simple and most commonly used approach to visualize the differentally abundant proteins between tow groups would be by using a volcano plot. We can use GLabR to easily plot a volcano plot as displayed in the code below.  


First we have to load our data, and normalize it appropriately. To do this we use the functions in GLabR that were previously described.  

```{r}
data = readr::read_delim("tests/testdata/combine_psm_fractions/PCB002_PSMs_Proteodiscover_output.txt") %>% 
  combine_psm_fractions() %>% 
  normalize_1plex() %>% 
  la_box_cox_norm() %>% 
  #Have to add the sample ID column manually by concating Sample and TMT.
  #Done outside of GLabR to allow the user freedom in determining naming convention and contents (. vs _ seperator etc) of what I would call Sample_ID
  dplyr::mutate(Sample_ID = paste0(Sample,".",TMT))
```

Next, we need to assign metadata that gives us information about what each of these samples are.

```{r}
# Loading metadata
md = readr::read_csv("tests/testdata/metadata.csv") %>% 
  #Removed redundant columns to prevent .x and .y columns in data_md
  dplyr::select(-Sample,-TMT)

# Appending md to data
data_md = dplyr::inner_join(data,md,by = "Sample_ID")
```

Now we can use the volcano_plot() function to plot the differences between two types of samples using conditions contained in the metadata. Here let's compare the greatest disease severity to healthy controls using the mayo score.

Note that the function has trouble dealing with weird column names, so to avoid this use column names without spaces, :s, or other characters that are dealt with differently in R

Also, note that you can change the p_theshold and log2fc threshold using arguments. Any proteins above these thresholds will be plotted. 

Note that the default output of the processing pipeline returns NA if there is no value associated. This can be problematic for downstream data processing. This is intentinal, as it requires dedicated thought by the user as to how best to deal with those values. Previous versions of the script always assigned them a value of 1. Here I have opted to use dplyr::mutate_ifand tidyr::replace(na) to transform these NAs to 0s.

```{r,fig.height =6, fig.width=10,eval = TRUE}
# Filtering our data to contain only the two conditions we want to test
# NEed to fix this later
f_data_md = data_md %>% 
  dplyr::filter(`Mayo_Endoscopic_Sub_Score` %in% c("Healthy_control","3: Severe disease (spontaneous bleeding, ulceration)")) %>% 
  dplyr::mutate_if(is.numeric,~tidyr::replace_na(.,0))

#

data = f_data_md
column_split_by = "Mayo_Endoscopic_Sub_Score"

 stat = data %>%
      dplyr::group_by(ProteinID)
      tryCatch({
      stat = rstatix::t_test(stat,as.formula(paste("box_cox_scaled_values",'~',column_split_by)))
      stat = rstatix::adjust_pvalue(stat,p.col = "p",output.col = "p.adj_fdr", method = "fdr")
      }, error = function(cond){
        return(NA)
        })
volcano_plot(f_data_md,"Mayo_Endoscopic_Sub_Score",p_threshold = 0.05,log2fc_threshold = 1)

```

Here we can see our volcano plot! We can note that there are 5 proteins that meet our criteria, and are called out by ProteinID on the plot. Note that this function uses the t_test function from the excellent rstatix package to determine statistical significance. Importantly, the pvalue reported has been adjusted for multiple comparisons using fdr. 


### Protein Idenfification. 

The next step in the process is figuring out what these proteins actually are/ what they do. In order to do this, we can use a combination of extract_sig_proteins and our annotate proteins function. 


extract_sig_proteins is a function that is intended to extract the proteinIDs that are significant given the desired parameters. If the same parameters are passed to the function as volcano plot, the proteins will match between the two. We can see this below. 


```{r,eval = FALSE}
sig_proteins = extract_sig_proteins(f_data_md,column_split_by = "Mayo_Endoscopic_Sub_Score",p_threshold = 0.05,fc_threshold = 1)

sig_proteins
```


Now that we have a list of our proteins we can figure out what they are using the annotate_proteins function. 


This function uses Uniprot's API to return results pertaining to the proteins. There is a package called UniprotR that handles this, but it was way too slow for long lists of protein IDs, so I made this solution. We can annotate these proteins as follows.  



```{r,eval = FALSE}

annotated_proteins = annotate_proteins(sig_proteins)

annotated_proteins
```

Let's say you wanted to get different information about these proteins. To do that, you can specify what columns to return through the columns argument. This argument takes a string of the column names that you want to add separated by columns. The complete list of field names that are accepted can be found here: https://www.uniprot.org/help/return_fields. 

In this example, lets say we want to get the GO biological process terms for our proteins (column name accessed through api by "go_p". We can do this as follows. 

```{r,eval = FALSE}
# Result table will have accesssion and GO biological process info
annotated_proteins_GO_p = annotate_proteins(sig_proteins,columns ="accession,go_p")

annotated_proteins_GO_p

```


Here we can see that information!  

### Phosphoproteomics


The above code is great for looking at proteomics, but sometimes we will want to look at phospo post translational modifications. 

The calc_phospho_ratio function allows us to do this! We can see an example below:


```{r}
proteomics_data = readr::read_delim("tests/testdata/psm_phospho_mod/PCB002_PSMs.txt")
phospho_data = readr::read_delim("tests/testdata/psm_phospho_mod/PCB001_PSM.txt")
metadata = readxl::read_excel("tests/testdata/psm_phospho_mod/metadata.xlsx")
#In this dataset, patient ID is the variable that we want to use to  combine the phospho and proteomic data. 
col_identifying_match = "PatientID"

phospho_ratio = calc_phospho_ratio(proteomics_data,phospho_data,metadata,col_identifying_match,1)

head(phospho_ratio)
```


In short, this code pairs up corresponding proteomics and phosphoproteomics data sets and then calculates the ratio of phosphorylation to using the box cox normalized data.


Sometimes, it will be useful to compare how many unique peptide-sequences-PTM are present in a phosphoenriched experiment, compared to a normal proteomics experiment. Not that in order to do this, you must  have your proteome discoverer analysis set up to report infomration about phosphorylated peptides. You can do this by setting dynamic settings to accomodate phospho mods in the sequest node, and enabling the ptmRS node. 


Once you have these results, you can compare them using the following code. Note that for this to make any sense to do, the experiments should contain the same samples. 

```{r}
phospho_enriched = readr::read_delim("tests/testdata/psm_phospho_mod/PCB001_PSM.txt")

proteomics = readr::read_delim("tests/testdata/phospho_venn_diagram/PCB002_proteomics_ptmRS_data.txt")

phospho_venn_diagram(proteomics,phospho_enriched)
```
Here we can see that the phospho enriched experiment enriched for phospho peptides, as we would expect. 


### Designing Experiments

GLabR also has functions that are intended to help design experiments. All together, these functions are designed to take a data frame of metadata, assign a plate number, randomly assign a 96 well location, assign a TMT plex number with an incorporated bridge channel (allowing for the option to group samples within the same plex by a factor if desired), and then randomly assign a TMT label for each sample in the plex. 

Note that these functions are designed to be modular. If you already have metadata with the plate numbers that were generated by some other method that you want to use (ie manually in excel), than no need to randomize the plates and so on.  

Right now, the 10 and 16 plex TMT options are supported.   

We can see an example as follows for a case where we have 108 samples across and would like to pair by mouse.  

```{r}
metadata = data.frame(sample_id = as.character(1:108),
                      mouse_id = rep(paste0("Mouse",1:27),each = 4)
                    )

#define metadata when using 16 plexes
final_metadata = metadata %>% 
  randomize_plates() %>% 
  randomize_wells() %>% 
  assign_plex_num("mouse_id") %>% 
  randomize_tmt()

head(final_metadata)
```


What if for some reason we wanted to use 10 plexes instead. We can do that as follows.


- need to fix this before stable release


### Peptide Quantification


GLabR also supports the automated formatting and quantification of peptides using the pierce peptide quantification kit and the export from defined templates (peptide_quant_od480_single_plate.spr, or peptide_quant_od480_multi_plate.spr) for SoftMax software. 


Note that these protocols assume that you added samples in locations indicated by the plate maps in the software. For example, peptide_quant_od480_single_plate.spr assumes that the standard curve has three replicates that are in columns 1:3, and the concentration decreases by 1/2 rowwise until you get to row G. Row H is a concentration of 0. All of the other wells have unknown samples. 

In the case that there are multiple plates, the peptide_quant_od480_multi_plate.spr assumes that you have loaded standards on one plate (three replicates that are in columns 1:3, and the concentration decreases by 1/2 rowwise until you get to row G. Row H is a concentration of 0), and that the unknowns are on other plates. 


We can use these functions as shown below.   

```{r}
# Single plate
single = peptide_quant("tests/testdata/peptide_quant/single_plate_od480.txt")

# Multiple plates
multi = peptide_quant("tests/testdata/peptide_quant/multi_plate_od480.txt")


head(multi)
```

Say you want to plot the standard curve used to generate the known vaues. You can do that, or any other function you want easily as follows:

```{r}
stds = multi %>% 
  filter(Sample_Type == "standards") %>% 
  ggplot(aes(Conc_ug_mL,Blank_OD480))+
  geom_point()+
  geom_smooth(method = "lm")+
  ggprism::theme_prism()

stds
```

Say you want to look at the details about the linear model used to predict the values. You could just generate the same model and inspect as follows:

```{r}
stds = multi %>% 
  filter(Sample_Type == "standards")

summary(lm(data = stds, Blank_OD480 ~ Conc_ug_mL))
```


### Pubchem

We can also use the annotate pubchem function to search Smiles for synonyms. 

Note that this can only find exact SMILE matches, it will not report results for an isomer like the web portal will

Also note that this will return warnings when a search string cannot be found in pubchems database.

lastly please not that the way I am determining what name to return isn't perfect. There are tons of synonyms for each compound, so returning the proper name doesnt seem to be trivial. Here if the first name at the top of the returned list has only numbers, I am then going for the next entry in list. In the future I can patch this with a better fix but this should be sufficient for now. 

```{r}
test = read.delim("tests/testdata/annotate_pubchem/Random_metabolomics_calls.txt") %>%
  pull(Smiles)


t = annotate_pubchem(test)

head(t)

```

Note right now this doesn't do a great job as sometimes the synonym returned isn't the one you would want ie "58-73-1" isnt helpful. Also is rather slow. Will figure out how to fix and update in future. 

# annotate_megadb

mega db is a database designed by Leigh-Ana to be all encompassing. In other words, it intends to contain an as comprehensive list as possible of proteins that may be present in stool samples. As of now, we 


I have hosted this database as a postgres sql database on aws due to the large size (~5gb). I subsequently stopped hosting it because it cost money-> I thought it would be in the free tier but I was wrong.

As such, the annotate_megadb function that was designed to make it as easy as possible to interact with this database from within R no longer works. I may bring back support for this in the future if it comes up. 


When it was supported, this is what  an example looked like:

```{r,eval = FALSE}
proteinids = read_csv("tests/testdata/megadb/proteinids.csv") %>% 
  pull(ProteinID)

megadb_results = annotate_megadb(proteinids) 
head(megadb_results)
```

Note that if a protein ID is not found when searched against uniprot, it is dropped from the results. This isn't necessarily ideal but I will have to think about the best way to address this in the future. Ideally, we would have all the information on the proteins that are contained withing the megadb.fasta file so then there would be no reason to search uniprot each time. 

# assign_food 

Dealing with food is challenging. We need a convenient way to classify foods which isn't exactly trivial, Fortunately, Leigh-Ana already did this work, so we just need a way to interact with her annotations. 

In this workflow we can first assign common names to our proteins as we did before using the megadb function. After this, we can sue the annotate_foods function to interact with the postgres sql database containing the food annotations.

This function also is currently not supported as the underlying sql database is not still being run on AWS. 

```{r, eval=FALSE}
#Assigning common names 
proteinids = read_csv("tests/testdata/megadb/proteinids.csv") %>% 
  pull(ProteinID)

megadb_results = annotate_megadb(proteinids) 

#pulling common names
common_names = megadb_results %>% 
  pull(common_name) %>% 
  unique()

#Now we can figure out what these foods are
foods = annotate_food(common_names)

head(foods)
```


# Donut plot  

Let's say we wanted to view the data as a donut plot. We can do that as follows:

```{rm,eval=FALSE}
data = inner_join(megadb_results,foods,by = "common_name")

donut_plot(data,"sample_type_group_3")
```

Let's say we wanted to get more granular- we could do that as well by selecting a different column! 

```{r,eval = FALSE}
donut_plot(data,"sample_type_group_4")
```
As a side note, we can use this function for other parts of the workflow as well. Let's say we wanted to see which database most of the proteins on a list of IDs are mapping to. We could do that as follows:

```{r,eval=FALSE}

donut_plot(megadb_results,"database")

```

# OD600 data from stratus plate reader  

We purchased a plate reader from cerillo bio that has a really small footprint and can easily be used in anaerobic chambers (https://cerillo.bio/stratus/). This plate reader exports the data as a csv file, but it isn't very useful in the format they provide. This function is provided to make the files much easier to work with.  

```{r}
parsed_data = parse_stratus("tests/testdata/parse_stratus/stratus_data_export.csv")

head(parsed_data)
```
# parse_cfu

Sometimes it might be useful to parse data from an excel sheet that contains colony counts This excel template can be found in the repo templates/parse_cfu_template.xlsx.

It looks like this: 
![](man/parse_cfu.png)
The idea is to enter all of the info pertaining to your plate, and to record the colony count with a nice, user friendly plate layout. I only count 1 dilution, whichever is the "countable" dilution and leave all the others blank. 

The problem with this template is that data in this format isn't great for data analysis. That is where there parse_cfu function comes in. Parse_cfu can take this excel file as an input, and parse it to a data frame- see below for an example. Note that in this experiment the spots were plated with 5uL, so that is passed to the function to allow for calculation of the cfu_mL

```{r}
CFU_data_frame = parse_cfu("tests/testdata/parse_cfu/colony_counts.xlsx",uL_plated = 5)

CFU_data_frame
```

^ here we can see the nicely parsed data! Now it is easy to plot/ manipulate the data to do whatever you would like.

note that all empty counts (read as NAs) get dropped from the count. This is intended behavior to make the final data frame more compact.  

```{r}
library(ggplot2)

ggplot(CFU_data_frame,aes(sample_id,cfu_mL))+
  geom_point()+
  coord_flip()


```