Data organisation with Spreadsheets

Overview

Teaching: XX min
Exercises: XX min

Questions

• How to organise tabular data?

Objectives

• Learn about spreadsheets, their strenghts and weaknesses

• How do we format data in spreadsheets for effective data use?

• Learn about common spreadsheet errors and how to correct them.

• Organize your data according to tidy data principles.

• Learn about text-based spreadsheet formats such as the comma-separated (CSV) or tab-separated formats.

Spreadsheet programs

Question

• What are basic principles for using spreadsheets for good data organization?

Objective

• Describe best practices for organizing data so computers can make the best use of data sets.

Keypoint

• Good data organization is the foundation of any research project.

Good data organization is the foundation of your research project. Most researchers have data or do data entry in spreadsheets. Spreadsheet programs are very useful graphical interfaces for designing data tables and handling very basic data quality control functions. See also @Broman:2018.

Spreadsheet outline

Spreadsheets are good for data entry. Therefore we have a lot of data in spreadsheets. Much of your time as a researcher will be spent in this ‘data wrangling’ stage. It’s not the most fun, but it’s necessary. We’ll teach you how to think about data organization and some practices for more effective data wrangling.

What this lesson will not teach you

• How to do statistics in a spreadsheet
• How to do plotting in a spreadsheet
• How to write code in spreadsheet programs

If you’re looking to do this, a good reference is Head First Excel, published by O’Reilly.

Why aren’t we teaching data analysis in spreadsheets

• Data analysis in spreadsheets usually requires a lot of manual work. If you want to change a parameter or run an analysis with a new dataset, you usually have to redo everything by hand. (We do know that you can create macros, but see the next point.)

• It is also difficult to track or reproduce statistical or plotting analyses done in spreadsheet programs when you want to go back to your work or someone asks for details of your analysis.

Many spreadsheet programs are available. Since most participants utilise Excel as their primary spreadsheet program, this lesson will make use of Excel examples. A free spreadsheet program that can also be used is LibreOffice. Commands may differ a bit between programs, but the general idea is the same.

Spreadsheet programs encompass a lot of the things we need to be able to do as researchers. We can use them for:

• Data entry
• Organizing data
• Subsetting and sorting data
• Statistics
• Plotting

Spreadsheet programs use tables to represent and display data. Data formatted as tables is also the main theme of this chapter, and we will see how to organise data into tables in a standardised way to ensure efficient downstream analysis.

Problems with Spreadsheets

Spreadsheets are good for data entry, but in reality we tend to use spreadsheet programs for much more than data entry. We use them to create data tables for publications, to generate summary statistics, and make figures.

Generating tables for publications in a spreadsheet is not optimal - often, when formatting a data table for publication, we’re reporting key summary statistics in a way that is not really meant to be read as data, and often involves special formatting (merging cells, creating borders, making it pretty). We advise you to do this sort of operation within your document editing software.

The latter two applications, generating statistics and figures, should be used with caution: because of the graphical, drag and drop nature of spreadsheet programs, it can be very difficult, if not impossible, to replicate your steps (much less retrace anyone else’s), particularly if your stats or figures require you to do more complex calculations. Furthermore, in doing calculations in a spreadsheet, it’s easy to accidentally apply a slightly different formula to multiple adjacent cells. When using a command-line based statistics program like R or SAS, it’s practically impossible to apply a calculation to one observation in your dataset but not another unless you’re doing it on purpose.

Using Spreadsheets for Data Entry and Cleaning

However, there are circumstances where you might want to use a spreadsheet program to produce “quick and dirty” calculations or figures, and data cleaning will help you use some of these features. Data cleaning also puts your data in a better format prior to importation into a statistical analysis program. We will show you how to use some features of spreadsheet programs to check your data quality along the way and produce preliminary summary statistics.

In this lesson, we will assume that you are most likely using Excel as your primary spreadsheet program - there are others (gnumeric, Calc from OpenOffice), and their functionality is similar, but Excel seems to be the program most used by biologists and biomedical researchers.

In this lesson we’re going to talk about:

1. Formatting data tables in spreadsheets
2. Formatting problems
3. Dates as data
4. Quality control
5. Exporting data

Formatting data tables in spreadsheets

Questions

• How do we format data in spreadsheets for effective data use?

Objectives

• Describe best practices for data entry and formatting in spreadsheets.

• Apply best practices to arrange variables and observations in a spreadsheet.

Keypoints

• Never modify your raw data. Always make a copy before making any changes.

• Keep track of all of the steps you take to clean your data in a plain text file.

• Organize your data according to tidy data principles.

The most common mistake made is treating spreadsheet programs like lab notebooks, that is, relying on context, notes in the margin, spatial layout of data and fields to convey information. As humans, we can (usually) interpret these things, but computers don’t view information the same way, and unless we explain to the computer what every single thing means (and that can be hard!), it will not be able to see how our data fits together.

Using the power of computers, we can manage and analyze data in much more effective and faster ways, but to use that power, we have to set up our data for the computer to be able to understand it (and computers are very literal).

This is why it’s extremely important to set up well-formatted tables from the outset - before you even start entering data from your very first preliminary experiment. Data organization is the foundation of your research project. It can make it easier or harder to work with your data throughout your analysis, so it’s worth thinking about when you’re doing your data entry or setting up your experiment. You can set things up in different ways in spreadsheets, but some of these choices can limit your ability to work with the data in other programs or have the you-of-6-months-from-now or your collaborator work with the data.

Note: the best layouts/formats (as well as software and interfaces) for data entry and data analysis might be different. It is important to take this into account, and ideally automate the conversion from one to another.

Keeping track of your analyses

When you’re working with spreadsheets, during data clean up or analyses, it’s very easy to end up with a spreadsheet that looks very different from the one you started with. In order to be able to reproduce your analyses or figure out what you did when a reviewer or instructor asks for a different analysis, you should

• create a new file with your cleaned or analyzed data. Don’t modify the original dataset, or you will never know where you started!

• keep track of the steps you took in your clean up or analysis. You should track these steps as you would any step in an experiment. We recommend that you do this in a plain text file stored in the same folder as the data file.

This might be an example of a spreadsheet setup:

Put these principles in to practice today during your exercises.

While versioning is out of scope for this course, you can look at the Carpentries lesson on ‘Git’ to learn how to maintain version control over your data. See also this blog post for a quick tutorial or @Perez-Riverol:2016 for a more research-oriented use-case.

Structuring data in spreadsheets

The cardinal rules of using spreadsheet programs for data:

1. Put all your variables in columns - the thing you’re measuring, like ‘weight’ or ‘temperature’.
2. Put each observation in its own row.
3. Don’t combine multiple pieces of information in one cell. Sometimes it just seems like one thing, but think if that’s the only way you’ll want to be able to use or sort that data.
4. Leave the raw data raw - don’t change it!
5. Export the cleaned data to a text-based format like CSV (comma-separated values) format. This ensures that anyone can use the data, and is required by most data repositories.

For instance, we have data from patients that visited several hospitals from Brussels, Belgium. They recorded the date of the visit, the hospital, the patients’ gender, weight and blood group.

If we were to keep track of the data like this:

the problem is that the ABO and Rhesus groups are in the same Blood type column. So, if they wanted to look at all observations of the A group or look at weight distributions by ABO group, it would be tricky to do this using this data setup. If instead we put the ABO and Rhesus groups in different columns, you can see that it would be much easier.

An important rule when setting up a datasheet, is that columns are used for variables and rows are used for observations:

• columns are variables
• rows are observations
• cells are individual values

Challenge: We’re going to take a messy data and describe how we would clean it up.

1. Download a messy data by clicking here.

2. Open up the data in a spreadsheet program.

3. You can see that there are two tabs. The data contains various clinical variables recorded in various hospitals in Brussels during the first and second COVID-19 waves in 2020. As you can see, the data have been recorded differently during the March and November waves. Now you’re the person in charge of this project and you want to be able to start analyzing the data.

4. With the person next to you, identify what is wrong with this spreadsheet. Also discuss the steps you would need to take to clean up first and second wave tabs, and to put them all together in one spreadsheet.

Important: Do not forget our first piece of advice: to create a new file (or tab) for the cleaned data, never modify your original (raw) data. However in this case we will have taken care of this for you!

After you go through this exercise, we’ll discuss as a group what was wrong with this data and how you would fix it.

• Take about 7 minutes to work on this exercise.
• All the mistakes in the common mistakes section below are present in the messy dataset. If the exercise is done during a workshop, ask people what they saw as wrong with the data. As they bring up different points, you can refer to the common mistakes or expand a bit on the point they brought up.
• If you get a response where they’ve fixed the date, you can pause and go to the dates lesson. Or you can say you’ll come back to dates at the end.

Challenge: Once you have tidied up the data, answer the following questions:

• How many men and women took part in the study?
• How many A, AB, and B types have been tested?
• As above, but disregarding the contaminated samples?
• How many Rhesus + and - have been tested?
• How many universal donors (0-) have been tested?
• What is the average weight of AB men?
• How many samples have been tested in the different hospitals?

An excellent reference, in particular with regard to R scripting is the Tidy Data paper @Wickham:2014.

Common Spreadsheet Errors

Questions

• What are some common challenges with formatting data in spreadsheets and how can we avoid them?

Objectives

• Recognize and resolve common spreadsheet formatting problems.

Keypoints

• Avoid using multiple tables within one spreadsheet.
• Avoid spreading data across multiple tabs.
• Record zeros as zeros.
• Use an appropriate null value to record missing data.
• Don’t use formatting to convey information or to make your spreadsheet look pretty.
• Place comments in a separate column.
• Record units in column headers.
• Include only one piece of information in a cell.
• Avoid spaces, numbers and special characters in column headers.
• Avoid special characters in your data.
• Record metadata in a separate plain text file.

There are a few potential errors to be on the lookout for in your own data as well as data from collaborators or the Internet. If you are aware of the errors and the possible negative effect on downstream data analysis and result interpretation, it might motivate yourself and your project members to try and avoid them. Making small changes to the way you format your data in spreadsheets, can have a great impact on efficiency and reliability when it comes to data cleaning and analysis.

• Using multiple tables
• Using multiple tabs
• Not filling in zeros
• Using problematic null values
• Using formatting to convey information
• Using formatting to make the data sheet look pretty
• Placing comments or units in cells
• Entering more than one piece of information in a cell
• Using problematic field names
• Using special characters in data
• Inclusion of metadata in data table

Using multiple tables

A common strategy is creating multiple data tables within one spreadsheet. This confuses the computer, so don’t do this! When you create multiple tables within one spreadsheet, you’re drawing false associations between things for the computer, which sees each row as an observation. You’re also potentially using the same field name in multiple places, which will make it harder to clean your data up into a usable form. The example below depicts the problem:

In the example above, the computer will see (for example) row 4 and assume that all columns A-AF refer to the same sample. This row actually represents four distinct samples (sample 1 for each of four different collection dates - May 29th, June 12th, June 19th, and June 26th), as well as some calculated summary statistics (an average (avr) and standard error of measurement (SEM)) for two of those samples. Other rows are similarly problematic.

Using multiple tabs

But what about workbook tabs? That seems like an easy way to organize data, right? Well, yes and no. When you create extra tabs, you fail to allow the computer to see connections in the data that are there (you have to introduce spreadsheet application-specific functions or scripting to ensure this connection). Say, for instance, you make a separate tab for each day you take a measurement.

This isn’t good practice for two reasons:

1. you are more likely to accidentally add inconsistencies to your data if each time you take a measurement, you start recording data in a new tab, and

2. even if you manage to prevent all inconsistencies from creeping in, you will add an extra step for yourself before you analyze the data because you will have to combine these data into a single datatable. You will have to explicitly tell the computer how to combine tabs - and if the tabs are inconsistently formatted, you might even have to do it manually.

The next time you’re entering data, and you go to create another tab or table, ask yourself if you could avoid adding this tab by adding another column to your original spreadsheet. We used multiple tabs in our example of a messy data file, but now you’ve seen how you can reorganize your data to consolidate across tabs.

Your data sheet might get very long over the course of the experiment. This makes it harder to enter data if you can’t see your headers at the top of the spreadsheet. But don’t repeat your header row. These can easily get mixed into the data, leading to problems down the road. Instead you can freeze the column headers so that they remain visible even when you have a spreadsheet with many rows.

Not filling in zeros

It might be that when you’re measuring something, it’s usually a zero, say the number of times a rabbit is observed in the survey. Why bother writing in the number zero in that column, when it’s mostly zeros?

However, there’s a difference between a zero and a blank cell in a spreadsheet. To the computer, a zero is actually data. You measured or counted it. A blank cell means that it wasn’t measured and the computer will interpret it as an unknown value (also known as a null or missing value).

The spreadsheets or statistical programs will likely mis-interpret blank cells that you intend to be zeros. By not entering the value of your observation, you are telling your computer to represent that data as unknown or missing (null). This can cause problems with subsequent calculations or analyses. For example, the average of a set of numbers which includes a single null value is always null (because the computer can’t guess the value of the missing observations). Because of this, it’s very important to record zeros as zeros and truly missing data as nulls.

Using problematic null values

Example: using -999 or other numerical values (or zero) to represent missing data.

Solutions:

There are a few reasons why null values get represented differently within a dataset. Sometimes confusing null values are automatically recorded from the measuring device. If that’s the case, there’s not much you can do, but it can be addressed in data cleaning with a tool like OpenRefine before analysis. Other times different null values are used to convey different reasons why the data isn’t there. This is important information to capture, but is in effect using one column to capture two pieces of information. Like for [using formatting to convey information]((#formatting) it would be good here to create a new column like ‘data_missing’ and use that column to capture the different reasons.

Whatever the reason, it’s a problem if unknown or missing data is recorded as -999, 999, or 0.

Many statistical programs will not recognize that these are intended to represent missing (null) values. How these values are interpreted will depend on the software you use to analyze your data. It is essential to use a clearly defined and consistent null indicator.

Blanks (most applications) and NA (for R) are good choices. @White:2013 explain good choices for indicating null values for different software applications in their article:

Using formatting to convey information

Example: highlighting cells, rows or columns that should be excluded from an analysis, leaving blank rows to indicate separations in data.

Solution: create a new field to encode which data should be excluded.

Using formatting to make the data sheet look pretty

Example: merging cells.

Solution: If you’re not careful, formatting a worksheet to be more aesthetically pleasing can compromise your computer’s ability to see associations in the data. Merged cells will make your data unreadable by statistics software. Consider restructuring your data in such a way that you will not need to merge cells to organize your data.

Placing comments or units in cells

Most analysis software can’t see Excel or LibreOffice comments, and would be confused by comments placed within your data cells. As described above for formatting, create another field if you need to add notes to cells. Similarly, don’t include units in cells: ideally, all the measurements you place in one column should be in the same unit, but if for some reason they aren’t, create another field and specify the units the cell is in.

Entering more than one piece of information in a cell

Example: Recording ABO and Rhesus groups in one cell, such as A+, B+, A-, …

Solution: Don’t include more than one piece of information in a cell. This will limit the ways in which you can analyze your data. If you need both these measurements, design your data sheet to include this information. For example, include one column the ABO group and one for the Rhesus group.

Using problematic field names

Choose descriptive field names, but be careful not to include spaces, numbers, or special characters of any kind. Spaces can be misinterpreted by parsers that use whitespace as delimiters and some programs don’t like field names that are text strings that start with numbers.

Underscores (_) are a good alternative to spaces. Consider writing names in camel case (like this: ExampleFileName) to improve readability. Remember that abbreviations that make sense at the moment may not be so obvious in 6 months, but don’t overdo it with names that are excessively long. Including the units in the field names avoids confusion and enables others to readily interpret your fields.

Examples

Using special characters in data

Example: You treat your spreadsheet program as a word processor when writing notes, for example copying data directly from Word or other applications.

Solution: This is a common strategy. For example, when writing longer text in a cell, people often include line breaks, em-dashes, etc in their spreadsheet. Also, when copying data in from applications such as Word, formatting and fancy non-standard characters (such as left- and right-aligned quotation marks) are included. When exporting this data into a coding/statistical environment or into a relational database, dangerous things may occur, such as lines being cut in half and encoding errors being thrown.

General best practice is to avoid adding characters such as newlines, tabs, and vertical tabs. In other words, treat a text cell as if it were a simple web form that can only contain text and spaces.

Inclusion of metadata in data table

Example: You add a legend at the top or bottom of your data table explaining column meaning, units, exceptions, etc.

Solution: Recording data about your data (“metadata”) is essential. You may be on intimate terms with your dataset while you are collecting and analysing it, but the chances that you will still remember that the variable “sglmemgp” means single member of group, for example, or the exact algorithm you used to transform a variable or create a derived one, after a few months, a year, or more are slim.

As well, there are many reasons other people may want to examine or use your data - to understand your findings, to verify your findings, to review your submitted publication, to replicate your results, to design a similar study, or even to archive your data for access and re-use by others. While digital data by definition are machine-readable, understanding their meaning is a job for human beings. The importance of documenting your data during the collection and analysis phase of your research cannot be overestimated, especially if your research is going to be part of the scholarly record.

However, metadata should not be contained in the data file itself. Unlike a table in a paper or a supplemental file, metadata (in the form of legends) should not be included in a data file since this information is not data, and including it can disrupt how computer programs interpret your data file. Rather, metadata should be stored as a separate file in the same directory as your data file, preferably in plain text format with a name that clearly associates it with your data file. Because metadata files are free text format, they also allow you to encode comments, units, information about how null values are encoded, etc. that are important to document but can disrupt the formatting of your data file.

Additionally, file or database level metadata describes how files that make up the dataset relate to each other; what format they are in; and whether they supercede or are superceded by previous files. A folder-level readme.txt file is the classic way of accounting for all the files and folders in a project.

(Text on metadata adapted from the online course Research Data MANTRA by EDINA and Data Library, University of Edinburgh. MANTRA is licensed under a Creative Commons Attribution 4.0 International License.)

Exporting data

Question

• How can we export data from spreadsheets in a way that is useful for downstream applications?

Objectives

• Store spreadsheet data in universal file formats.
• Export data from a spreadsheet to a CSV file.

Keypoints

• Data stored in common spreadsheet formats will often not be read correctly into data analysis software, introducing errors into your data.

• Exporting data from spreadsheets to formats like CSV or TSV puts it in a format that can be used consistently by most programs.

Storing the data you’re going to work with for your analyses in Excel default file format (*.xls or *.xlsx - depending on the Excel version) isn’t a good idea. Why?

• Because it is a proprietary format, and it is possible that in the future, technology won’t exist (or will become sufficiently rare) to make it inconvenient, if not impossible, to open the file.

• Other spreadsheet software may not be able to open files saved in a proprietary Excel format.

• Different versions of Excel may handle data differently, leading to inconsistencies. Dates is a well-documented example of inconsistencies in data storage.

• Finally, more journals and grant agencies are requiring you to deposit your data in a data repository, and most of them don’t accept Excel format. It needs to be in one of the formats discussed below.

• The above points also apply to other formats such as open data formats used by LibreOffice / Open Office. These formats are not static and do not get parsed the same way by different software packages.

Storing data in a universal, open, and static format will help deal with this problem. Try tab-delimited (tab separated values or TSV) or comma-delimited (comma separated values or CSV). CSV files are plain text files where the columns are separated by commas, hence ‘comma separated values’ or CSV. The advantage of a CSV file over an Excel/SPSS/etc. file is that we can open and read a CSV file using just about any software, including plain text editors like TextEdit or NotePad. Data in a CSV file can also be easily imported into other formats and environments, such as SQLite and R. We’re not tied to a certain version of a certain expensive program when we work with CSV files, so it’s a good format to work with for maximum portability and endurance. Most spreadsheet programs can save to delimited text formats like CSV easily, although they may give you a warning during the file export.

To save a file you have opened in Excel in CSV format:

1. From the top menu select ‘File’ and ‘Save as’.
2. In the ‘Format’ field, from the list, select ‘Comma Separated Values’ (*.csv).
3. Double check the file name and the location where you want to save it and hit ‘Save’.

An important note for backwards compatibility: you can open CSV files in Excel!

A note on R and xls: There are R packages that can read xls files (as well as Google spreadsheets). It is even possible to access different worksheets in the xls documents.

But

• some of these only work on Windows
• this equates to replacing a (simple but manual) export to csv with additional complexity/dependencies in the data analysis R code
• data formatting best practice still apply
• Is there really a good reason why csv (or similar) is not adequate?

Caveats on commas

In some datasets, the data values themselves may include commas (,). In that case, the software which you use (including Excel) will most likely incorrectly display the data in columns. This is because the commas which are a part of the data values will be interpreted as delimiters.

For example, our data might look like this:

species_id,genus,species,taxa
AB,Amphispiza,bilineata,Bird
AH,Ammospermophilus,harrisi,Rodent, not censused
AS,Ammodramus,savannarum,Bird
BA,Baiomys,taylori,Rodent

In the record AH,Ammospermophilus,harrisi,Rodent, not censused the value for taxa includes a comma (Rodent, not censused). If we try to read the above into Excel (or other spreadsheet program), we will get something like this:

The value for taxa was split into two columns (instead of being put in one column D). This can propagate to a number of further errors. For example, the extra column will be interpreted as a column with many missing values (and without a proper header). In addition to that, the value in column D for the record in row 3 (so the one where the value for ‘taxa’ contained the comma) is now incorrect.

If you want to store your data in csv format and expect that your data values may contain commas, you can avoid the problem discussed above by putting the values in quotes (“”). Applying this rule, our data might look like this:

species_id,genus,species,taxa
"AB","Amphispiza","bilineata","Bird"
"AH","Ammospermophilus","harrisi","Rodent, not censused"
"AS","Ammodramus","savannarum","Bird"
"BA","Baiomys","taylori","Rodent"

Now opening this file as a csv in Excel will not lead to an extra column, because Excel will only use commas that fall outside of quotation marks as delimiting characters.

Alternatively, if you are working with data that contains commas, you likely will need to use another delimiter when working in a spreadsheet¹. In this case, consider using tabs as your delimiter and working with TSV files. TSV files can be exported from spreadsheet programs in the same way as CSV files.

If you are working with an already existing dataset in which the data values are not included in “” but which have commas as both delimiters and parts of data values, you are potentially facing a major problem with data cleaning. If the dataset you’re dealing with contains hundreds or thousands of records, cleaning them up manually (by either removing commas from the data values or putting the values into quotes - “”) is not only going to take hours and hours but may potentially end up with you accidentally introducing many errors.

Cleaning up datasets is one of the major problems in many scientific disciplines. The approach almost always depends on the particular context. However, it is a good practice to clean the data in an automated fashion, for example by writing and running a script. The Python and R lessons will give you the basis for developing skills to build relevant scripts.

Summary

A typical data analysis worflow is illustrated in the figure above, where data is repeatedly tranformed, visualised, modelled. This iteration is repeated multiple times until the data is understood. In many real-life cases, however, most time is spent cleaning up and preparing the data, rather than actually analysing and understanding it.

An agile data analysis workflow, with several fast iterations of the transform/visualise/model cycle is only feasible is the data is formatted in a predictable way and one can reason about the data without having to look at it and/or fix it.

Key Points

• Good data organization is the foundation of any research project.

R and RStudio

Overview

Teaching: XX min
Exercises: XX min

Questions

• What are R and RStudio?

Objectives

• Describe the purpose of the RStudio Script, Console, Environment, and Plots panes.

• Organize files and directories for a set of analyses as an R project, and understand the purpose of the working directory.

• Use the built-in RStudio help interface to search for more information on R functions.

• Demonstrate how to provide sufficient information for troubleshooting with the R user community.

What is R? What is RStudio?

The term R is used to refer to both the programming language, the environment for statistical computing and the software that interprets the scripts written using it.

RStudio is currently a very popular way to not only write your R scripts but also to interact with the R software¹. To function correctly, RStudio needs R and therefore both need to be installed on your computer.

The RStudio IDE Cheat Sheet provides much more information that will be covered here, but can be useful to learn keyboard shortcuts and discover new features.

Why learn R?

R does not involve lots of pointing and clicking, and that’s a good thing

The learning curve might be steeper than with other software, but with R, the results of your analysis do not rely on remembering a succession of pointing and clicking, but instead on a series of written commands, and that’s a good thing! So, if you want to redo your analysis because you collected more data, you don’t have to remember which button you clicked in which order to obtain your results; you just have to run your script again.

Working with scripts makes the steps you used in your analysis clear, and the code you write can be inspected by someone else who can give you feedback and spot mistakes.

Working with scripts forces you to have a deeper understanding of what you are doing, and facilitates your learning and comprehension of the methods you use.

R code is great for reproducibility

Reproducibility means that someone else (including your future self) can obtain the same results from the same dataset when using the same analysis code.

R integrates with other tools to generate manuscripts or reports from your code. If you collect more data, or fix a mistake in your dataset, the figures and the statistical tests in your manuscript or report are updated automatically.

An increasing number of journals and funding agencies expect analyses to be reproducible, so knowing R will give you an edge with these requirements.

R is interdisciplinary and extensible

With 10000+ packages² that can be installed to extend its capabilities, R provides a framework that allows you to combine statistical approaches from many scientific disciplines to best suit the analytical framework you need to analyze your data. For instance, R has packages for image analysis, GIS, time series, population genetics, and a lot more.

R works on data of all shapes and sizes

The skills you learn with R scale easily with the size of your dataset. Whether your dataset has hundreds or millions of lines, it won’t make much difference to you.

R is designed for data analysis. It comes with special data structures and data types that make handling of missing data and statistical factors convenient.

R can connect to spreadsheets, databases, and many other data formats, on your computer or on the web.

R produces high-quality graphics

The plotting functionalities in R are extensive, and allow you to adjust any aspect of your graph to convey most effectively the message from your data.

R has a large and welcoming community

Thousands of people use R daily. Many of them are willing to help you through mailing lists and websites such as Stack Overflow, or on the RStudio community. These broad user community extends to specialised areas such as bioinformatics.

Not only is R free, but it is also open-source and cross-platform

Anyone can inspect the source code to see how R works. Because of this transparency, there is less chance for mistakes, and if you (or someone else) find some, you can report and fix bugs.

Knowing your way around RStudio

Let’s start by learning about RStudio, which is an Integrated Development Environment (IDE) for working with R.

The RStudio IDE open-source product is free under the Affero General Public License (AGPL) v3. The RStudio IDE is also available with a commercial license and priority email support from RStudio, Inc.

We will use the RStudio IDE to write code, navigate the files on our computer, inspect the variables we are going to create, and visualize the plots we will generate. RStudio can also be used for other things (e.g., version control, developing packages, writing Shiny apps) that we will not cover during the workshop.

The RStudio window is divided into 4 “Panes”:

• the Source for your scripts and documents (top-left, in the default layout)
• your Environment/History (top-right),
• your Files/Plots/Packages/Help/Viewer (bottom-right), and
• the R Console (bottom-left).

The placement of these panes and their content can be customized (see menu, Tools -> Global Options -> Pane Layout).

One of the advantages of using RStudio is that all the information you need to write code is available in a single window. Additionally, with many shortcuts, autocompletion, and highlighting for the major file types you use while developing in R, RStudio will make typing easier and less error-prone.

Getting set up

It is good practice to keep a set of related data, analyses, and text self-contained in a single folder, called the working directory. All of the scripts within this folder can then use relative paths to files that indicate where inside the project a file is located (as opposed to absolute paths, which point to where a file is on a specific computer). Working this way makes it a lot easier to move your project around on your computer and share it with others without worrying about whether or not the underlying scripts will still work.

RStudio provides a helpful set of tools to do this through its “Projects” interface, which not only creates a working directory for you, but also remembers its location (allowing you to quickly navigate to it) and optionally preserves custom settings and open files to make it easier to resume work after a break. Go through the steps for creating an “R Project” for this tutorial below.

1. Start RStudio.
2. Under the File menu, click on New project. Choose New directory, then New project.
3. Enter a name for this new folder (or “directory”), and choose a convenient location for it. This will be your working directory for this session (or whole course) (e.g., bioc-intro).
4. Click on Create project.
5. (Optional) Set Preferences to ‘Never’ save workspace in RStudio.

RStudio’s default preferences generally work well, but saving a workspace to .RData can be cumbersome, especially if you are working with larger datasets. To turn that off, go to Tools –> ‘Global Options’ and select the ‘Never’ option for ‘Save workspace to .RData’ on exit.’

To avoid character encoding issue between Windows and other operating systems, we are going to set UTF-8 by default:

Organizing your working directory

Using a consistent folder structure across your projects will help keep things organized, and will also make it easy to find/file things in the future. This can be especially helpful when you have multiple projects. In general, you may create directories (folders) for scripts, data, and documents.

• data/ Use this folder to store your raw data and intermediate datasets you may create for the need of a particular analysis. For the sake of transparency and provenance, you should always keep a copy of your raw data accessible and do as much of your data cleanup and preprocessing programmatically (i.e., with scripts, rather than manually) as possible. Separating raw data from processed data is also a good idea. For example, you could have files data/raw/tree_survey.plot1.txt and ...plot2.txt kept separate from a data/processed/tree.survey.csv file generated by the scripts/01.preprocess.tree_survey.R script.
• documents/ This would be a place to keep outlines, drafts, and other text.
• scripts/ (or src) This would be the location to keep your R scripts for different analyses or plotting, and potentially a separate folder for your functions (more on that later).

You may want additional directories or subdirectories depending on your project needs, but these should form the backbone of your working directory.

For this course, we will need a data/ folder to store our raw data, and we will use data_output/ for when we learn how to export data as CSV files, and fig_output/ folder for the figures that we will save.

Challenge: create your project directory structure

Create an RStudio project for this workshop called e.g. bioc-intro, as described above under Getting set up, if you haven’t already done that. Under the Files tab on the right of the screen, click on New Folder and create a folder named data within your newly created working directory (e.g., ~/bioc-intro/data). (Alternatively, type dir.create("data") at your R console.) Repeat these operations to create a data_output/ and a fig_output folders.

Your working directory should now look like this:

The working directory

The working directory is an important concept to understand. It is the place from where R will be looking for and saving the files. When you write code for your project, it should refer to files in relation to the root of your working directory and only need files within this structure.

Using RStudio projects makes this easy and ensures that your working directory is set properly. If you need to check it, you can use getwd(). If for some reason your working directory is not what it should be, you can change it in the RStudio interface by navigating in the file browser where your working directory should be, and clicking on the blue gear icon More, and select Set As Working Directory. Alternatively you can use setwd("/path/to/working/directory") to reset your working directory. However, your scripts should not include this line because it will fail on someone else’s computer.

Example

The schema below represents the working directory bioc-intro with the data and fig_output sub-directories, and 2 files in the latter:

bioc-intro/data/
          /fig_output/fig1.pdf
          /fig_output/fig2.png

If we were in the working directory, we could refer to the fig1.pdf file using the relative path fig_output/fig1.pdf or the absolute path /home/user/bioc-intro/fig_output/fig1.pdf.

If we were in the data directory, we would use the relative path ../fig_output/fig1.pdf or the same absolute path /home/user/bioc-intro/fig_output/fig1.pdf.

Interacting with R

The basis of programming is that we write down instructions for the computer to follow, and then we tell the computer to follow those instructions. We write, or code, instructions in R because it is a common language that both the computer and we can understand. We call the instructions commands and we tell the computer to follow the instructions by executing (also called running) those commands.

There are two main ways of interacting with R: by using the console or by using scripts (plain text files that contain your code). The console pane (in RStudio, the bottom left panel) is the place where commands written in the R language can be typed and executed immediately by the computer. It is also where the results will be shown for commands that have been executed. You can type commands directly into the console and press Enter to execute those commands, but they will be forgotten when you close the session.

Because we want our code and workflow to be reproducible, it is better to type the commands we want in the script editor, and save the script. This way, there is a complete record of what we did, and anyone (including our future selves!) can easily replicate the results on their computer. Note, however, that merely typing the commands in the script does not automatically run them - they still need to be sent to the console for execution.

RStudio allows you to execute commands directly from the script editor by using the Ctrl + Enter shortcut (on Macs, Cmd + Return will work, too). The command on the current line in the script (indicated by the cursor) or all of the commands in the currently selected text will be sent to the console and executed when you press Ctrl + Enter. You can find other keyboard shortcuts in this RStudio cheatsheet about the RStudio IDE.

At some point in your analysis you may want to check the content of a variable or the structure of an object, without necessarily keeping a record of it in your script. You can type these commands and execute them directly in the console. RStudio provides the Ctrl + 1 and Ctrl + 2 shortcuts allow you to jump between the script and the console panes.

If R is ready to accept commands, the R console shows a > prompt. If it receives a command (by typing, copy-pasting or sent from the script editor using Ctrl + Enter), R will try to execute it, and when ready, will show the results and come back with a new > prompt to wait for new commands.

If R is still waiting for you to enter more data because it isn’t complete yet, the console will show a + prompt. It means that you haven’t finished entering a complete command. This is because you have not ‘closed’ a parenthesis or quotation, i.e. you don’t have the same number of left-parentheses as right-parentheses, or the same number of opening and closing quotation marks. When this happens, and you thought you finished typing your command, click inside the console window and press Esc; this will cancel the incomplete command and return you to the > prompt.

How to learn more during and after the course?

The material we cover during this course will give you an initial taste of how you can use R to analyse data for your own research. However, you will need to learn more to do advanced operations such as cleaning your dataset, using statistical methods, or creating beautiful graphics³. The best way to become proficient and efficient at R, as with any other tool, is to use it to address your actual research questions. As a beginner, it can feel daunting to have to write a script from scratch, and given that many people make their code available online, modifying existing code to suit your purpose might make it easier for you to get started.

Seeking help

Use the built-in RStudio help interface to search for more information on R functions

One of the fastest ways to get help, is to use the RStudio help interface. This panel by default can be found at the lower right hand panel of RStudio. As seen in the screenshot, by typing the word “Mean”, RStudio tries to also give a number of suggestions that you might be interested in. The description is then shown in the display window.

I know the name of the function I want to use, but I’m not sure how to use it

If you need help with a specific function, let’s say barplot(), you can type:

?barplot

If you just need to remind yourself of the names of the arguments, you can use:

args(lm)

I want to use a function that does X, there must be a function for it but I don’t know which one…

If you are looking for a function to do a particular task, you can use the help.search() function, which is called by the double question mark ??. However, this only looks through the installed packages for help pages with a match to your search request

??kruskal

If you can’t find what you are looking for, you can use the rdocumentation.org website that searches through the help files across all packages available.

Finally, a generic Google or internet search “R <task>” will often either send you to the appropriate package documentation or a helpful forum where someone else has already asked your question.

I am stuck… I get an error message that I don’t understand

Start by googling the error message. However, this doesn’t always work very well because often, package developers rely on the error catching provided by R. You end up with general error messages that might not be very helpful to diagnose a problem (e.g. “subscript out of bounds”). If the message is very generic, you might also include the name of the function or package you’re using in your query.

However, you should check Stack Overflow. Search using the [r] tag. Most questions have already been answered, but the challenge is to use the right words in the search to find the answers:

http://stackoverflow.com/questions/tagged/r

The Introduction to R can also be dense for people with little programming experience but it is a good place to understand the underpinnings of the R language.

The R FAQ is dense and technical but it is full of useful information.

Asking for help

The key to receiving help from someone is for them to rapidly grasp your problem. You should make it as easy as possible to pinpoint where the issue might be.

Try to use the correct words to describe your problem. For instance, a package is not the same thing as a library. Most people will understand what you meant, but others have really strong feelings about the difference in meaning. The key point is that it can make things confusing for people trying to help you. Be as precise as possible when describing your problem.

If possible, try to reduce what doesn’t work to a simple reproducible example. If you can reproduce the problem using a very small data frame instead of your 50000 rows and 10000 columns one, provide the small one with the description of your problem. When appropriate, try to generalize what you are doing so even people who are not in your field can understand the question. For instance instead of using a subset of your real dataset, create a small (3 columns, 5 rows) generic one. For more information on how to write a reproducible example see this article by Hadley Wickham.

To share an object with someone else, if it’s relatively small, you can use the function dput(). It will output R code that can be used to recreate the exact same object as the one in memory:

## iris is an example data frame that comes with R and head() is a
## function that returns the first part of the data frame
dput(head(iris))

structure(list(Sepal.Length = c(5.1, 4.9, 4.7, 4.6, 5, 5.4), 
    Sepal.Width = c(3.5, 3, 3.2, 3.1, 3.6, 3.9), Petal.Length = c(1.4, 
    1.4, 1.3, 1.5, 1.4, 1.7), Petal.Width = c(0.2, 0.2, 0.2, 
    0.2, 0.2, 0.4), Species = structure(c(1L, 1L, 1L, 1L, 1L, 
    1L), .Label = c("setosa", "versicolor", "virginica"), class = "factor")), row.names = c(NA, 
6L), class = "data.frame")

If the object is larger, provide either the raw file (i.e., your CSV file) with your script up to the point of the error (and after removing everything that is not relevant to your issue). Alternatively, in particular if your question is not related to a data frame, you can save any R object to a file[^export]:

saveRDS(iris, file="/tmp/iris.rds")

The content of this file is however not human readable and cannot be posted directly on Stack Overflow. Instead, it can be sent to someone by email who can read it with the readRDS() command (here it is assumed that the downloaded file is in a Downloads folder in the user’s home directory):

some_data <- readRDS(file="~/Downloads/iris.rds")

Last, but certainly not least, always include the output of sessionInfo() as it provides critical information about your platform, the versions of R and the packages that you are using, and other information that can be very helpful to understand your problem.

sessionInfo()

R version 4.1.3 (2022-03-10)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.2.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] knitr_1.39

loaded via a namespace (and not attached):
[1] compiler_4.1.3 magrittr_2.0.3 tools_4.1.3    stringi_1.7.6  highr_0.9     
[6] stringr_1.4.0  xfun_0.31      evaluate_0.15 

Where to ask for help?

• The person sitting next to you during the course. Don’t hesitate to talk to your neighbour during the workshop, compare your answers, and ask for help.
• Your friendly colleagues: if you know someone with more experience than you, they might be able and willing to help you.
• Stack Overflow: if your question hasn’t been answered before and is well crafted, chances are you will get an answer in less than 5 min. Remember to follow their guidelines on how to ask a good question.
• The R-help mailing list: it is read by a lot of people (including most of the R core team), a lot of people post to it, but the tone can be pretty dry, and it is not always very welcoming to new users. If your question is valid, you are likely to get an answer very fast but don’t expect that it will come with smiley faces. Also, here more than anywhere else, be sure to use correct vocabulary (otherwise you might get an answer pointing to the misuse of your words rather than answering your question). You will also have more success if your question is about a base function rather than a specific package.
• If your question is about a specific package, see if there is a mailing list for it. Usually it’s included in the DESCRIPTION file of the package that can be accessed using packageDescription("name-of-package"). You may also want to try to email the author of the package directly, or open an issue on the code repository (e.g., GitHub).
• There are also some topic-specific mailing lists (GIS, phylogenetics, etc…), the complete list is here.

More resources

• The Posting Guide for the R mailing lists.

• How to ask for R help useful guidelines

• This blog post by Jon Skeet has quite comprehensive advice on how to ask programming questions.

• The reprex package is very helpful to create reproducible examples when asking for help. The rOpenSci community call “How to ask questions so they get answered” (Github link and video recording) includes a presentation of the reprex package and of its philosophy.

R packages

Loading packages

As we have seen above, R packages play a fundamental role in R. The make use of a package’s functionality, assuming it is installed, we first need to load it to be able to use it. This is done with the library() function. Below, we load ggplot2.

library("ggplot2")

Installing packages

The default package repository is The Comprehensive R Archive Network (CRAN), and any package that is available on CRAN can be installed with the install.packages() function. Below, for example, we install the dplyr package that we will learn about late.

install.packages("dplyr")

This command will install the dplyr package as well as all its dependencies, i.e. all the packages that it relies on to function.

Bioconductor packages are managed installed using a dedicated package, namely BiocManager, that can be installed from CRAN with

install.packages("BiocManager")

Individual packages such as SummarizedExperiment (we will use it later), DESeq2 (for RNA-Seq analysis), and many more can then be installed with BiocManager::install.

BiocManager::install("SummarizedExperiment")
BiocManager::install("DESeq2")

Key Points

• Start using R and RStudio

Introduction to R

Overview

Teaching: XX min
Exercises: XX min

Questions

• First commands in R

Objectives

• Define the following terms as they relate to R: object, assign, call, function, arguments, options.

• Assign values to objects in R.

• Learn how to name objects

• Use comments to inform script.

• Solve simple arithmetic operations in R.

• Call functions and use arguments to change their default options.

• Inspect the content of vectors and manipulate their content.

• Subset and extract values from vectors.

• Analyze vectors with missing data.

Creating objects in R

You can get output from R simply by typing math in the console:

3 + 5

[1] 8

12 / 7

[1] 1.714286

However, to do useful and interesting things, we need to assign values to objects. To create an object, we need to give it a name followed by the assignment operator <-, and the value we want to give it:

weight_kg <- 55

<- is the assignment operator. It assigns values on the right to objects on the left. So, after executing x <- 3, the value of x is 3. The arrow can be read as 3 goes into x. For historical reasons, you can also use = for assignments, but not in every context. Because of the slight differences in syntax, it is good practice to always use <- for assignments.

In RStudio, typing Alt + - (push Alt at the same time as the - key) will write <- in a single keystroke in a PC, while typing Option + - (push Option at the same time as the - key) does the same in a Mac.

Naming variables {-}

Objects can be given any name such as x, current_temperature, or subject_id. You want your object names to be explicit and not too long. They cannot start with a number (2x is not valid, but x2 is). R is case sensitive (e.g., weight_kg is different from Weight_kg). There are some names that cannot be used because they are the names of fundamental functions in R (e.g., if, else, for, see here for a complete list). In general, even if it’s allowed, it’s best to not use other function names (e.g., c, T, mean, data, df, weights). If in doubt, check the help to see if the name is already in use. It’s also best to avoid dots (.) within an object name as in my.dataset. There are many functions in R with dots in their names for historical reasons, but because dots have a special meaning in R (for methods) and other programming languages, it’s best to avoid them. It is also recommended to use nouns for object names, and verbs for function names. It’s important to be consistent in the styling of your code (where you put spaces, how you name objects, etc.). Using a consistent coding style makes your code clearer to read for your future self and your collaborators. In R, some popular style guides are Google’s, the tidyverse’s style and the Bioconductor style guide. The tidyverse’s is very comprehensive and may seem overwhelming at first. You can install the lintr package to automatically check for issues in the styling of your code.

Objects vs. variables What are known as objects in R are known as variables in many other programming languages. Depending on the context, object and variable can have drastically different meanings. However, in this lesson, the two words are used synonymously. For more information see: https://cran.r-project.org/doc/manuals/r-release/R-lang.html#Objects

When assigning a value to an object, R does not print anything. You can force R to print the value by using parentheses or by typing the object name:

weight_kg <- 55    # doesn't print anything
(weight_kg <- 55)  # but putting parenthesis around the call prints the value of `weight_kg`

[1] 55

weight_kg          # and so does typing the name of the object

[1] 55

Now that R has weight_kg in memory, we can do arithmetic with it. For instance, we may want to convert this weight into pounds (weight in pounds is 2.2 times the weight in kg):

2.2 * weight_kg

[1] 121

We can also change an object’s value by assigning it a new one:

weight_kg <- 57.5
2.2 * weight_kg

[1] 126.5

This means that assigning a value to one object does not change the values of other objects For example, let’s store the animal’s weight in pounds in a new object, weight_lb:

weight_lb <- 2.2 * weight_kg

and then change weight_kg to 100.

weight_kg <- 100

Challenge:

What do you think is the current content of the object weight_lb? 126.5 or 220?

Comments

The comment character in R is #, anything to the right of a # in a script will be ignored by R. It is useful to leave notes, and explanations in your scripts.

RStudio makes it easy to comment or uncomment a paragraph: after selecting the lines you want to comment, press at the same time on your keyboard Ctrl + Shift + C. If you only want to comment out one line, you can put the cursor at any location of that line (i.e. no need to select the whole line), then press Ctrl + Shift + C.

Challenge

What are the values after each statement in the following?

mass <- 47.5            # mass?
age  <- 122             # age?
mass <- mass * 2.0      # mass?
age  <- age - 20        # age?
mass_index <- mass/age  # mass_index?

Functions and their arguments

Functions are “canned scripts” that automate more complicated sets of commands including operations assignments, etc. Many functions are predefined, or can be made available by importing R packages (more on that later). A function usually gets one or more inputs called arguments. Functions often (but not always) return a value. A typical example would be the function sqrt(). The input (the argument) must be a number, and the return value (in fact, the output) is the square root of that number. Executing a function (‘running it’) is called calling the function. An example of a function call is:

b <- sqrt(a)

Here, the value of a is given to the sqrt() function, the sqrt() function calculates the square root, and returns the value which is then assigned to the object b. This function is very simple, because it takes just one argument.

The return ‘value’ of a function need not be numerical (like that of sqrt()), and it also does not need to be a single item: it can be a set of things, or even a dataset. We’ll see that when we read data files into R.

Arguments can be anything, not only numbers or filenames, but also other objects. Exactly what each argument means differs per function, and must be looked up in the documentation (see below). Some functions take arguments which may either be specified by the user, or, if left out, take on a default value: these are called options. Options are typically used to alter the way the function operates, such as whether it ignores ‘bad values’, or what symbol to use in a plot. However, if you want something specific, you can specify a value of your choice which will be used instead of the default.

Let’s try a function that can take multiple arguments: round().

round(3.14159)

[1] 3

Here, we’ve called round() with just one argument, 3.14159, and it has returned the value 3. That’s because the default is to round to the nearest whole number. If we want more digits we can see how to do that by getting information about the round function. We can use args(round) or look at the help for this function using ?round.

args(round)

function (x, digits = 0) 
NULL

?round

We see that if we want a different number of digits, we can type digits=2 or however many we want.

round(3.14159, digits = 2)

[1] 3.14

If you provide the arguments in the exact same order as they are defined you don’t have to name them:

round(3.14159, 2)

[1] 3.14

And if you do name the arguments, you can switch their order:

round(digits = 2, x = 3.14159)

[1] 3.14

It’s good practice to put the non-optional arguments (like the number you’re rounding) first in your function call, and to specify the names of all optional arguments. If you don’t, someone reading your code might have to look up the definition of a function with unfamiliar arguments to understand what you’re doing. By specifying the name of the arguments you are also safeguarding against possible future changes in the function interface, which may potentially add new arguments in between the existing ones.

Vectors and data types

A vector is the most common and basic data type in R, and is pretty much the workhorse of R. A vector is composed by a series of values, such as numbers or characters. We can assign a series of values to a vector using the c() function. For example we can create a vector of animal weights and assign it to a new object weight_g:

weight_g <- c(50, 60, 65, 82)
weight_g

[1] 50 60 65 82

A vector can also contain characters:

molecules <- c("dna", "rna", "protein")
molecules

[1] "dna"     "rna"     "protein"

The quotes around “dna”, “rna”, etc. are essential here. Without the quotes R will assume there are objects called dna, rna and protein. As these objects don’t exist in R’s memory, there will be an error message.

There are many functions that allow you to inspect the content of a vector. length() tells you how many elements are in a particular vector:

length(weight_g)

[1] 4

length(molecules)

[1] 3

An important feature of a vector, is that all of the elements are the same type of data. The function class() indicates the class (the type of element) of an object:

class(weight_g)

[1] "numeric"

class(molecules)

[1] "character"

The function str() provides an overview of the structure of an object and its elements. It is a useful function when working with large and complex objects:

str(weight_g)

 num [1:4] 50 60 65 82

str(molecules)

 chr [1:3] "dna" "rna" "protein"

You can use the c() function to add other elements to your vector:

weight_g <- c(weight_g, 90) # add to the end of the vector
weight_g <- c(30, weight_g) # add to the beginning of the vector
weight_g

[1] 30 50 60 65 82 90

In the first line, we take the original vector weight_g, add the value 90 to the end of it, and save the result back into weight_g. Then we add the value 30 to the beginning, again saving the result back into weight_g.

We can do this over and over again to grow a vector, or assemble a dataset. As we program, this may be useful to add results that we are collecting or calculating.

An atomic vector is the simplest R data type and is a linear vector of a single type. Above, we saw 2 of the 6 main atomic vector types that R uses: "character" and "numeric" (or "double"). These are the basic building blocks that all R objects are built from. The other 4 atomic vector types are:

• "logical" for TRUE and FALSE (the boolean data type)
• "integer" for integer numbers (e.g., 2L, the L indicates to R that it’s an integer)
• "complex" to represent complex numbers with real and imaginary parts (e.g., 1 + 4i) and that’s all we’re going to say about them
• "raw" for bitstreams that we won’t discuss further

You can check the type of your vector using the typeof() function and inputting your vector as the argument.

Vectors are one of the many data structures that R uses. Other important ones are lists (list), matrices (matrix), data frames (data.frame), factors (factor) and arrays (array).

Challenge:

We’ve seen that atomic vectors can be of type character, numeric (or double), integer, and logical. But what happens if we try to mix these types in a single vector?

Solution

R implicitly converts them to all be the same type

Challenge:

What will happen in each of these examples? (hint: use class() to check the data type of your objects):

num_char <- c(1, 2, 3, "a")
num_logical <- c(1, 2, 3, TRUE)
char_logical <- c("a", "b", "c", TRUE)
tricky <- c(1, 2, 3, "4")

Solution

class(num_char)

[1] "character"

class(num_logical)

[1] "numeric"

class(char_logical)

[1] "character"

class(tricky)

[1] "character"

Challenge:

Why do you think it happens?

Solution

Vectors can be of only one data type. R tries to convert (coerce) the content of this vector to find a common denominator that doesn’t lose any information.

Homework Challenge:

How many values in combined_logical are "TRUE" (as a character) in the following example:

num_logical <- c(1, 2, 3, TRUE)
char_logical <- c("a", "b", "c", TRUE)
combined_logical <- c(num_logical, char_logical)

Homework Challenge:

In R, we call converting objects from one class into another class coercion. These conversions happen according to a hierarchy, whereby some types get preferentially coerced into other types. Can you draw a diagram that represents the hierarchy of how these data types are coerced?

Solution

logical → numeric → character ← logical

Subsetting vectors

If we want to extract one or several values from a vector, we must provide one or several indices in square brackets. For instance:

molecules <- c("dna", "rna", "peptide", "protein")
molecules[2]

[1] "rna"

molecules[c(3, 2)]

[1] "peptide" "rna"    

We can also repeat the indices to create an object with more elements than the original one:

more_molecules <- molecules[c(1, 2, 3, 2, 1, 4)]
more_molecules

[1] "dna"     "rna"     "peptide" "rna"     "dna"     "protein"

R indices start at 1. Programming languages like Fortran, MATLAB, Julia, and R start counting at 1, because that’s what human beings typically do. Languages in the C family (including C++, Java, Perl, and Python) count from 0 because that’s simpler for computers to do.

Finally, it is also possible to get all the elements of a vector except some specified elements using negative indices:

molecules ## all molecules

[1] "dna"     "rna"     "peptide" "protein"

molecules[-1] ## all but the first one

[1] "rna"     "peptide" "protein"

molecules[-c(1, 3)] ## all but 1st/3rd ones

[1] "rna"     "protein"

molecules[c(-1, -3)] ## all but 1st/3rd ones

[1] "rna"     "protein"

Conditional subsetting

Another common way of subsetting is by using a logical vector. TRUE will select the element with the same index, while FALSE will not:

weight_g <- c(21, 34, 39, 54, 55)
weight_g[c(TRUE, FALSE, TRUE, TRUE, FALSE)]

[1] 21 39 54

Typically, these logical vectors are not typed by hand, but are the output of other functions or logical tests. For instance, if you wanted to select only the values above 50:

## will return logicals with TRUE for the indices that meet
## the condition
weight_g > 50

[1] FALSE FALSE FALSE  TRUE  TRUE

## so we can use this to select only the values above 50
weight_g[weight_g > 50]

[1] 54 55

You can combine multiple tests using & (both conditions are true, AND) or | (at least one of the conditions is true, OR):

weight_g[weight_g < 30 | weight_g > 50]

[1] 21 54 55

weight_g[weight_g >= 30 & weight_g == 21]

numeric(0)

Here, < stands for “less than”, > for “greater than”, >= for “greater than or equal to”, and == for “equal to”. The double equal sign == is a test for numerical equality between the left and right hand sides, and should not be confused with the single = sign, which performs variable assignment (similar to <-).

A common task is to search for certain strings in a vector. One could use the “or” operator | to test for equality to multiple values, but this can quickly become tedious. The function %in% allows you to test if any of the elements of a search vector are found:

molecules <- c("dna", "rna", "protein", "peptide")
molecules[molecules == "rna" | molecules == "dna"] # returns both rna and dna

[1] "dna" "rna"

molecules %in% c("rna", "dna", "metabolite", "peptide", "glycerol")

[1]  TRUE  TRUE FALSE  TRUE

molecules[molecules %in% c("rna", "dna", "metabolite", "peptide", "glycerol")]

[1] "dna"     "rna"     "peptide"

Homework Challenge:

Can you figure out why "four" > "five" returns TRUE?

Solution

"four" > "five"

[1] TRUE

When using > or < on strings, R compares their alphabetical order. Here "four" comes after "five", and therefore is greater than it.

Names

It is possible to name each element of a vector. The code chunk below show a initial vector without any names, how names are set, and retrieved.

x <- c(1, 5, 3, 5, 10)
names(x) ## no names

NULL

names(x) <- c("A", "B", "C", "D", "E")
names(x) ## now we have names

[1] "A" "B" "C" "D" "E"

When a vector has names, it is possible to access elements by their name, in addition to their index.

x[c(1, 3)]

A C 
1 3 

x[c("A", "C")]

A C 
1 3 

Missing data

As R was designed to analyze datasets, it includes the concept of missing data (which is uncommon in other programming languages). Missing data are represented in vectors as NA.

When doing operations on numbers, most functions will return NA if the data you are working with include missing values. This feature makes it harder to overlook the cases where you are dealing with missing data. You can add the argument na.rm = TRUE to calculate the result while ignoring the missing values.

heights <- c(2, 4, 4, NA, 6)
mean(heights)

[1] NA

max(heights)

[1] NA

mean(heights, na.rm = TRUE)

[1] 4

max(heights, na.rm = TRUE)

[1] 6

If your data include missing values, you may want to become familiar with the functions is.na(), na.omit(), and complete.cases(). See below for examples.

## Extract those elements which are not missing values.
heights[!is.na(heights)]

[1] 2 4 4 6

## Returns the object with incomplete cases removed.
## The returned object is an atomic vector of type `"numeric"`
## (or `"double"`).
na.omit(heights)

[1] 2 4 4 6
attr(,"na.action")
[1] 4
attr(,"class")
[1] "omit"

## Extract those elements which are complete cases.
## The returned object is an atomic vector of type `"numeric"`
## (or `"double"`).
heights[complete.cases(heights)]

[1] 2 4 4 6

Homework Challenge:

1. Using this vector of heights in inches, create a new vector with the NAs removed.

heights <- c(63, 69, 60, 65, NA, 68, 61, 70, 61, 59, 64, 69, 63, 63, NA, 72, 65, 64, 70, 63, 65)

1. Use the function median() to calculate the median of the heights vector.
2. Use R to figure out how many people in the set are taller than 67 inches.

Solution

heights_no_na <- heights[!is.na(heights)]
## or
heights_no_na <- na.omit(heights)

median(heights, na.rm = TRUE)

[1] 64

heights_above_67 <- heights_no_na[heights_no_na > 67]
length(heights_above_67)

[1] 6

Generating vectors

Constructors {-}

There exists some functions to generate vectors of different type. To generate a vector of numerics, one can use the numeric() constructor, providing the length of the output vector as parameter. The values will be initialised with 0.

numeric(3)

[1] 0 0 0

numeric(10)

 [1] 0 0 0 0 0 0 0 0 0 0

Note that if we ask for a vector of numerics of length 0, we obtain exactly that:

numeric(0)

numeric(0)

There are similar constructors for characters and logicals, named character() and logical() respectively.

Homework Challenge:

What are the defaults for character and logical vectors?

Solution

character(2) ## the empty charater

[1] "" ""

logical(2)   ## FALSE

[1] FALSE FALSE

Replicate elements

The rep function allow to repeat a value a certain number of times. If we want to initiate a vector of numerics of length 5 with the value -1, for example, we could do the following:

rep(-1, 5)

[1] -1 -1 -1 -1 -1

Similarly, to generate a vector populated with missing values, which is often a good way to start, without setting assumptions on the data to be collected:

rep(NA, 5)

[1] NA NA NA NA NA

rep can take vectors of any length as input (above, we used vectors of length 1) and any type. For example, if we want to repeat the values 1, 2 and 3 five times, we would do the following:

rep(c(1, 2, 3), 5)

 [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

Homework Challenge:

What if we wanted to repeat the values 1, 2 and 3 five times, but obtain five 1s, five 2s and five 3s in that order? There are two possibilities - see ?rep or ?sort for help.

Solution

rep(c(1, 2, 3), each = 5)

 [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

sort(rep(c(1, 2, 3), 5))

 [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

Sequence generation

Another very useful function is seq, to generate a sequence of numbers. For example, to generate a sequence of integers from 1 to 20 by steps of 2, one would use:

seq(from = 1, to = 20, by = 2)

 [1]  1  3  5  7  9 11 13 15 17 19

The default value of by is 1 and, given that the generate of a sequence of one value to another with steps of 1 is frequently used, there’s a shortcut:

seq(1, 5, 1)

[1] 1 2 3 4 5

seq(1, 5) ## default by

[1] 1 2 3 4 5

1:5

[1] 1 2 3 4 5

To generate a sequence of numbers from 1 to 20 of final length of 3, one would use:

seq(from = 1, to = 20, length.out = 3)

[1]  1.0 10.5 20.0

Random samples and permutations

A last group of useful functions are those that generate random data. The first one, sample, generates a random permutation of another vector. For example, to draw a random order to 10 students oral exam, I first assign each student a number from 1 to then (for instance based on the alphabetic order of their name) and then:

sample(1:10)

 [1]  9  4  7  1  2  5  3 10  6  8

Without further arguments, sample will return a permutation of all elements of the vector. If I want a random sample of a certain size, I would set this value as second argument. Below, I sample 5 random letters from the alphabet contained in the pre-defined letters vector:

sample(letters, 5)

[1] "s" "a" "u" "x" "j"

If I wanted an output larger than the input vector, or being able to draw some elements multiple times, I would need to set the replace argument to TRUE:

sample(1:5, 10, replace = TRUE)

 [1] 2 1 5 5 1 1 5 5 2 2

Homework Challenge:

When trying the functions above out, you will have realised that the samples are indeed random and that one doesn’t get the same permutation twice. To be able to reproduce these random draws, one can set the random number generation seed manually with set.seed() before drawing the random sample.

• Test this feature with your neighbour. First draw two random permutations of 1:10 independently and observe that you get different results.

• Now set the seed with, for example, set.seed(123) and repeat the random draw. Observe that you now get the same random draws.

• Repeat by setting a different seed.

Solution

Different permutations

sample(1:10)

 [1]  9  1  4  3  6  2  5  8 10  7

sample(1:10)

 [1]  4  9  7  6  1 10  8  3  2  5

Same permutations with seed 123

set.seed(123)
sample(1:10)

 [1]  3 10  2  8  6  9  1  7  5  4

set.seed(123)
sample(1:10)

 [1]  3 10  2  8  6  9  1  7  5  4

A different seed

set.seed(1)
sample(1:10)

 [1]  9  4  7  1  2  5  3 10  6  8

set.seed(1)
sample(1:10)

 [1]  9  4  7  1  2  5  3 10  6  8

Homework reading: Drawing samples from a normal distribution

The last function we are going to see is rnorm, that draws a random sample from a normal distribution. Two normal distributions of means 0 and 100 and standard deviations 1 and 5, noted noted N(0, 1) and N(100, 5), are shown below

The three arguments, n, mean and sd, define the size of the sample, and the parameters of the normal distribution, i.e the mean and its standard deviation. The defaults of the latter are 0 and 1.

rnorm(5)

[1]  0.69641761  0.05351568 -1.31028350 -2.12306606 -0.20807859

rnorm(5, 2, 2)

[1]  1.3744268 -0.1164714  2.8344472  1.3690969  3.6510983

rnorm(5, 100, 5)

[1] 106.45636  96.87448  95.62427 100.71678 107.12595

Now that we have learned how to write scripts, and the basics of R’s data structures, we are ready to start working with larger data, and learn about data frames.

Key Points

• How to interact with R

Starting with data

Overview

Teaching: XX min
Exercises: XX min

Questions

• First data analysis in R

Objectives

• Describe what a data frame is.

• Load external data from a .csv file into a data frame.

• Summarize the contents of a data frame.

• Describe what a factor is.

• Convert between strings and factors.

• Reorder and rename factors.

• Format dates.

• Export and save data.

Presentation of the gene expression data

We are going to use part of the data published by Blackmore et al. (2017), The effect of upper-respiratory infection on transcriptomic changes in the CNS. The goal of the study was to determine the effect of an upper-respiratory infection on changes in RNA transcription occuring in the cerebellum and spinal cord post infection. Gender matched eight week old C57BL/6 mice were inoculated with saline or with Influenza A by intranasal route and transcriptomic changes in the cerebellum and spinal cord tissues were evaluated by RNA-seq at days 0 (non-infected), 4 and 8.

The dataset is stored as a comma separated value (CSV) file. Each row holds information for a single RNA expression measurement, and the columns represent:

We have already downloaded the data and it can be found in the directory course-data/data/GSE96870/.

You are now ready to load the data:

rna <- read.csv("course-data/data/GSE96870/rnaseq.csv")

This statement doesn’t produce any output because, as you might recall, assignments don’t display anything. If we want to check that our data has been loaded, we can see the contents of the data frame by typing its name

rna

Wow… that was a lot of output. At least it means the data loaded properly. Let’s check the top (the first 6 lines) of this data frame using the function head():

head(rna)

     gene     sample expression     organism age    sex  infection  strain time
1     Asl GSM2545336       1170 Mus musculus   8 Female InfluenzaA C57BL/6    8
2    Apod GSM2545336      36194 Mus musculus   8 Female InfluenzaA C57BL/6    8
3 Cyp2d22 GSM2545336       4060 Mus musculus   8 Female InfluenzaA C57BL/6    8
4    Klk6 GSM2545336        287 Mus musculus   8 Female InfluenzaA C57BL/6    8
5   Fcrls GSM2545336         85 Mus musculus   8 Female InfluenzaA C57BL/6    8
6  Slc2a4 GSM2545336        782 Mus musculus   8 Female InfluenzaA C57BL/6    8
      tissue mouse ENTREZID
1 Cerebellum    14   109900
2 Cerebellum    14    11815
3 Cerebellum    14    56448
4 Cerebellum    14    19144
5 Cerebellum    14    80891
6 Cerebellum    14    20528
                                                                       product
1                               argininosuccinate lyase, transcript variant X1
2                                       apolipoprotein D, transcript variant 3
3 cytochrome P450, family 2, subfamily d, polypeptide 22, transcript variant 2
4                         kallikrein related-peptidase 6, transcript variant 2
5                Fc receptor-like S, scavenger receptor, transcript variant X1
6          solute carrier family 2 (facilitated glucose transporter), member 4
     ensembl_gene_id external_synonym chromosome_name   gene_biotype
1 ENSMUSG00000025533    2510006M18Rik               5 protein_coding
2 ENSMUSG00000022548             <NA>              16 protein_coding
3 ENSMUSG00000061740             2D22              15 protein_coding
4 ENSMUSG00000050063             Bssp               7 protein_coding
5 ENSMUSG00000015852    2810439C17Rik               3 protein_coding
6 ENSMUSG00000018566           Glut-4              11 protein_coding
                            phenotype_description
1           abnormal circulating amino acid level
2                      abnormal lipid homeostasis
3                        abnormal skin morphology
4                         abnormal cytokine level
5 decreased CD8-positive alpha-beta T cell number
6              abnormal circulating glucose level
  hsapiens_homolog_associated_gene_name
1                                   ASL
2                                  APOD
3                                CYP2D6
4                                  KLK6
5                                 FCRL2
6                                SLC2A4

## Try also
## View(rna)

Note

read.csv() assumes that fields are delineated by commas, however, in several countries, the comma is used as a decimal separator and the semicolon (;) is used as a field delineator. If you want to read in this type of files in R, you can use the read.csv2() function. It behaves exactly like read.csv() but uses different parameters for the decimal and the field separators. If you are working with another format, they can be both specified by the user. Check out the help for read.csv() by typing ?read.csv to learn more. There is also the read.delim() function for reading tab separated data files. It is important to note that all of these functions are actually wrapper functions for the main read.table() function with different arguments. As such, the data above could have also been loaded by using read.table() with the separation argument as ,. The code is as follows:

rna <- read.table(file = "course-data/data/GSE96870/rnaseq.csv",
                  sep = ",",
                  header = TRUE)

The header argument has to be set to TRUE to be able to read the headers as by default read.table() has the header argument set to FALSE.

What are data frames?

Data frames are the de facto data structure for most tabular data, and what we use for statistics and plotting.

A data frame can be created by hand, but most commonly they are generated by the functions read.csv() or read.table(); in other words, when importing spreadsheets from your hard drive (or the web).

A data frame is the representation of data in the format of a table where the columns are vectors that all have the same length. Because columns are vectors, each column must contain a single type of data (e.g., characters, integers, factors). For example, here is a figure depicting a data frame comprising a numeric, a character, and a logical vector.

We can see this when inspecting the structure of a data frame with the function str():

str(rna)

'data.frame':	32428 obs. of  19 variables:
 $ gene                                 : chr  "Asl" "Apod" "Cyp2d22" "Klk6" ...
 $ sample                               : chr  "GSM2545336" "GSM2545336" "GSM2545336" "GSM2545336" ...
 $ expression                           : int  1170 36194 4060 287 85 782 1619 288 43217 1071 ...
 $ organism                             : chr  "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" ...
 $ age                                  : int  8 8 8 8 8 8 8 8 8 8 ...
 $ sex                                  : chr  "Female" "Female" "Female" "Female" ...
 $ infection                            : chr  "InfluenzaA" "InfluenzaA" "InfluenzaA" "InfluenzaA" ...
 $ strain                               : chr  "C57BL/6" "C57BL/6" "C57BL/6" "C57BL/6" ...
 $ time                                 : int  8 8 8 8 8 8 8 8 8 8 ...
 $ tissue                               : chr  "Cerebellum" "Cerebellum" "Cerebellum" "Cerebellum" ...
 $ mouse                                : int  14 14 14 14 14 14 14 14 14 14 ...
 $ ENTREZID                             : int  109900 11815 56448 19144 80891 20528 97827 118454 18823 14696 ...
 $ product                              : chr  "argininosuccinate lyase, transcript variant X1" "apolipoprotein D, transcript variant 3" "cytochrome P450, family 2, subfamily d, polypeptide 22, transcript variant 2" "kallikrein related-peptidase 6, transcript variant 2" ...
 $ ensembl_gene_id                      : chr  "ENSMUSG00000025533" "ENSMUSG00000022548" "ENSMUSG00000061740" "ENSMUSG00000050063" ...
 $ external_synonym                     : chr  "2510006M18Rik" NA "2D22" "Bssp" ...
 $ chromosome_name                      : chr  "5" "16" "15" "7" ...
 $ gene_biotype                         : chr  "protein_coding" "protein_coding" "protein_coding" "protein_coding" ...
 $ phenotype_description                : chr  "abnormal circulating amino acid level" "abnormal lipid homeostasis" "abnormal skin morphology" "abnormal cytokine level" ...
 $ hsapiens_homolog_associated_gene_name: chr  "ASL" "APOD" "CYP2D6" "KLK6" ...

Inspecting data.frame Objects

We already saw how the functions head() and str() can be useful to check the content and the structure of a data frame. Here is a non-exhaustive list of functions to get a sense of the content/structure of the data. Let’s try them out!

Size:

• dim(rna) - returns a vector with the number of rows in the first element, and the number of columns as the second element (the dimensions of the object)
• nrow(rna) - returns the number of rows
• ncol(rna) - returns the number of columns

Content:

• head(rna) - shows the first 6 rows
• tail(rna) - shows the last 6 rows

Names:

• names(rna) - returns the column names (synonym of colnames() for data.frame objects)
• rownames(rna) - returns the row names

Summary:

• str(rna) - structure of the object and information about the class, length and content of each column
• summary(rna) - summary statistics for each column

Note: most of these functions are “generic”, they can be used on other types of objects besides data.frame.

Challenge:

Based on the output of str(rna), can you answer the following questions?

• What is the class of the object rna?
• How many rows and how many columns are in this object?
• How many genes (as defined by the gene variable) have been measured in this experiment?

Solution

• class: data frame
• how many rows: 32428, how many columns: 19
• how many genes: 1477

Indexing and subsetting data frames

Our rna data frame has rows and columns (it has 2 dimensions), if we want to extract some specific data from it, we need to specify the “coordinates” we want from it. Row numbers come first, followed by column numbers. However, note that different ways of specifying these coordinates lead to results with different classes.

# first element in the first column of the data frame (as a vector)
rna[1, 1]
# first element in the 6th column (as a vector)
rna[1, 6]
# first column of the data frame (as a vector)
rna[, 1]
# first column of the data frame (as a data.frame)
rna[1]
# first three elements in the 7th column (as a vector)
rna[1:3, 7]
# the 3rd row of the data frame (as a data.frame)
rna[3, ]
# equivalent to head_rna <- head(rna)
head_rna <- rna[1:6, ]
head_rna

: is a special function that creates numeric vectors of integers in increasing or decreasing order, test 1:10 and 10:1 for instance. See section \@ref(sec:genvec) for details.

You can also exclude certain indices of a data frame using the “-” sign:

rna[, -1]          ## The whole data frame, except the first column
rna[-c(7:32428), ] ## Equivalent to head(rna)

Data frames can be subset by calling indices (as shown previously), but also by calling their column names directly:

rna["gene"]       # Result is a data.frame
rna[, "gene"]     # Result is a vector
rna[["gene"]]     # Result is a vector
rna$gene          # Result is a vector

In RStudio, you can use the autocompletion feature to get the full and correct names of the columns.

Challenge

1. Create a data.frame (rna_200) containing only the data in row 200 of the rna dataset.

2. Notice how nrow() gave you the number of rows in a data.frame?

• Use that number to pull out just that last row in the inital rna data frame.

• Compare that with what you see as the last row using tail() to make sure it’s meeting expectations.

• Pull out that last row using nrow() instead of the row number.

• Create a new data frame (rna_last) from that last row.

1. Use nrow() to extract the row that is in the middle of the rna dataframe. Store the content of this row in an object named rna_middle.

2. Combine nrow() with the - notation above to reproduce the behavior of head(rna), keeping just the first through 6th rows of the rna dataset.

Solution

## 1.
rna_200 <- rna[200, ]
## 2.
## Saving `n_rows` to improve readability and reduce duplication
n_rows <- nrow(rna)
rna_last <- rna[n_rows, ]
## 3.
rna_middle <- rna[n_rows / 2, ]
## 4.
rna_head <- rna[-(7:n_rows), ]

Factors

Factors represent categorical data. They are stored as integers associated with labels and they can be ordered or unordered. While factors look (and often behave) like character vectors, they are actually treated as integer vectors by R. So you need to be very careful when treating them as strings.

Once created, factors can only contain a pre-defined set of values, known as levels. By default, R always sorts levels in alphabetical order. For instance, if you have a factor with 2 levels:

sex <- factor(c("male", "female", "female", "male", "female"))

R will assign 1 to the level "female" and 2 to the level "male" (because f comes before m, even though the first element in this vector is "male"). You can see this by using the function levels() and you can find the number of levels using nlevels():

levels(sex)

[1] "female" "male"  

nlevels(sex)

[1] 2

Sometimes, the order of the factors does not matter, other times you might want to specify the order because it is meaningful (e.g., “low”, “medium”, “high”), it improves your visualization, or it is required by a particular type of analysis. Here, one way to reorder our levels in the sex vector would be:

sex ## current order

[1] male   female female male   female
Levels: female male

sex <- factor(sex, levels = c("male", "female"))
sex ## after re-ordering

[1] male   female female male   female
Levels: male female

In R’s memory, these factors are represented by integers (1, 2, 3), but are more informative than integers because factors are self describing: "female", "male" is more descriptive than 1, 2. Which one is “male”? You wouldn’t be able to tell just from the integer data. Factors, on the other hand, have this information built in. It is particularly helpful when there are many levels (like the species names in our example dataset).

Converting to factors {-}

If you need to convert a factor to a character vector, you use as.character(x).

as.character(sex)

[1] "male"   "female" "female" "male"   "female"

Renaming factors {-}

When your data is stored as a factor, you can use the plot() function to get a quick glance at the number of observations represented by each factor level. Let’s look at the number of males and females in our data.

plot(sex)

If we want to rename these factor, it is sufficient to change its levels:

levels(sex)

[1] "male"   "female"

levels(sex) <- c("M", "F")
sex

[1] M F F M F
Levels: M F

plot(sex)

Homework Challenge:

• Rename “female” and “male” to “Female” and “Male” respectively.

Challenge:

We have seen how data frames are created when using read.csv(), but they can also be created by hand with the data.frame() function. There are a few mistakes in this hand-crafted data.frame. Can you spot and fix them? Don’t hesitate to experiment!

animal_data <- data.frame(
       animal = c(dog, cat, sea cucumber, sea urchin),
       feel = c("furry", "squishy", "spiny"),
       weight = c(45, 8 1.1, 0.8))

Solution

• missing quotations around the names of the animals
• missing one entry in the “feel” column (probably for one of the furry animals)
• missing one comma in the weight column

Homework Challenge:

Can you predict the class for each of the columns in the following example?

Check your guesses using str(country_climate):

• Are they what you expected? Why? Why not?

• Try again by adding stringsAsFactors = TRUE after the last variable when creating the data frame? What is happening now? stringsAsFactors can also be set when reading text-based spreadsheets into R using read.csv().

country_climate <- data.frame(
       country = c("Canada", "Panama", "South Africa", "Australia"),
       climate = c("cold", "hot", "temperate", "hot/temperate"),
       temperature = c(10, 30, 18, "15"),
       northern_hemisphere = c(TRUE, TRUE, FALSE, "FALSE"),
       has_kangaroo = c(FALSE, FALSE, FALSE, 1)
       )

The automatic conversion of data type is sometimes a blessing, sometimes an annoyance. Be aware that it exists, learn the rules, and double check that data you import in R are of the correct type within your data frame. If not, use it to your advantage to detect mistakes that might have been introduced during data entry (a letter in a column that should only contain numbers for instance).

Learn more in this RStudio tutorial

Matrices

Before proceeding, now that we have learnt about dataframes, let’s recap package installation and learn about a new data type, namely the matrix. Like a data.frame, a matrix has two dimensions, rows and columns. But the major difference is that all cells in a matrix must be of the same type: numeric, character, logical, … In that respect, matrices are closer to a vector than a data.frame.

The default constructor for a matrix is matrix. It takes a vector of values to populate the matrix and the number of row and/or columns¹. The values are sorted along the columns, as illustrated below.

m <- matrix(1:9, ncol = 3, nrow = 3)
m

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

Homework Challenge:

Using the function installed.packages(), create a character matrix containing the information about all packages currently installed on your computer. Explore it.

Solution:

## create the matrix
ip <- installed.packages()
head(ip)
## try also View(ip)
## number of package
nrow(ip)
## names of all installed packages
rownames(ip)
## type of information we have about each package
colnames(ip)

It is often useful to create large random data matrices as test data. The exercise below asks you to create such a matrix with random data drawn from a normal distribution of mean 0 and standard deviation 1, which can be done with the rnorm() function.

Homework Challenge:

Construct a matrix of dimension 1000 by 3 of normally distributed data (mean 0, standard deviation 1)

Solution

set.seed(123)
m <- matrix(rnorm(3000), ncol = 3)
dim(m)

[1] 1000    3

head(m)

            [,1]        [,2]       [,3]
[1,] -0.56047565 -0.99579872 -0.5116037
[2,] -0.23017749 -1.03995504  0.2369379
[3,]  1.55870831 -0.01798024 -0.5415892
[4,]  0.07050839 -0.13217513  1.2192276
[5,]  0.12928774 -2.54934277  0.1741359
[6,]  1.71506499  1.04057346 -0.6152683

Summary of R objects

So far, we have seen several types of R object varying in the number of dimensions and whether they could store a single of multiple data types:

• vector: one dimension (they have a length), single type of data.
• matrix: two dimensions, single type of data.
• data.frame: two dimensions, one type per column.

Lists

A data type that we haven’t seen yet, but that is useful to know, and follows from the summary that we have just seen are lists:

• list: one dimension, every item can be of a different data type.

Below, let’s create a list containing a vector of numbers, characters, a matrix, a dataframe and another list:

l <- list(1:10, ## numeric
          letters, ## character
          installed.packages(), ## a matrix
          cars, ## a data.frame
          list(1, 2, 3)) ## a list
length(l)

[1] 5

str(l)

List of 5
 $ : int [1:10] 1 2 3 4 5 6 7 8 9 10
 $ : chr [1:26] "a" "b" "c" "d" ...
 $ : chr [1:399, 1:16] "abind" "annotate" "AnnotationDbi" "AnnotationHub" ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:399] "abind" "annotate" "AnnotationDbi" "AnnotationHub" ...
  .. ..$ : chr [1:16] "Package" "LibPath" "Version" "Priority" ...
 $ :'data.frame':	50 obs. of  2 variables:
  ..$ speed: num [1:50] 4 4 7 7 8 9 10 10 10 11 ...
  ..$ dist : num [1:50] 2 10 4 22 16 10 18 26 34 17 ...
 $ :List of 3
  ..$ : num 1
  ..$ : num 2
  ..$ : num 3

List subsetting is done using [] to subset a new sub-list or [[]] to extract a single element of that list (using indices or names, of the list is named).

l[[1]] ## first element

 [1]  1  2  3  4  5  6  7  8  9 10

l[1:2] ## a list of length 2

[[1]]
 [1]  1  2  3  4  5  6  7  8  9 10

[[2]]
 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
[20] "t" "u" "v" "w" "x" "y" "z"

l[1]   ## a list of length 1

[[1]]
 [1]  1  2  3  4  5  6  7  8  9 10

Exporting and saving data

Exporting tabular data {-}

We have seen how to read a text-based spreadsheet into R using the read.table family of functions. To export a data.frame to a text-based spreadsheet, we can use the write.table set of functions (write.csv, write.delim, …). They all take the variable to be exported and the file to be exported to. For example, to export the rna data to an rna.csv file in the data_output directory, we would execute:

write.csv(rna, file = "data_output/rna.csv")

This new csv file can now be shared with other collaborators who aren’t familiar with R.

Saving data {-}

Exporting data to a spreadsheet has several limitations, such as those described in the first chapter such as possible inconsistencies with , and . for decimal separators and lack of variable type definitions. Furthermore, exporting data to a spreadsheet is only relevant for rectangular data such as dataframes and matrices.

A more general way to save data, that is specific to R and is guaranteed to work on any operating system, is to use the save function. Saving objects will generate a binary representation of the object on disk, a R Data file (rda extension) that guarantees to produce the same object once loaded back into R using the load function.

save(rna, file = "data_output/rna.rda")
rm(rna)
load("data_output/rna.rda")
head(rna)

Note about how the function load loads the object in the file directly in the global environment.

There is also the saveRDS and readRDS functions that save R objects to binary files (using the rds extension here) and read these back into R. From a user’s perspective, main difference is that, load loads an object in the global environment while readRDS reads the data from disk and returns it. It is thus necessary to store the output of readRDS in a variable:

saveRDS(rna, file = "data_output/rna.rds")
rm(rna)
rna <- readRDS("data_output/rna.rds")
head(rna)

To conclude, when it comes to saving data from R that will be loaded again in R, saving and loading is the preferred approach. If tabular data need to be shared with somebody that is not using R, then exporting to a text-based spreadsheet is a good alternative.

Key Points

• Tabular data in R

Manipulating and analyzing data with dplyr

Overview

Teaching: XX min
Exercises: XX min

Questions

• Data analysis in R using the tidyverse meta-package

Objectives

• Describe the purpose of the dplyr and tidyr packages.

• Describe several of their functions that are extremely useful to manipulate data.

• Describe the concept of a wide and a long table format, and see how to reshape a data frame from one format to the other one.

Data Manipulation using dplyr and tidyr

Bracket subsetting is handy, but it can be cumbersome and difficult to read, especially for complicated operations.

Some packages can greatly facilitate our task when we manipulate data. Packages in R are basically sets of additional functions that let you do more stuff. The functions we’ve been using so far, like str() or data.frame(), come built into R; Loading packages can give you access to other specific functions. Before you use a package for the first time you need to install it on your machine, and then you should import it in every subsequent R session when you need it.

• The package dplyr provides powerful tools for data manipulation tasks. It is built to work directly with data frames, with many manipulation tasks optimized.

• As we will see latter on, sometimes we want a data frame to be reshaped to be able to do some specific analyses or for visualization. The package tidyr addresses this common problem of reshaping data and provides tools for manipulating data in a tidy way.

To learn more about dplyr and tidyr after the workshop, you may want to check out this handy data transformation with dplyr cheatsheet and this one about tidyr.

• The tidyverse package is an “umbrella-package” that installs several useful packages for data analysis which work well together, such as tidyr, dplyr, ggplot2, tibble, etc. These packages help us to work and interact with the data. They allow us to do many things with your data, such as subsetting, transforming, visualizing, etc.

To install and load the tidyverse package type:

# no need to run this command as it's already been installed for you
# BiocManager::install("tidyverse")

## load the tidyverse packages, incl. dplyr
library("tidyverse")

Loading data with tidyverse

Instead of read.csv(), we will read in our data using the read_csv() function, from the tidyverse package readr, .

rna <- read_csv("course-data/data/GSE96870/rnaseq.csv")

## view the data
rna

# A tibble: 32,428 × 19
   gene    sample  expression organism   age sex   infection strain  time tissue
   <chr>   <chr>        <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Asl     GSM254…       1170 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 2 Apod    GSM254…      36194 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 3 Cyp2d22 GSM254…       4060 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 4 Klk6    GSM254…        287 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 5 Fcrls   GSM254…         85 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 6 Slc2a4  GSM254…        782 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 7 Exd2    GSM254…       1619 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 8 Gjc2    GSM254…        288 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 9 Plp1    GSM254…      43217 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
10 Gnb4    GSM254…       1071 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
# … with 32,418 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

Notice that the class of the data is now referred to as a “tibble”.

Tibbles tweak some of the behaviors of the data frame objects we introduced in the previously. The data structure is very similar to a data frame. For our purposes the only differences are that:

1. It displays the data type of each column under its name. Note that <dbl> is a data type defined to hold numeric values with decimal points.

2. It only prints the first few rows of data and only as many columns as fit on one screen.

We are now going to learn some of the most common dplyr functions:

• select(): subset columns
• filter(): subset rows on conditions
• mutate(): create new columns by using information from other columns
• group_by() and summarize(): create summary statistics on grouped data
• arrange(): sort results
• count(): count discrete values

Selecting columns and filtering rows

To select columns of a data frame, use select(). The first argument to this function is the data frame (rna), and the subsequent arguments are the columns to keep.

select(rna, gene, sample, tissue, expression)

# A tibble: 32,428 × 4
   gene    sample     tissue     expression
   <chr>   <chr>      <chr>           <dbl>
 1 Asl     GSM2545336 Cerebellum       1170
 2 Apod    GSM2545336 Cerebellum      36194
 3 Cyp2d22 GSM2545336 Cerebellum       4060
 4 Klk6    GSM2545336 Cerebellum        287
 5 Fcrls   GSM2545336 Cerebellum         85
 6 Slc2a4  GSM2545336 Cerebellum        782
 7 Exd2    GSM2545336 Cerebellum       1619
 8 Gjc2    GSM2545336 Cerebellum        288
 9 Plp1    GSM2545336 Cerebellum      43217
10 Gnb4    GSM2545336 Cerebellum       1071
# … with 32,418 more rows

To select all columns except certain ones, put a “-“ in front of the variable to exclude it.

select(rna, -tissue, -organism)

# A tibble: 32,428 × 17
   gene    sample   expression   age sex   infection strain  time mouse ENTREZID
   <chr>   <chr>         <dbl> <dbl> <chr> <chr>     <chr>  <dbl> <dbl>    <dbl>
 1 Asl     GSM2545…       1170     8 Fema… Influenz… C57BL…     8    14   109900
 2 Apod    GSM2545…      36194     8 Fema… Influenz… C57BL…     8    14    11815
 3 Cyp2d22 GSM2545…       4060     8 Fema… Influenz… C57BL…     8    14    56448
 4 Klk6    GSM2545…        287     8 Fema… Influenz… C57BL…     8    14    19144
 5 Fcrls   GSM2545…         85     8 Fema… Influenz… C57BL…     8    14    80891
 6 Slc2a4  GSM2545…        782     8 Fema… Influenz… C57BL…     8    14    20528
 7 Exd2    GSM2545…       1619     8 Fema… Influenz… C57BL…     8    14    97827
 8 Gjc2    GSM2545…        288     8 Fema… Influenz… C57BL…     8    14   118454
 9 Plp1    GSM2545…      43217     8 Fema… Influenz… C57BL…     8    14    18823
10 Gnb4    GSM2545…       1071     8 Fema… Influenz… C57BL…     8    14    14696
# … with 32,418 more rows, and 7 more variables: product <chr>,
#   ensembl_gene_id <chr>, external_synonym <chr>, chromosome_name <chr>,
#   gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

This will select all the variables in rna except tissue and organism.

To choose rows based on a specific criteria, use filter():

filter(rna, sex == "Male")

# A tibble: 14,740 × 19
   gene    sample  expression organism   age sex   infection strain  time tissue
   <chr>   <chr>        <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Asl     GSM254…        626 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 2 Apod    GSM254…      13021 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 3 Cyp2d22 GSM254…       2171 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 4 Klk6    GSM254…        448 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 5 Fcrls   GSM254…        180 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 6 Slc2a4  GSM254…        313 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 7 Exd2    GSM254…       2366 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 8 Gjc2    GSM254…        310 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 9 Plp1    GSM254…      53126 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
10 Gnb4    GSM254…       1355 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
# … with 14,730 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

filter(rna, sex == "Male" & infection == "NonInfected")

# A tibble: 4,422 × 19
   gene    sample  expression organism   age sex   infection strain  time tissue
   <chr>   <chr>        <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Asl     GSM254…        535 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 2 Apod    GSM254…      13668 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 3 Cyp2d22 GSM254…       2008 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 4 Klk6    GSM254…       1101 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 5 Fcrls   GSM254…        375 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 6 Slc2a4  GSM254…        249 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 7 Exd2    GSM254…       3126 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 8 Gjc2    GSM254…        791 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 9 Plp1    GSM254…      98658 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
10 Gnb4    GSM254…       2437 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
# … with 4,412 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

Now let’s imagine we are interested in the human homologs of the mouse genes analysed in this dataset. This information can be found in the last column of the rna tibble, named hsapiens_homolog_associated_gene_name. Some mouse gene have no human homologs. These can be retrieved using a filter() in the chain, and the is.na() function that determines whether something is an NA.

rna_NA <- filter(rna, is.na(hsapiens_homolog_associated_gene_name))
select(rna_NA, gene, hsapiens_homolog_associated_gene_name)

# A tibble: 4,290 × 2
   gene     hsapiens_homolog_associated_gene_name
   <chr>    <chr>                                
 1 Prodh    <NA>                                 
 2 Tssk5    <NA>                                 
 3 Vmn2r1   <NA>                                 
 4 Gm10654  <NA>                                 
 5 Hexa     <NA>                                 
 6 Sult1a1  <NA>                                 
 7 Gm6277   <NA>                                 
 8 Tmem198b <NA>                                 
 9 Adam1a   <NA>                                 
10 Ebp      <NA>                                 
# … with 4,280 more rows

If we want to keep only mouse gene that have a human homolog, we can insert a “!” symbol that negates the result, so we’re asking for every row where hsapiens_homolog_associated_gene_name is not an NA.

rna_no_NA <- filter(rna, !is.na(hsapiens_homolog_associated_gene_name))
select(rna_no_NA, gene, hsapiens_homolog_associated_gene_name)

# A tibble: 28,138 × 2
   gene    hsapiens_homolog_associated_gene_name
   <chr>   <chr>                                
 1 Asl     ASL                                  
 2 Apod    APOD                                 
 3 Cyp2d22 CYP2D6                               
 4 Klk6    KLK6                                 
 5 Fcrls   FCRL2                                
 6 Slc2a4  SLC2A4                               
 7 Exd2    EXD2                                 
 8 Gjc2    GJC2                                 
 9 Plp1    PLP1                                 
10 Gnb4    GNB4                                 
# … with 28,128 more rows

Pipes

What if you want to select and filter at the same time? There are three ways to do this: use intermediate steps, nested functions, or pipes.

With intermediate steps, you create a temporary data frame and use that as input to the next function, like this:

rna2 <- filter(rna, sex == "Male")
rna3 <- select(rna2, gene, sample, tissue, expression)
rna3

# A tibble: 14,740 × 4
   gene    sample     tissue     expression
   <chr>   <chr>      <chr>           <dbl>
 1 Asl     GSM2545340 Cerebellum        626
 2 Apod    GSM2545340 Cerebellum      13021
 3 Cyp2d22 GSM2545340 Cerebellum       2171
 4 Klk6    GSM2545340 Cerebellum        448
 5 Fcrls   GSM2545340 Cerebellum        180
 6 Slc2a4  GSM2545340 Cerebellum        313
 7 Exd2    GSM2545340 Cerebellum       2366
 8 Gjc2    GSM2545340 Cerebellum        310
 9 Plp1    GSM2545340 Cerebellum      53126
10 Gnb4    GSM2545340 Cerebellum       1355
# … with 14,730 more rows

This is readable, but can clutter up your workspace with lots of intermediate objects that you have to name individually. With multiple steps, that can be hard to keep track of.

You can also nest functions (i.e. one function inside of another), like this:

rna3 <- select(filter(rna, sex == "Male"), gene, sample, tissue, expression)
rna3

# A tibble: 14,740 × 4
   gene    sample     tissue     expression
   <chr>   <chr>      <chr>           <dbl>
 1 Asl     GSM2545340 Cerebellum        626
 2 Apod    GSM2545340 Cerebellum      13021
 3 Cyp2d22 GSM2545340 Cerebellum       2171
 4 Klk6    GSM2545340 Cerebellum        448
 5 Fcrls   GSM2545340 Cerebellum        180
 6 Slc2a4  GSM2545340 Cerebellum        313
 7 Exd2    GSM2545340 Cerebellum       2366
 8 Gjc2    GSM2545340 Cerebellum        310
 9 Plp1    GSM2545340 Cerebellum      53126
10 Gnb4    GSM2545340 Cerebellum       1355
# … with 14,730 more rows

This is handy, but can be difficult to read if too many functions are nested, as R evaluates the expression from the inside out (in this case, filtering, then selecting).

The last option, pipes, are a recent addition to R. Pipes let you take the output of one function and send it directly to the next, which is useful when you need to do many things to the same dataset.

Pipes in R look like %>% and are made available via the magrittr package, installed automatically with dplyr. If you use RStudio, you can type the pipe with Ctrl

• Shift + M if you have a PC or Cmd + Shift + M if you have a Mac.

In the above code, we use the pipe to send the rna dataset first through filter() to keep rows where sex is Male, then through select() to keep only the gene, sample, tissue, and expressioncolumns.

The pipe %>% takes the object on its left and passes it directly as the first argument to the function on its right, we don’t need to explicitly include the data frame as an argument to the filter() and select() functions any more.

rna %>%
  filter(sex == "Male") %>%
  select(gene, sample, tissue, expression)

# A tibble: 14,740 × 4
   gene    sample     tissue     expression
   <chr>   <chr>      <chr>           <dbl>
 1 Asl     GSM2545340 Cerebellum        626
 2 Apod    GSM2545340 Cerebellum      13021
 3 Cyp2d22 GSM2545340 Cerebellum       2171
 4 Klk6    GSM2545340 Cerebellum        448
 5 Fcrls   GSM2545340 Cerebellum        180
 6 Slc2a4  GSM2545340 Cerebellum        313
 7 Exd2    GSM2545340 Cerebellum       2366
 8 Gjc2    GSM2545340 Cerebellum        310
 9 Plp1    GSM2545340 Cerebellum      53126
10 Gnb4    GSM2545340 Cerebellum       1355
# … with 14,730 more rows

Some may find it helpful to read the pipe like the word “then”. For instance, in the above example, we took the data frame rna, then we filtered for rows with sex == "Male", then we selected columns gene, sample, tissue, and expression.

The dplyr functions by themselves are somewhat simple, but by combining them into linear workflows with the pipe, we can accomplish more complex manipulations of data frames.

If we want to create a new object with this smaller version of the data, we can assign it a new name:

rna3 <- rna %>%
  filter(sex == "Male") %>%
  select(gene, sample, tissue, expression)

rna3

# A tibble: 14,740 × 4
   gene    sample     tissue     expression
   <chr>   <chr>      <chr>           <dbl>
 1 Asl     GSM2545340 Cerebellum        626
 2 Apod    GSM2545340 Cerebellum      13021
 3 Cyp2d22 GSM2545340 Cerebellum       2171
 4 Klk6    GSM2545340 Cerebellum        448
 5 Fcrls   GSM2545340 Cerebellum        180
 6 Slc2a4  GSM2545340 Cerebellum        313
 7 Exd2    GSM2545340 Cerebellum       2366
 8 Gjc2    GSM2545340 Cerebellum        310
 9 Plp1    GSM2545340 Cerebellum      53126
10 Gnb4    GSM2545340 Cerebellum       1355
# … with 14,730 more rows

Homework Challenge:

Using pipes, subset the rna data to keep genes with an expression higher than 50000 in female mice at time 0, and retain only the columns gene, sample, time, expression and age.

Solution

rna %>%
  filter(expression > 50000,
         sex == "Female",
         time == 0 ) %>%
  select(gene, sample, time, expression, age)

# A tibble: 9 × 5
  gene   sample      time expression   age
  <chr>  <chr>      <dbl>      <dbl> <dbl>
1 Plp1   GSM2545337     0     101241     8
2 Atp1b1 GSM2545337     0      53260     8
3 Plp1   GSM2545338     0      96534     8
4 Atp1b1 GSM2545338     0      50614     8
5 Plp1   GSM2545348     0     102790     8
6 Atp1b1 GSM2545348     0      59544     8
7 Plp1   GSM2545353     0      71237     8
8 Glul   GSM2545353     0      52451     8
9 Atp1b1 GSM2545353     0      61451     8

Mutate

Frequently you’ll want to create new columns based on the values of existing columns, for example to do unit conversions, or to find the ratio of values in two columns. For this we’ll use mutate().

To create a new column of time in hours:

rna %>%
  mutate(time_hours = time * 24) %>%
  select(time, time_hours)

# A tibble: 32,428 × 2
    time time_hours
   <dbl>      <dbl>
 1     8        192
 2     8        192
 3     8        192
 4     8        192
 5     8        192
 6     8        192
 7     8        192
 8     8        192
 9     8        192
10     8        192
# … with 32,418 more rows

You can also create a second new column based on the first new column within the same call of mutate():

rna %>%
  mutate(time_hours = time * 24,
         time_mn = time_hours * 60) %>%
  select(time, time_hours, time_mn)

# A tibble: 32,428 × 3
    time time_hours time_mn
   <dbl>      <dbl>   <dbl>
 1     8        192   11520
 2     8        192   11520
 3     8        192   11520
 4     8        192   11520
 5     8        192   11520
 6     8        192   11520
 7     8        192   11520
 8     8        192   11520
 9     8        192   11520
10     8        192   11520
# … with 32,418 more rows

Homework Challenge

Create a new data frame from the rna data that meets the following criteria: contains only the gene, chromosome_name, phenotype_description, sample, and expression columns and a new column giving the log expression the gene. This data frame must only contain genes located on autosomes and associated with a phenotype_description.

Hint: think about how the commands should be ordered to produce this data frame!

Solution

rna %>%
  filter(chromosome_name != "X", chromosome_name != "Y") %>%
  mutate(log_expression = log(expression)) %>%
  select(gene, chromosome_name, phenotype_description, sample, log_expression) %>%
  filter(!is.na(phenotype_description))

# A tibble: 21,054 × 5
   gene    chromosome_name phenotype_description           sample log_expression
   <chr>   <chr>           <chr>                           <chr>           <dbl>
 1 Asl     5               abnormal circulating amino aci… GSM25…           7.06
 2 Apod    16              abnormal lipid homeostasis      GSM25…          10.5 
 3 Cyp2d22 15              abnormal skin morphology        GSM25…           8.31
 4 Klk6    7               abnormal cytokine level         GSM25…           5.66
 5 Fcrls   3               decreased CD8-positive alpha-b… GSM25…           4.44
 6 Slc2a4  11              abnormal circulating glucose l… GSM25…           6.66
 7 Gjc2    11              Purkinje cell degeneration      GSM25…           5.66
 8 Gnb4    3               decreased anxiety-related resp… GSM25…           6.98
 9 Tnc     4               abnormal CNS synaptic transmis… GSM25…           5.39
10 Trf     9               abnormal circulating phosphate… GSM25…           9.18
# … with 21,044 more rows

Split-apply-combine data analysis

Many data analysis tasks can be approached using the split-apply-combine paradigm: split the data into groups, apply some analysis to each group, and then combine the results. dplyr makes this very easy through the use of the group_by() function.

rna %>%
  group_by(gene)

# A tibble: 32,428 × 19
# Groups:   gene [1,474]
   gene    sample  expression organism   age sex   infection strain  time tissue
   <chr>   <chr>        <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Asl     GSM254…       1170 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 2 Apod    GSM254…      36194 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 3 Cyp2d22 GSM254…       4060 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 4 Klk6    GSM254…        287 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 5 Fcrls   GSM254…         85 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 6 Slc2a4  GSM254…        782 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 7 Exd2    GSM254…       1619 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 8 Gjc2    GSM254…        288 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 9 Plp1    GSM254…      43217 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
10 Gnb4    GSM254…       1071 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
# … with 32,418 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

The group_by() function doesn’t perform any data processing, it groups the data into subsets: in the example above, our initial tibble of 32428 observations is split into 1474 groups based on the gene variable.

We could similarly decide to group the tibble by the samples:

rna %>%
  group_by(sample)

# A tibble: 32,428 × 19
# Groups:   sample [22]
   gene    sample  expression organism   age sex   infection strain  time tissue
   <chr>   <chr>        <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Asl     GSM254…       1170 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 2 Apod    GSM254…      36194 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 3 Cyp2d22 GSM254…       4060 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 4 Klk6    GSM254…        287 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 5 Fcrls   GSM254…         85 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 6 Slc2a4  GSM254…        782 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 7 Exd2    GSM254…       1619 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 8 Gjc2    GSM254…        288 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 9 Plp1    GSM254…      43217 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
10 Gnb4    GSM254…       1071 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
# … with 32,418 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

Here our initial tibble of 32428 observations is split into 22 groups based on the sample variable.

Once the data have been combined, subsequent operations will be applied on each group independently.

The summarize() function

group_by() is often used together with summarize(), which collapses each group into a single-row summary of that group.

group_by() takes as arguments the column names that contain the categorical variables for which you want to calculate the summary statistics. So to compute the mean expression by gene:

rna %>%
  group_by(gene) %>%
  summarize(mean_expression = mean(expression))

# A tibble: 1,474 × 2
   gene    mean_expression
   <chr>             <dbl>
 1 Aamp            4751.  
 2 Abca12             4.55
 3 Abcc8           2498.  
 4 Abhd14a          525.  
 5 Abi2            4909.  
 6 Abi3bp          1002.  
 7 Abl2            2124.  
 8 Acadl           2053.  
 9 Acap3           3536.  
10 Acbd4           1431.  
# … with 1,464 more rows

We could also want to calculate the mean expression levels of all genes in each sample:

rna %>%
  group_by(sample) %>%
  summarize(mean_expression = mean(expression))

# A tibble: 22 × 2
   sample     mean_expression
   <chr>                <dbl>
 1 GSM2545336           2062.
 2 GSM2545337           1766.
 3 GSM2545338           1668.
 4 GSM2545339           1696.
 5 GSM2545340           1682.
 6 GSM2545341           1638.
 7 GSM2545342           1594.
 8 GSM2545343           2107.
 9 GSM2545344           1712.
10 GSM2545345           1700.
# … with 12 more rows

But we can can also group by multiple columns:

rna %>%
  group_by(gene, infection, time) %>%
  summarize(mean_expression = mean(expression))

`summarise()` has grouped output by 'gene', 'infection'. You can override using
the `.groups` argument.

# A tibble: 4,422 × 4
# Groups:   gene, infection [2,948]
   gene    infection    time mean_expression
   <chr>   <chr>       <dbl>           <dbl>
 1 Aamp    InfluenzaA      4         4870   
 2 Aamp    InfluenzaA      8         4763.  
 3 Aamp    NonInfected     0         4603.  
 4 Abca12  InfluenzaA      4            4.25
 5 Abca12  InfluenzaA      8            4.14
 6 Abca12  NonInfected     0            5.29
 7 Abcc8   InfluenzaA      4         2609.  
 8 Abcc8   InfluenzaA      8         2292.  
 9 Abcc8   NonInfected     0         2576.  
10 Abhd14a InfluenzaA      4          547.  
# … with 4,412 more rows

Once the data is grouped, you can also summarize multiple variables at the same time (and not necessarily on the same variable). For instance, we could add a column indicating the median expression by gene and by condition:

rna %>%
  group_by(gene, infection, time) %>%
  summarize(mean_expression = mean(expression),
            median_expression = median(expression))

`summarise()` has grouped output by 'gene', 'infection'. You can override using
the `.groups` argument.

# A tibble: 4,422 × 5
# Groups:   gene, infection [2,948]
   gene    infection    time mean_expression median_expression
   <chr>   <chr>       <dbl>           <dbl>             <dbl>
 1 Aamp    InfluenzaA      4         4870               4708  
 2 Aamp    InfluenzaA      8         4763.              4813  
 3 Aamp    NonInfected     0         4603.              4717  
 4 Abca12  InfluenzaA      4            4.25               4.5
 5 Abca12  InfluenzaA      8            4.14               4  
 6 Abca12  NonInfected     0            5.29               5  
 7 Abcc8   InfluenzaA      4         2609.              2424. 
 8 Abcc8   InfluenzaA      8         2292.              2224  
 9 Abcc8   NonInfected     0         2576.              2578  
10 Abhd14a InfluenzaA      4          547.               523  
# … with 4,412 more rows

Challenge

Calculate the mean expression level of gene “Dok3” by timepoints.

Solution

rna %>%
  filter(gene == "Dok3") %>%
  group_by(time) %>%
  summarize(mean = mean(expression))

# A tibble: 3 × 2
   time  mean
  <dbl> <dbl>
1     0  169 
2     4  156.
3     8   61 

Counting

When working with data, we often want to know the number of observations found for each factor or combination of factors. For this task, dplyr provides count(). For example, if we wanted to count the number of rows of data for each infected and non infected samples, we would do:

rna %>%
    count(infection)

# A tibble: 2 × 2
  infection       n
  <chr>       <int>
1 InfluenzaA  22110
2 NonInfected 10318

The count() function is shorthand for something we’ve already seen: grouping by a variable, and summarizing it by counting the number of observations in that group. In other words, rna %>% count() is equivalent to:

rna %>%
    group_by(infection) %>%
    summarise(n = n())

# A tibble: 2 × 2
  infection       n
  <chr>       <int>
1 InfluenzaA  22110
2 NonInfected 10318

Previous example shows the use of count() to count the number of rows/observations for one factor (i.e., infection). If we wanted to count combination of factors, such as infection and time, we would specify the first and the second factor as the arguments of count():

rna %>%
    count(infection, time)

# A tibble: 3 × 3
  infection    time     n
  <chr>       <dbl> <int>
1 InfluenzaA      4 11792
2 InfluenzaA      8 10318
3 NonInfected     0 10318

which is equivalent to this:

rna %>%
  group_by(infection, time) %>%
  summarize(n = n())

`summarise()` has grouped output by 'infection'. You can override using the
`.groups` argument.

# A tibble: 3 × 3
# Groups:   infection [2]
  infection    time     n
  <chr>       <dbl> <int>
1 InfluenzaA      4 11792
2 InfluenzaA      8 10318
3 NonInfected     0 10318

It is sometimes useful to sort the result to facilitate the comparisons. We can use arrange() to sort the table. For instance, we might want to arrange the table above by time:

rna %>%
  count(infection, time) %>%
  arrange(time)

# A tibble: 3 × 3
  infection    time     n
  <chr>       <dbl> <int>
1 NonInfected     0 10318
2 InfluenzaA      4 11792
3 InfluenzaA      8 10318

or by counts:

rna %>%
  count(infection, time) %>%
  arrange(n)

# A tibble: 3 × 3
  infection    time     n
  <chr>       <dbl> <int>
1 InfluenzaA      8 10318
2 NonInfected     0 10318
3 InfluenzaA      4 11792

To sort in descending order, we need to add the desc() function:

rna %>%
  count(infection, time) %>%
  arrange(desc(n))

# A tibble: 3 × 3
  infection    time     n
  <chr>       <dbl> <int>
1 InfluenzaA      4 11792
2 InfluenzaA      8 10318
3 NonInfected     0 10318

Challenge

1. How many genes were analysed in each sample?
2. Use group_by() and summarize() to evaluate the sequencing depth (the sum of all counts) in each sample. Which sample has the highest sequencing depth?
3. Pick one sample and evaluate the number of genes by biotype
4. Identify genes associated with “abnormal DNA methylation” phenotype description, and calculate their mean expression (in log) at time 0, time 4 and time 8.

Solution

## 1.
rna %>%
  count(sample)

# A tibble: 22 × 2
   sample         n
   <chr>      <int>
 1 GSM2545336  1474
 2 GSM2545337  1474
 3 GSM2545338  1474
 4 GSM2545339  1474
 5 GSM2545340  1474
 6 GSM2545341  1474
 7 GSM2545342  1474
 8 GSM2545343  1474
 9 GSM2545344  1474
10 GSM2545345  1474
# … with 12 more rows

## 2.
rna %>%
  group_by(sample) %>%
  summarize(seq_depth = sum(expression)) %>%
  arrange(desc(seq_depth))

# A tibble: 22 × 2
   sample     seq_depth
   <chr>          <dbl>
 1 GSM2545350   3255566
 2 GSM2545352   3216163
 3 GSM2545343   3105652
 4 GSM2545336   3039671
 5 GSM2545380   3036098
 6 GSM2545353   2953249
 7 GSM2545348   2913678
 8 GSM2545362   2913517
 9 GSM2545351   2782464
10 GSM2545349   2758006
# … with 12 more rows

## 3.
rna %>%
  filter(sample == "GSM2545336") %>%
  count(gene_biotype) %>%
  arrange(desc(n))

# A tibble: 13 × 2
   gene_biotype                           n
   <chr>                              <int>
 1 protein_coding                      1321
 2 lncRNA                                69
 3 processed_pseudogene                  59
 4 miRNA                                  7
 5 snoRNA                                 5
 6 TEC                                    4
 7 polymorphic_pseudogene                 2
 8 unprocessed_pseudogene                 2
 9 IG_C_gene                              1
10 scaRNA                                 1
11 transcribed_processed_pseudogene       1
12 transcribed_unitary_pseudogene         1
13 transcribed_unprocessed_pseudogene     1

## 4.
rna %>%
  filter(phenotype_description == "abnormal DNA methylation") %>%
  group_by(gene, time) %>%
  summarize(mean_expression = mean(log(expression))) %>%
  arrange()

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

# A tibble: 6 × 3
# Groups:   gene [2]
  gene   time mean_expression
  <chr> <dbl>           <dbl>
1 Xist      0            6.95
2 Xist      4            6.34
3 Xist      8            7.13
4 Zdbf2     0            6.27
5 Zdbf2     4            6.27
6 Zdbf2     8            6.19

Reshaping data

In the rna tibble, the rows contain expression values (the unit) that are associated with a combination of 2 other variables: gene and sample.

All the other columns correspond to variables describing either the sample (organism, age, sex,…) or the gene (gene_biotype, ENTREZ_ID, product…). The variables that don’t change with genes or with samples will have the same value in all the rows.

rna %>%
  arrange(gene)

# A tibble: 32,428 × 19
   gene  sample    expression organism   age sex   infection strain  time tissue
   <chr> <chr>          <dbl> <chr>    <dbl> <chr> <chr>     <chr>  <dbl> <chr> 
 1 Aamp  GSM25453…       5621 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 2 Aamp  GSM25453…       4049 Mus mus…     8 Fema… NonInfec… C57BL…     0 Cereb…
 3 Aamp  GSM25453…       3797 Mus mus…     8 Fema… NonInfec… C57BL…     0 Cereb…
 4 Aamp  GSM25453…       4375 Mus mus…     8 Fema… Influenz… C57BL…     4 Cereb…
 5 Aamp  GSM25453…       4095 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
 6 Aamp  GSM25453…       3867 Mus mus…     8 Male  Influenz… C57BL…     8 Cereb…
 7 Aamp  GSM25453…       3578 Mus mus…     8 Fema… Influenz… C57BL…     8 Cereb…
 8 Aamp  GSM25453…       5097 Mus mus…     8 Male  NonInfec… C57BL…     0 Cereb…
 9 Aamp  GSM25453…       4202 Mus mus…     8 Fema… Influenz… C57BL…     4 Cereb…
10 Aamp  GSM25453…       4701 Mus mus…     8 Male  Influenz… C57BL…     4 Cereb…
# … with 32,418 more rows, and 9 more variables: mouse <dbl>, ENTREZID <dbl>,
#   product <chr>, ensembl_gene_id <chr>, external_synonym <chr>,
#   chromosome_name <chr>, gene_biotype <chr>, phenotype_description <chr>,
#   hsapiens_homolog_associated_gene_name <chr>

This structure is called a long-format, as one column contains all the values, and other column(s) list(s) the context of the value.

In certain cases, the long-format is not really “human-readable”, and another format, a wide-format is preferred, as a more compact way of representing the data. This is typically the case with gene expression values that scientists are used to look as matrices, were rows represent genes and columns represent samples.

In this format, it would become therefore straightforward to explore the relationship between the gene expression levels within, and between, the samples.

# A tibble: 1,474 × 23
   gene    GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
   <chr>        <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
 1 Asl           1170        361        400        586        626        988
 2 Apod         36194      10347       9173      10620      13021      29594
 3 Cyp2d22       4060       1616       1603       1901       2171       3349
 4 Klk6           287        629        641        578        448        195
 5 Fcrls           85        233        244        237        180         38
 6 Slc2a4         782        231        248        265        313        786
 7 Exd2          1619       2288       2235       2513       2366       1359
 8 Gjc2           288        595        568        551        310        146
 9 Plp1         43217     101241      96534      58354      53126      27173
10 Gnb4          1071       1791       1867       1430       1355        798
# … with 1,464 more rows, and 16 more variables: GSM2545342 <dbl>,
#   GSM2545343 <dbl>, GSM2545344 <dbl>, GSM2545345 <dbl>, GSM2545346 <dbl>,
#   GSM2545347 <dbl>, GSM2545348 <dbl>, GSM2545349 <dbl>, GSM2545350 <dbl>,
#   GSM2545351 <dbl>, GSM2545352 <dbl>, GSM2545353 <dbl>, GSM2545354 <dbl>,
#   GSM2545362 <dbl>, GSM2545363 <dbl>, GSM2545380 <dbl>

To convert the gene expression values from rna into a wide-format, we need to create a new table where the values of the sample column would become the names of column variables.

The key point here is that we are still following a tidy data structure, but we have reshaped the data according to the observations of interest: expression levels per gene instead of recording them per gene and per sample.

The opposite transformation would be to transform column names into values of a new variable.

We can do both these of transformations with two tidyr functions, pivot_longer() and pivot_wider() (see here for details).

Pivoting the data into a wider format

Let’s first select the 3 first columns of rna and use pivot_wider() to transform data in a wide-format.

rna_exp <- rna %>%
  select(gene, sample, expression)
rna_exp

# A tibble: 32,428 × 3
   gene    sample     expression
   <chr>   <chr>           <dbl>
 1 Asl     GSM2545336       1170
 2 Apod    GSM2545336      36194
 3 Cyp2d22 GSM2545336       4060
 4 Klk6    GSM2545336        287
 5 Fcrls   GSM2545336         85
 6 Slc2a4  GSM2545336        782
 7 Exd2    GSM2545336       1619
 8 Gjc2    GSM2545336        288
 9 Plp1    GSM2545336      43217
10 Gnb4    GSM2545336       1071
# … with 32,418 more rows

pivot_wider takes three main arguments:

1. the data to be transformed;
2. the names_from : the column whose values will become new column names;
3. the values_from: the column whose values will fill the new columns.

rna_wide <- rna_exp %>%
  pivot_wider(names_from = sample,
              values_from = expression)
rna_wide

# A tibble: 1,474 × 23
   gene    GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
   <chr>        <dbl>      <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
 1 Asl           1170        361        400        586        626        988
 2 Apod         36194      10347       9173      10620      13021      29594
 3 Cyp2d22       4060       1616       1603       1901       2171       3349
 4 Klk6           287        629        641        578        448        195
 5 Fcrls           85        233        244        237        180         38
 6 Slc2a4         782        231        248        265        313        786
 7 Exd2          1619       2288       2235       2513       2366       1359
 8 Gjc2           288        595        568        551        310        146
 9 Plp1         43217     101241      96534      58354      53126      27173
10 Gnb4          1071       1791       1867       1430       1355        798
# … with 1,464 more rows, and 16 more variables: GSM2545342 <dbl>,
#   GSM2545343 <dbl>, GSM2545344 <dbl>, GSM2545345 <dbl>, GSM2545346 <dbl>,
#   GSM2545347 <dbl>, GSM2545348 <dbl>, GSM2545349 <dbl>, GSM2545350 <dbl>,
#   GSM2545351 <dbl>, GSM2545352 <dbl>, GSM2545353 <dbl>, GSM2545354 <dbl>,
#   GSM2545362 <dbl>, GSM2545363 <dbl>, GSM2545380 <dbl>

Pivoting data into a longer format

In the opposite situation we are using the column names and turning them into a pair of new variables. One variable represents the column names as values, and the other variable contains the values previously associated with the column names.

pivot_longer() takes four main arguments:

1. the data to be transformed;
2. the names_to: the new column name we wish to create and populate with the current column names;
3. the values_to: the new column name we wish to create and populate with current values;
4. the names of the columns to be used to populate the names_to and values_to variables (or to drop).

To recreate rna_long from rna_long we would create a key called sample and value called expression and use all columns except gene for the key variable. Here we drop gene column with a minus sign.

Notice how the new variable names are to be quoted here.

rna_long <- rna_wide %>%
    pivot_longer(names_to = "sample",
                 values_to = "expression",
                 -gene)
rna_long

# A tibble: 32,428 × 3
   gene  sample     expression
   <chr> <chr>           <dbl>
 1 Asl   GSM2545336       1170
 2 Asl   GSM2545337        361
 3 Asl   GSM2545338        400
 4 Asl   GSM2545339        586
 5 Asl   GSM2545340        626
 6 Asl   GSM2545341        988
 7 Asl   GSM2545342        836
 8 Asl   GSM2545343        535
 9 Asl   GSM2545344        586
10 Asl   GSM2545345        597
# … with 32,418 more rows

We could also have used a specification for what columns to include. This can be useful if you have a large number of identifying columns, and it’s easier to specify what to gather than what to leave alone. Here the starts_with() function can help to retrieve sample names without having to list them all! Another possibility would be to use the : operator!

rna_wide %>%
    pivot_longer(names_to = "sample",
                 values_to = "expression",
                 cols = starts_with("GSM"))

# A tibble: 32,428 × 3
   gene  sample     expression
   <chr> <chr>           <dbl>
 1 Asl   GSM2545336       1170
 2 Asl   GSM2545337        361
 3 Asl   GSM2545338        400
 4 Asl   GSM2545339        586
 5 Asl   GSM2545340        626
 6 Asl   GSM2545341        988
 7 Asl   GSM2545342        836
 8 Asl   GSM2545343        535
 9 Asl   GSM2545344        586
10 Asl   GSM2545345        597
# … with 32,418 more rows

rna_wide %>%
    pivot_longer(names_to = "sample",
                 values_to = "expression",
                 GSM2545336:GSM2545380)

# A tibble: 32,428 × 3
   gene  sample     expression
   <chr> <chr>           <dbl>
 1 Asl   GSM2545336       1170
 2 Asl   GSM2545337        361
 3 Asl   GSM2545338        400
 4 Asl   GSM2545339        586
 5 Asl   GSM2545340        626
 6 Asl   GSM2545341        988
 7 Asl   GSM2545342        836
 8 Asl   GSM2545343        535
 9 Asl   GSM2545344        586
10 Asl   GSM2545345        597
# … with 32,418 more rows

Challenge

Subset genes located on X and Y chromosomes from the rna data frame and spread the data frame with sex as columns, chromosome_name as rows, and the mean expression of genes located in each chromosome as the values, as in the following tibble:

You will need to summarize before reshaping!

Let’s first calculate the mean expression level of X and Y linked genes from male and female samples…

 rna %>%
  filter(chromosome_name == "Y" | chromosome_name == "X") %>%
  group_by(sex, chromosome_name) %>%
  summarize(mean = mean(expression))

`summarise()` has grouped output by 'sex'. You can override using the `.groups`
argument.

# A tibble: 4 × 3
# Groups:   sex [2]
  sex    chromosome_name  mean
  <chr>  <chr>           <dbl>
1 Female X               3504.
2 Female Y                  3 
3 Male   X               2497.
4 Male   Y               2117.

And pivot the table to wide format

rna_1 <- rna %>%
  filter(chromosome_name == "Y" | chromosome_name == "X") %>%
  group_by(sex, chromosome_name) %>%
  summarize(mean = mean(expression)) %>%
  pivot_wider(names_from = sex,
              values_from = mean)

`summarise()` has grouped output by 'sex'. You can override using the `.groups`
argument.

rna_1

# A tibble: 2 × 3
  chromosome_name Female  Male
  <chr>            <dbl> <dbl>
1 X                3504. 2497.
2 Y                   3  2117.

Now take that data frame and transform it with pivot_longer() so each row is a unique chromosome_name by gender combination.

Solution

rna_1 %>%
  pivot_longer(names_to = "gender",
               values_to = "mean",
               - chromosome_name)

# A tibble: 4 × 3
  chromosome_name gender  mean
  <chr>           <chr>  <dbl>
1 X               Female 3504.
2 X               Male   2497.
3 Y               Female    3 
4 Y               Male   2117.

Homework Challenge

1. Use the rna dataset to create an expression table where each row represents the mean expression levels of genes and columns represent the different timepoints.
2. Use the previous table containing mean expression levels per timepoint and create a new column containing fold-changes between timepoint 8 and timepoint 0, and fold-changes between timepoint 8 and timepoint 4. Convert this table into a long-format table gathering the foldchanges calculated.

Solution

Let’s first calculate the mean expression by gene and by time

 rna %>%
   group_by(gene, time) %>%
 summarize(mean_exp = mean(expression))

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

# A tibble: 4,422 × 3
# Groups:   gene [1,474]
   gene     time mean_exp
   <chr>   <dbl>    <dbl>
 1 Aamp        0  4603.  
 2 Aamp        4  4870   
 3 Aamp        8  4763.  
 4 Abca12      0     5.29
 5 Abca12      4     4.25
 6 Abca12      8     4.14
 7 Abcc8       0  2576.  
 8 Abcc8       4  2609.  
 9 Abcc8       8  2292.  
10 Abhd14a     0   591.  
# … with 4,412 more rows

before using the pivot_wider() function

rna_time <- rna %>%
 group_by(gene, time) %>%
 summarize(mean_exp = mean(expression)) %>%
 pivot_wider(names_from = time,
          values_from = mean_exp)

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

rna_time

# A tibble: 1,474 × 4
# Groups:   gene [1,474]
   gene        `0`     `4`     `8`
   <chr>     <dbl>   <dbl>   <dbl>
 1 Aamp    4603.   4870    4763.  
 2 Abca12     5.29    4.25    4.14
 3 Abcc8   2576.   2609.   2292.  
 4 Abhd14a  591.    547.    432.  
 5 Abi2    4881.   4903.   4945.  
 6 Abi3bp  1175.   1061.    762.  
 7 Abl2    2170.   2078.   2131.  
 8 Acadl   2059.   2099    1995.  
 9 Acap3   3745    3446.   3431.  
10 Acbd4   1219.   1410.   1668.  
# … with 1,464 more rows

Notice that this generates a tibble with some column names starting by a number. If we wanted to select the column corresponding to the timepoints, we could not use the column names directly… What happens when we select the colum 4?

 rna %>%
  group_by(gene, time) %>%
  summarize(mean_exp = mean(expression)) %>%
  pivot_wider(names_from = time,
              values_from = mean_exp) %>%
   select(gene, 4)

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

# A tibble: 1,474 × 2
# Groups:   gene [1,474]
   gene        `8`
   <chr>     <dbl>
 1 Aamp    4763.  
 2 Abca12     4.14
 3 Abcc8   2292.  
 4 Abhd14a  432.  
 5 Abi2    4945.  
 6 Abi3bp   762.  
 7 Abl2    2131.  
 8 Acadl   1995.  
 9 Acap3   3431.  
10 Acbd4   1668.  
# … with 1,464 more rows

To select the timepoint 4, we would have to quote the column name, with backticks “`”

 rna %>%
  group_by(gene, time) %>%
  summarize(mean_exp = mean(expression)) %>%
  pivot_wider(names_from = time,
               values_from = mean_exp) %>%
  select(gene, `4`)

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

# A tibble: 1,474 × 2
# Groups:   gene [1,474]
   gene        `4`
   <chr>     <dbl>
 1 Aamp    4870   
 2 Abca12     4.25
 3 Abcc8   2609.  
 4 Abhd14a  547.  
 5 Abi2    4903.  
 6 Abi3bp  1061.  
 7 Abl2    2078.  
 8 Acadl   2099   
 9 Acap3   3446.  
10 Acbd4   1410.  
# … with 1,464 more rows

Another possibility would be to rename the column, choosing a name that doesn’t start by a number :

rna_time <- rna %>%
  group_by(gene, time) %>%
  summarize(mean_exp = mean(expression)) %>%
  pivot_wider(names_from = time,
               values_from = mean_exp) %>%
  rename("time0" = `0`, "time4" = `4`, "time8" = `8`)

`summarise()` has grouped output by 'gene'. You can override using the
`.groups` argument.

Calculate FoldChanges:

 rna_time %>%
  mutate(time_8_vs_0 = time8 / time0, time_8_vs_4 = time8 / time4)

# A tibble: 1,474 × 6
# Groups:   gene [1,474]
   gene      time0   time4   time8 time_8_vs_0 time_8_vs_4
   <chr>     <dbl>   <dbl>   <dbl>       <dbl>       <dbl>
 1 Aamp    4603.   4870    4763.         1.03        0.978
 2 Abca12     5.29    4.25    4.14       0.784       0.975
 3 Abcc8   2576.   2609.   2292.         0.889       0.878
 4 Abhd14a  591.    547.    432.         0.731       0.791
 5 Abi2    4881.   4903.   4945.         1.01        1.01 
 6 Abi3bp  1175.   1061.    762.         0.649       0.719
 7 Abl2    2170.   2078.   2131.         0.982       1.03 
 8 Acadl   2059.   2099    1995.         0.969       0.950
 9 Acap3   3745    3446.   3431.         0.916       0.996
10 Acbd4   1219.   1410.   1668.         1.37        1.18 
# … with 1,464 more rows

And use the pivot_longer() function:

 rna_time %>%
   mutate(time_8_vs_0 = time8 / time0, time_8_vs_4 = time8 / time4) %>%
  pivot_longer(names_to = "comparisons",
                values_to = "Fold_changes",
               time_8_vs_0:time_8_vs_4)

# A tibble: 2,948 × 6
# Groups:   gene [1,474]
   gene      time0   time4   time8 comparisons Fold_changes
   <chr>     <dbl>   <dbl>   <dbl> <chr>              <dbl>
 1 Aamp    4603.   4870    4763.   time_8_vs_0        1.03 
 2 Aamp    4603.   4870    4763.   time_8_vs_4        0.978
 3 Abca12     5.29    4.25    4.14 time_8_vs_0        0.784
 4 Abca12     5.29    4.25    4.14 time_8_vs_4        0.975
 5 Abcc8   2576.   2609.   2292.   time_8_vs_0        0.889
 6 Abcc8   2576.   2609.   2292.   time_8_vs_4        0.878
 7 Abhd14a  591.    547.    432.   time_8_vs_0        0.731
 8 Abhd14a  591.    547.    432.   time_8_vs_4        0.791
 9 Abi2    4881.   4903.   4945.   time_8_vs_0        1.01 
10 Abi2    4881.   4903.   4945.   time_8_vs_4        1.01 
# … with 2,938 more rows

Exporting data

Now that you have learned how to use dplyr to extract information from or summarize your raw data, you may want to export these new data sets to share them with your collaborators or for archival.

Similar to the read_csv() function used for reading CSV files into R, there is a write_csv() function that generates CSV files from data frames.

Before using write_csv(), we are going to create a new folder, data_output, in our working directory (if we haven’t already created it) that will store this generated dataset. We don’t want to write generated datasets in the same directory as our raw data. It’s good practice to keep them separate. The data folder should only contain the raw, unaltered data, and should be left alone to make sure we don’t delete or modify it. In contrast, our script will generate the contents of the data_output directory, so even if the files it contains are deleted, we can always re-generate them.

Let’s use write_csv() to save the rna_wide table that we have created previously.

write_csv(rna_wide, file = "data_output/rna_wide.csv")

Key Points

• Tabular data in R using the tidyverse meta-package

Data visualization

Overview

Teaching: XX min
Exercises: XX min

Questions

• Visualization in R

Objectives

• Produce scatter plots, boxplots, line plots, etc. using ggplot.

• Set universal plot settings.

• Describe what faceting is and apply faceting in ggplot.

• Modify the aesthetics of an existing ggplot plot (including axis labels and color).

• Build complex and customized plots from data in a data frame.

Data Visualization

We start by loading the required packages. ggplot2 is included in the tidyverse package.

library("tidyverse")

If not still in the workspace, load the data we saved in the previous lesson.

rna <- read.csv("course-data/data/GSE96870/rnaseq.csv")

The Data Visualization Cheat Sheet will cover the basics and more advanced features of ggplot2 and will help, in addition to serve as a reminder, getting an overview of the many data representations available in the package. The following video tutorials (part 1 and 2) by Thomas Lin Pedersen are also very instructive.

Plotting with ggplot2

ggplot2 is a plotting package that makes it simple to create complex plots from data in a data frame. It provides a more programmatic interface for specifying what variables to plot, how they are displayed, and general visual properties. The theoretical foundation that supports the ggplot2 is the Grammar of Graphics (@Wilkinson:2005). Using this approach, we only need minimal changes if the underlying data change or if we decide to change from a bar plot to a scatterplot. This helps in creating publication quality plots with minimal amounts of adjustments and tweaking.

There is a book about ggplot2 (@ggplot2book) that provides a good overview, but it is outdated. The 3rd edition is in preparation and will be freely available online. The ggplot2 webpage (https://ggplot2.tidyverse.org) provides ample documentation.

ggplot2 functions like data in the ‘long’ format, i.e., a column for every dimension, and a row for every observation. Well-structured data will save you lots of time when making figures with ggplot2.

ggplot graphics are built step by step by adding new elements. Adding layers in this fashion allows for extensive flexibility and customization of plots.

The idea behind the Grammar of Graphics it is that you can build every graph from the same 3 components: (1) a data set, (2) a coordinate system, and (3) geoms — i.e. visual marks that represent data points ¹.

To build a ggplot, we will use the following basic template that can be used for different types of plots:

ggplot(data = <DATA>, mapping = aes(<MAPPINGS>)) +  <GEOM_FUNCTION>()

• use the ggplot() function and bind the plot to a specific data frame using the data argument

ggplot(data = rna)

• define a mapping (using the aesthetic (aes) function), by selecting the variables to be plotted and specifying how to present them in the graph, e.g. as x/y positions or characteristics such as size, shape, color, etc.

ggplot(data = rna, mapping = aes(x = expression))

• add ‘geoms’ – geometries, or graphical representations of the data in the plot (points, lines, bars). ggplot2 offers many different geoms; we will use some common ones today, including:

  * `geom_point()` for scatter plots, dot plots, etc.
  * `geom_histogram()` for histograms
  * `geom_boxplot()` for, well, boxplots!
  * `geom_line()` for trend lines, time series, etc.
  

To add a geom(etry) to the plot use the + operator. Let’s use geom_histogram() first:

ggplot(data = rna, mapping = aes(x = expression)) +
  geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The + in the ggplot2 package is particularly useful because it allows you to modify existing ggplot objects. This means you can easily set up plot templates and conveniently explore different types of plots, so the above plot can also be generated with code like this:

# Assign plot to a variable
rna_plot <- ggplot(data = rna,
                   mapping = aes(x = expression))

# Draw the plot
rna_plot + geom_histogram()

Homework Challenge

You have probably noticed an automatic message that appears when drawing the histogram:

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Change the arguments bins or binwidth of geom_histogram() to change the number or width of the bins.

Solution

# change bins
ggplot(rna, aes(x = expression)) +
    geom_histogram(bins = 15)

# change binwidth
ggplot(rna, aes(x = expression)) +
    geom_histogram(binwidth = 2000)

We can observe here that the data are skewed to the right. We can apply log2 transformation to have a more symmetric distribution. Note that we add here a small constant value (+1) to avoid having -Inf values returned for expression values equal to 0.

rna <- rna %>%
  mutate(expression_log = log2(expression + 1))

If we now draw the histogram of the log2-transformed expressions, the distribution is indeed closer to a normal distribution.

ggplot(rna, aes(x = expression_log)) + geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

From now on we will work on the log-transformed expression values.

Homework Challenge

Another way to visualize this transformation is to consider the scale of the observations. For example, it may be worth changing the scale of the axis to better distribute the observations in the space of the plot. Changing the scale of the axes is done similarly to adding/modifying other components (i.e., by incrementally adding commands). Try making this modification:

• Represent the un-transformed expression on the log10 scale; see scale_x_log10(). Compare it with the previous graph. Why do you now have warning messages appearing?

Solution

ggplot(data = rna,mapping = aes(x = expression))+
  geom_histogram() +
  scale_x_log10()

Warning: Transformation introduced infinite values in continuous x-axis

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 507 rows containing non-finite values (stat_bin).

Notes

• Anything you put in the ggplot() function can be seen by any geom layers that you add (i.e., these are global plot settings). This includes the x- and y-axis mapping you set up in aes().
• You can also specify mappings for a given geom independently of the mappings defined globally in the ggplot() function.
• The + sign used to add new layers must be placed at the end of the line containing the previous layer. If, instead, the + sign is added at the beginning of the line containing the new layer, ggplot2 will not add the new layer and will return an error message.

# This is the correct syntax for adding layers
rna_plot +
  geom_histogram()

# This will not add the new layer and will return an error message
rna_plot
  + geom_histogram()

Building your plots iteratively

We will now draw a scatter plot with two continuous variables and the geom_point() function. This graph will represent the log2 fold changes of expression comparing time 8 versus time 0 and time 4 versus time 0. To this end, we first need to compute the means of the log-transformed expression values by gene and time then the log fold changes by subtracting the mean log expressions between time 8 and time 0 and between time 4 and time 0. Note that we also include here the gene biotype that we will use later on to represent the genes. We will save the fold changes in a new data frame called rna_fc.

rna_fc <- rna %>% select(gene, time,
                         gene_biotype, expression_log) %>%
  group_by(gene, time, gene_biotype) %>%
  summarize(mean_exp = mean(expression_log)) %>%
  pivot_wider(names_from = time,
              values_from = mean_exp) %>%
  mutate(time_8_vs_0 = `8` - `0`, time_4_vs_0 = `4` - `0`)

`summarise()` has grouped output by 'gene', 'time'. You can override using the
`.groups` argument.

We can then build a ggplot with the newly created dataset rna_fc. Building plots with ggplot2 is typically an iterative process. We start by defining the dataset we’ll use, lay out the axes, and choose a geom:

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
  geom_point()

Then, we start modifying this plot to extract more information from it. For instance, we can add transparency (alpha) to avoid overplotting:

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
  geom_point(alpha = 0.3)

We can also add colors for all the points:

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
  geom_point(alpha = 0.3, color = "blue")

Or to color each gene in the plot differently, you could use a vector as an input to the argument color. ggplot2 will provide a different color corresponding to different values in the vector. Here is an example where we color with gene_biotype:

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0)) +
  geom_point(alpha = 0.3, aes(color = gene_biotype))

We can also specify the colors directly inside the mapping provided in the ggplot() function. This will be seen by any geom layers and the mapping will be determined by the x- and y-axis set up in aes().

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
                                color = gene_biotype)) +
  geom_point(alpha = 0.3)

Finally, we could also add a diagonal line with the geom_abline() function:

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
                                color = gene_biotype)) +
  geom_point(alpha = 0.3) +
  geom_abline(intercept = 0)

Notice that we can change the geom layer from geom_point to geom_jitter and colors will be still determined by gene_biotype.

ggplot(data = rna_fc, mapping = aes(x = time_4_vs_0, y = time_8_vs_0,
                                color = gene_biotype)) +
  geom_jitter(alpha = 0.3) +
  geom_abline(intercept = 0)

Challenge

Use what you just learned to create a scatter plot of expression_log over sample from the rna dataset with the time showing in different colors. Is this a good way to show this type of data?

Solution

ggplot(data = rna, mapping = aes(y = expression_log, x = sample)) +
    geom_point(aes(color = time))

Boxplot

We can use boxplots to visualize the distribution of gene expressions within each sample:

ggplot(data = rna,
         mapping = aes(y = expression_log, x = sample)) +
  geom_boxplot()

By adding points to boxplot, we can have a better idea of the number of measurements and of their distribution:

ggplot(data = rna,
         mapping = aes(y = expression_log, x = sample)) +
  geom_jitter(alpha = 0.2, color = "tomato") +
  geom_boxplot(alpha = 0)