deseq.R

#!/bin/env Rscript

library(readr)
library(DESeq2)
library(BiocParallel)
library(dplyr)
library(stringr)
library(tidyr)
library(cowplot)
library(ggrepel)
register(MulticoreParam(24))

datakey <- read_csv('datakey.csv')
df <- read_tsv('repCounts.tsv')  %>%
        setNames(str_replace(names(.),'.bam',''))
countMat <- df %>%
        select(grep('SAE',names(.)))
colData <- data_frame(sample = names(countMat)) %>%
        inner_join(datakey) %>%
        filter(treatment == 'non-frag') %>%
        mutate(Grade = factor(Grade)) %>%
        mutate(Grade = relevel(Grade, ref='L')) %>%
        data.frame()
rownames(colData) = colData$sample
deseqRes <- countMat %>%
        subset(select= names(.) %in% colData$sample) %>%
        DESeqDataSetFromMatrix(countData=.,
                                                colData = colData,
                                                design = ~Grade) %>%
        DESeq(parallel=T) %>%
        results(parallel=T) %>%
        data.frame %>%
        cbind(df[,1:4],.) %>%
        filter(baseMean!=0) %>%
        replace_na(list(pvalue=1,padj=1)) %>%
        mutate(signif = ifelse(padj < 0.05,'Significant','Not significant')) %>%
        tbl_df

p <- ggplot(data=deseqRes, aes(x= baseMean,y = log2FoldChange, color = signif)) +
        geom_point() +
        geom_text_repel(data=deseqRes %>% filter(signif == 'Significant'), aes(label = id)) +
        scale_color_manual(values=c('black','red')) +
        scale_x_log10()
ggsave(p,file='volcano.pdf')