Course Description:

This course will focus on utilizing the power of ggplot2 in R for analyzing genomes and comparing genomic data. Participants will learn how to create customized, publication-quality graphics for visualizing various genomic features, such as gene expression, sequence motifs, and genomic variations. Through hands-on exercises and real-world examples, participants will gain practical skills in data visualization and interpretation in the context of genomics.

Course Objectives:

• Understand the principles of ggplot2 and its grammar of graphics.
• Learn to create different types of plots for genome analysis using ggplot2.
• Gain proficiency in customizing plots to effectively visualize genomic data.
• Explore advanced techniques for comparing genomes and visualizing genomic variations.
• Develop skills in interpreting and presenting genomic data using ggplot2.

Prerequisites:

• Basic knowledge of R programming.
• Familiarity with genomic data formats (e.g., BED, VCF, FASTA) is helpful but not required.

Target Audience:

• Bioinformaticians, biologists, and researchers working with genomic data.
• Students and professionals interested in advanced data visualization and data manipulation techniques using ggplot2, dplyr, and the tidy data framework in R.

Introduction to ggplot2

ggplot2 is a popular R package for data visualization, known for its flexibility and powerful capabilities. It is particularly useful for genomic data analysis due to its ability to create complex and customized plots with relatively simple code. Here’s an overview of ggplot2 and its advantages for genomic data analysis:

1. Grammar of Graphics: ggplot2 is based on the Grammar of Graphics, which provides a consistent framework for constructing plots. This makes it easier to understand and modify plots, especially when dealing with complex datasets.
2. Layered Approach: ggplot2 uses a layered approach to plot construction, where different elements of a plot (such as data points, lines, and labels) are added as separate layers. This allows for easy customization and modification of plots.
3. Flexibility: ggplot2 offers a high level of flexibility, allowing users to create a wide variety of plots, including scatter plots, bar plots, histograms, and more. This flexibility is particularly useful for genomic data analysis, where different types of plots may be needed to visualize different aspects of the data.
4. Integration with tidyverse: ggplot2 is part of the tidyverse ecosystem, which includes a suite of packages for data manipulation and analysis. This integration makes it easy to use ggplot2 with other tidyverse packages, such as dplyr and tidyr, for seamless data analysis workflows.
5. Themes and Scales: ggplot2 provides a wide range of themes and scales for customizing the appearance of plots. This includes options for changing colors, fonts, and other visual elements, as well as scales for transforming data (e.g., log scales).
6. Publication-Quality Plots: ggplot2 is designed to produce high-quality plots suitable for publication. This includes options for exporting plots in various formats (e.g., PDF, PNG) at different resolutions.

Overall, ggplot2 is a powerful tool for genomic data analysis, thanks to its flexibility, ease of use, and ability to create publication-quality plots. It is widely used in the bioinformatics and genomics communities for visualizing complex datasets and exploring relationships within genomic data.

Basic principles of the grammar of graphics

The Grammar of Graphics is a framework for creating and understanding data visualizations, introduced by Leland Wilkinson in his book “The Grammar of Graphics.” It provides a structured way to think about how plots are constructed, allowing for a more systematic and intuitive approach to data visualization. The basic principles of the Grammar of Graphics include:

1. Data: The starting point for creating a plot is the data itself. The data is typically organized in a tabular format, with rows representing individual observations and columns representing variables.
2. Aesthetic Mapping: Aesthetic mappings define how variables in the data are mapped to visual properties of the plot, such as position, color, size, and shape. For example, mapping a variable to the x-axis position will create a scatter plot, while mapping it to the color will create a colored scatter plot.
3. Geometric Objects (Geoms): Geometric objects represent the actual graphical elements that are drawn on the plot, such as points, lines, bars, and polygons. Geoms are used to visually represent the data in different ways.
4. Scales: Scales control how data values are mapped to visual properties, such as mapping a range of data values to a range of colors. Scales can also transform data values (e.g., log scale) and provide legends and axis labels.
5. Facets: Faceting involves creating multiple plots (or panels) based on the values of one or more variables. This is useful for visualizing relationships within subgroups of the data.
6. Statistical Transformations (Stats): Statistical transformations are used to summarize or transform the data before plotting. For example, calculating a mean or median for a group of data points before plotting.
7. Coordinate System: The coordinate system defines the space in which the plot is drawn. Common coordinate systems include Cartesian coordinates (x and y axes) and polar coordinates (radius and angle).

By following the principles of the Grammar of Graphics, it becomes easier to create complex and customized plots, as well as to understand and modify existing plots. The ggplot2 package in R is based on the Grammar of Graphics and provides a powerful implementation of these principles for creating data visualizations.

Creating Basic Plots

Scatter plots, line plots, and bar plots for visualizing genomic data

Scatter plots, line plots, and bar plots are commonly used to visualize genomic data. Here’s how each type of plot can be used:

1. Scatter plots: Scatter plots are used to visualize the relationship between two continuous variables. In genomic data, scatter plots can be used to show the relationship between two different genomic features, such as gene expression levels or methylation levels at different genomic loci. Each point in the scatter plot represents a data point, and the position of the point corresponds to the values of the two variables being compared.
2. Line plots: Line plots are used to visualize changes in a variable over time or across different conditions. In genomic data, line plots can be used to show how gene expression levels change across different experimental conditions or how methylation levels change along a genomic region. Each line in the plot represents the trend for a specific gene or genomic region, with the x-axis representing the position along the genome and the y-axis representing the value of the variable being measured.
3. Bar plots: Bar plots are used to compare the values of a variable across different categories. In genomic data, bar plots can be used to compare the expression levels of different genes or the methylation levels of different genomic regions. Each bar in the plot represents the value of the variable for a specific category, and the height of the bar represents the magnitude of the value.

These types of plots can be created using the ggplot2 package in R, which provides a flexible and powerful framework for creating customized plots. By using these plots, researchers can gain insights into the patterns and trends present in genomic data, helping to further our understanding of the underlying biological processes.

Adding layers, annotations, and customizing aesthetics

In ggplot2, you can add layers, annotations, and customize aesthetics to enhance your plots. Here’s a brief overview of how to do this:

1. Adding Layers: You can add layers to a ggplot by using the + operator. Each layer can represent different aspects of your data or additional elements you want to include in the plot. For example, to add a scatter plot layer to a ggplot, you can use the geom_point() function:
   RCopy code

   ggplot(data = my_data, aes(x = x_variable, y = y_variable)) +
geom_point()

2. Annotations: Annotations are used to add text, labels, or geometric shapes to your plot to provide additional information. You can add annotations using functions like geom_text(), geom_label(), or geom_rect(). For example, to add a text annotation to a plot:
   RCopy code

   ggplot(data = my_data, aes(x = x_variable, y = y_variable)) +
geom_point() +
geom_text(aes(label = label_variable), vjust = -1)

3. Customizing Aesthetics: You can customize the aesthetics of your plot, such as colors, shapes, and sizes, using the aes() function within each geom_ layer. Additionally, you can use the theme() function to customize the overall appearance of the plot, including axis labels, titles, and grid lines. For example, to change the color and shape of points in a scatter plot:
   RCopy code

   ggplot(data = my_data, aes(x = x_variable, y = y_variable)) +
geom_point(aes(color = group_variable, shape = shape_variable), size = 3) +
scale_color_manual(values = c("red", "blue", "green")) +
theme_minimal()


These are just a few examples of how you can add layers, annotations, and customize aesthetics in ggplot2. The package offers a wide range of functions and options for creating highly customized and informative plots for your genomic data analysis.

Visualizing Gene Expression

Creating box plots, violin plots, and heatmaps for gene expression analysis

To create box plots, violin plots, and heatmaps for gene expression analysis in R using ggplot2 and other relevant packages, you can follow these examples:

1. Box plots:

2. Violin plots:

3. Heatmaps:

Make sure to replace 'data' with your actual data frame or matrix containing gene expression values. These examples should help you get started with creating these types of plots for gene expression analysis.

Using facetting and grouping to compare gene expression across samples

Facetting and grouping are powerful techniques in ggplot2 that allow you to compare gene expression across samples. Facetting creates small multiples of a plot, each showing a subset of the data, while grouping allows you to overlay multiple groups within the same plot. Here’s how you can use facetting and grouping to compare gene expression across samples:

1. Facetting:
   RCopy code

   ggplot(data, aes(x = gene, y = expression)) +
geom_boxplot() +
facet_wrap(~sample_group, scales = "free_x")


   This code creates a box plot for each gene, with the plots grouped by sample_group. The scales = "free_x" argument ensures that each facet has its own x-axis scale.

2. Grouping:
   RCopy code

   ggplot(data, aes(x = gene, y = expression, color = sample_group)) +
geom_point() +
scale_color_manual(values = c("red", "blue", "green")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))


   This code creates a scatter plot of gene expression, with each sample group represented by a different color. The scale_color_manual() function allows you to specify the colors for each group.

Combining facetting and grouping can provide a comprehensive view of gene expression patterns across different sample groups, helping you to identify trends and patterns in your data.

Sequence Motif Analysis

Plotting sequence logos and motif enrichment plots

To plot sequence logos and motif enrichment plots, you can use the ggseqlogo and ggplot2 packages in R, respectively. Here’s how you can create these plots:

1. Sequence Logos with ggseqlogo: First, install and load the ggseqlogo package:
   RCopy code

   install.packages("ggseqlogo")
library(ggseqlogo)


   Then, create a sequence logo plot:

   RCopy code

   # Assuming 'alignment' is a DNA sequence alignment object
ggseqlogo(logo_data = as.list(alignment))


   Replace alignment with your actual DNA sequence alignment object.

2. Motif Enrichment Plot with ggplot2: First, prepare your data. The motif column should contain the motif names, and the enrichment column should contain the enrichment scores.
   RCopy code

   # Example data
motif_data <- data.frame(
motif = c("motif1", "motif2", "motif3"),
enrichment = c(1.5, 2.0, 1.2)
)


   Then, create a motif enrichment plot:

   RCopy code

   library(ggplot2)

   ggplot(motif_data, aes(x = motif, y = enrichment)) +
geom_bar(stat = "identity") +
labs(x = "Motif", y = "Enrichment Score") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))


   Customize the plot as needed for your data and analysis.

These examples should help you get started with plotting sequence logos and motif enrichment plots in R.

Customizing plots for visualizing sequence motifs in genomic sequences

To customize plots for visualizing sequence motifs in genomic sequences, you can use tools like ggplot2 in R or matplotlib in Python. Here’s a basic example of how you can create a sequence motif plot and customize it:

R (using ggplot2)

This code creates a plot showing two sequence motifs (motif_A and motif_B) across different positions in the sequence. You can customize the plot further by adjusting colors, line types, axis labels, and other aesthetics to fit your specific needs.

Comparing Genomes

Creating genomic alignment plots using ggplot2

Creating genomic alignment plots using ggplot2 can be done by first preparing the alignment data and then using ggplot2 to visualize it. Here’s a basic example using synthetic data:

R (using ggplot2)

This code creates a genomic alignment plot where each nucleotide is represented by a colored tile. Each row represents a sequence, and each column represents a position in the alignment. You can customize the plot further by adjusting colors, labels, and other aesthetics to fit your specific needs.

For more complex genomic alignment plots, you may need to preprocess your alignment data to extract relevant information such as conservation scores, gaps, or additional annotations. Then, you can use ggplot2 to visualize the data in a meaningful way.

Visualizing synteny and genomic rearrangements

Visualizing synteny and genomic rearrangements can be done using tools like Circos, which is specifically designed for this purpose. Circos allows you to create circular plots that represent relationships between different genomic regions. Here’s a basic example of how you can use Circos to visualize synteny and genomic rearrangements:

1. Install Circos: First, you’ll need to download and install Circos from the official website: Circos – Circular Genome Data Visualization.
2. Prepare Data: Prepare your genomic data in a format that Circos can read. This typically involves creating configuration files that specify the positions of genes or genomic regions along with any synteny or rearrangement information.
3. Create Configuration Files: Circos uses configuration files to define the layout and appearance of the plot. You’ll need to create a configuration file that specifies the genomic regions to plot and any synteny or rearrangement information.
4. Run Circos: Once you have your data and configuration files ready, you can run Circos to generate the plot. Circos will create a circular plot that visualizes the synteny and genomic rearrangements in your data.

Here’s a basic example of a Circos configuration file for visualizing synteny:

In this example, links.txt contains the synteny or rearrangement information between genomic regions. You can customize the appearance of the plot by modifying the configuration file to suit your needs.

Please refer to the Circos documentation for more information on how to create and customize Circos plots for visualizing synteny and genomic rearrangements.

Genomic Variations

Plotting variant allele frequency and variant annotation plots

To plot variant allele frequency (VAF) and variant annotation plots, you can use R and ggplot2. Here’s a basic example for each:

Variant Allele Frequency Plot

This code creates a plot showing the VAF across different positions in the genome. You can customize the plot further by adjusting colors, line types, labels, and other aesthetics to fit your specific needs.

Variant Annotation Plot

For variant annotation plots, you can use tools like ggplot2 to visualize the annotations associated with each variant. Here’s a basic example:

This code creates a bar plot showing the count of variants for each annotation type across different positions in the genome. You can customize the plot further by adjusting colors, labels, and other aesthetics to fit your specific needs.

These examples should help you get started with plotting variant allele frequency and variant annotation plots in R using ggplot2.

Using gggenomes plot for visualizing genomic features across multiple genomes

To use the gggenomes package in R for visualizing genomic features across multiple genomes, you first need to install and load the package. Then, you can use the gggenome function to create the plot. Here’s a basic example:

1. Install and load gggenomes:

2. Prepare your data:
   • Create a data frame containing the genomic features of interest for each genome.
   • The data frame should include columns for genome, feature, start, and end coordinates of the feature.
3. Create the plot: Use the gggenome function to create the plot. Specify the genome, feature, start, and end columns from your data frame, along with any additional customization options.

Here’s a basic example using synthetic data:

This code creates a genomic feature plot using gggenomes, with different features (e.g., genes, exons, promoters) represented across multiple genomes. You can customize the plot further by adjusting colors, labels, and other aesthetics to fit your specific needs.

For more advanced usage and customization options, refer to the gggenomes documentation and examples.

Advanced Customization Techniques

Custom themes, scales, and color palettes for enhancing plot aesthetics

Customizing themes, scales, and color palettes can greatly enhance the aesthetics of your plots in ggplot2. Here’s how you can customize these elements:

Custom Theme

You can create a custom theme using the theme() function. For example, to create a minimalistic theme with a white background and no gridlines:

Custom Scales

You can create custom scales using the scale_ functions. For example, to customize the color scale:

Custom Color Palettes

You can create custom color palettes using the scale_fill_manual() or scale_color_manual() functions. For example, to use a custom color palette for fill colors:

These are just a few examples of how you can customize themes, scales, and color palettes in ggplot2. Experiment with different options to create plots that best suit your data and visualization needs.

Animating genomic data using gganimate

To animate genomic data using gganimate, you can create a series of plots that represent different frames of the animation, each showing a snapshot of the data at a specific time point. Here’s a basic example using synthetic data:

1. Install and load packages:
   RCopy code

   install.packages("ggplot2")
install.packages("gganimate")
library(ggplot2)
library(gganimate)

2. Prepare your data: Create a data frame containing the genomic data to be animated. The data frame should include columns for position, value, and any other variables needed for the plot.
3. Create the animated plot: Use ggplot to create the initial plot and transition_states or transition_time to create the animation. Here’s an example for animating genomic data across positions:
   RCopy code

   # Example data
   genomic_data <- data.frame(
   position = 1:10,
   value = runif(10)
   )

   # Create initial plot
   p <- ggplot(data = genomic_data, aes(x = position, y = value)) +
   geom_line() +
   labs(title = 'Genomic Data Animation') +
   transition_states(states = position, transition_length = 2, state_length = 1)

   # Animate the plot
animate(p, nframes = 100, fps = 10)


   In this example, transition_states is used to create an animation where each frame represents a different position in the genome.

4. Save the animation: You can save the animation as a GIF or video file using the anim_save function:
   RCopy code

   anim_save("genomic_animation.gif")


This is a basic example of how you can use gganimate to animate genomic data. You can customize the plot and animation settings to fit your specific needs and data format.

Overview of ggplot2, dplyr, and the tidy data framework

ggplot2, dplyr, and the tidy data framework are three key components of the tidyverse in R that work together to facilitate data manipulation and visualization. Here’s an overview of each:

1. ggplot2:
   • Purpose: ggplot2 is a powerful data visualization package in R that implements the grammar of graphics, allowing you to create a wide variety of plots with a consistent and intuitive syntax.
   • Key Concepts:
     • Aesthetics: Mapping variables in your data to visual properties of the plot (e.g., x and y axes, color, size).
     • Geoms: Geometric objects that represent the data in the plot (e.g., points, lines, bars).
     • Facets: Creating small multiples of a plot to show subsets of the data.
     • Themes: Customizing the appearance of the plot (e.g., axis labels, gridlines).
   • Example:
     RCopy code

     library(ggplot2)
ggplot(data = iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species)) +
geom_point() +
labs(title = "Scatter Plot of Iris Dataset", x = "Sepal Length", y = "Sepal Width")

2. dplyr:
   • Purpose: dplyr is a package for data manipulation that provides a set of functions for common data manipulation tasks, such as filtering, sorting, grouping, and summarizing data.
   • Key Functions:
     • filter(): Subset rows of a data frame based on conditions.
     • arrange(): Reorder rows of a data frame.
     • select(): Select columns of a data frame.
     • mutate(): Add new variables or modify existing variables.
     • group_by(): Group data by one or more variables.
     • summarize(): Summarize grouped data (e.g., calculate means, counts).
   • Example:
     RCopy code

     library(dplyr)
iris_filtered <- iris %>%
filter(Sepal.Length > 5.0) %>%
select(Species, Sepal.Length, Sepal.Width)

3. Tidy Data Framework:
   • Purpose: The tidy data framework is an approach to organizing data that makes it easy to work with in R. In tidy data:
     • Each variable is in a column.
     • Each observation is in a row.
     • Each type of observational unit forms a table.
   • Benefits:
     • Makes data manipulation and visualization easier.
     • Facilitates interoperability between tidyverse packages.
   • Example:
     mathematicaCopy code

     | Species | Sepal.Length | Sepal.Width | Petal.Length | Petal.Width |
|---------|--------------|-------------|--------------|-------------|
| setosa | 5.1 | 3.5 | 1.4 | 0.2 |
| setosa | 4.9 | 3.0 | 1.4 | 0.2 |
| ... | ... | ... | ... | ... |


Together, ggplot2, dplyr, and the tidy data framework provide a powerful toolkit for data manipulation and visualization in R, making it easier to explore and analyze data.

Data Manipulation with dplyr

o filter, sort, and summarize genomic data using dplyr, you can use the functions filter(), arrange(), group_by(), and summarize(). Here’s a basic example using synthetic genomic data:

Example Data

Let’s assume you have a data frame called genomic_data with columns gene, position, mutation, and effect.

Filtering

To filter the genomic data to include only mutations that have a certain effect:

This code filters the genomic_data data frame to include only rows where the effect column is equal to “missense_variant”.

Sorting

To sort the genomic data by position:

This code sorts the genomic_data data frame by the position column in ascending order.

Summarizing

To summarize the genomic data to calculate the average position of mutations for each gene:

This code groups the genomic_data data frame by the gene column and then calculates the average position of mutations for each gene.

These are just a few examples of how you can use dplyr to filter, sort, and summarize genomic data in R. Depending on your specific needs, you can combine these functions in different ways to perform more complex data manipulations.

Joining and merging genomic datasets

To join and merge genomic datasets in R, you can use the dplyr package along with its left_join(), right_join(), inner_join(), full_join(), and bind_rows() functions. Here’s a basic overview of how to do this:

Example Data

Let’s assume you have two data frames, genomic_data1 and genomic_data2, with columns gene, position, and mutation.

Joining Datasets

1. Inner Join: Keep only rows that have matching values in both datasets.
   RCopy code

   inner_merged_data <- inner_join(genomic_data1, genomic_data2, by = "gene")

2. Left Join: Keep all rows from the left dataset and matching rows from the right dataset.
   RCopy code

   left_merged_data <- left_join(genomic_data1, genomic_data2, by = "gene")

3. Right Join: Keep all rows from the right dataset and matching rows from the left dataset.
   RCopy code

   right_merged_data <- right_join(genomic_data1, genomic_data2, by = "gene")

4. Full Join: Keep all rows from both datasets, filling in missing values with NA.
   RCopy code

   full_merged_data <- full_join(genomic_data1, genomic_data2, by = "gene")


Merging Datasets

1. Bind Rows: Combine rows from two datasets into a single dataset.
   RCopy code

   merged_data <- bind_rows(genomic_data1, genomic_data2)


These functions allow you to combine genomic datasets based on a common column (e.g., gene) or combine all rows from multiple datasets. Adjust the by argument in the join functions to specify the column(s) to join on.

Tidy Data Organization

Principles of tidy data and organizing genomic data into a tidy format

The principles of tidy data, as introduced by Hadley Wickham, are designed to make data cleaning, manipulation, and analysis easier and more intuitive. These principles can be applied to organizing genomic data into a tidy format, which is particularly useful for data analysis using tools like dplyr and ggplot2 in R. Here are the key principles of tidy data and how they can be applied to genomic data:

1. Each variable forms a column: In tidy data, each variable is placed in its own column. For genomic data, this means that each attribute (e.g., gene name, position, mutation type) should be in a separate column.
2. Each observation forms a row: Each row in the dataset represents a single observation. In genomic data, each row could represent a specific gene or genomic region.
3. Each type of observational unit forms a table: Each dataset should represent a single type of observational unit. For genomic data, this could mean having separate tables for gene annotations, mutation data, and expression levels.

Applying these principles to genomic data can involve restructuring the data to ensure that each variable is in its own column, each observation is in its own row, and each type of genomic data (e.g., genes, mutations, expression levels) is organized into separate tables or data frames. Here’s an example of how genomic data could be organized into a tidy format:

By organizing genomic data into a tidy format, you can easily perform data manipulation, analysis, and visualization using tools like dplyr and ggplot2, making your analysis more efficient and reproducible.

Benefits of tidy data for data analysis and visualization

Tidy data offers several benefits for data analysis and visualization, especially when using tools like dplyr and ggplot2 in R. Some key benefits include:

1. Simplicity and readability: Tidy data makes datasets easier to understand and work with, as each variable is in a separate column and each observation is in a separate row. This makes it easier to interpret data and share findings with others.
2. Ease of data manipulation: Tidy data simplifies data manipulation tasks such as filtering, sorting, and summarizing. Tools like dplyr are designed to work with tidy data, allowing for more efficient and intuitive data manipulation.
3. Compatibility with tidyverse tools: Tidy data is compatible with the tidyverse ecosystem of packages in R, which includes tools like dplyr, ggplot2, and tidyr. These packages are designed to work together seamlessly, making it easier to perform complex data analyses and visualizations.
4. Facilitates data visualization: Tidy data is well-suited for visualization using ggplot2, which is designed to work with tidy data. This allows for the creation of complex and informative plots with minimal code.
5. Improved reproducibility: Tidy data and tidyverse tools promote a workflow that is more reproducible. By organizing data in a consistent and standardized format, it is easier to reproduce analyses and share code with others.

Overall, tidy data offers a more organized and efficient approach to data analysis and visualization, making it easier to explore, understand, and communicate insights from your data.

Case Studies and Applications

Cleaning and visualizing genomic data: a case study in tidy analysis

Assuming that the reader is familiar with the absolute basics of the dplyr and ggplot2 packages, especially the

operator. (If not, here’s a good set of tutorials).

%>%

Background: gene expression in starvation

Through the process of gene regulation, a cell can control which genes are transcribed from DNA to RNA- what we call being “expressed”. (If a gene is never turned into RNA, it may as well not be there at all). This provides a sort of “cellular switchboard” that can activate some systems and deactivate others, which can speed up or slow down growth, switch what nutrients are transported into or out of the cell, and respond to other stimuli. A gene expression microarray lets us measure how much of each gene is expressed in a particular condition. We can use this to figure out the function of a specific gene (based on when it turns on and off), or to get an overall picture of the cell’s activity.

Brauer 2008 used microarrays to test the effect of starvation and growth rate on baker’s yeast (S. cerevisiae, a popular model organism for studying molecular genomics because of its simplicity)1. Basically, if you give yeast plenty of nutrients (a rich media), except that you sharply restrict its supply of one nutrient, you can

control the growth rate to whatever level you desire (we do this with a tool called a chemostat). For example, you could limit the yeast’s supply of glucose (sugar, which the cell metabolizes to get energy and carbon), of leucine (an essential amino acid), or of ammonium (a source of nitrogen).

“Starving” the yeast of these nutrients lets us find genes that:

Raise or lower their activity in response to growth rate. Growth-rate dependent expression patterns can tell us a lot about cell cycle control, and how the cell responds to stress.

Respond differently when different nutrients are being limited. These

genes may be involved in the transport or metabolism of those nutrients.

Sounds pretty cool, right? So let’s get started!

Tidying data with dplyr and tidyr

You can download the original gene expression dataset here, or just run the following R code (you’ll need the readr package installed)2:

delim(“http://varianceexplained.org/files/Brauer2008_DataSet1.tds”, delim = “\t”)

This is the exact data that was published with the paper (though for some reason the link on the journal’s page is now broken). It thus serves as a good example of tidying a biological dataset “found in the wild.”

Let’s take a look at the data:

View(original_data)

dim(original_data)

## [1] 5537 40

Each of those columns like , and so on represents gene expression

G0.05

N0.3

values for that sample, as measured by the microarray. The column titles show the

condition: , for instance, means the limiting nutrient was glucose and the

G0.05

growth rate was .05. A higher value means the gene was more expressed in that sample, lower means the gene was less expressed. In total the yeast was grown with six limiting nutrients and six growth rates, which makes 36 samples, and therefore 36 columns, of gene expression data.

Before we do any analysis, we’re going to need to tidy this data, as described in Wickham’s “Tidy Data” paper (especially p5-13). Tidy data follows these rules:

1. Each variable forms a column.
2. Each observation forms a row.
3. Each type of observational unit forms a table.

What is “untidy” about this dataset?

Column headers are values, not variable names. Our column names contain the values of two variables: nutrient (G, N, P, etc) and growth rate (0.05-0.3). For this reason, we end up with not one observation per row, but 36! This is a very common issue in biological datasets: you often see one-row- per-gene and one-column-per-sample, rather than one-row-per-gene-per- sample.

Multiple variables are stored in one column. The column contains

NAME

||

lots of information, split up by like

’s. If we examine one of the names, it looks

SFB2 || ER to Golgi transport || molecular function unknown || YNL049C || 1082129

which seems to have both some systematic IDs and some biological information about the gene. If we’re going to use this programmatically, we need to split up the information into multiple columns.

The more effort you put up front into tidying your data, the easier it will be to explore interactively. Since the analysis steps are where you’ll actually be answering questions about your data, it’s worth putting up this effort!

Multiple variables are stored in one column

Let’s take a closer look at that NAME column.

original_data$NAME[1:3]

## [1] “SFB2 || ER to Golgi transport || molecular function unknown || YNL04 ## [2] “|| biological process unknown || molecular function unknown || YNL095C ||

## [3] “QRI7 || proteolysis and peptidolysis || metalloendopeptidase activit

The details of each of these fields isn’t annotated in the paper, but we can figure out

m

ost of it. It contains:

Gene name e.g. SFB2. Note that not all genes have a name. Biological process e.g. “proteolysis and peptidolysis” Molecular function e.g. “metalloendopeptidase activity”

Systematic ID e.g. YNL049C. Unlike a gene name, every gene in this dataset has a systematic ID.3

Another ID number e.g. 1082129. I don’t know what this number means, and it’s not annotated in the paper. Oh, well.

Having all give of these in the same column is very inconvenient. For example, if I have another dataset with information about each gene, I can’t merge the two.

Luckily, the tidyr package provides the function for exactly this case.

separate

library(dplyr) library(tidyr)

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\|

We simply told what column we’re splitting, what the new column names

separate

should be, and what we were separating by. What does the data look like now?

View(cleaned_data)

Just like that, we’ve separated one column into five.

Two more things. First, when we split by , we ended up with whitespace at

||

the start and end of some of the columns, which is inconvenient:

head(cleaned_data$BP)

## [1] ” ER to Golgi transport ” ” biological process unknown ” ## [3] ” proteolysis and peptidolysis ” ” mRNA polyadenylylation* ”

## [5] ” vesicle fusion* ” ” biological process unknown ”

We’ll solve that with dplyr’s whitespace”) function.

, along with the built-in

mutate_each

(“trim

trimws

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\| mutate_each(funs(trimws), name:systematic_name)

Finally, we don’t even know what the column represents (if you can figure

number

it out, let me know!) And while we’re at it, we’re not going to use the , or

GID

YORF

columns in this analysis either. We may as well drop them, which we can

GWEIGHT

do with dplyr’s .

select

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\| mutate_each(funs(trimws), name:systematic_name) %>%

select(–number, –GID, –YORF, –GWEIGHT)

Column headers are values, not variable names

Let’s take a closer look at all those column headers like , and .

G0.05

N0.2

P.15

Limiting nutrient. This has six possible values: glucose (G), ammonium (N), sulfate (S), phosphate (P), uracil (U) or leucine (L).

Growth rate: A number, ranging from .05 to .3. .05 means slow growth (the yeast were being starved hard of that nutrient) while .3 means fast growth. (Technically, this value measures the dilution rate from the chemostat).

Expression level. These are the values currently stored in those columns, as measured by the microarray. (Note that the paper already did some pre- processing and normalization on these values, which we’re ignoring here).

The rules of tidy data specify that each variable forms one column, and this is not even remotely the case- we have 36 columns when we should have 3. That means our data is trapped in our column names. If you don’t see why this is a problem, consider: how would you put growth rate on the x-axis of a graph? How would you filter to look only the glucose condition?

Luckily, the tidyr package has a solution ready. The documentation for notes (emphasis mine):

gather

Gather takes multiple columns and collapses into key-value pairs, duplicating all other columns as needed. You use

gather()

when you notice that you have

columns that are not variables.

Hey, that’s us! So let’s apply as our next step:

gather

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\| mutate_each(funs(trimws), name:systematic_name) %>%

select(–number, –GID, –YORF, –GWEIGHT) %>%

gather(sample, expression, G0.05:U0.3)

Now when we view it, it looks like:

View(cleaned_data)

Notice that the dataset no longer consists of one-row-per-gene: it’s one-row- per-gene-per-sample. This has previously been called “melting” a dataset, or turning it into “long” format. But I like the term “gather”: it shows that we’re taking these 36 columns and pulling them together.

One last problem. That column really contains two variables,

sample

nutrient

and . We already learned what to do when we have two variables in one

rate

column: use :

separate

library(dplyr) library(tidyr) library(stringr)

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\| mutate_each(funs(trimws), name:systematic_name) %>%

select(–number, –GID, –YORF, –GWEIGHT) %>%

gather(sample, expression, G0.05:U0.3) %>%

separate(sample, c(“nutrient”, “rate”), sep = 1, convert = TRUE)

This time, instead of telling to split the strings based on a particular

separate

delimiter, we told it to separate it after the first character (that is, after

G/P/S/N/L/U ). We also told it convert = TRUE to tell it that it should notice the

0

05/0.1/etc value is a number and convert it.

.

Take a look at those six lines of code, a mini-sonnet of data cleaning. Doesn’t it read less like code and more like instructions? (“First we separated the NAME column into its five parts, and trimmed each. We selected out columns we didn’t

need…”) That’s the beauty of the operator and the dplyr/tidyr verbs.

%>%

Why keep this data tidy? Visualization with ggplot2

So we went through this effort to get this dataset into this structure, and you’re probably wondering why. In particular, why did we have to bother gathering those expression columns into one-row-per-gene-per-sample?

Well, suppose we have a single yeast gene we’re interested in. Let’s say LEU1, a gene involved in the leucine synthesis pathway.

cleaned_data %>%

filter(name == “LEU1”) %>%

View()

We now have 36 data points (six conditions, six growth rates), and for each we have a limiting nutrient, a growth rate, and the resulting expresion. We’re probably interested in how both the growth rate and the limiting nutrient affect the gene’s expression.

36 points is too many to look at manually. So it’s time to bring in some visualization. To do that, we simply pipe the results of our filtering right into ggplot2:

library(ggplot2) cleaned_data %>%

filter(name == “LEU1”) %>%

ggplot(aes(rate, expression, color = nutrient)) +

geom_line()

What a story this single gene tells! The gene’s expression is far higher (more “turned on”) when the cell is being starved of leucine than in any other condition, because in that case the cell has to synthesize its own leucine. And as the amount of leucine in the environment (the growth rate) increases, the cell can focus less on leucine production, and the expression of those genes go down. We’ve just gotten one snapshot of our gene’s regulatory network, and how it responds to external stimuli.

We don’t have to choose one gene to visualize- LEU1 is just one gene in the leucine biosynthesis process. Recall that we have that information in the column, so we can filter for all genes in that process, and then facet to create sub- plots for each.

BP

cleaned_data %>%

filter(BP == “leucine biosynthesis”) %>% ggplot(aes(rate, expression, color = nutrient)) + geom_line() +

facet_wrap(~name)

LEU1, LEU2, and LEU4 all show a similar pattern, where starvation of leucine causes higher gene expression. (Interestingly, LEU4 responds to glucose starvation as well. Any geneticists have some idea why?). LEU9 is a little more ambiguous but is still highest expressed under leucine starvation. We already know what these genes do, but this hints at how we might be able to find other genes that are involved in leucine synthesis, including ones we don’t yet know.

Let’s play with graph a little more. These trends look vaguely linear. Maybe we should show best points with best fit lines instead:

cleaned_data %>%

filter(BP == “leucine biosynthesis”) %>% ggplot(aes(rate, expression, color = nutrient)) + geom_point() +

geom_smooth(method = “lm”, se = FALSE) +

facet_wrap(~name)

The options for exploratory data analysis are endless. We could instead look at sulfur metabolism:

cleaned_data %>%

filter(BP == “sulfur metabolism”) %>% ggplot(aes(rate, expression, color = nutrient)) + geom_point() +

geom_smooth(method = “lm”, se = FALSE) +

facet_wrap(~name + systematic_name, scales = “free_y”)

(Notice that we have to facet by the this process have traditional names).

systematic_name

here, since not all genes in

If you were an interested molecular biologist, you could go a long way just by examining various biological processes, or looking at lists of genes you were interested in. This is a great way to explore a dataset interactively, and to hint at methods of further analysis. For example, we notice that we can fit linear models to many of these expression profiles- that will come in handy in the next post:

Conclusion: specific data, generalized tools

I earlier pointed to a list of available workflows from Bioconductor, which teach ways to analyze many kinds of genomic data using packages specialized for those purposes. These kinds of guides are incredibly valuable, and Bioconductor has built up an excellent set of tools for analyzing genomic data, many of which come with their own data processing and visualization methods.

So why bother teaching a dplyr/ggplot2 approach? Because these tools are

useful everywhere. Consider the dozen lines of code we used to clean and visualize our data:

library(dplyr) library(tidyr) library(ggplot2)

cleaned_data <- original_data %>%

separate(NAME, c(“name”, “BP”, “MF”, “systematic_name”, “number”), sep = “\\|\\| mutate_each(funs(trimws), name:systematic_name) %>%

select(–number, –GID, –YORF, –GWEIGHT) %>%

gather(sample, expression, G0.05:U0.3) %>%

separate(sample, c(“nutrient”, “rate”), sep = 1, convert = TRUE)

cleaned_data %>%

filter(BP == “leucine biosynthesis”) %>% ggplot(aes(rate, expression, color = nutrient)) + geom_point() +

geom_smooth(method = “lm”, se = FALSE) +

facet_wrap(~name)

With this code we’re able to go from published data to a new conclusion about leucine biosynthesis (or sulfate metabolism, or whatever biological process or genes we’re interested in). But the functions we use aren’t specific to our dataset or data

fo

rmat: in fact, they’re not related to biology at all. Instead, these tools are

building blocks, or “atoms,” in a grammar of data manipulation and visualization that applies to nearly every kind of data.

In June 2015 I finished my PhD in computational biology and joined Stack Overflow as a data scientist. People sometimes ask me whether it was difficult to switch my focus away from biology. But I found that my experience was almost equally applicable to data science on a website, because I’d been working with dplyr, ggplot2 and other tidy tools and performing analyses like the one in this post. If I’d instead limited myself to specialized tools that looked something like:

library(brauerAnalyze)

cleaned_data <- read_brauer_data(“http://varianceexplained.org/files/Brauer2008_Da plot_biological_process(cleaned_data, “leucine biosynthesis”)

Then I might have gotten that task done in fewer lines of code, but I wouldn’t have learned to approach new datasets in the future.

This isn’t meant to disparage Bioconductor in any way: scientists should use whatever tool gets their job done. But educators can and should focus on teaching tools that can be universally applied. In turn, students can take these tools and build new packages, for Bioconductor and elsewhere, that analyze novel forms of data.

Next time

You’ll probably agree that the above analysis gives a compelling case for keeping your data in a tidy format (one-row-per-gene-per-sample) when it comes to making plots. (Just try creating the above “leucine biosynthesis” plot without first tidying the data- it is a hassle, and especially impossible in ggplot2).

But graphing is nothing special. Almost everything– from exploratory analysis to modeling to parallelization- is easier when your data is structured according to tidy data principles. In future posts, I’ll show how to perform per-gene linear regressions

with broom, use GO term enrichment to find functional groups of genes that show differential expression, and many other tidy analyses.

Tips and best practices for effective data visualization and data manipulation in genomics

Effective data visualization and manipulation are crucial in genomics for gaining insights from complex biological data. Here are some tips and best practices:

Data Visualization:

1. Choose the right plot: Select a plot type that best represents your data and highlights the patterns or trends you want to show (e.g., scatter plots, bar plots, heatmaps).
2. Use informative labels: Ensure that axis labels, legends, and titles are clear and descriptive, making it easy for viewers to understand the plot.
3. Avoid clutter: Avoid overcrowding your plots with too much information. Use facets or multiple plots to show different aspects of the data.
4. Color wisely: Use colors purposefully to highlight different categories or groups in your data. Consider colorblind-friendly palettes for accessibility.
5. Interactivity: Use interactive plots when appropriate to allow viewers to explore the data themselves (e.g., using plotly in R).
6. Consistency: Maintain a consistent style and color scheme across multiple plots to make them easier to compare.

Data Manipulation:

1. Tidy data format: Organize your genomic data into a tidy format (as described earlier) to simplify data manipulation tasks using packages like dplyr and tidyr.
2. Filtering and subsetting: Use functions like filter() and slice() in dplyr to extract specific subsets of your data based on conditions.
3. Grouping and summarizing: Use group_by() and summarize() in dplyr to group your data by variables of interest and calculate summary statistics.
4. Joining datasets: Use left_join(), right_join(), inner_join(), or full_join() in dplyr to combine datasets based on common variables.
5. Avoiding unnecessary duplication: Remove duplicate rows or columns in your data using distinct() or select() in dplyr.
6. Handling missing data: Use na.omit() or complete.cases() to remove rows with missing data, or use mutate() to fill missing values with appropriate replacements.
7. Efficient coding: Use pipelines (%>%) in dplyr to chain together multiple data manipulation steps, making your code more readable and efficient.

By following these tips and best practices, you can effectively visualize and manipulate genomic data to uncover meaningful insights and drive biological discoveries.

Conclusion and Future Directions

Summary of Key Concepts and Skills Learned:

• ggplot2: Learned how to create a wide range of plots for visualizing genomic data, including scatter plots, bar plots, and heatmaps. Explored customization options for themes, scales, and color palettes to enhance plot aesthetics.
• dplyr: Learned how to manipulate genomic data, including filtering, sorting, grouping, and summarizing data. Used dplyr functions to perform these operations efficiently.
• Tidy Data Framework: Understood the principles of tidy data and applied them to organize genomic data into a tidy format. Learned how tidy data simplifies data manipulation and analysis.

Resources for Further Learning:

Books:

• • R for Data Science by Hadley Wickham and Garrett Grolemund: A comprehensive guide to data science using R, including chapters on ggplot2 and dplyr.
  • Bioinformatics Data Skills by Vince Buffalo: Covers a wide range of bioinformatics data analysis techniques, including using R for data manipulation and visualization.

In this project, we’ve learned how to effectively analyze and visualize genomic data using ggplot2, dplyr, and the tidy data framework in R. These skills are valuable for exploring complex biological datasets and generating insights that can contribute to scientific discoveries in genomics.

For future directions, consider expanding your skills by exploring advanced topics in data visualization (e.g., interactive plots, animation), mastering additional dplyr functions for complex data manipulations, and exploring other packages in the tidyverse and Bioconductor ecosystems for genomic data analysis. Continuing to stay updated with new developments in R and bioinformatics will also be beneficial for enhancing your skills and staying current in the field.