---
title: Analysis of NSD3 Isoform Expression from TCGA-LUSC Data 
author: David Dilworth and Dalia Barsyte-Lovejoy
date: '2018-08-14'
slug: NSD3-isoform
output:
  pdf_document: default
  fig_caption: yes
html_document: default
word_document: default
tags:
  thumbnail: ''
categories:
- NSD3
- Open NB
description: 'exp-028'
---
  
```{r, echo = FALSE}
library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=FALSE, echo = TRUE, message = FALSE, warning = FALSE)
```
  
# Objective

There are two major NSD3 isoforms expressed, long (aa 1-1437) and short (aa 1-645, differing in sequence from 620-645). Importantly, the short isoform was shown to be required for the maintenance of acute myeloid leukemia (AML) [1]. This isoform lacks a SET domain and thus methyltransferase activity. Additionally, my preliminary experiments suggest it is the short isoform that mediates NSD3 involvement in observed EMT phenotypes. It is not clear if these the two isoforms are simply co-expressed or independently regulated. As a first attempt to study the differential regulation of the two isoforms, I have used the TCGA-LUSC data set to analyze relative expression levels in the context of squamous cell lung cancer (LUSC) [2]. 

# Code & Results

This analysis will rely on three R packages; tidyverse to clean, tidy, and plot the data; ggbeeswarm to extend the functionality of ggplot to include beeswarm plots; and ggfortify for auto-ploting linear regression models. These packages are first loaded below.


```{r}
library(tidyverse)
library(ggbeeswarm)
library(ggfortify)
```


To analyze the isoform data downloaded from Firebrowse (http://firebrowse.org/), we need to also download and load the correct NSD3 UCSC isoform ids used for expression analysis by Firebrowse, as they likely differ from the most current. I was able to download the correct annotations from NIH NCI Genomic Data Commons (https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files).  


For reference, I've included a schematic representation of NSD3 isoforms here.

![Schematic Representation of NSD3 Isoforms. Intron lengths have been truncated where indicated with a double-slash.](./NSD3_isoforms.png)


```{r}
# The following code will filter the annotation file and pull all UCSC isoform ids 
# linked to NSD3 (WHSC1L1).

NSD3.isoforms <- 
  read_tsv(file = "TCGA.hg19.June2011.gaf") %>%   
  filter(FeatureDBSource == "UCSCgene" &  
         str_detect(Gene, "WHSC1L1")   &  
         CompositeType == "gene") %>%  
  pull(FeatureID)
```

  
Next, normalized isoform expression data is loaded into R and filtered by NSD3 isoform ids identified in the previous step. I've downloaded this data from the Firebrowse website. File name: illuminahiseq_rnaseqv2-RSEM_isoforms_normalized  (MD5) (http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/LUSC/20160128/gdac.broadinstitute.org_LUSC.Merge_rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_isoforms_normalized__data.Level_3.2016012800.0.0.tar.gz)


```{r}

# Read in isoform expression data, filter for NSD3 isoforms, and prepare for plotting 
# with ggplot.

iso <- read_tsv("LUSC.isoform.txt") %>%  
        filter(`Hybridization REF` %in% NSD3.isoforms) %>%  
        gather(Patient, RSEM, 2:553) %>%  
        mutate(Isoform = `Hybridization REF`) %>%  
        select(Isoform, Patient, RSEM) %>%  
        mutate(RSEM = as.numeric(RSEM))
```


The first plot we will look at is a box and point plot of normalized RSEM values for each NSD3 isoform. Because of the large variance in the data, I also plotted with a log transformed y-axis.


```{r, fig.width=4,fig.height=3, fig.cap="Boxplot of isoform expression from TCGA-LUSC. We have cut off the y-axis at 1500 RSEM for easier visualization. This will remove highly expressed NSD3 short outliers."}  

# Boxplot of isoform expression from TCGA-LUSC. I've cut off the y-axis at 1500 RSEM 
# for easier visualization, removing highly expressed NSD3 short outliers (n=174).

iso %>% 
    ggplot(aes(Isoform, RSEM, fill = Isoform)) +  
    geom_boxplot() +  
    geom_quasirandom(pch= 21, size = 0.05, alpha =0.2) +  
    scale_y_continuous("NSD3 Isoform Expression (RSEM)",  
                       limits = c(0,1500)) +  
    theme_classic() +  
    theme(axis.text.x = element_blank()) +  
    xlab("")

```
\

Next, the same data is plotted with a log transformed y-axis.


```{r, fig.width=4, fig.height=3, fig.cap="Boxplot of isoform expression from TCGA-LUSC. We have cut off the y-axis at 1500 RSEM for easier visualization. This will remove highly expressed NSD3 short outliers."}

# Boxplot with log transformed y axis to visualize all data points.

iso %>% 
    ggplot(aes(Isoform, RSEM, fill = Isoform)) +  
    geom_boxplot() +  
    geom_quasirandom(pch= 21, size = 0.01, alpha =0.2) +  
    scale_y_continuous("NSD3 Isoform Expression (RSEM)",  
                       labels = scales::comma, trans = "log10") +  
    theme_classic() +  
    theme(axis.text.x = element_blank()) +  
    xlab("")

```
\


Next, I plotted the expression of short versus long isoforms to determine if there is any trend. Because there are three long isoforms, which seem to have similar expression levels, I grouped all three for the below comparison. First, linear regression is used to model the data. I found a moderate positive correlation between long and short expression with an R-squared value of 0.492. Plotting several measures of the linear regression model shows that a linear model is generally acceptable to map the relationship with few outliers. 


```{r, fig.width=6, fig.height=6, fig.cap = "Several plots to access whether a linear regression model is appropriate to model the data"}

# Perfrom linear regression analysis on Short vs Long isoform expression analysis. 
# Long isoform expression values have been grouped. 

fit <- iso %>% spread(Isoform, RSEM) %>%  
  mutate(Long = uc003xli.2 + uc010lwe.2 + uc011lbm.1) %>%  
  lm(uc003xlj.2~Long, .)  

summary(fit)  

autoplot(fit)  

```
\

```{r, fig.width=4, fig.height=3, fig.cap="Scattter plot of NSD3 Short Vs Long Expression. R-squared value is shown as measure of linear regression model fit"}
# Scatter plot including R-squared value.

iso %>% spread(Isoform, RSEM) %>%   
  mutate(Long = uc003xli.2 + uc010lwe.2 + uc011lbm.1) %>%   
  ggplot(aes(Long, uc003xlj.2)) +  
  geom_point() +  
  geom_smooth(method='lm',formula=y~x) +  
  annotate("text", label = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5)),  
           x = 12000, y =35000) +  
  scale_x_log10("RSEM - NSD3 Long", limits = c(100, 50000)) +  
  scale_y_log10("RSEM - NSD3 Short", limits = c(100, 50000)) +  
  theme_classic() 

```
\

Lastly, as another visualization I've plotted the distribution of the ratio between long and short NSD3 expression, which we can see is relatively uniform with a positive peak, indicating a similar ratio across samples with higher short form expression relative to the long isoforms.


```{r, fig.width=4, fig.height=3, fig.cap="Distribution of the Ratio of Long vs Short Isoform Expression"}

# Distribution of Long vs Short Isoform Ratio

iso %>% spread(Isoform, RSEM) %>%   
  mutate(Long = uc003xli.2 + uc010lwe.2 + uc011lbm.1) %>%   
  mutate(Ratio = uc003xlj.2 / Long) %>%   
  ggplot(aes(Ratio)) +  
  geom_density(fill = "grey") +  
  theme_classic() +  
  xlab("Ratio of NSD3-Short to NSD3-Long") +  
  ylab("Count") +  
  scale_x_continuous(limits= c(-10, 70))  
  
```
\


\pagebreak

# Observations

Here I observed that expression of long NSD3 isoforms vs short isoform is moderately correlated in the context of the TCGA-LUSC data-set. This observation suggests that the isoforms are likely co-regulated, as one increases so does the other. The analysis will be helpful in framing shared and separate functions of the two isoforms. 
  
# Acknowledgement

The results shown here are in whole or part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov


# References

1. Shen C, Ipsaro JJ, Shi J, et al. NSD3-short is an adaptor protein that couples BRD4 to the CHD8 chromatin remodeler. Molecular cell. 2015;60(6):847-859. doi:10.1016/j.molcel.2015.10.033.  
2. Weinstein JN, Collisson EA, Mills GB, et al. The Cancer Genome Atlas Pan-Cancer Analysis Project. Nature genetics. 2013;45(10):1113-1120. doi:10.1038/ng.2764.

****
ExpID-028

