.Rhistory

10000/4
50*8
15163/(17000+22200)
15163/((17000+22200)-(15163))
counts <- table(mtcars$gear)
barplot(counts, main="Car Distribution",
xlab="Number of Gears")
counts
class(counts)
mtcars$gear
-log10(0.99)
-log10(0.99)*-10
-10*log10(0.99)
log10(0.99)*10
5000*0.2
66000*1.0275
67815*1.15
67815*1.15*1.15
67815*1.3
67815*1.12
66000*0.7
66000*0.65
31+43
74-24
24/31
24/43
24/34
24/36
19.26
19/26
38/26
19/38
22/31
22/34
30/43
20/34
26/36
18/26
27/38
18/31
25/34
66*1.4
66*1.5
930/58
58*8
58*10
1400/58
14*58
15*58
60/5
300/60
library(Vennerable)
diff_dir <- "diff-binding/"
cont_vs_24 <- read.csv(paste0(diff_dir, "diff-binding-res-annot-FDR-5%-24hr.csv"), stringsAsFactors = F)
cont_vs_48 <- read.csv(paste0(diff_dir, "diff-binding-res-annot-FDR-5%-48hr.csv"), stringsAsFactors = F)
cont_vs_72 <- read.csv(paste0(diff_dir, "diff-binding-res-annot-FDR-5%-72hr.csv"), stringsAsFactors = F)
cont_vs_24 <- read.csv(paste0("~/Desktop/", "diff-binding-res-annot-FDR-5%-24hr.csv"), stringsAsFactors = F)
cont_vs_48 <- read.csv(paste0("~/Desktop/", "diff-binding-res-annot-FDR-5%-48hr.csv"), stringsAsFactors = F)
cont_vs_72 <- read.csv(paste0("~/Desktop/", "diff-binding-res-annot-FDR-5%-72hr.csv"), stringsAsFactors = F)
list1 <- list(cont_vs_24$peak_ID, cont_vs_48$peak_ID, cont_vs_72$peak_ID)
names(list1) <- c("Cont vs 24", "Cont vs 48", "Cont vs 72")
vlist1 <- Venn(list1)
vlist1 <- vlist1[, c("Cont vs 24", "Cont vs 48", "Cont vs 72")]
ppi=300
png(paste0(diff_dir, "venn_diagram_diff_bound.png"), width=7*ppi, height=7*ppi, res=ppi)
plot(vlist1, doWeights = TRUE)
dev.off()
66000*2.75
66000*0.0275
67815*0.15
67815+10172.25
80000-67815
12185/67815
11/5.530089
11/0.468566
58*3
1400/3
1400/58
80000/12
80000/12*0.86
(80000/12*)0.86
(80000/12)*0.86
500000*4000
7702477/12307899
100000/5e6
1/62.581500
3*58
140/58
1400/58
58*24
58*27
1400+175
64/10
64/16
4*10
8*56
setwd('~/Documents/GitHub/RNA-seq-Differential-Expression-workshop-June-2021/'),
dds <- readRDS("DESeq2.rdata")
setwd('~/Documents/GitHub/RNA-seq-Differential-Expression-workshop-June-2021/')
dds <- readRDS("DESeq2.rdata")
setwd('/Users/omw/Documents/GitHub/Bulk-RNA-seq_workshop_Part-2_June-2021-/')
dds <- readRDS("DESeq2.rdata")
list.files()
list.files("data")
dds <- readRDS("data/DESeq2.rdata")
dds <- readRDS("data/DESeq2.rdata")
library(DESeq2)
# Read in the count data:
cts <- as.matrix(read.table("Day-2/all_counts.txt",
sep="\t",
header = TRUE,
row.names=1,
stringsAsFactors = F))
list.files()
# Read in the count data:
cts <- as.matrix(read.table("data/all_counts.txt",
sep="\t",
header = TRUE,
row.names=1,
stringsAsFactors = F))
# Read in the metadata:
colData <- read.csv("data/sample_metadata.csv", row.names=1)
# Construct a DESeq2 dataset from the raw counts, the metadata, and the desired design variable to be tested for differential expression.
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = colData,
design = ~ group)
group
colData
# Construct a DESeq2 dataset from the raw counts, the metadata, and the desired design variable to be tested for differential expression.
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = colData,
design = ~ group)
# Apply the DESeq2 analysis pipeline:
dds <- DESeq(dds)
# extract results
res <- results(dds,
name = "group_Dex_vs_untreated",
alpha = 0.05,
lfcThreshold = 0)
resultsNames(dds)
colData$tx.group <- factor(colData$tx.group, levels=c("untreated", "Dex", "Alb", "Alb_Dex"))
# Construct a DESeq2 dataset from the raw counts, the metadata, and the desired design variable to be tested for differential expression.
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = colData,
design = ~ tx.group)
# Apply the DESeq2 analysis pipeline:
dds <- DESeq(dds)
# extract results
res <- results(dds,
name = "group_Dex_vs_untreated",
alpha = 0.05,
lfcThreshold = 0)
resultsNames()
resultsNames(dds)
# extract results
res <- results(dds,
name = "tx.group_Dex_vs_untreated",
alpha = 0.05,
lfcThreshold = 0)
saveRDS(dds, file = "data/DESeq2.rds")
dds <- readRDS("data/DESeq2.rdata")
dds <- readRDS("data/DESeq2.rds")
# read in the flat file we downloaded and have a look at it
anno <- read.delim("data/GRCh38.p12_ensembl-97.txt", stringsAsFactors = T, header = T)
anno <- anno[order(anno$Chromosome.scaffold.name),]
dim(anno)
# have a look at the first few rows
head(anno)
tab1 <- table(anno$Chromosome.scaffold.name)
tab1[1:22]
any(duplicated(anno$Gene.stable.ID))
# use match() to find corresponding indicies (rows) for each ENSG ID
mat1 <- match(rownames(res_ord), anno$Gene.stable.ID)
table(is.na(mat1))
# order by adj Pval
res_ord <- res[order(res$padj),]
# use match() to find corresponding indicies (rows) for each ENSG ID
mat1 <- match(rownames(res_ord), anno$Gene.stable.ID)
table(is.na(mat1))
# add gene names to results as a new column
res_ord$gene <- as.character(anno$Gene.name[mat1])
head(res_ord, 20)
res_ord$chr <- as.character(anno$Chromosome.scaffold.name[mat1])
res_ord$start <- as.character(anno$Gene.start..bp.[mat1])
res_ord$end <- as.character(anno$Gene.end..bp.[mat1])
res_ord$strand <- as.character(anno$Strand[mat1])
### volcano plot
plot(res$log2FoldChange, -log10(res$pvalue),
main = "Volcano plot",
las = 1, col = "indianred",
ylab = "- log10 P-value", xlab = "log2 Fold change")
# add horizontal lines to help guide interpretation
abline(h=-log10(0.05/nrow(res)), lty = 1, col = "black") # Bonferonni
abline(h=-log10(0.05), lty = 2, col = "black") # nominal P-value
## beautified volcano plot
# save a dataframe from the results() output
res_tmp <- as.data.frame(res_ord)
# add a column that will be used to save the colors we want to plot
res_tmp$cols <- c()
# set the significance cut off (alpha) and fold change threshold to be used for coloring of genes
alpha <- 0.05/nrow(res)
fc_cutoff <- 2
# loop through our dataframe and add values to the color column based on magnitude of alpha and LFCs
res_tmp$cols <- NA
for(i in 1:nrow(res_tmp)){
# if pvalue is NA don't assign a color - no plotting
if(is.na(res_tmp$pvalue[i])){
res_tmp$cols[i] <- NA
}
# if pvalue is < alpha AND LFC is > FC cutoff color red
else if(res_tmp$pvalue[i]<=alpha & res_tmp$log2FoldChange[i] > fc_cutoff){
res_tmp$cols[i] <- "indianred"
}
# if pvalue is < alpha AND LFC is < -FC cutoff color red
else if(res_tmp$pvalue[i]<=alpha & res_tmp$log2FoldChange[i] < -fc_cutoff){
res_tmp$cols[i] <- "indianred"
}
# if pvalue is < alpha AND LFC is not within cut off value color blue
else if(res_tmp$pvalue[i]<=alpha & res_tmp$log2FoldChange[i]>-fc_cutoff & res_tmp$log2FoldChange[i]<fc_cutoff){
res_tmp$cols[i] <- "cornflowerblue"
}
# if pvalue is > alpha AND LFC is > cut off value color gray
else if(res_tmp$pvalue[i]>alpha & res_tmp$log2FoldChange[i] > fc_cutoff){
res_tmp$cols[i] <- "gray47"
}
# if pvalue is > alpha and LFC is < -cut off value color gray
else if(res_tmp$pvalue[i]>alpha & res_tmp$log2FoldChange[i] < -fc_cutoff){
res_tmp$cols[i] <- "gray47"
}
# if pvalue is > alpha and LFC is not within cutoff values color light gray
else if(res_tmp$pvalue[i]>alpha & res_tmp$log2FoldChange[i] < fc_cutoff){
res_tmp$cols[i] <- "gray10"
}
}
res_tmp$ENSG <- rownames(res_tmp)
# generate the plot
p = ggplot(res_tmp, aes(log2FoldChange, -log10(pvalue))) +
geom_point(aes(col=col), alpha = 0.5, size =2.5, colour = res_tmp$cols, fill = res_tmp$cols)  +
xlab("Log2 fold change") + ylab("-log10 Q-value") +
ylim(0, 9) +
xlim(-5, 11) +
geom_hline(yintercept = -log10(alpha), color = "black", linetype = "dashed", size = 0.4) +
# add vertical fold change lines
geom_vline(xintercept = fc_cutoff, colour = "black", linetype="dotted") +
geom_vline(xintercept = -fc_cutoff, colour = "black", linetype="dotted")
theme(legend.key = element_blank()) +
ggtitle("Control vs Dex")
# print the plot
print(p)
library(ggplot)
library(ggplot2)
library(vsn)
library(dplyr)
library(pheatmap)
library(gplots)
library(RColorBrewer)
library(ComplexHeatmap)
library(readr)
library(circlize)
library(EnhancedVolcano)
library(apeglm)
library(ggplot2)
# generate the plot
p = ggplot(res_tmp, aes(log2FoldChange, -log10(pvalue))) +
geom_point(aes(col=col), alpha = 0.5, size =2.5, colour = res_tmp$cols, fill = res_tmp$cols)  +
xlab("Log2 fold change") + ylab("-log10 Q-value") +
ylim(0, 9) +
xlim(-5, 11) +
geom_hline(yintercept = -log10(alpha), color = "black", linetype = "dashed", size = 0.4) +
# add vertical fold change lines
geom_vline(xintercept = fc_cutoff, colour = "black", linetype="dotted") +
geom_vline(xintercept = -fc_cutoff, colour = "black", linetype="dotted")
theme(legend.key = element_blank()) +
ggtitle("Control vs Dex")
p = ggplot(res_tmp, aes(log2FoldChange, -log10(pvalue))) +
geom_point(aes(col=col), alpha = 0.5, size =2.5, colour = res_tmp$cols, fill = res_tmp$cols)  +
xlab("Log2 fold change") + ylab("-log10 Q-value") +
ylim(0, 9) +
xlim(-5, 11) +
geom_hline(yintercept = -log10(alpha), color = "black", linetype = "dashed", size = 0.4) +
# add vertical fold change lines
geom_vline(xintercept = fc_cutoff, colour = "black", linetype="dotted") +
geom_vline(xintercept = -fc_cutoff, colour = "black", linetype="dotted")
# print the plot
print(p)
p2 <- p +
# add labels to genes w/ LFC > 2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange > 2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = 0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# add labels to genes w/ LFC < -2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange < -2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = -0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# print the plot
print(p2)
library(ggrepel)
p2 <- p +
# add labels to genes w/ LFC > 2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange > 2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = 0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# add labels to genes w/ LFC < -2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange < -2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = -0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# print the plot
print(p2)
p2 <- p +
# add labels to genes w/ LFC > 2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange > 2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = 0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# add labels to genes w/ LFC < -2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange < -2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = -0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# print the plot
print(p2)
p2 <- p +
# add labels to genes w/ LFC > 2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange > 2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = 0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3) +
# add labels to genes w/ LFC < -2 and above alpha threshold
geom_label_repel(data = subset(res_tmp, log2FoldChange < -2 & pvalue < alpha), aes(label = gene),
box.padding   = 0.35,
nudge_x = -0.1,
nudge_y = 0.1,
point.padding = 1,
label.size = 0.1,
segment.size = 0.3,
segment.color = 'grey50', size = 3)
# print the plot
print(p2)
# subset @ 5% adjusted pval sig. level
res_order_FDR_05 <- res_ord[res_ord$padj<0.05,]
nrow(res_order_FDR_05)
res_ord <- res_ord[!is.na(res_ord$padj),]
# subset @ 5% adjusted pval sig. level
res_order_FDR_05 <- res_ord[res_ord$padj<0.05,]
nrow(res_order_FDR_05)
# write both to csv files
write.csv(as.data.frame(res_ord), file= "DE_results.csv")
write.csv(as.data.frame(res_order_FDR_05), file="DE_results.FDR.0.05.csv")
plotMA(res_ord, ylim=c(-6,6), main = "Raw Log2 Fold change")
# calculate shrunken fold change estimate
res_shrink <- lfcShrink(dds,
coef=paste0(resultsNames(dds)[which(resultsNames(dds)=="tx.group_Dex_vs_untreated")]),
type="apeglm")
BiocManager::install("apeglm")
library(apeglm)
plotMA(res_ord, ylim=c(-6,6), main = "Raw Log2 Fold change")
# calculate shrunken fold change estimate
res_shrink <- lfcShrink(dds,
coef=paste0(resultsNames(dds)[which(resultsNames(dds)=="tx.group_Dex_vs_untreated")]),
type="apeglm")
par(mfrow=c(2,1))
plotMA(res_ord, ylim=c(-6,6), main = "Raw Log2 Fold change")
plotMA(res_shrink, ylim=c(-6,6), main = "Shrunken Log2 Fold change")
rld <- rlog(dds, blind = FALSE)
ind_to_keep <- c(which(colData(rld)$group=="untreated"), which(colData(rld)$group=="Dex"))
# set up gene expression matrix
mat1 <- assay(rld)[rownames(res_order_FDR_05), ind_to_keep]
# scale matrix by each col. values
mat_scaled = t(apply(mat1, 1, scale))
# set up colors for heatmap
col = colorRamp2(c(-3, 0, 3), c("blue", "white", "red"))
cols1 <- brewer.pal(11, "Paired")
cols2 <- brewer.pal(9, "Greens")
# subset coldata for samples in untx and ex groups
colData_sub <- colData(dds)[ind_to_keep, ]
# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$group,
col = list(Group = c("untreated" = cols1[1], "Dex" = cols1[2])),
show_legend = TRUE)
# se up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
which="column", rot = 45,
gp = gpar(fontsize = 10)))
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE)
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")
mat_scaled=mat1
# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$group,
col = list(Group = c("untreated" = cols1[1], "Dex" = cols1[2])),
show_legend = TRUE)
# se up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
which="column", rot = 45,
gp = gpar(fontsize = 10)))
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE)
# plot the heatmap
par(mfrow=c(1,1))
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")
head(mat1)
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE, show_column_names = FALSE)
# plot the heatmap
par(mfrow=c(1,1))
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")
max(mat1)
# set up colors for heatmap
col = colorRamp2(c(0, 9, 18), c("blue", "white", "red"))
cols1 <- brewer.pal(11, "Paired")
cols2 <- brewer.pal(9, "Greens")
# subset coldata for samples in untx and ex groups
colData_sub <- colData(dds)[ind_to_keep, ]
# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$group,
col = list(Group = c("untreated" = cols1[1], "Dex" = cols1[2])),
show_legend = TRUE)
# se up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
which="column", rot = 45,
gp = gpar(fontsize = 10)))
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE, show_column_names = FALSE)
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")
# set up colors for heatmap
col = colorRamp2(c(0, 9, 18), c("white", "blue", "red"))
cols1 <- brewer.pal(11, "Paired")
cols2 <- brewer.pal(9, "Greens")
# subset coldata for samples in untx and ex groups
colData_sub <- colData(dds)[ind_to_keep, ]
# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$group,
col = list(Group = c("untreated" = cols1[1], "Dex" = cols1[2])),
show_legend = TRUE)
# se up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
which="column", rot = 45,
gp = gpar(fontsize = 10)))
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE, show_column_names = FALSE)
# plot the heatmap
par(mfrow=c(1,1))
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")
# set up colors for heatmap
col = colorRamp2(c(0, 18), c("white", "red"))
cols1 <- brewer.pal(11, "Paired")
cols2 <- brewer.pal(9, "Greens")
# subset coldata for samples in untx and ex groups
colData_sub <- colData(dds)[ind_to_keep, ]
# set up annotation bar for samples
ha1 = HeatmapAnnotation(Group = colData_sub$group,
col = list(Group = c("untreated" = cols1[1], "Dex" = cols1[2])),
show_legend = TRUE)
# se up column annotation labels (samples)
ha = columnAnnotation(x = anno_text(colData_sub$SRR,
which="column", rot = 45,
gp = gpar(fontsize = 10)))
# generate heatmap object
ht1 = Heatmap(mat_scaled, name = "Expression", col = col,
top_annotation = c(ha1),
bottom_annotation = c(ha),
show_row_names = FALSE, show_column_names = FALSE)
# plot the heatmap
par(mfrow=c(1,1))
draw(ht1, row_title = "Genes", column_title = "Hierachical clustering of DEGs (padj<0.05)")