In this document, we will explore a small table of data about the genomes of various organisms. I also include a file, scrape_ensembl.R which shows how I scraped this data from the Ensembl website.
We start by reading in the TSV (tab-separated value) data using the readr package:
read_tsv returns a tibble object, which is similar to a data.frame, being tabular data with different types of variables as columns (character, double, integer, etc.). One difference, aside from the presentation when you show the object, is that the characters are not automatically converted to factor, as they are with read.delim or the like.
Let’s convert the first column to a factor:
orgs$type <-factor(orgs$type)
Let’s make some plots of the number of basepairs in the genomes of these organisms, split by the type of organism. There’s one plant with a really large genome. We could use the stripplot function, which is in the lattice package: this draws points for each level of a factor, similar to a boxplot, but is preferred when there are not many observations (e.g. < 20).
library(lattice)with(orgs, stripplot(bp ~ type))
with(orgs, stripplot(log10(bp) ~ type))
Compare to boxplot. Again, given the small number of samples, stripplot is preferred.
with(orgs, boxplot(bp ~ type))
If you haven’t seen with(...) before, it’s a way to temporarily use the variables (columns) of a data.frame as if they were variables in the environment, that is, without using the $. As we will see below, ggplot2 also offers such a convenient syntax for working with column names.
Most of the organisms have around 1 billion basepairs in their genome.
That large genome is wheat.
orgs[which.max(orgs$bp),]
# A tibble: 1 × 8
type name species bp coding noncoding pseudo tx
<fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 plant Wheat Triticum_aestivum 13427354022 104272 10156 NA 166589
Let’s remove it for now to make plots where we can see the range of the other organisms.
orgs <- orgs[orgs$name !="Wheat",]
I plot the number of genes that code for proteins over the total number of basepairs:
# right click to Esc from identify()# with(orgs, identify(bp, coding, labels=orgs$name))
Another way to make a plot with labels is using ggplot2. In ggplot2, you build a plot by adding layers with +: a call to ggplot() sets up the coordinate system and maps columns to visual properties via aes(), and then geom_* functions draw the data. Common geoms include geom_point(), geom_line(), geom_histogram(), geom_density(), geom_boxplot(), and geom_smooth().
# A tibble: 5 × 2
name bp
<chr> <dbl>
1 Orangutan 3446771396
2 Chimpanzee 3309577922
3 Elephant 3196760833
4 Human 3096649726
5 Gorilla 3040677044
Another way to get quick summaries of a dataset is skim from the skimr package:
library(skimr)skim(orgs)
Data summary
Name
orgs
Number of rows
37
Number of columns
8
_______________________
Column type frequency:
character
2
factor
1
numeric
5
________________________
Group variables
None
Variable type: character
skim_variable
n_missing
complete_rate
min
max
empty
n_unique
whitespace
name
0
1
3
19
0
37
0
species
0
1
8
24
0
37
0
Variable type: factor
skim_variable
n_missing
complete_rate
ordered
n_unique
top_counts
type
0
1
FALSE
6
mam: 14, fis: 8, inv: 7, bir: 4
Variable type: numeric
skim_variable
n_missing
complete_rate
mean
sd
p0
p25
p50
p75
p100
hist
bp
0
1.00
1.574812e+09
1.111087e+09
12157105
729664433.00
1233186341
2.521924e+09
3446771396
▆▇▁▆▅
coding
0
1.00
1.981016e+04
5.999450e+03
6692
17488.00
20033
2.143700e+04
39498
▂▅▇▁▁
noncoding
0
1.00
5.999380e+03
1.008250e+04
340
1042.00
2644
6.492000e+03
55401
▇▁▁▁▁
pseudo
3
0.92
1.276790e+03
3.056010e+03
21
69.25
382
8.727500e+02
14589
▇▁▁▁▁
tx
0
1.00
4.099754e+04
3.876149e+04
7126
21326.00
29196
3.811800e+04
199234
▇▁▁▁▁
For many operations, it’s useful to use the dplyr package to manipulate data.frames or tibbles.
Before I show dplyr operations, let me introduce the pipe operator |>, which sends the left-hand side to the first argument of the function on the right-hand side.
A simple example of the pipe:
x <-c(1:9,1000)x |>mean(trim=.1)
[1] 5.5
# equivalent tomean(x, trim=.1)
[1] 5.5
Using dplyr, let’s add a new column tallying up coding and non-coding genes:
dplyr is really useful for grouping and summarization tasks. The following line groups the organisms by type and summarizes the number for each type, the mean ratio of genes which are coding and the standard deviation of the ratio. One type has an SD of NA, why?
We can use summarization as a mid-point and then send those summaries to a plotting function. For example, points and line ranges:
# here we have to name the median something other than `bp`orgs |>group_by(type) |>summarize(basepairs=median(bp),min=min(bp),max=max(bp)) |>ggplot(aes(type, basepairs, ymin=min, ymax=max)) +geom_pointrange()
ggplot2 is great for grouping points together with summarizing graphical elements, such as lines through sets of points:
Finally, I like to show a useful plot constructed using ggplot2, just with simulated data. It is common that wants to show overlapping densities or histograms, which is not a simple task with base R graphics. It’s fairly simple with ggplot2.
We create some simulated data drawing three sets of samples from three distributions.
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
Appendix: |> vs %>% and the placeholder
R 4.1 introduced a native pipe |> that covers the most common use case: passing the left-hand side to the first argument of the right-hand side. For almost everything in this course, |> and %>% (from the magrittr package) are interchangeable.
The one place they differ is the placeholder, a way to pipe into an argument that is not the first, or that is nested inside another call.
Base R uses _ as a placeholder, but only as a direct named argument.
library(magrittr)orgs |>lm(coding ~ bp, data = _) |>coef() # base R
orgs %>%lm(coding ~ bp, data = .) %>% coef # magrittr can drop ()
(Intercept) bp
1.894803e+04 5.474499e-07
Magrittr’s . can appear anywhere on the right-hand side, including nested inside another call. Here we compute trim from the data itself:
x <-c(1:9, 1000)x |>mean(trim =0.1)
[1] 5.5
x %>%mean(trim =1/length(.)) # trim = 1/10 = 0.1
[1] 5.5
# x |> mean(trim = 1/length(_)) # syntax error
Magrittr also provides %<>%, a compound assignment pipe that applies a function and reassigns the result to the left-hand side. There is no base R equivalent; the plain <- form is always a safe replacement:
# these two are equivalentorgs$type %<>% factor # magrittr compound assignmentorgs$type <-factor(orgs$type) # base R
