差异片段分析overlappeak

2023-03-22  本文已影响0人  pudding815

!!无重复样本

####H3K4me3差异分析

library(stringr)

gr1<-toGRanges(data = "H3K4me3_hCG_0h_peaks.narrowPeak",

format="narrowPeak",header=FALSE)

gr2<-toGRanges(data = "H3K4me3_hCG_4h_peaks.narrowPeak",

format="narrowPeak",header=FALSE)

class(gr1$score)

class(gr2$score)

#[1] "integer"

ol <- findOverlapsOfPeaks(gr1, gr2)

ol <- addMetadata(ol, colNames="score", FUN=mean)

# gr1和gr2都有的peaks数量

length(ol$peaklist$`gr1///gr2`)

[1] 12423

# 仅gr2有的peaks数量

length(ol$peaklist$gr2)

[1] 25005

# 仅gr1有的peaks数量

length(ol$peaklist$gr1)

[1] 429

ol$venn_cnt

makeVennDiagram(ol, fill=c("#009E73", "#F0E442"), # circle fill color

col=c("#D55E00", "#0072B2"), #circle border color

cat.col=c("#D55E00", "#0072B2")) # label color, keep same as circle border colors

#data--save

unique4h<-as.data.frame(ol$peaklist$gr2)

export(unique4h,"unique4h_peak,xlsx")

unique0h<-as.data.frame(ol$peaklist$gr1)

export(unique0h,"unique0h_peak.xlsx")

me3_commonpeak<-as.data.frame(ol$peaklist$`gr1///gr2`)

export(me3_commonpeak,"me3_commonpeak.xlsx")

####H3K27ac差异分析

library(stringr)

gr3<-toGRanges(data = "H3K27ac_hCG_0h_peaks.narrowPeak",

format="narrowPeak",header=FALSE)

gr4<-toGRanges(data = "H3K27ac_hCG_4h_peaks.narrowPeak",

format="narrowPeak",header=FALSE)

class(gr3$score)

class(gr4$score)

#[1] "integer"

oll <- findOverlapsOfPeaks(gr3, gr4)

oll <- addMetadata(oll, colNames="score", FUN=mean)

# gr3和gr4都有的peaks数量

length(oll$peaklist$`gr3///gr4`)

[1] 1502

# 仅gr3有的peaks数量

length(oll$peaklist$gr3)

[1] 3770

# 仅gr4有的peaks数量

length(oll$peaklist$gr4)

[1] 8299

oll$venn_cnt

makeVennDiagram(oll, fill=c("#009E73", "#F0E442"), # circle fill color

col=c("#D55E00", "#0072B2"), #circle border color

cat.col=c("#D55E00", "#0072B2")) # label color, keep same as circle border colors

#data--save

ac_unique4h<-as.data.frame(oll$peaklist$gr4)

export(ac_unique4h,"ac_unique4h_peak.xlsx")

ac_unique0h<-as.data.frame(oll$peaklist$gr3)

export(ac_unique0h,"ac_unique0h_peak.xlsx")

ac_commonpeak<-as.data.frame(oll$peaklist$`gr3///gr4`)

export(ac_commonpeak,"ac_commonpeak.xlsx")

上一篇下一篇

猜你喜欢

热点阅读