-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
557 lines (329 loc) · 22.4 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
---
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()
```