Week 2 Project

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## Loading required package: ggplot2

Read in the ENCFF166SFX expression data file from ENCODE

expressionData <- read.table("ENCFF166SFX.tsv", header = TRUE)

Determine the distribution of gene lengths in this sample

When first looking at the data, the most basic information we have is gene lengths. Transcriptomic data was collected and mapped to genes, creating a basis for gene expression analysis. Let’s plot the distribution of gene lengths in this sample. The majority of the gene lengths in this sample are less than 25000 bp. However, since the scale of the histogram extends to 200000 bp, it is likely that there are a few genes that are larger than 25000 bp and reach a length that is a little over 200000 bp. We will investigate this further when we create the next plot.

ggplot(expressionData, aes(x = length)) +
  geom_histogram(binwidth = 100, color = "blue") +
  labs(
    title = "Distribution of Gene Length",                     
    x = "Length (bp)",
    y = "Count"
  )         

Analyze the relationship between gene length and the effective length for each gene ID

Effective length is the actual mappable length of a read in comparison to its mapped gene. It is likely that the longer the gene length, the longer the effective length is because having a longer gene length creates a larger window for a continuous mappable read to map to. This plot does indeed show that gene length and effective length have a strong correlation. Also, we can now see that there are several genes that greater than 25000 bp in length. There is a gene that is actually greater than 200000 bp long, as was seen in the previous histogram. Effective length for each gene allows us to calculate normalized RNA-seq metrics such as TPM and FPKM, which then further tell us relative gene expression levels for specific genes.

ggplot(data = expressionData, mapping = aes(x = length, y = effective_length)) +
  geom_point(color = "blue", size = 1, alpha = 1.0) + 
  geom_smooth(method = "lm", se = TRUE, linewidth = 0.25, lty = "dashed", color = "purple") +
  labs(
    title = "The Relationship of Gene Length and Effective Length",                     
    x = "Length (bp)",
    y = "Effective Length (bp)"
  )         

## `geom_smooth()` using formula = 'y ~ x'

Identify high levels of gene expression

The goal of RNA-seq is to determine what relative levels of gene expression exist in a sample, especially a sample that has a disease and/or been administered a treatment. To take a snapshot of what could be happening within this patient, let’s visualize high levels of expression (high FPKM) since this is more apparent than downregulated gene expression. RN7SL1 and RN7SL2 have the highest levels of gene expression in this filtered data frame. These genes encode for proteins that are part of the signal recognition particle complex (https://pmc.ncbi.nlm.nih.gov/articles/PMC5859127/).

highExpression <- expressionData |>
  filter(FPKM > 5000)
ggplot(highExpression, aes(x = gene_id, y = FPKM, fill = FPKM)) +
  geom_col() +
  scale_fill_viridis_c(option = "turbo") +
  scale_x_discrete(name = "Gene Name", labels = c("ENSG00000198804.2" = "MT-CO1", "ENSG00000198888.2" =
  "MT-ND1", "ENSG00000198899.2" = "MT-ATP6", "ENSG00000198938.2" = "MT-CO3", "ENSG00000222328.1" =
  "RNU2-2","ENSG00000263934.4" = "SNORD3A", "ENSG00000274012.1" = "RN7SL2", "ENSG00000276168.1" = "RN7SL1",
  "ENSG00000277027.1" = "RMRP", "ENSG00000277209.1" = "RPPH1", "ENSG00000283293.1" = "RN7SK")) +
  labs(
    title = "Highly Expressed Genes",
    subtitle = "Genes with FPKM > 5000"
  ) +
  theme(axis.text.x = element_text(angle = 50, vjust = 1, hjust = 1))