增强子被定义为增强子标记的峰区，如H3K27ac和H3K4me1，它们位于距注释转录起始位点（TSS）至少2 kb的位置。

将12.5 kb内的增强子峰缝合为一簇，计算每个增强子的H3K27ac结合强度，用于进一步的信号排序，以使用ROSE算法（版本2）识别SEs。每个增强子根据基因组区域中H3K27ac或H3K4me1的信号进行排序，包括缝合增强子和12.5 kb内没有相邻增强子的单个增强子。H3K27ac或H3K4me1信号超过位于H3K27ac或H3K4me1芯片序列强度分布图中信号曲线斜率为1的增强子信号的缝合增强子或单个增强子被视为SEs。信号较弱的增强子被认为是TEs。

TEs和SEs的注释是使用homer（v4.11.1）和默认参数执行的。R软件包ChIPseeker（v1.20.0）用于分析SEs和TE在整个基因组中的分布。

# 组成型增强子排序

bedtools sort -i mm10.enhancer.bed >mm10.enhancer.sort.bed

# 距离在12.5 kb内的增强子归为一簇

bedtools cluster -d 125000 -i mm10.enhancer.sort.bed >mm10.cluster.enhancer.bed

# 计算每个增强子的H3K27ac结合强度

bedtools coverage -c -a mm10.cluster.enhancer.bed -b MESC_H3K27ac/MESC_H3K27ac.rmdup.bam \>mm10.cluster.enhancer.H3K27ac.profile_tmp

# 对每簇增强子的结合强度做簇内加和
# 注意：-g 指定以那一列分组，指定的应该是标记分簇的数字所在的列;
# -c 表示对coverage所在的列计算加和 (-o sum)，注意列需要根据实际指定

bedtools groupby -i mm10.cluster.enhancer.H3K27ac.profile_tmp -g 5 -c 6 -o sum \>mm10.cluster.enhancer.H3K27ac.profile

以下为R代码请在R中运行
#####以下为R代码

enhancer = read.table("mm10.cluster.enhancer.H3K27ac.profile",header=F, row.names=NULL, sep="\t")
head(enhancer)

# 注意查看丰度信息是否在第二列，若不在，则需做相应修改

H3K27ac = sort(enhancer$V2)plot(H3K27ac, col=2, type="l")

# 计算拐点, 代码取自ROSE

numPts_below_line <- function(myVector,slope,x){yPt <- myVector[x]b <- yPt-(slope*x)xPts <- 1:length(myVector)return(sum(myVector<=(xPts*slope+b)))
}inputVector <- H3K27ac

#set those regions with more control than ranking equal to zero（将控制力大于排名的区域设置为零）

inputVector[inputVector<0]<-0

# This is the slope of the line we want to slide. This is the diagonal.（这是我们想要滑动的线的斜率。这是对角线。）

slope <- (max(inputVector)-min(inputVector))/length(inputVector)

# Find the x-axis point where a line passing through that point has the minimum number of points below it. (ie. tangent)。（查找 x 轴点，其中通过该点的直线在其下方具有最小点数。（即切线）。）
# 该点就是切点

xPt <- floor(optimize(numPts_below_line, lower=1, upper=length(inputVector),myVector= inputVector,slope=slope)$minimum)y_cutoff <- inputVector[xPt] #The y-value at this x point. This is our cutoff.b <- y_cutoff-(slope* xPt)
abline(v= xPt,h= y_cutoff,lty=2,col=8)
points(xPt,y_cutoff,pch=16,cex=0.9,col=2)
abline(coef=c(b,slope),col=2)
title(paste("x=",xPt,"\ny=",signif(y_cutoff,3),"\nFold over Median=",signif(y_cutoff/median(inputVector),3),"x\nFold over Mean=",signif(y_cutoff/mean(inputVector),3),"x",sep=""))

#Number of regions with zero signal（零信号区域数）

axis(1,sum(inputVector==0),sum(inputVector==0),col.axis="pink",col="pink")

## 超级增强子cluster

enhancer[enhancer$V2>=y_cutoff,1]