[求助]HOMER寻找超级增强子后作ROSE类的图

kylinson

先说说为什么要用homer寻找超级增强子因为ROSE要求输入bam文件。而我拿不到bam文件，只能拿到bed、bedGraph之类的文件，通常情况下bed文件或者bedGraph文件转换为bam文件是不可能的。无奈之下只有先用bed文件在homer环境下生成tagDirectory。后续继续使用homer寻找超级增强子。
homer的output只有一个表格。然后我看看这个表格怎么作图呢

Making your own Super Enhancer Plot

Find super enhancers like you normally would, but add the option "-superSlope -1000" - the idea is to include ALL potential peaks as 'super enhancers' so that we can plot them together.  Open the resulting peak file in Excel.  The 6th column ("Normalized Tag Count") contains the super enhancer score for each region.  Simply ploting this column as a line plot will give you a sense of what your plot will look like.  To get an official "Young-lab style" plot you'll have to do some Excel algebra to normalize score by the total.

嗯，含糊其辞。不知道这个normalization方法到底是什么。

然后看看ROSE的CalllSuper.R

对我来说有点难度。。。。

所以问题来了，请问有人能根据homer生成的superEnhancers.txt画出ROSE的那种图嘛

附上表格header及前两行数据，以及SE图的示例。

#]#]# Peak finding parameters:

#]#]#

#]#]# super enhancer stitching window = 12500

#]#]# super enhancer slope moving window = 10

#

<span]#

<span]# total peaks = 34910

#]#]# peaks found using tags on both strands

#]#]# fragment length = 262

#]#]# Total tags = 19369571.0

#]#]# Approximate IP efficiency = 55.94%

#]#]# expected tags per peak = 2.847

#]#]# effective number of tags used for normalization = 10000000.0

#]#]# FDR rate threshold = 0.001000000

#]#]# FDR tag threshold = 14.0

#]#]#

#]#]# Fold over local region required = 4.00

#]#]# Putative peaks filtered by local signal = 485963

#

<span]#

<span]# Maximum fold under expected unique positions for tags = 2.00

#]#]#

#]#]#

#]#]#PeakID        chr        start        end        strand        Normalized Tag Count        superEnhancer slope        findPeaks Score        Fold Change vs Local        p-value vs Local        Clonal Fold Change

chr3-27789<span]chr3-27789<span]chr6-2384        chr6        30710774        30769390        +        4826.6        58616.000        151.000000        4.47        5.64e-51        0.78

chr19-5188<span]chr19-5188<span]chr5-9537        chr5        111272574        111360908        +        4083.7        88334.000        75.000000        4.30        7.80e-24        0.89

chr1-45761        chr1        235091137        235134056        +        3987.7        42919.000        28.000000        7.21        8.45e-18        0.95

kylinson

没人回复啊。那就献上我自己的代码了。

args <- commandArgs(TRUE)
input<-args[1]
output<-args[2]
inputfile<-paste(input,"superEnhancers.tmp",sep="/")
outputfile<-paste(output, "SuperEnhancersPlot.pdf",sep=".")
print(inputfile)
print(outputfile)
a<-read.csv(inputfile,sep = "\t",header = T)
if(dim(a)[1]<2){print("number of identified super enhancers below 2,quiting...");q()}
for (i in 1:dim(a)[1]){a$rank[i]<-1-i/dim(a)[1]}
library("ggplot2")
ggplot()+geom_line(aes(x=a$rank,y=a$Normalized.Tag.Count/max(a$Normalized.Tag.Count)),color="#FF0000",size=3)+labs(x="Ranked Peaks",y="Super Enhancer Score",title="Super Enhancer plot")
ggsave(outputfile,dpi = 600)