6.5 – Discrete probability distributions

Binomial distribution

Discrete refers to particular outcomes — data that take on only specific separate values. Discrete data types include all of the categorical types we introduced earlier, including binary, ordinal, and nominal.

The binomial probability distribution is a discrete distribution for the number of successes, k, in a sequence of n independent trials, where the outcome of each trial can take on only one of two possible outcomes. For cases of 0 or 1, yes or no, “heads” or “tails,” male or female, pass or fail, we talk about the binomial distribution, because the outcomes are discrete and there can be only two possible (binary) outcomes.

The mathematical function of the binomial is written as

where the binomial coefficient is given by

and X refers to the number of ways to choose “success” from n observations.

Consider an example.

We have to define what we mean by success. For coin toss, this might be the number of heads.

The mean for the binomial this is given simply as

where X is “Heads” (the category of successes for our example), and p corresponds to the probability the selected event occurs, in this case, “Heads.”

The variance of the binomial distribution is given by

Here’s a density plot of two trials with success 2% with n (x) equal to 20 (Fig 1).

Figure 1. Plot generated with KMggplot2 Rcmdr plugin.

Here’s the R code

Create the trials, 1 through 20, then create an object to hold the number of trials

nSize=1:20
Size <- length(nSize); Size

R returns

[1] 20

Assign the probability value to an object

prob <- 0.02

Next, calculate the mean, mu, and the variance, var, for the binomial with prob = 0.02 and number of trials Size = 20

mu <- Size*prob
var <- Size*prob*(1-prob)

Print the mean and variance; lets assign them to an object then print the object

stats <- c(mu, var); stats

And R returns

[1] 0.400 0.392

And here’s a real-world example. Twinning in humans is rare. In Hawaiʻi in the 1990’s the rate of twin births (monozygotic and dizygotic) was about 20 for every 1000 births or 2%. “Success” here then is twin births.

Figure 2. Example of binomial-like distribution: reported twins born in Hawaiʻi.

Interestingly, rates of twins have since increased in Hawaiʻi (31 out of 1000 births) and in the United States overall (33 out of 1000 births) (Table 2, NCHS Data Brief No. 80, 2012). Data for Fig 2 were for year 2009.

Out of ten births, what is the probability of two twin births in Hawaiʻi?

You can solve this with your calculator (yikes!), or take advantage of online calculators (GraphPad QuickCalcs), or use R and Rcmdr.

In R, simply type at the prompt

dbinom(2,10,0.031)
[1] 0.03361446

Try in R Commander.

Rcmdr → Distributions → Discrete distributions → Binomial distribution → Binomial probabilities …

Figure 3. Rcmdr menu to get binomial probability.

Note I used p = 0.033 the rate for entire USA. Here’s the output

> .Table <- data.frame(Pr=dbinom(0:10, size=10, prob=0.033))
> rownames(.Table) <- 0:10
> .Table
Pr
0  7.149320e-01
1  2.439789e-01
2  3.746728e-02 ← Answer, 0.0375 or 3.75%
3  3.409639e-03
4  2.036263e-04
5  8.338782e-06
6  2.371422e-07
7  4.624430e-09
8  5.918028e-11
9  4.487991e-13
10 1.531579e-15

And here is the output for our example from Hawaiʻi (p = 0.031).

> .Table <- data.frame(Pr=dbinom(0:10, size=10, prob=0.031))
> rownames(.Table) <- 0:10
> .Table
Pr
0  7.298570e-01
1  2.334940e-01
2  3.361446e-02 ← Answer, 0.0336 or 3.36%
3  2.867694e-03
4  1.605494e-04
5  6.163507e-06
6  1.643178e-07
7  3.003893e-09
8  3.603741e-11
9  2.561999e-13
10 8.196283e-16

We use the binomial distribution as the foundation for the binomial test, ie, the test of an observed proportion against an expected population level proportion in a Bernoulli trial.

Hypergeometric distribution

The binomial distribution is used for cases of sampling with replacement from a population. When sampling without replacement is done, the hypergeometric distribution is used. It is the number of successes, k, in a sequence of n independent trials drawn from a fixed population.

This sampling scheme means that each draw is no longer independent — with each draw you decrease the remaining number of observations and thus change the proportion.

The P-value from the hypergeometric distribution is often used to study if genes from a list of Gene Ontology (eg, biological process functional terms like “cell cycle,” or “apoptosis”) are enriched in the sets of genes with expression differences across the treatment groups (cf, discussion in Cao and Zhang 2014). The R package gprofiler2, which links to a Cloud service at g:Profiler, the  can be used to do gene enrichment analyses with stunning visualizations (eg, Manhattan plot) suitable for publication.

The mathematical function of the hypergeometric is written as

where N is the population size, K is the number of successes in that population, and n and k are defined as above. Lets look apply this to the twinning problem.

In 2009, 2200 women gave birth in Hawaiʻi County, Hawaiʻi. Out of 10 births, what is the probability of 2 twin births in Hawaiʻi?

Assuming “risk” of twinning is the same rate as in rest of USA, then we have expected 72 successes in this population (0.033 * 2200).

Here’s the graph (Fig 4).

Figure 4. Plot of hypergeometric distribution twinning Hawaiʻi.

where the X axis values shows the number of events with successes (twin births). Taking the bin 2 (we wanted to know about probability of 2 out of ten), we can draw a line back to the Y-axis to get our probability — looks like about 5% roughly. Plot drawn with KMggplot2

To get the actual probability,

Rcmdr → Distributions → Discrete distributions → Hypergeometric distribution → Hypergeometric probabilities …

Figure 5. Rcmdr menu to get hypergeometric probability.

where m is the number of successes, n is the number of “failures,” and k is the number of trials.

> .Table
Pr
0 0.716453457
1 0.243438645
2 0.036688041 ← Answer, 0.0367 or 3.67%
3 0.003228871

The reference to white and black balls and urns is a device described by Bernoulli himself and has been used by others ever since to discuss probability problems (called the urn problem) and so I apply it here to be consistent. The urn contains a number of white (x) and a number of black (y) balls mixed together. One ball is drawn randomly from the urn — what color is it? The ball is then is either returned into the urn (replacement) or it is left out (without replacement) as in the hypergeometric problem, and the selection process is repeated.

Besides applications in gambling and balls-in-urns problems, this distribution is the basis for many tests of gene enrichment from microarray analyses. The hypergeometric forms the basis of the Fisher Exact test (see Chapter 9.5).

Discrete uniform distribution

For discrete cases of “1,” “2,” “3,” “4,” “5,” or “6,” on the single toss of fair dice, we can talk about the discrete uniform distribution because all possible outcomes are equally likely. If you are branded as a “card-counter” in Las Vegas, all you’ve done is reached an understanding of the uniform distribution of card suits!

One biological example would be the fate of a random primary oocyte in the human (mammal) female — three out of four will become polar bodies, eventually reabsorbed, whereas one in four will develop into a secondary oocyte (egg); the uniform distribution has to do with the counts of the products — each of the four primary oocytes has the same (apparently) chance (25%) of becoming the egg (Edwards and Beard 1997).

Analyses in biological research often use a uniform distribution as a baseline or null model (Gotelli and Ulrich 2011). For DNA sequencing, the uniform distribution provides a baseline for assessing coverage bias, eg, GC bias (Ross et al 2013). The uniform distribution underlies the Jukes–Cantor evolutionary model. The Jukes–Cantor model assumes equal substitution rates among all four nucleotides (A, T, C, G) and equal base frequencies (Felsenstein 2003).

The uniform distribution exists also for continuous data types, where every value within a specified interval has equal probability density. For example, DNA shearing — by endonucleases, proteases, or tansposases — during initial library prep for sequencing reactions, breakpoints may be modeled as uniformly distributed along fragment length. Non-uniform distributions reveal cleavage bias. Sequence-dependent cleavage bias is concluded when certain DNA fragments or sequences are generated or detected more (or less) frequently than would be expected by chance (Meyer and Liu 2014). This generates many false positives — that is, identifying a biological signal or molecule that isn’t actually there (Nix et a; 2008).

Poisson distribution

An extension from the binomial case is that, rather than following success or failure, you may have the following scenario. Consider a wind-dispersed seed released from a plant. If we markup the area around the plant in grids, we could then count the number of seeds within each grid. Most grids will have no seeds, some grids will have one seed, a few grids may have two seeds, etc. Multiple seeds in grids is a rare event. The graph might look like

Figure 6. Example, poisson-like graph: the number of wind-dispersed seeds within each grid.

In molecular biology, use of Poisson distribution is common. In digital PCR for example, a DNA sample is diluted and divided into thousands of small reaction volumes or partitions. Thus, each partition contains either zero molecules, one molecule, or more than one molecule, where the loading of molecules into partitions is assumed to be independent, random, and rare per partition. In sequencing, Poisson models are used to describe random sampling of molecules and reads, especially when coverage is generated by non-deterministic processes like shotgun sequencing.

The Poisson has interesting properties, one being that the expected mean, λ, is equal to the variance. Lambda is, with respect to the Poisson distribution, the average number of times an event happens over a specific period or area.

An equation is

where μ is the mean (or could substitute with variance!), e is the natural logarithm, and X is number of successes you are interested in. For example, if μ = 1, what is the probability of observing a grid with five seeds? Simple enough to do this by hand, but let’s use Rcmdr instead. Here’s the graph (Fig 7) from Rcmdr (KMggplot2 plugin).

Figure 7. The probability of observing a grid with five seeds, poisson μ = 1 (ggplot2).

and for the actual probability we have from R

dpois(5, lambda = 1)
[1] 0.003065662

Rcmdr → Distributions → Discrete distributions → Poisson distribution → Poisson probabilities … (Fig 8)

The only thing to enter is the mean (some call μ lambda with symbol λ).

Figure 8. Rcmdr menu, poisson probability.

Here’s the output from R. For intervals 0, 1, 2, 3, …, 6 (Rcmdr just enters this range for you!

> .Table <- data.fram(Pr=dpois(0:6, lambda = 1))
> rownames(.Table) <- 0:6
> .Table
Pr
0 0.3678794412
1 0.3678794412
2 0.1839397206
3 0.0613132402
4 0.0153283100
5 0.0030656620 ← Answer, 0.0307 or 3.07%
6 0.0005109437

Next — Continuous distributions

And finally, for ratio (continuous) scale data, which can take on any value, we can express the chance that probability of a given point as a continuous function, with the normal distribution being one of the most important examples (there are others, like the F-distribution). Many statistical procedures assume that the data we use can be viewed as having come from a “normally distributed population.” See Chapter 6.6.

Questions

1.2. Quarterback sacks by game for the NFL team Seahawks, years 2011 through 2022 are summarized below (data extracted from https://www.pro-football-reference.com/ ).

a) Assuming a Poisson distribution, what are the mean (lambda) and variance?

b) The table covers a total of 112 games. How many sacks (events) were observed?

c) What is the probability of the Seahawks getting zero sacks in a game (in 2022, a season was 17 games; prior years a season was 16 games)?

Quiz Chapter 6.5

Discrete probability

Start Quiz

Question

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 6 contents

• Introduction
• Some preliminaries
• Ratios and proportions
• Combinations and permutations
• Types of probability
• Discrete probability distributions
• Continuous distributions
• Normal distribution and the normal deviate
• Moments
• Chi-square distribution
• t distribution
• F distribution
• References and suggested readings

4.5 – Scatter plots

Introduction

Scatter plots, also called scatter diagrams, scatterplots, or XY plots, display associations between two quantitative, ratio-scaled variables. Each point in the graph is identified by two values: its X value and its Y value. The horizontal axis is used to display the dispersion of the X variable, while the vertical axis displays the dispersion of the Y variable.

The graphs we just looked at with Tufte’s examples Anscombe’s quartet data were scatter plots (Chapter 4 – How to report statistics).

Here’s another example of a scatter plot, data from Francis Galton, data contained in the R package HistData.

Figure 1. Scatterplot of mid-parent (vertical axis) and their adult children’s (horizontal axis) height, in inches. data from Galton’s 1885 paper, “Regression towards mediocrity in hereditary stature.” The red line is the linear regression fitted line, or “trend” line, which is interpreted in this case as the heritability of height.

The commands I used to make the plot in Figure 1 were

library(HistData)
data(GaltonFamilies, package="HistData")
attach(GaltonFamilies)

plot(childHeight~midparentHeight, xlab="Mid-parent height", ylab="Mean adult children's height", xlim=c(55, 80), ylim=c(55,80), cex=0.8, cex.axis=1.2, cex.lab=1.3, pch=c(19), data=GaltonFamilies)
abline(lm(childHeight~midparentHeight), col="red", lwd=2)

I forced the plot function to use the same range of values, set by providing values for xlim and ylim; the default values of the plot command picks a range of data that fits each variable independently. Thus, the default X axis values ranged from 64 to 76 and the Y variable values ranged from 55 to 80. This has the effect of shifting the data, reducing the amount of white space, which a naïve reading of Tufte would suggest is a good idea, but at the expense of allowing the reader to see what would be the main point of the graph: that the children are, on average, shorter than the parents, mean height = 67 vs. 69 inches, respectively. Therefore, Galton’s title begins with the word “regression,” as in the definition of regression as a “return to a former … state” (Oxford Dictionary).

For completeness, cex sets the size of the points (default = 1), and therefore cex.axis and cex.lab apply size changes to the axes and labels, respectively; pch refers to the graph elements or plotting characters, further discussed below (see Fig 8); lm() is a call to the linear model function; col refers to color.

Figure 2 shows the same plot, but without attention to the axis scales, and, more in keeping with Tufte’s principle of maximize data, minimize white space.

Figure 2. Same plot as Figure 1, but with default settings for axis scales.

Take a moment to compare the graphs in Figure 1 and 2. Setting the scales equal allows you to see that the mid-parent heights were less variable, between 65 and 75 inches, than the mean children height, which ranged from 55 to 80 inches.

And another example, Figure 3. This plot is from the ggplot2() function and was generated from within R Commander’s KMggplot2 plug-in.

Figure 3. Finishing times in minutes of 1278 runners by age and gender at the 2013 Jamba Juice Banana 5K in Honolulu, Hawaiʻi. Loess smoothing functions by groups of female (red) and male (blue) runners are plotted along with 95% confidence intervals.

Figure 3 is a busy plot. Because there were so many data points, it is challenging to view any discernible pattern, unlike Figure 1 and 2 plots, which featured less data. Use of the loess smoothing function, a transformation of the data to reduce data “noise” to reveal a continuous function, helps reveal patterns in the data:

1. across most ages, men completed 5K faster than did females and
2. there was an inverse, nonlinear association between runner’s age and time to complete the 5K race.

Take a look at the X-axis. Some runners ages were reported as less than 5 years old (trace the points down to the axis to confirm), and yet many of these youngsters were completing the 5K race in less than 30 minutes. That’s 6-minute mile pace! What might be some explanations for how pre-schoolers could be running so fast?

Design criteria

As in all plotting, maximize information to background. Keep white space minimal and avoid distorting relationships. Some things to consider:

1. keep axes same length
2. do not connect the dots UNLESS you have a continuous function
3. do not draw a trend line UNLESS you are implying causation

Scatter plots in R

We have many options in R to generate scatter plots. We have already demonstrated use of plot() to make scatter plots. Here we introduce how to generate the plot in R Commander.

Rcmdr: Graphs → Scatterplot…

Rcmdr uses the scatterplot function from the car package. In recent versions of R Commander the available options for the scatterplot command are divided into two menu tabs, Data and Options, shown in Figure 4 and Figure 5.

Figure 4. First menu popup in R Commander Scatterplot command, Rcmdr ver. 2.2-3.

Select X and Y variables, choose Plot by groups if multiple grounds are included, eg, male, female, then click Options tab to complete.

Figure 5. Second menu popup in R Commander scatterplot command., Rcmdr ver. 2.2-3

Set graph options including axes labels and size of the points.

Here is another example using the default settings in scatterplot() function in the car package, now the default scatter plot command via R Commander (Fig. 4), along with the same graph, but modified to improve the look and usefulness of the graph (Fig. 6). The data set was Puromycin in the package datasets.

Figure 6. Default scatterplot, package car, from R Commander, version 2.2-4.

Grid lines in graphs should be avoided unless you intend to draw attention to values of particular data points. I prefer to position the figure legend within the frame of the graph, eg, the open are at the bottom right of the graph. Modified graph shown in Figure 7.

Figure 7. Modified scatterplot, same data from Figure 6

R commands used to make the scatter plot in Figure 7 were

scatterplot(rate~conc|state, col=c("blue", "red"), cex=1.5, pch=c(19,19), 
bty="n", reg=FALSE, grid=FALSE, legend.coords="bottomright")

A comment about graph elements in R

In some ways R is too rich in options for making graphs. There are the plot functions in the base package, there’s lattice and ggplot2 which provide many options for graphics, and more. The advice is to start slowly, for example taking advantage of the Figure 8 displays R’s plotting characters and the number you would invoke to retrieve that plotting character.

Figure 8. R plotting characters pch = 1 – 25 along with examples of color.

colors()

  [1] "white"                "aliceblue"            "antiquewhite" 

[655] "yellow3"              "yellow4"              "yellowgreen"  

There’s a lot more to R plotting. For example, you are not limited to just 25 possible characters. R can print any of the ASCII characters 32:127 or from the extended ASCII code 128:255. See Wikipedia to see the listing of ASCII characters.

In addition to pch, you can change the size of the plotting character with cex, color of the character with col for lines or outlines, and bg to change background color. lwd is used to set the character’s border thickness.

Here’s the R code used to generate the graph in Figure 8. Remember, any line beginning with # is a comment line, not an R command.

#create a vector with 26 numbers, from 0 to 25
stuff <-c(0:25)
plot(stuff, pch=c(32:58), cex = 2.5, col = c(1:26), 'xlab' = "pch number", 'ylab' = "Plotting character number")

Is it “scatter plot” or “scatterplot”?

Spelling matters, of course, and yet there are many words for which the correct spelling seems to be like “beauty,” it is “in the eye of the beholder” (Molly Bawn, 1878, by Margaret Hungerford). Scatter plot is one of these — is it one word or two, or is it something else entirely?

Scatter plot is one of these terms: you’ll find it spelled as “scatterplot” or as “scatter plot,” in the dictionary (eg, Oxford English dictionary), with no guidance to choose between them. And I’m not just talking about the differences between British and American English spelling traditions. 

Thus, in these situations as an author, you can turn to which of the spellings is in common use. I first looked at some of the statistics books on my shelves. I selected 14 (bio)statistics textbooks and checked the index and if present, chapters on graphics for term usage.

Table 1. Frequency of use of different terms for scatter plot in 14 (bio)statistics books currently on Mike’s shelves.

Not much help, basically, it is a tie between “scatter plot” and “scatterplot.”

Next, I searched six journals for the interval 1990 – 2016 for use of these terms. Results presented in Table 6 along with journal impact factor for 2014 and number of issues (Table 2).

Table 2. Impact factor and number of issues 1990 – 2016 for six science journals.

My methods? I used the journal’s online search functions for the various usage for scatter plot and the results are shown in Figure 9.

Figure 9. Usage of terms for X Y plots in research articles normalized to number of issues in six journals between 1990 and 2016.

The journals have different numbers of articles; I partially corrected for this by calculating the ratio number of articles with one of the terms divided by the number of issues for the interval 1990 – 2016. It would have been better to count all of the articles, but even I found that to be an excessive effort given the point I’m trying to make here.

Not much help there, although we can see a trend favoring “scatter plot” over any of the other options.

And finally, to completely work over the issue I present results from use of Google’s Ngram Viewer. Ngram Viewer allows you to search words in all of the texts that Google’s folks have scanned into digital form. I searched on the terms in texts between 1950 and 2015, and results are displayed in Figure 10 and Figure 11.

Figure 10. Results from Ngram Viewer for American English, “scatterplot” (blue), “scatter plot” (red), “scatter diagram” (green), “scattergram” (orange), and “XY plot” (purple).

And the same plot, but this time for British sources

Figure 11. Results from Ngram Viewer for British English. See Figure 10 for key.

Conclusion? It looks like “scatterplot” (blue line) is the preferred usage, but it is close. Except for “scattergram” and “XY plot,” which, apparently, are rarely used. After all of this, it looks like you’re free to make your choice between “scatterplot” or “scatter plot.” I will continue to use “scatter plot.”

Bland-Altman plot

Also known as Tukey mean-difference plot, the Bland-Altman plot is used to describe agreement between two ratio scale variables (Bland and Altman 1986, Giavarina 2015), for example agreement between two different methods used to measure the same samples.

Consider use of imageJ by two different observers to record number of pixels of a unit measure (1 cm) on a series of digital images — there’s subjectivity in drawing the lines (where to start, where to end) — do they agree?  Data set below. blandr package, blandr.draw function.

Figure 12. Bland-Altman plot of 1 cm unit measure in pixel number by imageJ from digital images by two independent observers. Purple central region is 95% CI. Lower and upper dashed horizontal lines represent bounds of “acceptable” agreement.

blandr.draw(Obs1, Obs2)

The plot makes it easy to identify questionable points, for example, the one point in upper right quadrant looks suspect.

Visualizing large data sets

Since the early 2000s, size of data sets in biology have increased dramatically, leaving many statistical graphics solutions ill-equipped to tell stories. The data scientists involved with the genomics and health care communities have introduced a number of useful visualizations. We introduce just two, the volcano plot (Fig 13) and the Manhattan-style plot (Fig 14).

Volcano plot

The volcano plot (Fig 13) is a common way to show events that differ between two groups of subjects (eg, p-values), and is common often provides to help visualize gene expression studies of an exposed group vs a control group (eg, fold changes). For the volcano plot, fold changes are reflected on horizontal axis, indicating how much the gene expression level has increased or decreased. The vertical axis typically represents the negative logarithm of the p-value (-logP), which indicates the statistical significance of the change.

Below is R code — with assistance from Gemini AI (Fall 2026) — for a set of five replicate simulated gene expression fold change data sets on stress test (NaCl, H₂O₂) on yeast grown on YPDA plates. Gene expression profiles were simulated for a set of  100 genes (oxidative pathway, DNA repair pathway); as noted, such studies typically include many more genes. 

Figure 13. Volcano plot, gene expression fold change on simulated stress test (NaCl, H₂O₂) yeast data sets, 100 genes (oxidative pathway, DNA repair pathway).

Inspecting the plots (Fig 13), we see a symmetric distribution: The overall distribution of genes is cleanly centered around a Log2 Fold Change (LFC) of zero. This mirrors a realistic transcriptomic profile where the vast majority of the background genome (shown in grey) remains stable, displaying only technical and biological noise. There is a distinct, asymmetrical shift toward positive fold changes among specific gene pathways (targeted activation). Genes belonging to the Oxidative Stress pathway (red) and DNA Repair pathway (blue) consistently surpass both the statistical threshold () and the biological threshold (). Additionally, we see evidence of synergistic effects: The simultaneous application of both stressors (“Both”) drives these pathways to their highest expression levels, pushing them into the top-right quadrant of the plots as highly significant, upregulated features.

R code — with assistance from Gemini AI (Fall 2026) — for yeast gene expression and volcano plots; graphs by base package only.

# ==============================================================================
# 1. SIMULATE 5 REASONABLE DATASETS WITH BIOLOGICAL & TECHNICAL NOISE
# ==============================================================================
set.seed(123) # New seed for realistic distribution

n_genes <- 100
n_datasets <- 5
groups <- c("Control", "Salt", "H2O2", "Both")
n_groups <- length(groups)

# Real-world experiment: each dataset has 3 biological replicates per group
n_reps_per_group <- 3

# Define gene pathways (15 oxidative stress, 15 DNA repair, 70 background/unchanged)
gene_ids <- paste0("Gene_", 1:n_genes)
pathways <- c(rep("Oxidative_Stress", 15), rep("DNA_Repair", 15), rep("Background", 70))

# Initialize a list to hold the 5 realistic datasets
datasets <- list()

for (d in 1:n_datasets) {

# 1. Establish a realistic base expression distribution across 100 genes
# Real data has high-expressing and low-expressing genes
base_expression <- rnorm(n_genes, mean = 7, sd = 2)

# Create a blank sample matrix: 100 genes x 12 samples (4 groups x 3 replicates)
sample_names <- paste0(rep(groups, each = n_reps_per_group), "_R", 1:n_reps_per_group)
expr_matrix <- matrix(0, nrow = n_genes, ncol = length(sample_names))
rownames(expr_matrix) <- gene_ids
colnames(expr_matrix) <- sample_names

# 2. Add realistic baseline biological variation to every sample
for (col_idx in 1:ncol(expr_matrix)) {
expr_matrix[, col_idx] <- base_expression + rnorm(n_genes, mean = 0, sd = 0.4)
}

# 3. Inject subtle, realistic treatment effects (not over-exaggerated)
for (i in 1:n_genes) {
# Every gene gets a tiny bit of random directional drift across groups (real-world noise)
drift <- rnorm(n_groups, mean = 0, sd = 0.15)

if (pathways[i] == "Oxidative_Stress") {
# Subtle upregulation in H2O2 and Both (mean fold change ~ 1.0 to 1.3)
drift[3] <- drift[3] + rnorm(1, mean = 1.1, sd = 0.2) # H2O2
drift[4] <- drift[4] + rnorm(1, mean = 1.3, sd = 0.2) # Both
} else if (pathways[i] == "DNA_Repair") {
# Subtle upregulation in Salt and Both (mean fold change ~ 0.8 to 1.2)
drift[2] <- drift[2] + rnorm(1, mean = 0.9, sd = 0.2) # Salt
drift[4] <- drift[4] + rnorm(1, mean = 1.2, sd = 0.2) # Both
}

# Apply the group drifts to the respective sample columns
for (g_idx in 1:n_groups) {
start_col <- (g_idx - 1) * n_reps_per_group + 1
end_col <- g_idx * n_reps_per_group
expr_matrix[i, start_col:end_col] <- expr_matrix[i, start_col:end_col] + drift[g_idx]
}
}

datasets[[d]] <- expr_matrix
}

# ==============================================================================
# 2. STATISTICAL ANALYSIS (Calculated directly from replicates)
# ==============================================================================
# We will compare the double stressor 'Both' vs 'Control' using real t-tests.
all_lfc <- matrix(0, nrow = n_genes, ncol = n_datasets)
all_pvals <- matrix(0, nrow = n_genes, ncol = n_datasets)
rownames(all_lfc) <- rownames(all_pvals) <- gene_ids

for (d in 1:n_datasets) {
mat <- datasets[[d]]

# Identify the columns for each group
control_cols <- grep("^Control_", colnames(mat))
both_cols <- grep("^Both_", colnames(mat))

for (i in 1:n_genes) {
control_vals <- mat[i, control_cols]
both_vals <- mat[i, both_cols]

# Real Log2 Fold Change (Difference of means since data is on log scale)
all_lfc[i, d] <- mean(both_vals) - mean(control_vals)

# Real Student's t-test
ttest_res <- t.test(both_vals, control_vals, var.equal = TRUE)
all_pvals[i, d] <- ttest_res$p.value
}
}

# ==============================================================================
# 3. GENERATE REALISTIC VOLCANO PLOTS
# ==============================================================================
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 2.5, 1.5))

lfc_threshold <- 1.0 
p_threshold <- 0.05

for (d in 1:n_datasets) {
lfc <- all_lfc[, d]
neg_log10_p <- -log10(all_pvals[, d])

col_vector <- rep("grey75", n_genes)

is_sig_up <- (lfc >= lfc_threshold) & (all_pvals[, d] <= p_threshold)
is_sig_down <- (lfc <= -lfc_threshold) & (all_pvals[, d] <= p_threshold)

col_vector[is_sig_up & pathways == "Oxidative_Stress"] <- "firebrick2"
col_vector[is_sig_up & pathways == "DNA_Repair"] <- "dodgerblue2"
col_vector[(is_sig_up | is_sig_down) & col_vector == "grey75"] <- "grey30"

plot(lfc, neg_log10_p,
col = col_vector,
pch = 19,
cex = 0.9,
xlab = expression(Log ~ Fold ~ Change),
ylab = expression(-Log ~ p - value),
main = paste("Dataset Replicate", d),
xlim = c(-2.5, 2.5), 
ylim = c(0, 5), 
font.main = 1)

abline(v = c(-lfc_threshold, lfc_threshold), col = "grey50", lty = 3)
abline(h = -log10(p_threshold), col = "grey50", lty = 3)

# Set the legend.
if (d == 1) {
legend("topleft", 
legend = c("Oxidative Stress (Up)", "DNA Repair (Up)", "Other Significant", "Not Significant"),
col = c("firebrick2", "dodgerblue2", "grey30", "grey75"), 
pch = 19, 
cex = 0.6,
box.lty = 0,
bg = "transparent")
}
}

# Reset layout parameters
par(mfrow = c(1, 1))

End R code.

About the simulation

1. Symmetric, Balanced Center: Most of the 100 genes (the 70 background genes) are strictly driven by standard biological and technical noise. This distributes them evenly around 0 on the X-axis, meaning you will naturally see negative fold changes alongside positive ones.
2. Realistic Effect Scales: True upregulations are dialed down to max out around to Fold Changes, mirroring standard biological responses on YPDA plates.
3. True T-Tests: Instead of faking p-values, each simulated dataset now contains 3 biological samples per group. The p-values are calculated dynamically using t.test(), generating a mathematically sound distribution pattern on the Y-axis.

Manhattan plot

A Manhattan plot (Fig 14) is a type of scatter plot used to display the statistical significance of large sets of data points simultaneously. Some feature of the experiment, for example location in the genome or sequential test identifiers, is displayed on the horizontal axis. Like the volcano plot (Fig 13), the vertical axis displays the statistical significance (eg, -log₁₀(p-value)).

Figure 14. Manhattan plot, gene expression fold change on simulated stress test (NaCl, H2O2) yeast data sets, 100 genes (oxidative pathway, DNA repair pathway).

Looking over the plots (Fig 14) we see evidence of pathway clustering: By arranging the 100 genes sequentially across three chromosomes, we see distinct localized peaks rather than a random distribution of significant points. We also see chromosome hotspots: Chromosome I (Oxidative Stress genes) and Chromosome II (DNA Repair genes) display pronounced vertical clusters that reliably break through the red significance line across all five experimental replicates. This simulates a real-world genomic architecture where functionally related gene clusters or operon-like structures respond in unison to environmental shifts.

Conversely, Chromosome III (Background genes, ie, baseline Control) remained entirely flat and well below the significance threshold, validating that the observed stress responses are highly specific and not an artifact of global transcriptional runaway.

R code — with assistance from Gemini AI (Fall 2026) —  for Manhattan plot.

# ==============================================================================
# 1. ASSIGN GEOMETRIC CHROMOSOME POSITIONS TO SIMULATED GENES
# ==============================================================================
# We will use the exact same 5 datasets generated in the previous step.
# To plot a Manhattan plot, every gene needs a Chromosome and a Base Pair (BP) position.

# Assign chromosomes: 
# Chromosome 1: Genes 1-40 (Oxidative Stress pathway)
# Chromosome 2: Genes 41-80 (DNA Repair pathway)
# Chromosome 3: Genes 81-100 (Background/Control pathway)
chroms <- c(rep(1, 40), rep(2, 40), rep(3, 20))

# Assign incremental Base Pair positions along each chromosome
bp_positions <- c(seq(1000, 40000, length.out = 40), 
seq(1000, 40000, length.out = 40), 
seq(1000, 20000, length.out = 20))

# To make a continuous X-axis across multiple chromosomes, we calculate cumulative positions
chr_lengths <- c(45000, 45000, 25000)
chr_offsets <- c(0, cumsum(chr_lengths)[1:2])
coord_x <- bp_positions + chr_offsets[chroms]

# Calculate where the chromosome tick marks should sit on the X-axis (the midpoints)
chr_midpoints <- chr_offsets + (chr_lengths / 2)

# ==============================================================================
# 2. GENERATE MANHATTAN PLOTS FOR ALL 5 REPLICATES
# ==============================================================================
# Set up a clean 2x3 plotting grid to display all 5 replicates together
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 2.5, 1.5))

# Standard scientific significance line (-log10(0.05) is ~1.3)
p_threshold <- 0.05
suggestive_line <- -log10(p_threshold)

for (d in 1:n_datasets) {
neg_log10_p <- -log10(all_pvals[, d])

# Alternating color scheme per chromosome to differentiate boundaries
# Chromosome 1 = Slate Blue, Chromosome 2 = Deep Blue, Chromosome 3 = Dark Grey
chr_colors <- c("lightslateblue", "steelblue4", "grey40")
col_vector <- chr_colors[chroms]

# Create the Manhattan plot grid
plot(coord_x, neg_log10_p,
col = col_vector,
pch = 19,
cex = 1.0,
xaxt = "n", # Suppress the standard numeric X-axis to draw custom Chromosome ticks
xlab = "Chromosome",
ylab = expression(-Log[10] ~ p - value),
main = paste("Dataset Replicate", d),
xlim = c(0, sum(chr_lengths)),
ylim = c(0, 5),
font.main = 1)

# Add the standard genome-wide significance horizontal threshold line
abline(h = suggestive_line, col = "firebrick3", lty = 2, lwd = 1.2)

# Draw the custom structural X-axis labeling the specific chromosomes
axis(1, at = chr_midpoints, labels = c("Chr I", "Chr II", "Chr III"))

# Add an explanatory legend to the first plot panel to map pathways to locations
if (d == 1) {
legend("topright", 
legend = c("Chr I (Oxidative Stress)", "Chr II (DNA Repair)", "Chr III (Background)"),
col = c("lightslateblue", "steelblue4", "grey40"), 
pch = 19, 
cex = 0.65,
box.lty = 0,
bg = "transparent")
}
}

# Reset layout parameters back to default single screen layout
par(mfrow = c(1, 1))

About the R code.

The Genomic Peaks: Because these simulated Oxidative Stress genes sit together on Chromosome I and DNA Repair genes sit on Chromosome II,  we see pronounced, beautiful vertical clusters extending upwards over the red significance line.

The Noise Flatline: Chromosome III hosts the background genes, creating a realistically flat landscape where points rarely cross the significance threshold.

Alternating Chromosome Tracks: Using the standard coordinate offset method, the data points stay organized left-to-right without overlapping.

Questions

1. Using our Comet assay data set (Table 1, Chapter 4.2), create scatter plots to show associations between tail length, tail percent, and olive moment. 
2. Explore different settings including size of points, amount of white area, and scale of the axes. Evaluate how these changes change the “story” told by the graph.

Quiz Chapter 4.5

Scatter plots

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Data sets

Number of pixels of 1 cm length from unit measure on ten digital images. Recorded by two independent observers, Obs1 and Obs2.

 Obs1, Obs2
171.026, 171.105
136.528, 138.521
148.084, 144.222
142.014, 140.057
150.213, 153.118
187.011, 195.092
168.760, 168.668
154.302 ,160.381
209.022, 209.876
240.067, 161.805

Gene expression

ccc

Chapter 4 contents

• Graphs and tables (How to report statistics)
• Bar (column) charts
• Histograms
• Box plot
• Mosaic plots
• Scatter plots
• Add a second Y axis
• Q-Q plot
• Ternary plots
• Heat maps
• Graph software
• References and suggested readings