RNAi and ku217 DEseq2 and correlation to ChIP peaks

First step is to load DESeq2 and the necessary packages

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
library(DESeq2)

OK great had to do a few things to get it to load: Do not compile from source (when it asks say n). I think this fixed problems.

Next step: I need to load the data. In this case I am loading the HTSeq count tables that I downloaded from trimming, mapping and counting on galaxy. First I have to load the sample table so that the object knows the conditions, then I can make the HTSeq count object.

sample_table_RNAi_ku217 = read.table("sample_table_RNAi_ku217.txt", header = TRUE, sep = "\t")
HTSeq_RNAi_ku217 <- DESeqDataSetFromHTSeqCount(sampleTable = sample_table_RNAi_ku217, directory = "HTseq_count_galaxy/", design = ~ Condition)

Next step: I wanter filter reads with low counts. I am going to filter out those reads that have less than 2 for total FPKM across all conditions. This is a proxy for saying that across biological replicates, there should be close to 1 FPKM in each condition. This is slightly more stringent than I did before, but I think it is good.

HTSeq_RNAi_ku217 <- HTSeq_RNAi_ku217[rowSums(counts(HTSeq_RNAi_ku217))>2]

Next for the intial analyses, I want to perform the read normalization and variance stabilization. Just so that I can take a look at the data. I am concerned right now it did not read in correctly

HTSeq_RNAi_ku217_vsd <- vst(HTSeq_RNAi_ku217, blind = FALSE)
HTSeq_RNAi_ku217_rlog <- rlog(HTSeq_RNAi_ku217, blind = FALSE)
head(assay(HTSeq_RNAi_ku217_vsd))
head(assay(HTSeq_RNAi_ku217_rlog))

OK I think it looks ok… itersting the values are so different depending on the assay type. Now I want to take a look at the PCA analysis.

plotPCA(HTSeq_RNAi_ku217_vsd, intgroup=c("Condition", "Replicate"))
plotPCA(HTSeq_RNAi_ku217_rlog, intgroup=c("Condition", "Replicate"))

OK few pieces here… there is still more varaince across the WTs, which is annoying, but PC2 clearly separates along the condition, so that is good. Also, looks like the file loaded correctly.

HTSeq_RNAi_ku217 <- DESeq(HTSeq_RNAi_ku217)
RNAi_N2_res <- results(HTSeq_RNAi_ku217, contrast = c("Condition", "nhr-25(RNAi)", "N2_L4440"))
ku217_N2_res <- results(HTSeq_RNAi_ku217, contrast = c("Condition", "nhr-25(ku217)", "N2_L4440"))
ku217_RNAi_res <- results(HTSeq_RNAi_ku217, contrast = c("Condition", "nhr-25(ku217)","nhr-25(RNAi)" ))
RNAi_N2_res
ku217_N2_res
ku217_RNAi_res

summary(RNAi_N2_res)
summary(ku217_N2_res)
summary(ku217_RNAi_res)

plotCounts(dds = HTSeq_RNAi_ku217, gene = "lin-14", intgroup = c("Condition"))

plotMA(RNAi_N2_res)
plotMA(ku217_N2_res)
plotMA(ku217_RNAi_res)

Ordering the results tables by the smallest pvalue, then subset for only those padj values that are less than .05, I could likely even go as low as .01, but I will start with .05.I am going to write the tables to files, then I can use these in a python analysis.

RNAi_N2_res_ordered <- RNAi_N2_res[order(RNAi_N2_res$pvalue),]
ku217_N2_res_ordered <- ku217_N2_res[order(ku217_N2_res$pvalue),]
ku217_RNAi_res_ordered <- ku217_RNAi_res[order(ku217_RNAi_res$pvalue),]
RNAi_N2_res_sig <- subset(RNAi_N2_res_ordered, padj < .05)
ku217_N2_res_sig <- subset(ku217_N2_res_ordered, padj < .05)
ku217_RNAi_res_sig <- subset(ku217_RNAi_res_ordered, padj < .05)
write.table(as.data.frame(RNAi_N2_res_sig), file="RNAi_N2_dea.txt", sep = "\t")
write.table(as.data.frame(ku217_N2_res_sig), file="ku217_N2_dea.txt", sep = "\t")
write.table(as.data.frame(ku217_RNAi_res_sig), file="ku217_RNAi_dea.txt", sep = "\t")

I wrote a python script that correlates the chip peaks with diff reg geens in the rna-seq against N2. I am going to read that file in here and then do some visualization.

peak_to_gene <- read.table('peak_to_drg.txt', header = TRUE, sep = '\t')

Try and plot all the distances as a histogram

hist(peak_to_gene$Peak.to.DRG.distance, breaks = 50, main = 'ChIP peak to diff reg gene', xlab = 'Peak to gene distance (basepairs)', xlim = c(-60000, 60000), col = '#2e6bf1')

Now I want to make a violin plot sep, by the types of peaks and then hopefully the types of chromatin. Question: do the distance change based on peak features?

peak_to_gene$motif.class <- as.factor(peak_to_gene$motif.class)
library(ggplot2)
violin <- ggplot(peak_to_gene, aes(x=peak_to_gene$motif.class, y=peak_to_gene$Peak.to.DRG.distance, fill = peak_to_gene$motif.class)) +  geom_violin(trim = FALSE) 
violin + stat_summary(fun.data=mean_sdl, mult = 1, geom="pointrange", width=.1)

peak_to_gene$chromatin.class <- as.factor(peak_to_gene$chromatin.class)
violin <- ggplot(peak_to_gene, aes(x=peak_to_gene$chromatin.class, y=peak_to_gene$Peak.to.DRG.distance, fill = peak_to_gene$chromatin.class)) +  geom_violin(trim = FALSE) 
violin + stat_summary(fun.data=mean_sdl, mult = 1, geom="pointrange", width=.1)