library(tidyverse)
library('DESeq2')
library(ggrepel)
library(dplyr)
library(scales)
library(viridis)
library(ggplot2)
library(RColorBrewer)
library(gplots)
library(EnhancedVolcano)
library(cowplot)
counts<- read.csv('tidied_counts.csv', header = T)
metadata <- read.csv('scleroderma_metadata.csv')
annotation <- read.csv("annotations.csv", header = T)
rownames(metadata) <- metadata$patient
rownames(counts) <- counts$precursors
counts <- counts %>% select(-X,-precursors,-total)
all(rownames(metadata))==colnames(counts)
## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
all(rownames(metadata) %in% colnames(counts))
## [1] TRUE
#metadata <- metadata %>% select()
counts <- counts[, rownames(metadata)]
all(rownames(metadata) == colnames(counts))
## [1] TRUE
####################################
#VEDOSS Vs Control
####################################
### Now we subset our data for the comparisons.
VEDOSS_Control_metadata <- filter(metadata, condition == "Control"| condition =="VEDOSS")
VEDOSS_Control_counts <- counts %>%
select(intersect(names(.), VEDOSS_Control_metadata$patient))
#this line defines our comparison order.
VEDOSS_Control_metadata$condition <- factor(VEDOSS_Control_metadata$condition, levels = c('Control','VEDOSS'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = VEDOSS_Control_counts,
colData = VEDOSS_Control_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 408
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_VEDOSS_vs_Control"
VEDOSS_VS_Control <-results(dds, alpha=0.1,name="condition_VEDOSS_vs_Control",cooksCutoff=FALSE)
VEDOSS_VS_Control.lfc <- lfcShrink(dds, coef=2,res = VEDOSS_VS_Control, type = "apeglm")
VEDOSS_VS_Control_cooks <-results(dds, alpha=0.1,name="condition_VEDOSS_vs_Control",cooksCutoff=T)
VEDOSS_VS_Control_cooks.lfc <- lfcShrink(dds, coef=2,res = VEDOSS_VS_Control_cooks, type = "apeglm")
####################################
#SSC_high Vs Control
####################################
### Now we subset our data for the comparisons.
high_control_metadata <- filter(metadata, condition == "Control"| condition =="SSC_high")
high_control_counts <- counts %>%
select(intersect(names(.), high_control_metadata$patient))
#this line defines our comparison order.
high_control_metadata$condition <- factor(high_control_metadata$condition, levels = c('Control','SSC_high'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = high_control_counts,
colData = high_control_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 424
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_SSC_high_vs_Control"
SSC_high_vs_Control <-results(dds, alpha=0.1,name="condition_SSC_high_vs_Control",cooksCutoff=FALSE)
SSC_high_vs_Control.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_Control, type = "apeglm")
SSC_high_vs_Control_cooks <-results(dds, alpha=0.1,name="condition_SSC_high_vs_Control",cooksCutoff=T)
SSC_high_vs_Control_cooks.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_Control_cooks, type = "apeglm")
####################################
#SSC_low Vs Control
####################################
### Now we subset our data for the comparisons.
low_control_metadata <- filter(metadata, condition == "Control"| condition =="SSC_low")
low_control_counts <- counts %>%
select(intersect(names(.), low_control_metadata$patient))
#this line defines our comparison order.
low_control_metadata$condition <- factor(low_control_metadata$condition, levels = c('Control','SSC_low'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = low_control_counts,
colData = low_control_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 386
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_SSC_low_vs_Control"
SSC_low_vs_Control <-results(dds,alpha=0.1,name="condition_SSC_low_vs_Control",cooksCutoff=FALSE)
SSC_low_vs_Control.lfc <- lfcShrink(dds, coef=2,res = SSC_low_vs_Control, type = "apeglm")
SSC_low_vs_Control_cooks <-results(dds,alpha=0.1,name="condition_SSC_low_vs_Control",cooksCutoff=T)
SSC_low_vs_Control_cooks.lfc <- lfcShrink(dds, coef=2,res = SSC_low_vs_Control_cooks, type = "apeglm")
####################################
#SSC_high vs SSC_low
####################################
high_low_metadata <- filter(metadata, condition == "SSC_high"| condition =="SSC_low")
high_low_counts <- counts %>%
select(intersect(names(.), high_low_metadata$patient))
#this line defines our comparison order.
high_low_metadata$condition <- factor(high_low_metadata$condition, levels = c('SSC_low','SSC_high'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = high_low_counts,
colData = high_low_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 328
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_SSC_high_vs_SSC_low"
SSC_high_vs_SSC_low <-results(dds, alpha=0.1,name="condition_SSC_high_vs_SSC_low",cooksCutoff=FALSE)
SSC_high_vs_SSC_low.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_SSC_low, type = "apeglm")
SSC_high_vs_SSC_low_cooks <-results(dds, alpha=0.1,name="condition_SSC_high_vs_SSC_low",cooksCutoff=T)
SSC_high_vs_SSC_low_cooks.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_SSC_low_cooks, type = "apeglm")
####################################
#SSC_low vs VEDOSS
####################################
low_VEDOSS_metadata <- filter(metadata, condition == "SSC_low"| condition =="VEDOSS")
low_VEDOS_counts <- counts %>%
select(intersect(names(.), low_VEDOSS_metadata$patient))
#this line defines our comparison order.
low_VEDOSS_metadata$condition <- factor(low_VEDOSS_metadata$condition, levels = c('VEDOSS','SSC_low'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = low_VEDOS_counts,
colData = low_VEDOSS_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 305
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_SSC_low_vs_VEDOSS"
SSC_low_vs_VEDOSS <-results(dds, alpha=0.1,name="condition_SSC_low_vs_VEDOSS",cooksCutoff=FALSE)
SSC_low_vs_VEDOSS.lfc <- lfcShrink(dds, coef=2,res = SSC_low_vs_VEDOSS, type = "apeglm")
SSC_low_vs_VEDOSS_cooks <-results(dds, alpha=0.1,name="condition_SSC_low_vs_VEDOSS",cooksCutoff=T)
SSC_low_vs_VEDOSS_cooks.lfc <- lfcShrink(dds, coef=2,res = SSC_low_vs_VEDOSS_cooks, type = "apeglm")
####################################
#SSC_high vs VEDOSS
####################################
high_VEDOSS_metadata <- filter(metadata, condition == "SSC_high"| condition =="VEDOSS")
high_VEDOSS_counts <- counts %>%
select(intersect(names(.), high_VEDOSS_metadata$patient))
#this line defines our comparison order.
high_VEDOSS_metadata$condition <- factor(high_VEDOSS_metadata$condition, levels = c('VEDOSS','SSC_high'))
###DESEQ
dds <- DESeqDataSetFromMatrix(countData = high_VEDOSS_counts,
colData = high_VEDOSS_metadata,
design = ~condition)
# Run DESeq analysis
keep <- rowSums(counts(dds)) >=10
count(keep == TRUE)
## [1] 359
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "condition_SSC_high_vs_VEDOSS"
SSC_high_vs_VEDOSS <-results(dds, alpha=0.1,name="condition_SSC_high_vs_VEDOSS",cooksCutoff=FALSE)
SSC_high_vs_VEDOSS.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_VEDOSS, type = "apeglm")
SSC_high_vs_VEDOSS_cooks <-results(dds, alpha=0.1,name="condition_SSC_high_vs_VEDOSS",cooksCutoff=T)
SSC_high_vs_VEDOSS_cooks.lfc <- lfcShrink(dds, coef=2,res = SSC_high_vs_VEDOSS_cooks, type = "apeglm")
####################################
#VEDOSS Vs Control
####################################
MA_plot <- as.data.frame(VEDOSS_VS_Control.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("VEDOSS vs Control")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
####################################
#SSC_high Vs Control
####################################
MA_plot <- as.data.frame(SSC_high_vs_Control.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("SSC high vs Control")
####################################
#SSC_low Vs Control
####################################
MA_plot <- as.data.frame(SSC_low_vs_Control.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("SSC low vs Control")
####################################
#SSC_high vs SSC_low
####################################
MA_plot <- as.data.frame(SSC_high_vs_SSC_low.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("SSC high vs SSC low")
####################################
#SSC_low vs VEDOSS
####################################
MA_plot <- as.data.frame(SSC_low_vs_VEDOSS.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("SSC low vs VEDOSS")
####################################
#SSC_high vs VEDOSS
####################################
MA_plot <- as.data.frame(SSC_high_vs_VEDOSS.lfc)
ggplot(MA_plot, aes(baseMean, log2FoldChange, colour=pvalue)) + geom_point(size=2) + scale_y_continuous(limits=c(-3, 3), oob=squish) + scale_x_log10() + geom_hline(yintercept = 0, colour="darkorchid4", size=1, linetype="longdash") + labs(x="mean of normalized counts", y="log fold change") + scale_colour_viridis(direction=-1, trans='sqrt') + theme_bw()+ggtitle("SSC high vs VEDOSS")
###create differential expression table
####################################
#VEDOSS Vs Control
####################################
Vedoss_vs_control.lfcDF <- data.frame(VEDOSS_VS_Control.lfc)
Vedoss_vs_control_sig <- Vedoss_vs_control.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(VEDOSS_vs_healthy = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(Vedoss_vs_control_sig, "VEDOSS_vs_healthy_significant_miRNAs.csv")
####################################
#SSC_high Vs Control
####################################
SSC_high_vs_Control.lfcDF <- data.frame(SSC_high_vs_Control.lfc)
SSC_high_vs_Control_sig <- SSC_high_vs_Control.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(high_vs_healthy = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(SSC_high_vs_Control_sig, "SSC_high_vs_Control_significant_miRNAs.csv")
####################################
#SSC_low Vs Control
####################################
SSC_low_vs_Control.lfcDF <- data.frame(SSC_low_vs_Control.lfc)
SSC_low_vs_Control_sig <- SSC_low_vs_Control.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(low_vs_healthy = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(SSC_low_vs_Control_sig, "SSC_low_vs_Control_significant_miRNAs.csv")
####################################
#SSC_high vs SSC_low
####################################
SSC_high_vs_SSC_low.lfcDF <- data.frame(SSC_high_vs_SSC_low.lfc)
SSC_high_vs_SSC_low_sig <- SSC_high_vs_SSC_low.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(high_vs_low = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(SSC_high_vs_SSC_low_sig, "SSC_high_vs_SSC_low_significant_miRNAs.csv")
####################################
#SSC_low vs VEDOSS
####################################
SSC_low_vs_VEDOSS.lfcDF <- data.frame(SSC_low_vs_VEDOSS.lfc)
SSC_low_vs_VEDOSS_sig <- SSC_low_vs_VEDOSS.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(low_vs_VEDOSS = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(SSC_low_vs_VEDOSS_sig, "SSC_low_vs_VEDOSS_significant_miRNAs.csv")
####################################
#SSC_high vs VEDOSS
####################################
SSC_high_vs_VEDOSS.lfcDF <- data.frame(SSC_high_vs_VEDOSS.lfc)
SSC_high_vs_VEDOSS_sig <- SSC_high_vs_VEDOSS.lfcDF %>% subset(pvalue<0.05) %>% filter(log2FoldChange > 0.585 | log2FoldChange < -0.585, baseMean >50) %>% mutate(high_vs_VEDOSS = ifelse(log2FoldChange > 0, "up", "down")) %>% rownames_to_column("precursor") %>% left_join(annotation,by="precursor") %>% select(-X)
write.csv(SSC_high_vs_VEDOSS_sig, "SSC_high_vs_VEDOSS_significant_miRNAs.csv")
significant_miRs <- Vedoss_vs_control_sig %>% full_join(SSC_high_vs_Control_sig, by="precursor") %>% full_join(SSC_low_vs_Control_sig, by="precursor") %>% full_join(SSC_high_vs_SSC_low_sig, by="precursor") %>% full_join(SSC_low_vs_VEDOSS_sig, by="precursor") %>% full_join(SSC_high_vs_VEDOSS_sig, by="precursor") %>% select(precursor,VEDOSS_vs_healthy,high_vs_healthy,low_vs_healthy,high_vs_low,low_vs_VEDOSS,high_vs_VEDOSS) %>% left_join(annotation,by='precursor') %>% select(-X,-Accession,-precursor) %>% unique()
write.csv(significant_miRs, "differentially_regulated_miRs.csv")
####################################
#VEDOSS Vs Control
####################################
Vedoss_vs_control.lfcDF <- as.data.frame(VEDOSS_VS_Control_cooks.lfc)
threshold_isf <- Vedoss_vs_control.lfcDF$pvalue < 0.05
Vedoss_vs_control.lfcDF$threshold <- threshold_isf
isf_table_ordered <- Vedoss_vs_control.lfcDF[order(Vedoss_vs_control.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
VEDOSS_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs VEDOSS vs Healthy',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
maxoverlapsConnectors = 20)
VEDOSS_volcano
####################################
#SSC_high Vs Control
####################################
SSC_high_vs_Control.lfcDF <- as.data.frame(SSC_high_vs_Control_cooks.lfc)
threshold_isf <- SSC_high_vs_Control.lfcDF$pvalue < 0.05
SSC_high_vs_Control.lfcDF$threshold <- threshold_isf
isf_table_ordered <- SSC_high_vs_Control.lfcDF[order(SSC_high_vs_Control.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
SSC_high_control_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs SSC high vs Healthy',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
max.overlaps = 100,
maxoverlapsConnectors = 20)
SSC_high_control_volcano
## Warning: ggrepel: 27 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
####################################
#SSC_low Vs Control
####################################
SSC_low_vs_Control.lfcDF <- as.data.frame(SSC_low_vs_Control_cooks.lfc)
threshold_isf <- SSC_low_vs_Control.lfcDF$pvalue < 0.05
SSC_low_vs_Control.lfcDF$threshold <- threshold_isf
isf_table_ordered <- SSC_low_vs_Control.lfcDF[order(SSC_low_vs_Control.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
SSC_low_control_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs SSC low vs Healthy',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
max.overlaps = 100,
maxoverlapsConnectors = 20)
SSC_low_control_volcano
####################################
#SSC_high vs SSC_low
####################################
SSC_high_vs_SSC_low.lfcDF <- as.data.frame(SSC_high_vs_SSC_low_cooks.lfc)
threshold_isf <- SSC_high_vs_SSC_low.lfcDF$pvalue < 0.05
SSC_high_vs_SSC_low.lfcDF$threshold <- threshold_isf
isf_table_ordered <- SSC_high_vs_SSC_low.lfcDF[order(SSC_high_vs_SSC_low.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
SSC_high_vs_low_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs SSC high vs SSC low',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
max.overlaps = 100,
maxoverlapsConnectors = 20)
SSC_high_vs_low_volcano
## Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
####################################
#SSC_low vs VEDOSS
####################################
SSC_low_vs_VEDOSS.lfcDF <- as.data.frame(SSC_low_vs_VEDOSS_cooks.lfc)
threshold_isf <- SSC_low_vs_VEDOSS.lfcDF$pvalue < 0.05
SSC_low_vs_VEDOSS.lfcDF$threshold <- threshold_isf
isf_table_ordered <- SSC_low_vs_VEDOSS.lfcDF[order(SSC_low_vs_VEDOSS.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
SSC_low_vs_VEDOSS_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs SSC low vs VEDOSS',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
max.overlaps = 100,
maxoverlapsConnectors = 20)
SSC_low_vs_VEDOSS_volcano
## Warning: Removed 1 rows containing missing values (`geom_vline()`).
####################################
#SSC_high vs VEDOSS
####################################
SSC_high_vs_VEDOSS.lfcDF <- as.data.frame(SSC_high_vs_VEDOSS_cooks.lfc)
threshold_isf <- SSC_high_vs_VEDOSS.lfcDF$pvalue < 0.05
SSC_high_vs_VEDOSS.lfcDF$threshold <- threshold_isf
isf_table_ordered <- SSC_high_vs_VEDOSS.lfcDF[order(SSC_high_vs_VEDOSS.lfcDF$padj),]
threshold_OE <- isf_table_ordered$threshold==TRUE
isf_table_ordered$threshold <- threshold_OE
isf_table_ordered <- isf_table_ordered[order(isf_table_ordered$padj),]
isf_table_ordered$precursor <- ""
isf_table_ordered$precursor <- rownames(isf_table_ordered)
isf_table_ordered <- left_join(isf_table_ordered,annotation, by="precursor")
genes_to_plot <- isf_table_ordered %>% filter(threshold =='TRUE' & log2FoldChange > 0.585|log2FoldChange < -0.585, baseMean >50) %>% select(miRNA)
genes_to_plot<- as.character(genes_to_plot[,1])
SSC_high_vs_VEDOSS_volcano <- EnhancedVolcano(isf_table_ordered, x = 'log2FoldChange', y = 'pvalue',lab=isf_table_ordered$miRNA,
title = 'ISF associated scleroderma microRNAs SSC high vs VEDOSS',
pCutoff = 0.05,
FCcutoff = 0.585,
xlim = c(min(isf_table_ordered[['log2FoldChange']], na.rm = TRUE), max(isf_table_ordered[['log2FoldChange']], na.rm = TRUE)),
ylim = c(0, max(-log10(isf_table_ordered[['pvalue']]), na.rm = TRUE)),
pointSize = 4.0,
labSize = 4.0,
labCol = 'black',
labFace = 'bold',
boxedLabels = TRUE,
colAlpha = 4/5,
legendPosition = 'None',
legendLabSize = 14,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
max.overlaps = 100,
maxoverlapsConnectors = 20)
SSC_high_vs_VEDOSS_volcano
