
###################################################################

# ChIPseeker for ChIP peak Annotation

setwd("F:/Chipseq-tutorial")

library(BiocManager)

BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")


library(ChIPseeker)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene) 
library(TxDb.Hsapiens.UCSC.hg38.knownGene) 


chipseq_file <- "SRR12331557_macs_output/SRR12331557.bed"
peak <- readPeakFile(chipseq_file)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
annotations <- annotatePeak(peak, tssRegion=c(-3000, 3000), TxDb=txdb)

annotations <- as.data.frame(annotations)
View(annotations)

library(org.Hs.eg.db)
library(AnnotationDbi)

keytypes(org.Hs.eg.db)

gene_symbol <- mapIds(org.Hs.eg.db, key = annotations$geneId , column = "SYMBOL" , keytype = "ENTREZID" ,multiVals = "first")

annotations <- cbind(annotations , gene_symbol)

View(annotations)

write.table(annotations, "bed annotation.txt", sep = "\t" , quote = F)


###################################################################

# Differential Binding Analysis of ChIP-seq peaksets

setwd("F:/Chipseq-tutorial")


library(BiocManager)

BiocManager::install("DiffBind")
library("DiffBind")


metadata <- read.csv("metadata.csv")
View(metadata)

dba_obj <- dba(sampleSheet = metadata)

dba.plotHeatmap(dba_obj)

dba_obj <- dba.count(dba_obj , summits = 250)

dba.plotHeatmap(dba_obj)

dba.plotPCA(dba_obj,  attributes=DBA_CONDITION, label=DBA_ID)

dba_obj <- dba.contrast(dba_obj, categories=DBA_CONDITION)

dba_obj <- dba.analyze(dba_obj) # with DESeq2

dba_obj <- dba.analyze(dba_obj, method=DBA_EDGER) #with edgeR

dba.show(dba_obj, bContrasts=T)

dba.plotVenn(dba_obj,contrast=1,method=DBA_ALL_METHODS)

dba.plotHeatmap(dba_obj, contrast=1)

dba.plotHeatmap(dba_obj, ColAttributes = DBA_CONDITION, contrast=1, correlations=FALSE)

dba.plotMA(dba_obj)

dba.plotVolcano(dba_obj)

dba.plotPCA(dba_obj, contrast = 1,label=DBA_ID)

dba.plotBox(dba_obj)

results <- dba.report(dba_obj , contrast = 1)

results <- as.data.frame(results)  

View(results)

write.table(results, "final result.txt", sep="\t", quote=F)


library(ChIPseeker)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg38.knownGene) 

View(results)

write.table(results[,1:3] , "result.bed" , sep = "\t" , quote = F , row.names = F ,col.names = F)

peak <- readPeakFile("result.bed")

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
annotations <- annotatePeak(peak, tssRegion=c(-3000, 3000), TxDb=txdb)

annotations <- as.data.frame(annotations)
View(annotations)

View(results)


annotations$start <- annotations$start - 1
annotations$end <- annotations$end -1

vec <- c()
for(i in 1:nrow(results)) {
  
  a <- results[i,2]
  b <- results[i,3]
  

  c <- annotations[annotations$start == a & annotations$end == b,c(6:12)]
  
  if (is.na(c[1,1]) == F) {
    
    vec <- rbind(vec , c)
    
  } else {
    
    vec <- rbind(vec , c("NA" , "NA"))
    
  }
  
}
  
nrow(vec) 
nrow(results)
 
results <- cbind(results , vec)

View(results)


library(org.Hs.eg.db)
library(AnnotationDbi)

keytypes(org.Hs.eg.db)

gene_symbol <- mapIds(org.Hs.eg.db, key = results$geneId , column = "SYMBOL" , keytype = "ENTREZID" ,multiVals = "first")

results <- cbind(results , gene_symbol)

results_sub <- subset(results , results$FDR < 0.05)
write.table(results_sub , "final_result.txt" , sep = "\t" , quote = F)

View(results_sub)

###################################################################





