setwd("C:/Users/S Thakur/OneDrive/Documents/microarray_tutorial/test") #set your working directory
library(GEOquery) #load libraries
library(limma)
library(clusterProfiler)
library(affy)
library(oligo)
library(clusterProfiler)
library(preprocessCore)
library(ggplot2)
library(splitstackshape)
library(dplyr)
library(AnnotationDbi)
library(org.Hs.eg.db) 
# volcano plot
dev.new()
volcano_plot <-volcanoplot(results, coef=1, main=sprintf("%d features pass our cutoffs",nrow(sig_genes1))) #results and sig_genes1 are coming from the differential gene exp analysis
points(sig_genes1[['logFC']],-log10(sig_genes1[['P.Value']]),col='red')
print(volcano_plot)
pdf("volcano_plot.pdf")
dev.off()
#heatmap 
data1 <- read.delim("ano.txt", check.names = FALSE) #data1 are coming from annotation
matrix = read.delim("exp.txt") #exp is coming from normalized exp values
mat_sym = merge(x=matrix,y=data1, by.x = "ID", by.y = "ID")
mat_sym = mat_sym[,c(1:5,20)]
mat_sym = na.omit(mat_sym)
mat_sym=mat_sym[!duplicated(mat_sym$Gene.Symbol_01),]

#row.names(mat_sym)=mat_sym$Gene.Symbol_01
#mat_sym = mat_sym[,-6]
#write.csv(mat_sym,"matrix.csv")
DEGs = read.csv("significant_genes.csv", header = T)
DEGs = DEGs[1:100,]

HP_significant = merge(mat_sym,DEGs, by.x = "ID", by.y = "ID")
row.names(HP_significant)= HP_significant$Gene.Symbol_01
HP_significant = HP_significant[,-c(1,6:12)]
HP_significant = scale(HP_significant)
library(pheatmap)
library(gplots)
fontsize_row = 10 - nrow(HP_significant) / 15
dev.new()
pheatmap(HP_significant, col=greenred(256), main="Heatmap", cluster_cols=T, 
         fontsize_row=fontsize_row, angle_col=45,  border_color=NA)
dev.copy(png, "p_heatmap.png")
dev.off()
#box plot
library(ggplot2)
args(ggplot)
gene_data <- data.frame(
  Gene = rep(c("DDR1")),  
  Condition = rep(c("drugresistant", "parent"), each = 5),  
  Expression = c(9.860,9.9168,9.6345,9.6378,9.9255,7.935,7.879,7.4567,7.4567,7.5788) # Replace with your expression values
)

# Create a gene expression box plot with color
ggplot(data = gene_data, aes(x = Condition, y = Expression, fill = Condition)) +
  geom_boxplot() +
  geom_boxplot(width = 0.5) +
  labs(
    title = "Box Plot of Normalized Gene Expression",
    x = "Condition",
    y = "Gene Expression"
  )