ENCODE Gene Expression Exploration
Christopher Isaac
2026-09-14
library(tidyverse)Introduction
RNA sequencing can be used to measure gene expression and identify which genes are active within a tissue sample. In this analysis, I explored RNA-seq gene quantification data from the ENCODE Project collected from dorsolateral prefrontal cortex tissue. The sample was obtained from a female participant over 90 years old from the RUSH Alzheimer’s disease study.
The goal of this exploratory analysis was to examine the overall distribution of gene expression, identify the most highly expressed genes, and investigate whether gene length is related to expression level.
Because this dataset represents a single sample, this analysis describes the gene-expression profile of this tissue but cannot determine which expression patterns are specifically caused by Alzheimer’s disease.
Data Import and Exploration
The gene expression data were downloaded from the ENCODE Project as a tab-separated values (TSV) file. The file contains gene-level RNA-seq quantification measurements.
expression_data <- read_tsv("ENCFF166SFX.tsv")## Rows: 59526 Columns: 17
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (2): gene_id, transcript_id(s)
## dbl (15): length, effective_length, expected_count, TPM, FPKM, posterior_mea...
##
## â„ą Use `spec()` to retrieve the full column specification for this data.
## â„ą Specify the column types or set `show_col_types = FALSE` to quiet this message.Before creating plots, I examined the structure of the dataset to identify the available variables and their data types.
expression_data## # A tibble: 59,526 Ă— 17
## gene_id `transcript_id(s)` length effective_length expected_count TPM FPKM
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 10904 10904 93 0 0 0 0
## 2 12954 12954 94 0 0 0 0
## 3 12956 12956 72 0 0 0 0
## 4 12958 12958 82 0 0 0 0
## 5 12960 12960 73 0 0 0 0
## 6 12962 12962 72 0 0 0 0
## 7 12964 12964 74 0 0 0 0
## 8 12965 12965 82 0 0 0 0
## 9 12967 12967 73 0 0 0 0
## 10 12969 12969 73 0 0 0 0
## # â„ą 59,516 more rows
## # â„ą 10 more variables: posterior_mean_count <dbl>,
## # posterior_standard_deviation_of_count <dbl>, pme_TPM <dbl>, pme_FPKM <dbl>,
## # TPM_ci_lower_bound <dbl>, TPM_ci_upper_bound <dbl>,
## # TPM_coefficient_of_quartile_variation <dbl>, FPKM_ci_lower_bound <dbl>,
## # FPKM_ci_upper_bound <dbl>, FPKM_coefficient_of_quartile_variation <dbl>Plot 1
Question
How does gene expression vary across genes in this brain sample?
glimpse(expression_data)## Rows: 59,526
## Columns: 17
## $ gene_id <chr> "10904", "12954", "12956", "129…
## $ `transcript_id(s)` <chr> "10904", "12954", "12956", "129…
## $ length <dbl> 93, 94, 72, 82, 73, 72, 74, 82,…
## $ effective_length <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ expected_count <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ TPM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ FPKM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ posterior_mean_count <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ posterior_standard_deviation_of_count <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ pme_TPM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ pme_FPKM <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ TPM_ci_lower_bound <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ TPM_ci_upper_bound <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ TPM_coefficient_of_quartile_variation <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ FPKM_ci_lower_bound <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ FPKM_ci_upper_bound <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ FPKM_coefficient_of_quartile_variation <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…names(expression_data)## [1] "gene_id"
## [2] "transcript_id(s)"
## [3] "length"
## [4] "effective_length"
## [5] "expected_count"
## [6] "TPM"
## [7] "FPKM"
## [8] "posterior_mean_count"
## [9] "posterior_standard_deviation_of_count"
## [10] "pme_TPM"
## [11] "pme_FPKM"
## [12] "TPM_ci_lower_bound"
## [13] "TPM_ci_upper_bound"
## [14] "TPM_coefficient_of_quartile_variation"
## [15] "FPKM_ci_lower_bound"
## [16] "FPKM_ci_upper_bound"
## [17] "FPKM_coefficient_of_quartile_variation"dim(expression_data)## [1] 59526 17summary(expression_data$TPM)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 0.00 16.80 1.52 267430.21ggplot(expression_data, aes(x = TPM)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
expression_plot1 <- expression_data %>%
mutate(log_TPM =log10(TPM + 1))ggplot(expression_plot1, aes(x = log_TPM)) +
geom_histogram(bins = 50)
ggplot(expression_plot1, aes(x = log_TPM)) +
geom_histogram(bins = 50) +
labs(
title = "Distribution of Gene Expression in Human Brain Tissue",
subtitle = "Dorsolateral prefrontal cortex RNA-seq sample",
x = "log10(TPM + 1)",
y = "Number of Genes"
) +
theme_minimal()
# Interpretation This histogram shows the distribution of gene
expression across the genes detected in the dorsolateral prefrontal
cortex sample. Most genes appear at relatively low expression levels,
while fewer genes show very high expression, demonstrating that gene
expression is unevenly distributed across the transcriptome.
The log transformation makes this pattern easier to visualize by compressing the range of very highly expressed genes without removing them from the dataset. Because this analysis represents a single sample, the plot describes the expression profile of this tissue but does not by itself show which genes are specifically altered by Alzheimer’s disease.
Plot 2
Question
Which genes are the most highly expressed in this brain sample?
To answer this question, the genes were ranked from highest to lowest TPM, and the ten genes with the highest expression values were selected.
expression_data %>%
arrange(desc(TPM)) %>%
select(gene_id, TPM) %>%
head(10)## # A tibble: 10 Ă— 2
## gene_id TPM
## <chr> <dbl>
## 1 ENSG00000276168.1 267430.
## 2 ENSG00000274012.1 227199.
## 3 ENSG00000283293.1 45111.
## 4 ENSG00000222328.1 26414.
## 5 ENSG00000277209.1 21961.
## 6 ENSG00000198804.2 9164.
## 7 ENSG00000277027.1 8207.
## 8 ENSG00000198899.2 5863.
## 9 ENSG00000198938.2 4786.
## 10 ENSG00000198888.2 4704.top10_genes <- expression_data %>%
arrange(desc(TPM)) %>%
select(gene_id, TPM) %>%
head(10)ggplot(top10_genes, aes(x = gene_id, y = TPM)) +
geom_col()
ggplot(top10_genes, aes(x = gene_id, y = TPM)) +
geom_col() +
coord_flip()
ggplot(top10_genes, aes(x = reorder(gene_id, TPM), y = TPM)) +
geom_col() +
coord_flip()
ggplot(top10_genes, aes(x = reorder(gene_id, TPM), y = TPM)) +
geom_col() +
coord_flip() +
labs(
title = "Top 10 Most Highly Expressed Genes",
subtitle = "Dorsolateral prefrontal cortex RNA-seq sample",
x = "Gene ID",
y = "TPM"
) +
theme_minimal()
top10_genes <- top10_genes %>%
mutate(
gene_symbol = case_when(
gene_id == "ENSG00000276168.1" ~ "RN7SL1",
gene_id == "ENSG00000274012.1" ~ "RN7SL2",
gene_id == "ENSG00000283293.1" ~ "RN7SK",
gene_id == "ENSG00000222328.1" ~ "RNU2-2",
gene_id == "ENSG00000277209.1" ~ "RPPH1",
gene_id == "ENSG00000198804.2" ~ "MT-CO1",
gene_id == "ENSG00000277027.1" ~ "RMRP",
gene_id == "ENSG00000198899.2" ~ "MT-ATP6",
gene_id == "ENSG00000198938.2" ~ "MT-CO3",
gene_id == "ENSG00000198888.2" ~ "MT-ND1"
)
)The Ensembl gene IDs were converted to more recognizable gene symbols to make the visualization easier to interpret.
ggplot(
top10_genes,
aes(x = reorder(gene_symbol, TPM), y = TPM)
) +
geom_col() +
coord_flip() +
scale_y_continuous(labels = scales::comma) +
labs(
title = "Top 10 Most Highly Expressed Genes",
subtitle = "Dorsolateral prefrontal cortex RNA-seq sample",
x = "Gene",
y = "TPM"
) +
theme_minimal()
## Interpretation This bar chart shows the ten genes with the highest
TPM values in the dorsolateral prefrontal cortex RNA-seq sample. RN7SL1
and RN7SL2 have substantially higher expression than the remaining
genes, showing that a small number of genes account for some of the
highest RNA abundance observed in this sample. Several other highly
expressed genes also show much lower TPM values by comparison,
reinforcing the uneven expression pattern seen in the first plot.
Because this analysis represents a single individual, these results describe which genes are most abundant in this tissue sample but do not indicate that their expression is specifically caused by Alzheimer’s disease.
Plot 3
Question
Is gene length related to expression level in this brain sample?
To examine the relationship between gene length and expression, genes with TPM values greater than zero were selected. Zero values were removed because a logarithmic scale cannot display zero.
summary(expression_data$length)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8 387 808 1339 1823 205012ggplot(expression_data, aes(x = length, y = TPM)) +
geom_point()
expression_plot3 <- expression_data %>%
filter(TPM > 0)ggplot(expression_plot3, aes(x = length, y = TPM)) +
geom_point(alpha = 0.2) +
scale_x_log10() +
scale_y_log10()
ggplot(expression_plot3, aes(x = length, y = TPM)) +
geom_point(alpha = 0.2) +
geom_smooth() +
scale_x_log10() +
scale_y_log10()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggplot(expression_plot3, aes(x = length, y = TPM)) +
geom_point(alpha = 0.2) +
geom_smooth() +
scale_x_log10() +
scale_y_log10(labels = scales::comma) +
labs(
title = "Relationship Between Gene Length and Expression",
subtitle = "Dorsolateral prefrontal cortex RNA-seq sample",
x = "Gene Length (log10 scale)",
y = "TPM (log10 scale)"
) +
theme_minimal()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Interpretation This scatterplot examines the relationship between
gene length and gene expression in the dorsolateral prefrontal cortex
RNA-seq sample. Gene expression varies widely across genes of similar
lengths, and the smoothed trend does not show a strong, consistent
relationship between increasing gene length and increasing TPM.
Both axes were displayed on logarithmic scales to make the wide ranges of gene length and expression easier to visualize. Because TPM is normalized in part for transcript length, it is reasonable that longer genes do not simply have higher TPM values. The greater uncertainty among the longest genes likely reflects the smaller number of observations in that range.
Conlcusion
This exploratory analysis identified several characteristics of gene expression within the dorsolateral prefrontal cortex RNA-seq sample. Most genes showed relatively low expression, while a much smaller number showed extremely high RNA abundance. The most highly expressed genes included several abundant non-coding RNAs and mitochondrial genes.
The analysis also showed that gene length does not have a simple relationship with normalized gene expression. Genes of similar lengths displayed a wide range of TPM values, suggesting that many other biological and regulatory factors influence gene expression.
Because this dataset contains a single sample, additional samples and appropriate comparison groups would be necessary to determine whether any of these expression patterns are specifically associated with Alzheimer’s disease.