dplyrHappy families are all alike;
every unhappy family is unhappy in its own way.Leo Tolstoy (first line of “Anna Karenina”)
Hadley Wickham points out the nice connection between the above quote and data formats 1. While data formatting might seem esoteric and boring, paying attention to it early on can pay off a great deal in the long term, which hopefully leads to a happy scientist. If one were to break data formatting into two phases, they might be these:
Working backward, you do the least work on data formatting if you choose to save your data initially in a format that is easiest to analyze with R. Since we will use the graphing package, ggplot2, by Hadley Wickham, we will use the data format, tidy data, that he advocates. Lest you worry that this is too Hadley specific, or even R specific, this kind of format is pretty popular now as a way to “layer” or “compose” plots that use different graphical elements to represent different variables in the data (e.g., plotly and seaborn in Python and AlgebraOfGraphics and TidierPlots.jl in Julia). To get a sense for what this format is and why it might be useful, consider some different ways to organize data on tuberculosis (TB) cases where you have three variables, country, year, and population, for each measurement of cases. First, each variable including cases could have its own column:
library(tidyverse)
table1# A tibble: 6 × 4
country year cases population
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583Here, each row represents one time you measured the TB cases. Alternatively, you could have each row represent a country with columns for different years. This means you need a table that measures the cases and one that measures the population:
table4a # cases# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 745 2666
2 Brazil 37737 80488
3 China 212258 213766table4b # population# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 19987071 20595360
2 Brazil 172006362 174504898
3 China 1272915272 1280428583To see why the this format might get cumbersome, you just have to suppose you want to add more variables: each new variable requires a new table. If you are interested in just one variable and how it changes across country and year, you are fine. But if you want to look across variables, say correlate population and caseload over time, then you have to slice two different tables. If you want to manipulate three variables, then you have to manipulate three different tables, etc. With the first format, we can include all the variables in a single table. The first format is the tidy format.
The tidy format has the following three rules2:
Visually, this looks like 
The primary benefit of tidy data is that every variable is a vector (column), which means that slicing data is just slicing columns. Though this may become complicated when the number of variables is large, there are some helper functions that will make slicing tidy data much easier. The slicing functions come from the package dplyr and the plotting package ggplot2 assumes tidy data.
To give a flavor of the power of tidy data and ggplot2, we’ll look at an example. The “GWAS Catalog” (https://www.ebi.ac.uk/gwas/home) keeps track of the SNPs associated with traits in genome-wide association studies (GWAS). We can download these data, take a reasonable subset to play with (cancer related traits), and plot some things. For example, here are SNPs associated with breast, colon, and prostate cancer:
# https://www.ebi.ac.uk/gwas/docs/file-downloads
# we'll load a subset of these data with traits related to cancer
# to clean things up a little, we filter out rows where the CHR_POS isn't a single number; this way, the x-axis label is nice automatically.
gwas_cancer = read_tsv("gwas-catalog-v1.0-associations-e115-r2026-01-19-cancer.tsv", col_types = cols(CHR_POS = col_number()))
filter(gwas_cancer, `DISEASE/TRAIT` == 'Prostate cancer' | `DISEASE/TRAIT` == "Colorectal cancer" | `DISEASE/TRAIT` == "Breast cancer") |>
ggplot(aes(x=CHR_POS)) +
geom_point(aes(y=PVALUE_MLOG, color=`DISEASE/TRAIT`))
Often, data will not be in the tidy format by default, so it will be necessary to format it. The first step is to figure out what the “variables” and “observations” are. Sometimes this may require carefully thinking about the experimental design used to create the data. Two common problems with data that are untidy are:
Typically, a dataset will only suffer one of the above problems. The first problem is dealt with using the pivot_longer function and the second with the pivot_wider function.
Recall that this table is tidy
table1# A tibble: 6 × 4
country year cases population
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583whereas these two are not
table4a# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 745 2666
2 Brazil 37737 80488
3 China 212258 213766table4b# A tibble: 3 × 3
country `1999` `2000`
<chr> <dbl> <dbl>
1 Afghanistan 19987071 20595360
2 Brazil 172006362 174504898
3 China 1272915272 1280428583The problem with table4a and table4b is that the variable year is spread across multiple columns. Thus, you need to gather it together into a single column, which makes the data frame longer. The columns 1999 and 2000 will be collected as values of the variable year, which we pass to the argument names_to. The first table contains the number of cases (which is our values_to argument):
# note that backticks `` here. we need them since the column names start with a number
# (i.e, this allows us to break normal R naming rulues!).
tidy4a = pivot_longer(table4a, cols = c(`1999`, `2000`), names_to = "year", values_to = "cases")
tidy4a# A tibble: 6 × 3
country year cases
<chr> <chr> <dbl>
1 Afghanistan 1999 745
2 Afghanistan 2000 2666
3 Brazil 1999 37737
4 Brazil 2000 80488
5 China 1999 212258
6 China 2000 213766and the second table contains the population size
tidy4b = pivot_longer(table4b, cols = c(`1999`, `2000`), names_to = "year", values_to = "population")
tidy4b# A tibble: 6 × 3
country year population
<chr> <chr> <dbl>
1 Afghanistan 1999 19987071
2 Afghanistan 2000 20595360
3 Brazil 1999 172006362
4 Brazil 2000 174504898
5 China 1999 1272915272
6 China 2000 1280428583To join these two results together, you use the function full_join
full_join(tidy4a, tidy4b)Joining with `by = join_by(country, year)`# A tibble: 6 × 4
country year cases population
<chr> <chr> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583which is the same as the first table
table1# A tibble: 6 × 4
country year cases population
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583These “joins” are related to joins you do with databases if you know SQL. We’ll talk more about this next week.
To get a little closer work at how this works conceptually, we’ll borrow an example from Hadley’s R for Data Science book3. Suppose we have the following blood pressure data from three individuals
bpdf = tribble(
~id, ~bp1, ~bp2,
"A", 100, 120,
"B", 140, 115,
"C", 120, 125
)where id refers to the individual and bp1 and bp2 are measurements of blood pressure at two different times. By the way, the tribble function is a handy way of creating a tibble data frame by hand.
We want to tidy table with three columns, id, which we already have, measurement, which will be either bp1 or bp2, and the value of the measurement. To do this we use the code
pivot_longer(bpdf, cols = c(bp1, bp2), names_to = "measurement", values_to = "value")# A tibble: 6 × 3
id measurement value
<chr> <chr> <dbl>
1 A bp1 100
2 A bp2 120
3 B bp1 140
4 B bp2 115
5 C bp1 120
6 C bp2 125To see better how this reshaping works, think about first about the id column. We are keeping this column and its values need to be repeated in new rows, one for each additional variable being pivoted, which is one additional variable here. 
The column names become a new column, measurement, which need to be repeated once for each original row of the dataset. 
The second new column captures the values, which are unwound row by row from the original table into the new table 
The following table is untidy
table2# A tibble: 12 × 4
country year type count
<chr> <dbl> <chr> <dbl>
1 Afghanistan 1999 cases 745
2 Afghanistan 1999 population 19987071
3 Afghanistan 2000 cases 2666
4 Afghanistan 2000 population 20595360
5 Brazil 1999 cases 37737
6 Brazil 1999 population 172006362
7 Brazil 2000 cases 80488
8 Brazil 2000 population 174504898
9 China 1999 cases 212258
10 China 1999 population 1272915272
11 China 2000 cases 213766
12 China 2000 population 1280428583because observations (i.e, every year) are spread across multiple rows. To tidy it, identify the column that names the variables, which in this case is type, and then the column with the values, which is count:
pivot_wider(table2, names_from = type, values_from = count)# A tibble: 6 × 4
country year cases population
<chr> <dbl> <dbl> <dbl>
1 Afghanistan 1999 745 19987071
2 Afghanistan 2000 2666 20595360
3 Brazil 1999 37737 172006362
4 Brazil 2000 80488 174504898
5 China 1999 212258 1272915272
6 China 2000 213766 1280428583Let’s look at the long alternate version of our blood pressure example.
long_bpdf = tribble(
~id, ~measurement, ~value,
"A", "bp1", 100,
"B", "bp1", 140,
"B", "bp2", 115,
"A", "bp2", 120,
"A", "bp3", 105
)
long_bpdf# A tibble: 5 × 3
id measurement value
<chr> <chr> <dbl>
1 A bp1 100
2 B bp1 140
3 B bp2 115
4 A bp2 120
5 A bp3 105To pivot this wider, we could do
pivot_wider(long_bpdf, names_from = measurement, values_from = value)# A tibble: 2 × 4
id bp1 bp2 bp3
<chr> <dbl> <dbl> <dbl>
1 A 100 120 105
2 B 140 115 NAThe pivot_wider function accomplishes this by looking for unique values in the measurement column since these will be the new columns. These values are
unique(long_bpdf$measurement)[1] "bp1" "bp2" "bp3"By default, the rows in the output are determined by all the variables that aren’t going into the new names or values. These are called the id_cols. Here, there is only one column, but in general there can be any number. The columns, bp1, bp2, and bp3, are now filled from the data in the original table and NA is inserted when there isn’t a corresponding row for that new columan row combination.
There are other ways your data might be untidy
separate or separate_wider_delim functions may be used.unite function is used.To read more about these, see the “vingnette” (or code tutorial) in R by typing vignette("pivot")
dplyrNow that you have made your data tidy, you probably want to slice and dice it. The dplyr package has handy functions for doing just this. These functions will making slicing MUCH EASIER than the base R way we’ve been doing it so far. Generally, dplyr is useful for doing the following five things
filter().arrange().select().mutate().summarize().Each of the functions above works in a similar way.
$.We’ll use the GWAS data you loaded above to demonstrate each of the five tasks.
filter()Let’s look at the GWAS cancer data with glimpse (which is from tidyverse).
glimpse(gwas_cancer)Rows: 25,272
Columns: 34
$ `DATE ADDED TO CATALOG` <date> 2019-03-27, 2019-03-27, 2019-03-27, 2020…
$ PUBMEDID <dbl> 28864454, 28864454, 28864454, 32117775, 3…
$ `FIRST AUTHOR` <chr> "Lacson JCA", "Lacson JCA", "Lacson JCA",…
$ DATE <date> 2017-09-01, 2017-09-01, 2017-09-01, 2020…
$ JOURNAL <chr> "Cancer Epidemiol Biomarkers Prev", "Canc…
$ LINK <chr> "www.ncbi.nlm.nih.gov/pubmed/28864454", "…
$ STUDY <chr> "Genome-Wide Testing of Exonic Variants a…
$ `DISEASE/TRAIT` <chr> "Breast cancer (estrogen-receptor positiv…
$ `INITIAL SAMPLE SIZE` <chr> "1,325 cases, 1,310 controls", "1,325 cas…
$ `REPLICATION SAMPLE SIZE` <chr> NA, NA, NA, "2,404 Han Chinese ancestry c…
$ REGION <chr> "10q26.13", "16q12.1", "10q26.13", "6p21.…
$ CHR_ID <chr> "10", "16", "10", "6", "8", "1", "1", "2"…
$ CHR_POS <dbl> 121577821, 52565276, 121577821, 34834333,…
$ `REPORTED GENE(S)` <chr> "FGFR2", "CASC16", "FGFR2", "UHRF1BP1", "…
$ MAPPED_GENE <chr> "FGFR2", "CASC16", "FGFR2", "BLTP3A", "JR…
$ UPSTREAM_GENE_ID <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ DOWNSTREAM_GENE_ID <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ SNP_GENE_IDS <chr> "ENSG00000066468", "ENSG00000249231", "EN…
$ UPSTREAM_GENE_DISTANCE <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ DOWNSTREAM_GENE_DISTANCE <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ `STRONGEST SNP-RISK ALLELE` <chr> "rs2981579-T", "rs4784227-T", "rs2981579-…
$ SNPS <chr> "rs2981579", "rs4784227", "rs2981579", "r…
$ MERGED <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ SNP_ID_CURRENT <dbl> 2981579, 4784227, 2981579, 35356162, 3736…
$ CONTEXT <chr> "intron_variant", "intron_variant", "intr…
$ INTERGENIC <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ `RISK ALLELE FREQUENCY` <chr> "NR", "NR", "NR", "0.0017", "0.0987", "NR…
$ `P-VALUE` <dbl> 4e-07, 7e-06, 5e-07, 4e-07, 5e-06, 3e-09,…
$ PVALUE_MLOG <dbl> 6.397940, 5.154902, 6.301030, 6.397940, 5…
$ `P-VALUE (TEXT)` <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ `OR or BETA` <dbl> 1.230, 1.230, 1.230, 4.332, 1.265, 5.933,…
$ `95% CI (TEXT)` <chr> NA, NA, NA, "[2.463–7.619]", NA, "NR z sc…
$ `PLATFORM [SNPS PASSING QC]` <chr> "Illumina [129771]", "Illumina [129771]",…
$ CNV <chr> "N", "N", "N", "N", "N", "N", "N", "N", "…Each observation (row) is a SNP from a GWAS study and we information on its location, the disease/trait associated, etc. Let’s filter just SNPs for colorectal cancer
filter(gwas_cancer, `DISEASE/TRAIT` == 'Colorectal cancer')# A tibble: 1,337 × 34
DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
<date> <dbl> <chr> <date> <chr> <chr> <chr>
1 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
2 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
3 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
4 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
5 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
6 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
7 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
8 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
9 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
10 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
# ℹ 1,327 more rows
# ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
# ℹ 27 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
# `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
# `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
# DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
# DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …We could apply another filter by just adding an argument, such as having a -log(p-value) (higher values of this mean more significant!) of greater than 10
filter(gwas_cancer, `DISEASE/TRAIT` == 'Colorectal cancer', PVALUE_MLOG > 10)# A tibble: 405 × 34
DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
<date> <dbl> <chr> <date> <chr> <chr> <chr>
1 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
2 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
3 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
4 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
5 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
6 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
7 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
8 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
9 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
10 2014-10-03 24448986 Zhang B 2014-01-21 Int J … www.… Geno…
# ℹ 395 more rows
# ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
# ℹ 27 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
# `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
# `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
# DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
# DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …Note that filter combines consecutive arguments by default using the “&”. You could equivalently give a single argument with the “&” to get the same slice
filter(gwas_cancer, `DISEASE/TRAIT` == 'Colorectal cancer' & PVALUE_MLOG > 10)# A tibble: 405 × 34
DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
<date> <dbl> <chr> <date> <chr> <chr> <chr>
1 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
2 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
3 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
4 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
5 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
6 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
7 2016-02-13 25990418 Al-Tassan NA 2015-05-20 Sci Rep www.… A ne…
8 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
9 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
10 2014-10-03 24448986 Zhang B 2014-01-21 Int J … www.… Geno…
# ℹ 395 more rows
# ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
# ℹ 27 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
# `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
# `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
# DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
# DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …arrange()The function arrange just sorts the rows based on the columns you specify. For example, to sort the colorectal cancer SNPs by DATE of the study,
colocancer = filter(gwas_cancer, `DISEASE/TRAIT` == 'Colorectal cancer')
arrange(colocancer, DATE)# A tibble: 1,337 × 34
DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
<date> <dbl> <chr> <date> <chr> <chr> <chr>
1 2008-06-16 17618283 Zanke BW 2007-07-08 Nat Ge… www.… Geno…
2 2008-06-16 17618284 Tomlinson I 2007-07-08 Nat Ge… www.… A ge…
3 2008-06-16 17934461 Broderick P 2007-10-14 Nat Ge… www.… A ge…
4 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
5 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
6 2008-06-16 18372901 Tenesa A 2008-03-30 Nat Ge… www.… Geno…
7 2008-06-16 18372905 Tomlinson IP 2008-03-30 Nat Ge… www.… A ge…
8 2008-06-16 18372905 Tomlinson IP 2008-03-30 Nat Ge… www.… A ge…
9 2008-06-16 18372905 Tomlinson IP 2008-03-30 Nat Ge… www.… A ge…
10 2008-06-16 18372905 Tomlinson IP 2008-03-30 Nat Ge… www.… A ge…
# ℹ 1,327 more rows
# ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
# ℹ 27 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
# `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
# `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
# DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
# DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …and use desc to make the sort a descending one,
arrange(colocancer, desc(DATE))# A tibble: 1,337 × 34
DATE ADDED TO CATALO…¹ PUBMEDID `FIRST AUTHOR` DATE JOURNAL LINK STUDY
<date> <dbl> <chr> <date> <chr> <chr> <chr>
1 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
2 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
3 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
4 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
5 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
6 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
7 2024-11-07 39392253 Bau DT 2024-10-11 Mol Ca… www.… Char…
8 2024-07-09 38872215 Tian J 2024-06-13 Genome… www.… Deve…
9 2024-07-09 38872215 Tian J 2024-06-13 Genome… www.… Deve…
10 2024-07-09 38872215 Tian J 2024-06-13 Genome… www.… Deve…
# ℹ 1,327 more rows
# ℹ abbreviated name: ¹`DATE ADDED TO CATALOG`
# ℹ 27 more variables: `DISEASE/TRAIT` <chr>, `INITIAL SAMPLE SIZE` <chr>,
# `REPLICATION SAMPLE SIZE` <chr>, REGION <chr>, CHR_ID <chr>, CHR_POS <dbl>,
# `REPORTED GENE(S)` <chr>, MAPPED_GENE <chr>, UPSTREAM_GENE_ID <chr>,
# DOWNSTREAM_GENE_ID <chr>, SNP_GENE_IDS <chr>, UPSTREAM_GENE_DISTANCE <dbl>,
# DOWNSTREAM_GENE_DISTANCE <dbl>, `STRONGEST SNP-RISK ALLELE` <chr>, …The select function simple selects specific columns (i.e., variables). To get only the disease, chromosome position, and -log(p-value) from the full data,
select(gwas_cancer, `DISEASE/TRAIT`, CHR_POS, PVALUE_MLOG)# A tibble: 25,272 × 3
`DISEASE/TRAIT` CHR_POS PVALUE_MLOG
<chr> <dbl> <dbl>
1 Breast cancer (estrogen-receptor positive) 121577821 6.40
2 Breast cancer (estrogen-receptor positive) 52565276 5.15
3 Breast cancer (progesterone-receptor positive) 121577821 6.30
4 Bladder cancer 34834333 6.40
5 Bladder cancer 142681389 5.30
6 Adult diffuse glioma (IDH mutation) 243651246 8.52
7 Adult diffuse glioma (IDH mutation) 243751309 9.70
8 Adult diffuse glioma (IDH mutation) 208188267 8.15
9 Adult diffuse glioma (IDH mutation) 241764203 9.52
10 Adult diffuse glioma (IDH mutation) 129503207 13.0
# ℹ 25,262 more rowsYou may want to add new variables that are functions of other variables. For example, you could create a column for the p-value from the PVALUE_MLOG (there is already one in the table, but we’ll make another)
sm_gwas_cancer = select(gwas_cancer, `DISEASE/TRAIT`, CHR_POS, PVALUE_MLOG, `P-VALUE`)
mutate(sm_gwas_cancer, new_p_value = 10^(-PVALUE_MLOG))# A tibble: 25,272 × 5
`DISEASE/TRAIT` CHR_POS PVALUE_MLOG `P-VALUE` new_p_value
<chr> <dbl> <dbl> <dbl> <dbl>
1 Breast cancer (estrogen-receptor p… 1.22e8 6.40 4 e- 7 4 e- 7
2 Breast cancer (estrogen-receptor p… 5.26e7 5.15 7 e- 6 7 e- 6
3 Breast cancer (progesterone-recept… 1.22e8 6.30 5 e- 7 5 e- 7
4 Bladder cancer 3.48e7 6.40 4 e- 7 4 e- 7
5 Bladder cancer 1.43e8 5.30 5 e- 6 5 e- 6
6 Adult diffuse glioma (IDH mutation) 2.44e8 8.52 3 e- 9 3 e- 9
7 Adult diffuse glioma (IDH mutation) 2.44e8 9.70 2 e-10 2.00e-10
8 Adult diffuse glioma (IDH mutation) 2.08e8 8.15 7 e- 9 7 e- 9
9 Adult diffuse glioma (IDH mutation) 2.42e8 9.52 3 e-10 3.00e-10
10 Adult diffuse glioma (IDH mutation) 1.30e8 13.0 9 e-14 9.00e-14
# ℹ 25,262 more rowssummarize()The function summarize (Hadley prefers the British spelling, summarise, which I refuse to acknowledge as an American 😤) collapses the data frame to a single row:
summarize(gwas_cancer, mean_mlog_p_value = mean(PVALUE_MLOG)) # A tibble: 1 × 1
mean_mlog_p_value
<dbl>
1 13.0This isn’t terribly useful until you use the group_by function to do the summarize action on data by “group”, which you specify according to values of variables. To see the average for each disease, group by DISEASE/TRAIT,
grouped_gwas_cancer = group_by(gwas_cancer, `DISEASE/TRAIT`)
summarize(grouped_gwas_cancer, mean_mlog_p_value = mean(PVALUE_MLOG))# A tibble: 929 × 2
`DISEASE/TRAIT` mean_mlog_p_value
<chr> <dbl>
1 (Z)-4-hydroxytamoxifen levels in tamoxifen-treated hormone… 8.05
2 (Z)-4-hydroxytamoxifen to tamoxifen ratio in tamoxifen-tre… 11.5
3 15-year breast cancer-specific survival (ER negative treat… 7.22
4 15-year breast cancer-specific survival (ER positive or PR… 8.52
5 15-year breast cancer-specific survival (ER positive treat… 6.70
6 15-year breast cancer-specific survival (ER positive, PR p… 8
7 15-year breast cancer-specific survival (grade 3 tumor) 7.35
8 15-year breast cancer-specific survival (treated with anth… 7.22
9 5-year breast cancer-specific survival (ER negative treate… 7.5
10 5-year breast cancer-specific survival (treated with tamox… 6.70
# ℹ 919 more rowsor add the journal the study was published in
grouped_gwas_cancer = group_by(gwas_cancer, `DISEASE/TRAIT`, JOURNAL)
meanp = summarize(grouped_gwas_cancer, mean_mlog_p_value = mean(PVALUE_MLOG))`summarise()` has grouped output by 'DISEASE/TRAIT'. You can override using the
`.groups` argument.meanp# A tibble: 1,298 × 3
# Groups: DISEASE/TRAIT [929]
`DISEASE/TRAIT` JOURNAL mean_mlog_p_value
<chr> <chr> <dbl>
1 (Z)-4-hydroxytamoxifen levels in tamoxifen-treated… Clin P… 8.05
2 (Z)-4-hydroxytamoxifen to tamoxifen ratio in tamox… Clin P… 11.5
3 15-year breast cancer-specific survival (ER negati… Breast… 7.22
4 15-year breast cancer-specific survival (ER positi… Breast… 8.52
5 15-year breast cancer-specific survival (ER positi… Breast… 6.70
6 15-year breast cancer-specific survival (ER positi… Breast… 8
7 15-year breast cancer-specific survival (grade 3 t… Breast… 7.35
8 15-year breast cancer-specific survival (treated w… Breast… 7.22
9 5-year breast cancer-specific survival (ER negativ… Breast… 7.5
10 5-year breast cancer-specific survival (treated wi… Breast… 6.70
# ℹ 1,288 more rowsFinally, let’s take the mean across all diseases for each journal and then sort by p-value to see the pattern with the journal a little better
arrange(summarize(group_by(gwas_cancer, JOURNAL), mean_mlog_p_value = mean(PVALUE_MLOG)), desc(mean_mlog_p_value))# A tibble: 154 × 2
JOURNAL mean_mlog_p_value
<chr> <dbl>
1 Circulation 154.
2 Nat Immunol 78.5
3 Kidney Int 58.5
4 Nat Metab 49.4
5 Science 35.5
6 Cancer Res Treat 26.8
7 Front Microbiol 22.3
8 Schizophr Bull 22.3
9 Commun Biol 20.7
10 Nat Genet 19.0
# ℹ 144 more rowsAnother nice function to use with summarize is the “counting” function n(), which can give the number of rows in each group. Since each row is a SNP, we can get the number of SNPs for each disease/trait this way:
arrange(summarize(group_by(gwas_cancer, `DISEASE/TRAIT`), n_SNPS = n()), desc(n_SNPS))# A tibble: 929 × 2
`DISEASE/TRAIT` n_SNPS
<chr> <int>
1 Prostate cancer 3206
2 Breast cancer 1860
3 Core binding factor acute myeloid leukemia 1744
4 Colorectal cancer 1337
5 Glioblastoma 1093
6 Oligodendroglioma 1089
7 Basal cell carcinoma 522
8 Lung cancer 431
9 Astrocytoma (high-grade) 425
10 Prostate cancer (late onset) 296
# ℹ 919 more rows|> operator to chain together functionsThe slicing and summarzing operations we’ve been doing can quickly get messy as we added more operations. Since each operation is applied to the results of the previous one, it ends up with a bunch of very nested parentheses, which can be hard to read. The pipe operator |> allows you to easily chain together operations so that its easier to read them from left to right. For example, the previous summarize and arrange operations could be put together like
gwas_cancer |>
group_by(`DISEASE/TRAIT`) |>
summarize(n_SNPS = n()) |>
arrange(desc(n_SNPS))# A tibble: 929 × 2
`DISEASE/TRAIT` n_SNPS
<chr> <int>
1 Prostate cancer 3206
2 Breast cancer 1860
3 Core binding factor acute myeloid leukemia 1744
4 Colorectal cancer 1337
5 Glioblastoma 1093
6 Oligodendroglioma 1089
7 Basal cell carcinoma 522
8 Lung cancer 431
9 Astrocytoma (high-grade) 425
10 Prostate cancer (late onset) 296
# ℹ 919 more rowsEffectively, the pipe takes what is on the left side of it and uses it as the first argument in the function on the right side. In the above example, gwas, our initial data table, is on the left and it gets used as the first argument in group_by. The group_by function returns a data.frame, which is then used as the argument to the function of the right of the next pipe, which is summarize. The rule of thumb then is that if you want to convert a command without a pipe to one where you use a pipe, take the first argument out of the function, put it on the left of the function separated by the pipe. Likewise, to get rid of the pipe, take the input on the left and put it in the first argument of the function on the right of the pipe and delete the pipe.
A very significant advantage of structuring your slicing / wrangling commands this way when you use RStudio is that RStudio will help you auto-complete column names. As you’ve probably seen already, typing these names can get cumbersome, so this is a big help. RStudio doesn’t do this when you put the data.frame object in the function directly (for some reason unknown to me).
For these problems, use the functions introduced above from the tidyverse packages like dplyr (e.g., filter, select, arrange, etc) as much as possible.
Using the GWAS data,
library(tidyverse)
gwas_cancer = read_tsv("gwas-catalog-v1.0-associations-e115-r2026-01-19-cancer.tsv", col_types = cols(CHR_POS = col_number()), na = "NR")produce a data table that shows the mean, maximum, and minimum of PVALUE_MLOG for each combination of disease and journal. Sort the table by the mean value of PVALUE_MLOG.
The RISK ALLELE FREQUENCY column in the GWAS data contains the frequency of the allele the increases the risk of the DISEASE/TRAIT in the population. For the three DISEASE/TRAITs “Prostate cancer”, “Lung cancer”, and “Breast cancer”, calculate the average value of the RISK ALLELE FREQUENCY for each disease. Which disease has the highest average RISK ALLELE FREQUENCY?
Recall the gene expression last that was briefly introduced last week:
library(readxl)
imprint = read_excel("babak-etal-2015_imprinted-mouse.xlsx", na = "NaN")Each row is a gene, each column is a tissue type, and each cell contains a gene expression measurement.
Make these data tidy (in one way) by collapsing the tissue columns and making the data table longer
You will need the pivot_longer function. The first trick here will be first to identify the “observations” and then the “names_to”, which is the variable that changes across each observation.
The second trick is specifying the columns across which you need to gather together into a single column. A hint for the second trick is that you can specify the columns you don’t want with the - operator or ! (NOT) operator.
The answer will be deceivingly simple! This is the elegance / frustration of R.
Using the tidy data from Problem 2, find the number of genes (across all tissue types) that have an expression value <= 10 and > 2. These genes are “paternally imprinted” in the Babak et al. dataset, which means they are only expressed from the maternal copy of the gene.
Hint: you will need to count each gene once even if it appears in multiple rows, which can be done with the distinct function.
We’ll use some the CDC data on COVID-19 deaths for this problem
library(socratadata)
covid19deaths = soc_read("https://data.cdc.gov/api/odata/v4/r8kw-7aab") |> as_tibble() |> filter(group == "By Week")These data are tidy. Make a wider version of this table by making columns for each state and using the covid_19_deaths column as the values. You can discard all the other columns except week_ending_date, covid_19_deaths, and the columns for each state. You should have only one row for each date.
Wickham, Hadley. 2014. J Stat Softw, 59:1–23. DOI: 10.18637/jss.v059.i10↩︎
Wickham, Hadley. 2014. J Stat Softw, 59:1–23. DOI: 10.18637/jss.v059.i10↩︎
Wickham, Hadley. 2023. R for Data Science (2e). O’Reilly Media. website↩︎