Home Next Lecture →

Lecture 01 — Introduction to R for Biostatistics

This lecture introduces the basic R workflow used in the course: working with R objects, reading data from files, visualizing data with ggplot2, and using a simple Bioconductor package for biological sequence analysis.

1. R Basics
2. Reading and Writing Data
3. Plotting with ggplot2
4. Monte Carlo Simulation
5. Introduction to Bioconductor

Learning objectives

By the end of this lecture, you should be able to:

1. R Basics

File: Rbasics.R

We begin with the basic syntax of R and the objects used throughout the rest of the course.

Running commands in RStudio

Run the current line or selected lines using:

Lines beginning with # are comments and are ignored by R.

For example:

# Store a value
x <- 10

# Use the stored value
x + 2

The assignment operator

<-

means store the value on the right in the object on the left.

Main topics in Rbasics.R

Basic mathematics

2 + 3
10 / 4
5^2
sqrt(25)

abs(-5)
log(10)
exp(1)
pi

Variables

x <- 10
y <- 4

x + y
x * y

Vectors

The function c() combines values into a vector.

x <- c(2, 4, 6, 8, 10)

x
x * 2
x^2

Descriptive statistics

mean(x)
median(x)
sd(x)
var(x)
min(x)
max(x)
summary(x)
quantile(x)

Selecting elements

R indexing begins at 1.

x <- c(5, 8, 12, 3, 15, 7)

x[1]
x[1:3]
x[c(2, 5)]
x[x > 7]

Data frames

A data frame is similar to a spreadsheet:

students <- data.frame(
  name  = c("Anita", "Rahul", "Meera", "Arjun", "Nisha"),
  age   = c(20, 21, 19, 22, 20),
  score = c(78, 85, 92, 67, 88)
)

Useful commands include:

students
str(students)
summary(students)
nrow(students)
ncol(students)
names(students)

Columns can be accessed using $:

students$score
mean(students$score)

Rows and columns can be selected using:

students[row, column]

For example:

students[1, ]
students[, 2]
students[students$score > 80, ]

Basic visualization

Rbasics.R introduces the base-R plotting system.

Scatter plot

plot(x, y)

Histogram

hist(heights)

A histogram shows the distribution of a numerical variable by grouping observations into intervals or bins.

Boxplot

boxplot(control, treated)

A boxplot gives a compact summary of a distribution, including:

Introductory statistical methods

The script also gives a first look at several statistical methods that will be discussed in more detail later.

Independent two-sample t-test

t.test(control, treated)

By default, R uses Welch’s two-sample t-test.

Paired t-test

t.test(before, after, paired = TRUE)

This is appropriate when the same subjects are measured twice.

Correlation

cor(height, weight)

Linear regression

model <- lm(weight ~ height)

summary(model)

Add the fitted line to a plot:

abline(model)

Make a prediction:

predict(
  model,
  newdata = data.frame(height = 172)
)

Random numbers and missing values

R can generate random data:

set.seed(123)

x <- rnorm(
  100,
  mean = 50,
  sd = 10
)

set.seed() makes the simulation reproducible.

Missing observations are represented by NA.

values <- c(10, 12, 15, NA, 18)

mean(values)
mean(values, na.rm = TRUE)
is.na(values)

2. Reading and Writing Data

File: ReadCSV.R

Most real analyses begin with data stored in a file rather than data entered directly into an R script.

The second part of the lecture therefore introduces file input/output and simple data manipulation.

Working directory

Check the current working directory:

getwd()

List the files in it:

list.files()

When possible, use relative paths rather than computer-specific absolute paths. This makes scripts easier to share through GitHub and run on another computer.

For example, if students.csv is in the current folder:

students <- read.csv("students.csv")

R can also open a file-selection window:

students <- read.csv(file.choose())

Examining imported data

After reading a dataset, first inspect it.

head(students)
tail(students)
str(students)
names(students)
dim(students)
summary(students)

These commands help answer basic questions such as:

Selecting data

Select rows using a logical condition:

students[students$score > 80, ]

Select columns:

students[, c("name", "score")]

Creating a new variable

A new column can be created directly.

students$result <- ifelse(
  students$score >= 75,
  "Pass",
  "Fail"
)

Sorting data

Sort by score:

students[order(students$score), ]

Sort from highest to lowest:

students[order(students$score, decreasing = TRUE), ]

Missing values

Check for missing observations:

is.na(students)

Count missing values in each column:

colSums(is.na(students))

Count all missing values:

sum(is.na(students))

Modifying data

Modify an existing row:

students[2, ] <- data.frame(
  name = "Rahul",
  age = 22,
  score = 90
)

Add a new row:

new_student <- data.frame(
  name = "Kiran",
  age = 21,
  score = 81
)

students <- rbind(students, new_student)

Writing a CSV file

write.csv(
  students,
  "students_modified.csv",
  row.names = FALSE
)

Using

row.names = FALSE

prevents R from writing an unwanted extra column containing row numbers.

Other file formats

Tab-separated files

data_tsv <- read.delim("students.tsv")

or:

data_tsv <- read.table(
  "students.tsv",
  header = TRUE,
  sep = "\t"
)

RDS

RDS is useful for storing an individual R object while preserving its R-specific structure.

saveRDS(students, "students.rds")

students2 <- readRDS("students.rds")

RData

Several R objects can be saved together:

save(
  x,
  y,
  students,
  file = "my_data.RData"
)

They can later be restored using:

load("my_data.RData")

3. Plotting with ggplot2

File: Plotggplot2.R

Base R provides functions such as plot(), hist(), and boxplot(). A widely used alternative is ggplot2, which provides a flexible layered system for constructing graphics.

Installing and loading ggplot2

A package only needs to be installed once:

install.packages("ggplot2")

It must be loaded in each new R session in which it is used:

library(ggplot2)

A useful distinction is therefore:

install.packages()   -> install once
library()            -> load for the current R session

Example dataset

students <- data.frame(
  name = c("Asha", "Ravi", "Meera", "Arun", "Neha"),
  hours_studied = c(2, 4, 5, 6, 8),
  score = c(55, 65, 72, 78, 90)
)

Base R version

plot(
  students$hours_studied,
  students$score,
  xlab = "Hours Studied",
  ylab = "Exam Score",
  main = "Study Time vs Exam Score"
)

The same plot with ggplot2

ggplot(
  students,
  aes(x = hours_studied, y = score)
) +
  geom_point(size = 3) +
  labs(
    title = "Study Time vs Exam Score",
    x = "Hours Studied",
    y = "Exam Score"
  ) +
  theme_minimal()

A useful way to read this code is:

data
  +
mapping of variables
  +
geometric objects
  +
labels
  +
theme

The main components are:

Adding a fitted line

A linear-model fit can be added using geom_smooth():

ggplot(
  students,
  aes(x = hours_studied, y = score)
) +
  geom_point(size = 3) +
  geom_smooth(
    method = "lm",
    se = FALSE
  ) +
  labs(
    title = "Study Time vs Exam Score",
    x = "Hours Studied",
    y = "Exam Score"
  ) +
  theme_minimal()

Here:

method = "lm"

requests a linear model, while

se = FALSE

suppresses the confidence band around the fitted line.

Common themes

theme_gray()
theme_minimal()
theme_classic()
theme_bw()

The important idea is that ggplot2 builds plots by adding layers with the + operator.

4. Monte Carlo Simulation

File: MonteCarlo.R

Monte Carlo methods use random sampling to estimate numerical quantities.
They are widely used in statistics, physics, biology, finance, and many other fields when an exact calculation is difficult or when we want to study the behavior of a random process.

In this introductory example, random points are used to estimate the value of \(\pi\).

Geometrical idea

Consider a square covering the region

\[ 0 \le x \le 1, \qquad 0 \le y \le 1. \]

Inside this square is a quarter-circle of radius 1.

A point \((x,y)\) lies inside the quarter-circle when

\[ x^2 + y^2 \le 1. \]

The area of the square is

\[ A_{\mathrm{square}} = 1, \]

while the area of the quarter-circle is

\[ A_{\mathrm{quarter\ circle}} = \frac{\pi r^2}{4} = \frac{\pi}{4}. \]

Therefore,

\[ \frac{ A_{\mathrm{quarter\ circle}} }{ A_{\mathrm{square}} } = \frac{\pi}{4}. \]

If points are generated uniformly inside the square, the fraction that fall inside the quarter-circle should approach \(\pi/4\) as the number of points becomes large.

Thus,

\[ \pi \approx 4 \times \frac{ \text{number of points inside the quarter-circle} }{ \text{total number of points} }. \]

Reproducible random numbers

The script begins with

set.seed(123)

Computer-generated random numbers are actually pseudorandom.
Using set.seed() fixes the starting point of the random-number generator so that the same script produces the same sequence of random values every time.

This is important for reproducible scientific analysis.

Generating random points

Choose the number of simulated points:

N <- 100000

Generate random \(x\) and \(y\) coordinates between 0 and 1:

x <- runif(N, 0, 1)
y <- runif(N, 0, 1)

The function

runif(N, min, max)

generates N random numbers from a uniform distribution between min and max.

Here, every location between 0 and 1 has equal probability of being sampled.

Testing whether points are inside the circle

The condition

inside <- x^2 + y^2 <= 1

creates a logical vector containing values such as

TRUE
FALSE
TRUE
TRUE
FALSE
...

TRUE means that the corresponding point lies inside the quarter-circle.

Because R treats

TRUE  -> 1
FALSE -> 0

when calculating a mean, the command

mean(inside)

directly gives the fraction of points inside the quarter-circle.

For example,

fraction_inside <- mean(inside)

estimates

\[ \frac{\pi}{4}. \]

Estimating pi

The estimate is therefore

pi_estimate <- 4 * fraction_inside

and can be compared with R’s built-in value:

pi_estimate
pi

With a sufficiently large value of N, the Monte Carlo estimate should be close to the true value of \(\pi\).

The result will not usually be exactly equal to \(\pi\) because the estimate is based on a finite random sample.

Complete numerical calculation

set.seed(123)

# Number of random points
N <- 100000

# Generate uniformly distributed coordinates
x <- runif(N, 0, 1)
y <- runif(N, 0, 1)

# TRUE if the point lies inside the quarter-circle
inside <- x^2 + y^2 <= 1

# Fraction of simulated points inside
fraction_inside <- mean(inside)

# Monte Carlo estimate of pi
pi_estimate <- 4 * fraction_inside

pi_estimate
pi

Visualizing the simulation

For visualization, a smaller number of points is convenient:

set.seed(123)

N <- 5000

x <- runif(N)
y <- runif(N)

inside <- x^2 + y^2 <= 1

Plot the simulated points:

plot(
  x, y,
  pch = 16,
  cex = 0.5,
  xlab = "x",
  ylab = "y",
  asp = 1,
  main = "Monte Carlo Estimation of Pi"
)

Important plotting arguments are:

The argument

asp = 1

is particularly important here. Without equal axis scaling, a geometrically circular curve may appear stretched or compressed on the screen.

Drawing the quarter-circle

The upper boundary of the quarter-circle follows from

\[ x^2 + y^2 = 1. \]

Solving for \(y\) gives

\[ y = \sqrt{1-x^2}. \]

It can be added to the existing plot using

curve(
  sqrt(1 - x^2),
  from = 0,
  to = 1,
  add = TRUE,
  lwd = 2
)

Here:

Final estimate from the plotted sample

The Monte Carlo estimate can again be calculated directly:

4 * mean(inside)

Since this visualization uses only 5000 points rather than 100000, its estimate may differ slightly more from the exact value of \(\pi\).

Effect of sample size

An important idea in Monte Carlo simulation is that estimates usually become more stable as the number of random samples increases.

For example, try:

N <- 100
N <- 1000
N <- 10000
N <- 100000

and compare the corresponding estimates of \(\pi\).

The exact sequence of estimates will depend on the random sample, but larger values of N generally give estimates closer to the true value.

This example therefore introduces several useful R ideas at once:

5. Introduction to Bioconductor

File: BioConductor.R

The final part of the lecture gives a short introduction to Bioconductor, an ecosystem of R packages designed for biological and genomic data analysis.

This example uses the Biostrings package to work with a DNA sequence.

Installing Bioconductor packages

Install BiocManager once:

install.packages("BiocManager")

Use it to install Bioconductor packages:

BiocManager::install("Biostrings")

Load the package:

library(Biostrings)

Creating a DNA sequence

dna <- DNAString(
  "ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"
)

dna

DNAString stores a DNA sequence as a biological sequence object rather than as ordinary text.

Sequence length

length(dna)

Nucleotide composition

Count each nucleotide:

letterFrequency(
  dna,
  letters = c("A", "C", "G", "T")
)

GC content

gc_count <- letterFrequency(
  dna,
  letters = c("G", "C")
)

gc_percent <- sum(gc_count) / length(dna) * 100

gc_percent

The GC percentage is

\[ \mathrm{GC\%} = \frac{G + C}{A + T + G + C} \times 100 \]

Reverse complement

reverseComplement(dna)

For double-stranded DNA, this gives the sequence of the complementary strand written in the standard 5' -> 3' direction.

Translation

Translate the DNA sequence into an amino-acid sequence:

protein <- translate(dna)

protein

This provides a simple example of how specialized R packages can represent and manipulate biological data directly.