Reading a Gene Counts Table

In RNA sequencing, a counts table records how many sequencing reads mapped to each gene in each sample. Each row is one gene and each column is one sample, and every value is a whole number called a count. We will load this table together with a metadata table that says which biological group each sample belongs to.

Loading The Tables

# read.csv reads a comma-separated file into a data frame (a table in R)
# row.names = 1 tells R to use the first column (the gene names) as row labels
counts <- read.csv("counts.csv", row.names = 1)

# head() prints just the first 6 rows so we can eyeball the data
head(counts)

       control_1 control_2 treated_1 treated_2
ACTB        1204      1350      1189      1275
GAPDH        980      1042       999      1010
TP53         145       160       210       198
MYC          320       298       540       512
GATA3         88        95        41        38
VIM          410       430       405       420

# load the sample information table the same way
meta <- read.csv("metadata.csv", row.names = 1)

# this table is small, so print the whole thing
meta

          condition replicate
control_1   control         1
control_2   control         2
treated_1   treated         1
treated_2   treated         2

Library Size Checks

# colSums() adds up every value in each column, one total per sample
# this per-sample total is the library size, a basic quality check
colSums(counts)

control_1 control_2 treated_1 treated_2 
 18452030  19113445  17984201  18650770 

Raw counts range from zero to tens of thousands, so a handful of very highly expressed genes dominate any plot. Taking a base-2 logarithm compresses that range and makes the samples easier to compare. We add 1 to every count first, a trick called a pseudocount, because log2(0) evaluates to -Inf (negative infinity) in R, which would distort the plot.

# add a pseudocount of 1, then take log base 2 of every value
log_counts <- log2(counts + 1)

# boxplot() draws one box per column so we can compare sample distributions
# las = 2 turns the sample names sideways so they all fit under the boxes
boxplot(log_counts,
        main = "Log2 counts per sample",
        ylab = "log2(count + 1)",
        las = 2)