fold_change

d=read.table("/Users/JieHuang1/Downloads/H3K27Ac_SAHA_uniquePeaks_SAHA_DMSO_500bp_summed_FoldChange.txt")
x=d[,1]
x=x-mean(x)
y=d[,2]
plot(x,y)
lm=lm(y~x+x^2)
summary(lm)

## 
## Call:
## lm(formula = y ~ x + x^2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6526 -0.1789 -0.0424  0.1024 12.3537 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.126e+00  1.712e-03  1826.2   <2e-16 ***
## x           4.967e-05  1.082e-07   459.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4007 on 54793 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7936 
## F-statistic: 2.107e+05 on 1 and 54793 DF,  p-value: < 2.2e-16

points(x,lm$fitted.values,col='red',type='l')

plot(lm$residuals,type='l');abline(0,0,col='red')

#stepwise difference of the fold_change
diff = matrix(0,length(x)-1)
for(i in 1: (length(x)-1)) diff[i]=y[i+1]-y[i]
 
mu=mean(diff)  #mean
s= sd(diff)   #standard deviation 



plot(diff,type='s');abline(mu+3*s,0,col='red');abline(0,0,col='green')

plot(54500:54794,diff[54500:54794],type='s');abline(mu+3*s,0,col='red');abline(0,0,col='green')

plot(54650:54794,diff[54650:54794],type='s');abline(mu+3*s,0,col='red');abline(0,0,col='green')

#which positions goes above 3*SD
which(diff>3*s)

##  [1]     2 54704 54711 54717 54719 54725 54726 54727 54731 54734 54736
## [12] 54739 54741 54744 54745 54748 54750 54751 54752 54754 54755 54756
## [23] 54757 54759 54760 54765 54775 54784 54786 54788 54791 54792 54794