Analysis of Gene Expression in Sample of Brain with Alzheimers
Isabel Di Donna
2026-09-10
Load in data
genexp <- read.table("ENCFF166SFX.tsv", header = TRUE)Investigate the data
dim(genexp)## [1] 59526 17names(genexp)## [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"str(genexp)## 'data.frame': 59526 obs. of 17 variables:
## $ gene_id : chr "10904" "12954" "12956" "12958" ...
## $ transcript_id.s. : chr "10904" "12954" "12956" "12958" ...
## $ length : num 93 94 72 82 73 72 74 82 73 73 ...
## $ effective_length : num 0 0 0 0 0 0 0 0 0 0 ...
## $ expected_count : num 0 0 0 0 0 0 0 0 0 0 ...
## $ TPM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ posterior_mean_count : num 0 0 0 0 0 0 0 0 0 0 ...
## $ posterior_standard_deviation_of_count : num 0 0 0 0 0 0 0 0 0 0 ...
## $ pme_TPM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ pme_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ TPM_ci_lower_bound : num 0 0 0 0 0 0 0 0 0 0 ...
## $ TPM_ci_upper_bound : num 0 0 0 0 0 0 0 0 0 0 ...
## $ TPM_coefficient_of_quartile_variation : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FPKM_ci_lower_bound : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FPKM_ci_upper_bound : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FPKM_coefficient_of_quartile_variation: num 0 0 0 0 0 0 0 0 0 0 ...summary(genexp)## gene_id transcript_id.s. length effective_length
## Length :59526 Length :59526 Min. : 8 Min. : 0
## N.unique :59526 N.unique :59526 1st Qu.: 387 1st Qu.: 287
## N.blank : 0 N.blank : 0 Median : 808 Median : 708
## Min.nchar: 4 Min.nchar: 4 Mean : 1339 Mean : 1240
## Max.nchar: 24 Max.nchar: 3455 3rd Qu.: 1823 3rd Qu.: 1723
## Max. :205012 Max. :204912
## expected_count TPM FPKM
## Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.0 1st Qu.: 0.00 1st Qu.: 0.00
## Median : 0.0 Median : 0.00 Median : 0.00
## Mean : 339.3 Mean : 16.80 Mean : 19.55
## 3rd Qu.: 46.0 3rd Qu.: 1.52 3rd Qu.: 1.77
## Max. :1251092.4 Max. :267430.21 Max. :311258.58
## posterior_mean_count posterior_standard_deviation_of_count pme_TPM
## Min. : 0.0 Min. : 0.000 Min. : 0.00
## 1st Qu.: 0.0 1st Qu.: 0.000 1st Qu.: 0.09
## Median : 0.0 Median : 0.000 Median : 0.30
## Mean : 339.3 Mean : 1.965 Mean : 16.80
## 3rd Qu.: 46.0 3rd Qu.: 0.000 3rd Qu.: 3.00
## Max. :1068373.8 Max. :22576.980 Max. :217493.25
## pme_FPKM TPM_ci_lower_bound TPM_ci_upper_bound
## Min. : 0.00 Min. :0.000e+00 Min. :0.000e+00
## 1st Qu.: 0.11 1st Qu.:0.000e+00 1st Qu.:2.530e-01
## Median : 0.36 Median :3.000e-03 Median :6.790e-01
## Mean : 20.22 Mean :1.502e+01 Mean :1.923e+01
## 3rd Qu.: 3.61 3rd Qu.:1.072e+00 3rd Qu.:4.465e+00
## Max. :261833.14 Max. :2.169e+05 Max. :2.181e+05
## TPM_coefficient_of_quartile_variation FPKM_ci_lower_bound FPKM_ci_upper_bound
## Min. :0.0000 Min. :0.000e+00 Min. :0.000e+00
## 1st Qu.:0.0759 1st Qu.:0.000e+00 1st Qu.:3.040e-01
## Median :0.3873 Median :4.000e-03 Median :8.180e-01
## Mean :0.3677 Mean :1.808e+01 Mean :2.316e+01
## 3rd Qu.:0.6553 3rd Qu.:1.292e+00 3rd Qu.:5.370e+00
## Max. :0.7098 Max. :2.612e+05 Max. :2.625e+05
## FPKM_coefficient_of_quartile_variation
## Min. :0.00000
## 1st Qu.:0.07591
## Median :0.38733
## Mean :0.36768
## 3rd Qu.:0.65535
## Max. :0.70981Load packages
library(ggplot2)
library(dplyr) ##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionPlot 1: A Distribution of Expressed Genes
The range for TPM values is too large to be directly plotted, so we will make a new column taking the log of TPM to compress the raw data. Then compress further so we get genes with a TPM > 0.
genexp$logTPM <- log10(genexp$TPM + 1)
expressed_genes <- genexp[genexp$TPM > 0, ]Now we make the histogram. This will show the distribution of genes with detectable expression in the brain tissue sample. For this example, the data is concentrated on the left side of the graph, which indicates many genes with low expression. The long right tail indicates fewer genes with higher levels of expression.
ggplot(expressed_genes, aes(x = log10(TPM + 1))) +
geom_histogram(
bins = 50,
fill = "skyblue",
color = "blue"
) +
labs(
title = "Distribution of Expressed Genes",
x = "Log of Gene Expression",
y = "Number of Genes"
) +
theme_minimal()
Plot 2: 10 Highest Expressed Genes
For this, we have to sort the data from highest TPM to lowest TPM
top_genes <- genexp %>%
arrange(desc(TPM))Then we want to keep only the first 10 rows and check that the code was successful
top10 <- head(top_genes, 10)
top10## gene_id transcript_id.s. length effective_length expected_count
## 1 ENSG00000276168.1 ENST00000618786.1 299 199 1251092.39
## 2 ENSG00000274012.1 ENST00000490232.3 300 200 1068223.36
## 3 ENSG00000283293.1 ENST00000636484.1 328 228 241793.72
## 4 ENSG00000222328.1 ENST00000410396.1 191 91 56506.22
## 5 ENSG00000277209.1 ENST00000516869.1 333 233 120292.74
## 6 ENSG00000198804.2 ENST00000361624.2 1542 1442 310661.74
## 7 ENSG00000277027.1 ENST00000363046.1 264 164 31642.18
## 8 ENSG00000198899.2 ENST00000361899.2 681 581 80074.14
## 9 ENSG00000198938.2 ENST00000362079.2 784 684 76961.80
## 10 ENSG00000198888.2 ENST00000361390.2 956 856 94658.58
## TPM FPKM posterior_mean_count
## 1 267430.21 311258.58 990565.47
## 2 227198.90 264433.89 1068373.81
## 3 45111.19 52504.34 237856.45
## 4 26413.69 30742.55 56506.16
## 5 21961.28 25560.45 120288.98
## 6 9164.24 10666.15 310659.39
## 7 8207.23 9552.29 31631.68
## 8 5862.60 6823.41 80075.12
## 9 4786.23 5570.63 76961.81
## 10 4703.92 5474.84 94657.66
## posterior_standard_deviation_of_count pme_TPM pme_FPKM TPM_ci_lower_bound
## 1 22576.98 202666.69 243984.09 194281.00
## 2 944.77 217493.25 261833.14 216914.00
## 3 1510.02 42475.06 51134.37 41894.50
## 4 4.04 25282.16 30436.38 25071.40
## 5 48.54 21019.66 25304.89 20898.80
## 6 18.26 8771.48 10559.70 8737.43
## 7 359.14 7853.16 9454.16 7652.58
## 8 31.45 5611.50 6755.51 5571.47
## 9 1.43 4581.18 5515.13 4547.77
## 10 6.86 4502.34 5420.23 4472.52
## TPM_ci_upper_bound TPM_coefficient_of_quartile_variation FPKM_ci_lower_bound
## 1 211101.00 0.017216000 233930.00
## 2 218085.00 0.000925635 261215.00
## 3 42958.50 0.004696260 50459.40
## 4 25487.00 0.002794370 30186.30
## 5 21140.80 0.002014630 25166.20
## 6 8803.78 0.001300840 10525.90
## 7 8046.70 0.008376890 9211.81
## 8 5650.71 0.002417070 6710.44
## 9 4613.26 0.002448150 5478.15
## 10 4531.38 0.002237810 5386.53
## FPKM_ci_upper_bound FPKM_coefficient_of_quartile_variation logTPM
## 1 254217.00 0.017235700 5.427212
## 2 262540.00 0.000869774 5.356408
## 3 51736.10 0.004722920 4.654294
## 4 30691.20 0.002822800 4.421846
## 5 25454.70 0.001968750 4.341677
## 6 10599.30 0.001191460 3.962144
## 7 9685.98 0.008379530 3.914250
## 8 6804.06 0.002367760 3.768164
## 9 5555.65 0.002424780 3.680084
## 10 5455.62 0.002191620 3.672552Then we code for the ggplot bar chart. This will showw the ten genes with the highest TPM values in the brain sample. They represent a small subset of genes expressed at high levels. For this example, there is a large difference between the most highly expressed gene and the tenth most highly expressed gene because, as we saw from the skewed bell curve from Plot 1, there are very few genes with high expression levels in the dataset.
ggplot(top10, aes(x = reorder(gene_id, TPM), y = TPM)) +
geom_col(fill = "pink", color = "hotpink") +
coord_flip() +
labs(
title = "10 Most Highly Expressed Genes",
x = "Gene ID",
y = "TPM"
) +
theme_minimal()
Plot 3: Gene lenth vs. Expression
A scatterplot to show the relationship between gene length and expression level. The expression is skewed, so we use the log instead of the raw data for gene expression. For this example, there appears to be no correlation between gene expression and gene length, suggesting influence from other biological factors.
ggplot(genexp, aes(x = length, y = logTPM)) +
geom_point(alpha = 0.2) +
labs(
title = "Gene Length and Expression",
x = "Gene Length",
y = "Log of Gene Expression"
) +
theme_minimal()