setwd("/Users/matteo/Desktop/") library(dplyr) library(gplots) library(DESeq2) library(ggplot2) #Read gene counts table/ normalize gene counts with DESeq2 diffexpressed<-read.table("Brains_combined_genecounts.txt",header=TRUE, row.names = 1) (conditionAdult <- factor(c(rep("H.melpomene", 12), rep("H.cydno", 11), rep("F1 Hybrid", 4),rep("H.melpomene", 5),rep("H.cydno", 5),rep("F1 Hybrid", 12) ))) (coldataAdult <- data.frame(row.names=colnames(diffexpressed), conditionAdult)) ddsAdult<- DESeqDataSetFromMatrix(countData=diffexpressed, colData=coldataAdult, design=~conditionAdult) transformedAdult<- rlog(ddsAdult) dds <- estimateSizeFactors(ddsAdult) normalized<-as.data.frame(counts(dds, normalized=TRUE)) dds2<-dds dds2$conditionAdult <- factor(dds2$conditionAdult, levels = c("H.cydno", "F1 Hybrid", "H.melpomene")) #Plot expression levels in melpomene, cydno, F1 hybrids #Example gene with statistically intermediate expression (in F1 hybrids) geneCounts <- plotCounts(dds2, gene = "HMEL017564g1", intgroup = "conditionAdult",returnData = TRUE) geneCounts <- plotCounts(dds2, gene = "HMEL034829g1", intgroup = "conditionAdult", returnData = TRUE) p <- ggplot(geneCounts, aes(x = conditionAdult, y = count)) + geom_dotplot(binaxis='y', stackdir='center',fill = c(rep("#2E86C1", 16 ), rep("#95A5A6", 16), rep("#CB4335", 17))) p + stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), geom="crossbar", width=0.5)+ theme(panel.background = element_blank(),axis.line = element_line(colour = "black")) + theme(axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_blank(), axis.text.y = element_blank(), axis.ticks.x=element_blank()) #Do the same for the other gene categories #Examples: mel-like = HMEL002663g3, cyd-like = HMEL021601g1, transgressive = HMEL031270g1