diamond plots
2026-09-14
#Step 1: load in gene expression data
# myGex = my gene expression
myGex <- read.table("ENCFF166SFX.tsv", header = TRUE)
head(myGex)## gene_id transcript_id.s. length effective_length expected_count TPM FPKM
## 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
## posterior_mean_count posterior_standard_deviation_of_count pme_TPM pme_FPKM
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## TPM_ci_lower_bound TPM_ci_upper_bound TPM_coefficient_of_quartile_variation
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## FPKM_ci_lower_bound FPKM_ci_upper_bound
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0
## FPKM_coefficient_of_quartile_variation
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0#Step 2: Plot 1
# plot 1: Log10 of TPM Distribution
plot1 <- ggplot(myGex, aes(x = log10(TPM + 1))) +
geom_histogram(bins = 50, fill = "darkslateblue", color = "white") +
# logg10(TPM + 1) compresses numbers so low and high expression cna be compared in the same axis; bins splits data into a different number of bars
labs(
title = "Distribution of Gene Expression in 90 yo Female Brain",
subtitle = "ENCODE Accession: ENCSR562BUN",
x = "Log10 (TPM +1)",
y = "Number of Genes"
) +
theme_minimal()
print(plot1)
Step 3: sort the top 15 highest expressed genes
top_15 <- myGex %>%
arrange(desc(TPM)) %>%
# arranges TPM from highest to lowest
slice(1:15)Step 4: Plot 2
# plot of top expressed genes
plot2 <- ggplot(top_15, aes(x = reorder(gene_id, TPM), y = TPM)) +
geom_col(fill = "lightblue") +
coord_flip() +
# turns vertical bars into horizontal bars to be able to read the gene names
labs(
title = "Top 15 Most Highly Expressed Genes",
x = "gene ID",
y = "Expression (TPM)"
) +
theme_light()
print(plot2)
Step 5: plot 3 for gene length vs expression
plot3 <- ggplot(myGex, aes(x = length, y = log10(TPM + 1))) +
# draws a dot for every gene, alpha denotes transparency so places with more genes (aka dots) looks darker
geom_point(alpha = 0.2, color = "lightcoral") +
# fits a trend line through the dots to show if genes have lower or higher expression
geom_smooth(method = "lm", color = "black", se = FALSE) +
scale_x_log10() +
labs(title = "Gene Length vs Gene Expression Level",
x = "Gene Length (bp)",
y = "Expression (Log10)") +
theme_bw()
print(plot3)## `geom_smooth()` using formula = 'y ~ x'