Step One: Setting up packages and loading in data

Included is a list of packages that were utilized for this project.

Since biomaRt is from the Bioconductor repository and not CRAN, install.packages() does not work with biomaRt. In order to download the biomaRt package, a different script is needed to set up the package.

Reading in the data is straightforward. This data was obtained from the ENCODE database from the RUSH Alzheimer’s study. (https://www.encodeproject.org/experiments/ENCSR562BUN/)

Gex <- read.table('ENCFF166SFX.tsv', header = TRUE)

The Gex dataframe contains a list of genes using the Ensembl gene ID. To make data visualization easier, biomaRt was utilized to download a list of the more recognizable HGNC gene symbols.

mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = "https://www.ensembl.org")
GeneID <- sub("\\..*", "", Gex$gene_id)
genelist <- getBM(filters = "ensembl_gene_id",
                  attributes = c("ensembl_gene_id", "hgnc_symbol"),
                  values = GeneID,
                  mart = mart)
genelist <- genelist %>%
  rename(gene_id = ensembl_gene_id)

The list produced containing the Ensembl gene IDs and their associate HGNC symbols were merged together. This can be achieved my removing the Ensembl gene variant suffixes and merging the two dataframes.

Step Two: Finding the highest expressed genes in the dataset

This list will contain genes where the FPKM is above 4000. This does not mean that genes under this threshold do not hold any significance, but this will create a short comprehensive list of the most expressed genes from the dataset.

topscripts <- Gex %>%
  filter(FPKM > 4000)
topscripts <- topscripts %>%
  mutate(gene_id = str_replace(gene_id, "\\..*", ""))

A dataframe containing the genes with over 4000 FPKM and their associate HGNC gene symbols was formed to create a plot.

plotgenes <- merge(topscripts, genelist, by = "gene_id", all.x = TRUE)

ggplot(data = plotgenes,
      aes(x = reorder(hgnc_symbol, FPKM), y = FPKM,
          fill = FPKM)) +
  labs(y = 'FPKM',
       x = 'Gene',
       fill = 'FPKM',
       title = 'Top Represented Genes by FPKM') +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90),
        panel.background = element_rect(fill = 'grey60')) +
  scale_fill_gradientn("FPKM", colors = rev(brewer.pal(9, name = "Oranges"))) +
  geom_col()

This plot mapping FPKM data shows genes found to have a high number of transcripts for this gene. As gene length can be a factor for how many reads map back (as in a longer gene sequence will likely have more transcripts mapped to it), FPKM attempts to normalize that data so longer genes are not over represented. This data shows the 12 highest genes expressed to the brain tissue and could suggest genes that are overexpressed due to Alzheimer’s.

Step Three: Plotting the relationship of FPKM vs. Expected Counts

The relationship between expected count and FPKM can show potential outliers in the expression data. The expected count is an RSEM output that takes all expression levels into consideration to estimate the amount of reads each gene show have. By plotting this against FPKM, outliers that have either higher or lower expression than expected may reveal genes that have a relationship to the disease in a specific patient.

highlight <- data.frame(xmin = 1800, xmax = 15000,
                        ymin = 200000, ymax = 350000)
# the outliers were determined visually after the creation of the initial plot, and this highlight dataframe was then created and added on to the initial plot

ggplot(data = Gex,
      aes(x = FPKM, y = expected_count)) +
  labs(x = 'FPKM',
       y = 'Expected Counts',
       title = 'FPKM vs Expected Counts') +
  stat_smooth(method = 'lm', formula = y ~ x, geom = 'smooth') +
  theme_bw() +
  theme(panel.background = element_rect(fill = "white")) +
  geom_point(shape = 1) +
  geom_rect(data = highlight,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            color = 'green', fill = NA, linewidth = 0.5, inherit.aes = FALSE)

# the visually confirmed outliers can then be plotted into their own graph and labelled for identification
outliers <- Gex %>%
  filter(expected_count > 250000, expected_count < 320000)
outliers <- outliers %>%
  mutate(gene_id = str_replace(gene_id, "\\..*", ""))
plotoutliers <- merge(outliers, genelist, by = "gene_id", all.x = TRUE)


ggplot(data = plotoutliers,
      aes(x = FPKM, y = expected_count)) +
  labs(x = 'FPKM',
       y = 'Expected Counts',
       title = 'FPKM vs Expected Counts Outliers') +
  theme_bw() +
  theme(panel.background = element_rect(fill = "white")) +
  geom_point(shape = 16) +
  geom_text_repel(label = plotoutliers$hgnc_symbol)

This plot highlights two genes that could have a connection to disease. Seeing as they sit above the regression line on the first plot, this plot suggests that higher expression of these genes was expected the pattern of their estimated count was lower than the number of transcripts mapped, which could suggest under expression of MALAT1 and MT-C01 in this patient’s brain tissue.

Step Four: Looking at variation within the FPKM column

To find genes with high expression, a boxplot can be made to find the outliers the higher or lower expression than the rest of the dataset. To achieve this, the Gex dataset was filtered so genes with no expected counts AND no FPKM values were removed too look at only genes with measured expression levels OR genes that were expected to be expressed. The new list of FPKM values were increased by 1 for log10 transformation so FPKM values of 0 would not be removed during plotting.

Z_filtered <- Gex[!(Gex$expected_count == 0 & Gex$FPKM == 0), ]

Z_filtered$FPKM <- log10(Z_filtered$FPKM + 1)

Z_filtered <- Z_filtered %>%
  mutate(gene_id = str_replace(gene_id, "\\..*", ""))

Z_filtered <- merge(Z_filtered, genelist, by = "gene_id", all.x = TRUE)

The following boxplot was created with the filtered and log transformed dataset to look at the overall variability in expression levels.

# this dataframe was created with the 5 highest transformed FPKM values to show genes with high expression.
boxoutliers <- Z_filtered %>%
  filter(hgnc_symbol %in% c('RN7SL1', 'RN7SL2', 'RN7SK', 'RNU2-2', 'RPPH1'))

ggplot(data = Z_filtered,
       aes(x = "",
           y = FPKM)) +
  labs(x = NULL,
       y = "log10(FPKM + 1)",
       title = "Variability of FPKM Values") +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white")) +
  geom_boxplot(width = 0.4, fill = 'cyan') +
  geom_text_repel(data = boxoutliers,
                  aes(x = "", y = FPKM, label = hgnc_symbol),
                  size = 4) +
  stat_summary(fun = median, geom = "point", color = 'red', size = 2) +
  stat_summary(fun = median, geom = "text", color = 'red', size = 3, vjust = -1, aes(label = round(after_stat(y), 1)))

This plot helps show that there is high variability in the FPKM values. While FPKM attempts to normalize the transcript count by gene length and sequencing depth, this plot helps highlight genes with high FPKM counts from the tissue showing potential genes of interest when studying this patient and the relationship of these genes to the disease.