Skip to content

Latest commit

 

History

History
62 lines (48 loc) · 2.08 KB

002.Figure 1E (Number of promoter-localized ACRs per gene in each tissue).md

File metadata and controls

62 lines (48 loc) · 2.08 KB
Raw

Number of promoter-localized ACRs per gene in each tissue

library(ggplot2)
library(ggpointdensity) 
library(viridis)
library("ChIPseeker")
library("GenomicFeatures")
options(stringsAsFactors = FALSE)
poppy_TxDb <- makeTxDbFromGFF("D:/001poppy_atac_new_genome/Papaver_somniferum.gene.gff3")
all_peak_result <- rbind()
kb_3_peak_result<- rbind()
for (tissue in c("capsule","petal","stem","leaf","taproot","fineroot")) {
  X = annotatePeak(paste0("ATAC_seq_",tissue,".narrowPeak"),tssRegion = c(-3000,3000),TxDb = poppy_TxDb)
  perk_related_to_gene <- as.data.frame(as.GRanges(X))
  all_peak <- as.data.frame(table(perk_related_to_gene[,19]))
  all_peak <- as.data.frame(table(all_peak[,2]))
  
  all_peak <- as.matrix(all_peak)
  all_peak <- as.data.frame(all_peak)
  all_peak[,1] <- as.numeric(all_peak[,1])
  all_peak[,2] <- as.numeric(all_peak[,2])
  
  all_peak[,2] <- all_peak[,2]/sum(all_peak[,2])
  all_peak[,3] <- tissue
  perk_related_to_gene  <- perk_related_to_gene[abs(perk_related_to_gene[,21]) <= 3000,]
  kb_3_peak <- as.data.frame(table(perk_related_to_gene[,19]))
  kb_3_peak <- as.data.frame(table(kb_3_peak[,2]))
  
  kb_3_peak <- as.matrix(kb_3_peak)
  kb_3_peak <- as.data.frame(kb_3_peak)
  kb_3_peak[,1] <- as.numeric(kb_3_peak[,1])
  kb_3_peak[,2] <- as.numeric(kb_3_peak[,2])
  
  kb_3_peak[,2] <- kb_3_peak[,2]/sum(kb_3_peak[,2])
  kb_3_peak[,3] <- tissue
  all_peak_result <- rbind(all_peak_result,all_peak)
  kb_3_peak_result <- rbind(kb_3_peak_result,kb_3_peak)
  
}

all_peak_result <- all_peak_result[all_peak_result[,1] <= 6,]
kb_3_peak_result <- kb_3_peak_result[kb_3_peak_result[,1] <= 6,]
colnames(kb_3_peak_result) <- c("Number","Frequence","Tissue")
kb_3_peak_result_plot <- kb_3_peak_result
kb_3_peak_result_plot$Number <- factor(kb_3_peak_result_plot$Number)


library(RColorBrewer)
display.brewer.all()
color_name <- brewer.pal(6,"Accent")

library(ggplot2)
ggplot(data = kb_3_peak_result_plot, mapping = aes(
  x = Tissue, fill = Number, y = Frequence
)) + geom_col()+
  scale_fill_manual(values=color_name)