The following script was done using a yeast transcriptomics data set from Kaggle (https://www.kaggle.com/datasets/costalaether/yeast-transcriptomics).

Step One: Loading in required packages

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tibble)
library(RColorBrewer)

Step Two: Reading in, transforming, and filtering data

The following code reads in the dataframe containing the conditions of the yeast culture with their associated treatment ID. The conditions and expression dataframes are obtainable from the data set. The GO dataframe is a merged dataframe for the three different types of GO terms (MF, CC, and BP) for each gene that is also found in the above data set.

conditions <- read.csv(file = "conditions_annotation.csv", header = TRUE, stringsAsFactors = FALSE)

GO <- read.csv(file = "mergedLabels.csv", header = TRUE, stringsAsFactors = FALSE)

expression <- read.csv(file = "SC_expression.csv", header = TRUE, stringsAsFactors = FALSE)

The goal of this project is to look at expression data under a specific condition in the data set, which be reported using the following code. There are two different conditions for each sample set.

unique(conditions$primary)
##  [1] "wildtype"                         "itc1"                            
##  [3] "swr1"                             "tet-STH1"                        
##  [5] "tet-INO80"                        "tet-control strain"              
##  [7] "Strain5"                          "Strain6"                         
##  [9] "Strain7"                          "anMX4"                           
## [11] "rrp6"                             "MATa MAL2-8c SUC2"               
## [13] "MATa ura3-52 his3-1 leu2-3"       "hybrid  cerevisiae x  paradoxus "
## [15] "hybrid  cerevisiae x  paradoxus"  "DDY3630"                         
## [17] "DDY4300"                          "DDY4301"                         
## [19] "E1A1"                             "E1A2"                            
## [21] "E1B1"                             "E1B2"                            
## [23] "E2A1"                             "E2A2"                            
## [25] "E2B1"                             "E2B2"                            
## [27] "G1A1"                             "G1A2"                            
## [29] "G1B1"                             "G1B2"                            
## [31] "G2A1"                             "G2A2"                            
## [33] "G2B1"                             "G2B2"                            
## [35] "YB210"                            "YB211"                           
## [37] "YB212"                            "YB213"                           
## [39] "wild type"                        "pbs2"                            
## [41] "ste12"                            "fy3"                             
## [43] "15 deg"                           "37 deg"                          
## [45] "30 deg"                           "37 deg "                         
## [47] "undocumented"
unique(conditions$secondary)
##  [1] "wildtype 1"         "wildtype 2"         "itc1-1_dUTP"       
##  [4] "swr1 mutant"        "<not provided>"     "tet-ino80"         
##  [7] "tet-control"        "phenol lysis"       "Mat a"             
## [10] "BY4741 background"  "rrp6"               "CEN.PK113-7D"      
## [13] " hxk1::KlLEU2"      "Citrinin"           "Untreated"         
## [16] "strain W303"        "ethanol"            "glucose"           
## [19] "biofuel generation" "none"               "galactose"         
## [22] "tunicamycin"        "salt stress"        "chemostat"         
## [25] "temperature"        "undocumented"

For this workflow, the expression within biofuel producing yeast strains will be compared. The conditions dataframe contains two identical treatment IDs (INICIA) but denoted separately in the expression data as one is followed by a ‘.1’. This change will need to be made in the filtered treatment dataframe.

biofuel_strain <- conditions[grepl("biofuel generation", conditions$secondary),]

biofuel <- biofuel_strain %>%
  mutate(ID = make.unique(as.character(ID)))
# Will see that INICIA is repeated and add a .1 to the second one. The treatment ID of this dataframe will now match that of the expression dataframe.

Step Three: Merging datasets

The biofuel dataframe will now be merged with the expression data based on the treatment IDs found within both dataframes. To achieve this, the treatment IDs from the biofuel dataframe will be turned into a list to be used as a filter to specify the data wanted from the expression dataframe.

names(expression)[1] <- "Gene"

expression_filter <- biofuel$ID
expression_filter
##  [1] "INICIR"   "INICIA"   "INICIA.1" "INICIC"   "INICIB"   "INICIF"  
##  [7] "INICIQ"   "INICII"   "INICIS"   "INICIN"   "INICSR"   "INICSA"
biofuel_expression <- expression %>%
  select(Gene, all_of(expression_filter))

The biofuel expression dataframe now only contains the gene expression data of the treatment groups where biofuel generation was studied.

Another step can be taken to merge the GO terms of each gene into the biofuel expression data.

biofuel_go <- left_join(biofuel_expression, GO, by = "Gene")

Using the pivot_longer() command allows us to now look at gene exression of individual genes between the treatment groups. Setting names_to = “treatment’ creates a column of the treatment groups for each gene that was observed. The values_to =”count” will create a new column where the TPM of each gene by treatment group is displayed.

biofuel_long <- biofuel_expression %>%
  pivot_longer(
    cols = c(2:13),
    names_to = "treatment",
    values_to = "count"
  )

Step Four: Summary Statistics and Data Visualization

To further prepare the data for visualization, we can check if there are any NAs in the count to column to see if they need to be accounted for. A tibble can then be built to include some of the summary statistics for each treatment group. In this case, the mean and median of the TPMs for each treatment group can be determined and displayed.

anyNA(biofuel_long$count)
## [1] FALSE
biofuel_tibble <- biofuel_long %>%
  group_by(treatment) %>%
  summarize(
    mean(count),
    median(count),
    n = n()
  )

biofuel_tibble
## # A tibble: 12 × 4
##    treatment `mean(count)` `median(count)`     n
##    <chr>             <dbl>           <dbl> <int>
##  1 INICIA             165.            37.1  6071
##  2 INICIA.1           165.            20.7  6071
##  3 INICIB             165.            43.2  6071
##  4 INICIC             165.            21.8  6071
##  5 INICIF             165.            50.6  6071
##  6 INICII             165.            20.3  6071
##  7 INICIN             165.            35.6  6071
##  8 INICIQ             165.            38.8  6071
##  9 INICIR             165.            37.2  6071
## 10 INICIS             165.            23.0  6071
## 11 INICSA             165.            41.8  6071
## 12 INICSR             165.            38.3  6071

For this plot, a cutoff of 3 SDs from the mean will be applied to remove highly expressed outliers. This will include approximately 99.7% of the original data, cutting off abnormal outliers.

biofuel_plot <- biofuel_long %>%
  filter(count <= (mean(count)) + 3 * sd(count))

A violin plot can then be put together to compare the median TPM of each group and the overall distribution of the data.

ggplot(biofuel_plot,
       aes(x = treatment, y = count)) +
  labs(x = 'Treatment Group',
       y = 'Count (TPM)',
       title = 'RNA Sequencing Counts by Treatment Groups',
       subtitle = 'RNA Expresion Levels of Biofuel Producing Yeast') +
  theme_bw() +
  geom_violin()