---
title: "Wapl Nipbl"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: embed
    
---

```{r setup, include=FALSE}

library(tidyverse)
library(readxl)
library(ggplot2)
library(plotly)
library(tableone)
knitr::opts_chunk$set(message = FALSE)

```

Weight Data {vertical_layout=fill data-icon="fa-archive" data-orientation=rows}
===================================== 

Column {data-width=10000}
------------------------------------------------------------------------------

### A look at the data

```{r raw data}

R_Wapl_and_Nipbl <- read_excel("R_Wapl_and_Nipbl.xlsx", 
                               sheet = "R") %>% mutate(Genotype = as.factor(Genotype))

colnames(R_Wapl_and_Nipbl)[colnames(R_Wapl_and_Nipbl)=="Embryo-Placenta"] <- "Embryo_Placenta"

#AD Look at tableone::createTableOne for some better options than `summary`

vars <- c("Mother", "Embryo", "Placenta", "Genotype", "Sex", "Embryo_Placenta")
tableOne <- CreateTableOne(vars = vars, strata = c("Genotype"), data = R_Wapl_and_Nipbl)
tableOne


```

Statistics {vertical_layout=scroll data-icon="fa-asterisk"}
===================================== 

Row {.tabset .tabset-fade data-height=10000}
------------------------------------------------------------------------------

### Normality Tests

***
Embryo Data
```{r embryo}

shapiro.test(R_Wapl_and_Nipbl$Embryo)

```

***
Placenta Data
```{r Placenta}
shapiro.test(R_Wapl_and_Nipbl$Placenta)

```

***
Embryo/Placenta Data
```{r EP}
shapiro.test(R_Wapl_and_Nipbl$Embryo_Placenta)

```

***

All variables passed normality test


### T-test

*** 
WAPL vs WT Genotype Comparison
```{r ttest}

t.test(Embryo-Placenta ~ Genotype, data = R_Wapl_and_Nipbl, subset = Genotype %in% c("WT", "WAPL"))


```


### One-way Anova

***
Embryo/Placenta vs Genotype Anova
```{r Anova}

anova <- aov(Embryo_Placenta ~ Genotype, data = R_Wapl_and_Nipbl)
summary(anova)

```

***
Post-hoc Testing
```{r posthoc a}
TukeyHSD(anova)

```
***
Overall Anova is significant indicating difference in genotypes 

We also can look at individual relationships using the posthoc testing

### Two-way Anova

***
Embryo/Placenta vs Genotype * Mother Two-way Anova
```{r 2Anova}

litter_genotype_anova <- aov(Embryo_Placenta ~ Genotype * Mother, data = R_Wapl_and_Nipbl) #star tests if these two categoricals interact and affect outcome
summary(litter_genotype_anova) #output will show what categorical variables have an affect on embryo_placenta weight
TukeyHSD(litter_genotype_anova) #posthoc testing for normally distruted data

```
***
We see that there is significant effects of litter on embryo/placenta and an effect of genotype on embryo/placenta but there is no interaction between the litter and gemotype. Thus, each of these effects on the data is not caused by their relationship.

Visualizations {data-orientation=columns data-icon="fa-chart-bar"}
===================================== 
Row {.tabset data-height=400}
------------------------------------------------------------------------------

### Visual Data 

```{r Violin Plot}

jn_theme <- function(){
  theme_bw() +
    theme(axis.text = element_text(size = 8, color = "Black"),
          axis.title = element_text(size = 14),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size=10),
          strip.background = element_blank(),
          plot.title = element_text(size = 20, hjust = 0.5),
          panel.grid.major.x = element_blank())
}

ggplot(R_Wapl_and_Nipbl, aes(x=Genotype, y=Embryo_Placenta)) + 
  geom_violin() + 
  geom_point(size = .5, height = 0, width = 0.05) + 
  stat_summary(fun.y = mean, geom="point", shape = 18, size  = 3, color = "red") + 
  stat_summary(fun.y = median,  geom = "point", shape = 3, size = 3, color = "blue") +
  jn_theme() + 
  scale_x_discrete(limits = c("WT", "WAPL", "NIPBL", "COMPOUND", "WAPL_NIPBL_FLOX")) 
  


```



### Interactive

```{r interactive}

jn_theme <- function(){
  theme_bw() +
    theme(axis.text = element_text(size = 14, color = "Black"),
          axis.title = element_text(size = 16),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size=14),
          strip.background = element_blank(),
          plot.title = element_text(size = 20, hjust = 0.5),
          panel.grid.major.x = element_blank())
}

plot1 <- ggplot(R_Wapl_and_Nipbl, aes(x=Genotype, y=Embryo_Placenta, label = Mother)) + 
  geom_violin() + 
  geom_point(size = .5, height = 0, width = 0.05) + 
  stat_summary(fun.y = mean, geom="point", shape = 18, size  = 3, color = "red") + 
  stat_summary(fun.y = median,  geom = "point", shape = 3, size = 3, color = "blue") +
  jn_theme() + 
  scale_x_discrete(limits = c("WT", "WAPL", "NIPBL", "COMPOUND", "WAPL_NIPBL_FLOX")) 

plot1





library(plotly)

ggplotly(plot1, tooltip = c("Embryo_Placenta", "Mother"))%>%layout(violinmode = 'group', violingap = 1, violingroupgap = 1)

```


RNA-Seq {data-orientation=columns data-icon="fa-sitemap"}
===================================== 
Row {.tabset data-height=400}
------------------------------------------------------------------------------

### WT vs WAPL

```{r RNA-Seq}

#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
 #      dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
#       dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_WAPL","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_WAPL","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
waplvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
waplvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")




```

### WT vs Wapl No Outlier

```{r no_out1}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
 #      dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
 #      dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
 #      dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_no_outlier.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_WAPL","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_WAPL","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
waplvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
waplvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")



```


### WT vs NIPBL

```{r RNA Seq2}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
#       dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
#       dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_NIPBL","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_NIPBL","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
nipblvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
nipblvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")



```

### WT vs Nipbl No Outlier

```{r no_out2}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
#       dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
#       dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_no_outlier.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_NIPBL","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_NIPBL","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
nipblvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
nipblvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")



```

### WT vs Compound

```{r RNA Seq3}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
#       dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
#       dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_COMPOUND","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_COMPOUND","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
compoundvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
compoundvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")



```

### WT vs Compound No Outlier

```{r no_out3}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
 #      dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
 #      dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_no_outlier.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_COMPOUND","WT")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_COMPOUND","WT"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
compoundvwt_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
compoundvwt_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")


```

### NIPBL vs WAPL

```{r RNA Seq4}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
 #      dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
 #      dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
 #      dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_NIPBL","DD_WAPL")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_NIPBL","DD_WAPL"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
nipblvwapl_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
nipblvwapl_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")


```


### Nipbl vs Compound

```{r RNA Seq6}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
#       dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
#       dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
 #      dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_NIPBL","DD_COMPOUND")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_NIPBL","DD_COMPOUND"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
nipblvcompound_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
nipblvcompound_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")


```



### Wapl vs Compound

```{r RNA Seq7}
#Jacob_Noeker 
#RNA Seq Analysis
#WAPL NIPBL Project


#install.packages("BiocManager")
library(DESeq2)


# Set up variables, output folders, etc
#
# This section does the following:

# 1. Checks for output folders and creates them if absent
# 2. Defines input read count files
# 3. Strips annotation columns and saves annotation info
# 4. Defines FDR threshold

# output directory for deseq2 files 
# create these if not present

deseq.folder <- "output/deseq2/"
norm.folder <- "output/vstnorm/"
degenes.folder <- "output/degenes/"

# check for subfolders and if not present, create folder
#ifelse(!dir.exists(file.path('.',deseq.folder)), 
#       dir.create(file.path('.',deseq.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',norm.folder)), 
 #      dir.create(file.path('.',norm.folder), recursive = TRUE), NA)
#ifelse(!dir.exists(file.path('.',degenes.folder)), 
 #      dir.create(file.path('.',degenes.folder), recursive = TRUE), NA)

# hardcoded column names from featureCounts output
# only 'EnsID' column name is dependent on GTF file
annot.cols <- c("Chr","Start","End","Strand","Length")
geneid.colname <- "Geneid"

# Define FDR threshold and FDR column name
fdr.thres <- 0.1

# Define FC threshold
fc.thres <- 0

# Define limits for MA and volcano plots
#fc.lim <- c(-2, 2)

# define size of pdf for MA and volcano plots (in inches)
#ma.width <- 10
#ma.height <- 8

# define size of PDF file for degPattern plots (in inches)
#pdf.width <- 20
#pdf.height <- 20


# read count tables
# count tables from featureCounts
# 1. The first line with the commant is commented out and not read
# 2. the second line wth the file names can be parsed within R using 
#    grep() to match specified sample names
counts <- read.table("counts_CT_simp.txt",
                     sep="\t",header=TRUE, row.names=1)
exptname <- 'wapl_nipbl'

# strip annotation columns but save annot data
annot.data <- cbind(rownames(counts), counts[,colnames(counts) %in% annot.cols])
colnames(annot.data)[1] <- geneid.colname

counts <- counts[,!colnames(counts) %in% annot.cols]

# sample list
# converted to input variable
all.samples <- c("DD_WAPL","DD_COMPOUND")

# reformat column names matching sample names
# TODO: throw error if sample name does not match column names
sample.vec <- NULL
for(i in 1:length(all.samples)){
  # which column names match the sample name
  idx <- grep(all.samples[i], colnames(counts))
  
  # create named vector to save metadata
  # this maps column names to sample names
  temp.vec <- rep(all.samples[i], length(idx))
  names(temp.vec) <- colnames(counts)[idx]
  sample.vec <- c(sample.vec, temp.vec)
}


# DESeq2 analysis
#
# Input: read counts in TSV format
# Output: differential expression results from DESeq2 in TSV format
#
# This section takes the cleaned count files and performs the following steps:
#
# 1. creates DESeq data object
# 2. creates sample lists for individual comparisons
# 3. performs DESeq2 comparisons between samples separately for plus & minus BMP2
#
# Also, obtains a list of genes to be used as background for downstream analysis
#

# create a list for read counts
all.reads <- list(counts)
names(all.reads) <- exptname

# create sample lists for individual comparisons
# create input variable with this info
all.comps <- list(a = c("DD_WAPL","DD_COMPOUND"))

# list of genes to be used as background
background.genes <- NULL

# list of DE results
res.list <- list()

# list of DE genes
deg.list <- list()

# list of upregulated genes
up.list <- list()

for(k in 1:length(all.reads)){
  read.mat <- all.reads[[k]]
  for(i in 1:length(all.comps)){
    # get treatment and control names
    treatment <- all.comps[[i]][1]
    control <- all.comps[[i]][2]
    
    # get column indices corresponding to samples in comparison
    sample.idx <- NULL
    for(j in all.comps[[i]]){
      sample.idx <- c(sample.idx, grep(j, colnames(read.mat)))
    }
    
    # extract reads for specific sample columns
    reads.sub <- read.mat[,sample.idx]
    
    # create design matrix by selection from sample.vec
    idx <- which(sample.vec %in% all.comps[[i]])
    design <- data.frame(row.names = names(sample.vec[idx]),
                         condition = sample.vec[idx])
    
    design$sampleid <- sub("_.+", "", rownames(design))
    
    # DESeq2 object
    dds <- DESeqDataSetFromMatrix(countData = reads.sub,
                                  colData = design,
                                  design = ~condition)                                         
    
    # run DESeq2
    dds <- DESeq(dds, betaPrior = FALSE)
    #res <- results(dds, contrast=c("condition", treatment, control),
     #              cooksCutoff=FALSE)
    
    # get results with lfcShrink
    res <- lfcShrink(dds, contrast=c("condition", treatment, control))
    res.list[[length(res.list) + 1]] <- res
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DESeq2 results to file
    write.table(cbind(geneid=rownames(res),res),
                paste0(deseq.folder, names(all.reads)[k], "_", treatment, "-vs-", control,".txt"),
                sep="\t",quote=FALSE, row.names=FALSE)
    
    # keep list of degenes
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   abs(res$log2FoldChange) >= fc.thres)
    deg.list[[length(deg.list) + 1]] <- rownames(res)[idx]
    names(deg.list)[length(deg.list)] <- names(all.comps)[i]
    
    # output DE genes to file
    write.table(deg.list[[i]], 
                paste0(degenes.folder, names(all.reads)[k], "_fdr", fdr.thres, "_fc", fc.thres, ".txt"),
                sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
    
    # list of upregulated genes (for easy reference)
    idx <- which(res$padj < fdr.thres & 
                   !is.na(res$padj) & 
                   res$log2FoldChange >= fc.thres)
    up.list[[length(up.list) + 1]] <- rownames(res)[idx]
    names(up.list)[length(up.list)] <- names(all.comps)[i]
    
    # get normalized data
    norm.vst <- vst(dds)
    
    # write normalized data to file
    write.table(cbind(geneid=rownames(assay(norm.vst), assay(norm.vst))),
                paste0(norm.folder, "vstnorm_", treatment, "-", control,".txt"),
                sep="\t", quote=FALSE, row.names=FALSE)
    
  }
  # collect genes that have baseMean != 0 in at least one comparison
  background.genes <- c(background.genes, rownames(res)[res$baseMean != 0])
}

# keep unique bg genes and sort by gene names
background.genes <- unique(background.genes)
background.genes <- background.genes[order(background.genes)]

# extract annotation data for bg genes and output Ensembl ID and gene symbols
background.annot <- annot.data[annot.data[,geneid.colname] %in% background.genes, 
                               geneid.colname]

# write background list to file
write.table(background.annot, paste0(deseq.folder, "background-genes_", ".txt"), 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)
waplvcompound_pca <- plotPCA(norm.vst, intgroup="condition")
library(ggrepel)
waplvcompound_pca + geom_text_repel(aes_string(x = "PC1", y = "PC2", label = "name"), color = "black")



```



About {data-orientation=rows data-icon="fa-comment-alt"}
===================================== 

Row {data-height=1000}
------------------------------------------------------------------------------

Dashboard made by Jacob Noeker