RChIP-Tech.Rmd

---
title: "R-ChIP MC Paper"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Common R Function
```{r Common R Function}
# a) Working directory
setwd("~/Desktop/R-ChIP-Tech/Bin")
# b) my own theme
library(ggplot2)
my_clear_theme = theme_bw() + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA),
        legend.background = element_blank())
# c) BinnedErrorBarPlot Function
#    if x2/y2 are given, then normalize data[,x1/y1] by dividing data[,x2/y2]
BinnedErrorBarPlot <- function(data,x1,y1,bin=15,x2=0,y2=0,xlab="x",ylab="y"){
  if(x2 == 0){
    bk <- as.double(quantile(data[,x1],probs = seq(0,1,1/bin)))
    status <- .bincode(data[,x1],bk,include.lowest = T)
    x <- data.frame(data[,x1],status)
  }else{
    bk <- as.double(quantile(data[,x1]/data[,x2],probs = seq(0,1,1/bin)))
    status <- .bincode(data[,x1]/data[,x2],bk,include.lowest = T)
    x <- data.frame(data[,x1]/data[,x2],status)
  }
  
  if(y2 == 0){
    y <- data.frame(data[,y1],status)
  }else{
    y <- data.frame(data[,y1]/data[,y2],status)
  }
  
  me_x <- c(); q25_x <- c(); q75_x <- c();
  me_y <- c(); q25_y <- c(); q75_y <- c();
  for(i in 1:bin){
    me_x[i]  <- median(x[status==i,1]);
    q25_x[i]  <- quantile(x[status==i,1],probs = 0.25);
    q75_x[i] <- quantile(x[status==i,1],probs = 0.75)
    me_y[i]  <- median(y[status==i,1])
    q25_y[i]  <- quantile(y[status==i,1],probs = 0.25)
    q75_y[i] <- quantile(y[status==i,1],probs = 0.75)
  }
  mat <- as.data.frame(cbind(me_x,q25_x,q75_x,me_y,q25_y,q75_y))
  
  library(ggplot2)
  p <- ggplot(mat,aes(x=me_x,y=me_y)) +
    geom_errorbar(aes(ymin = q25_y,ymax = q75_y),col="blue") + 
    geom_errorbarh(aes(xmin = q25_x,xmax = q75_x),col="blue") +
    geom_point(col="red") +
    xlab(xlab) +
    ylab(ylab) +
    theme(panel.grid.major = element_line(colour="NA"),
          panel.grid.minor = element_line(colour="NA"),
          panel.background = element_rect(fill="NA"),
          panel.border = element_rect(colour="black", fill=NA))
  print(p)
  cor.test(x[,1],y[,1],method = "spearman")
}
# d) BinnedHeatmapPlot Function
#    if x2/y2/z2 are given, then normalize data[,x1/y1/z1] by dividing data[,x2/y2/z2]
BinnedHeatmapPlot <- function(data,z1,x1,y1,bin=6,x2=0,y2=0,z2=0){
  if(x2 == 0){
    bk1 <- as.double(quantile(data[,x1],probs = seq(0,1,1/bin)))
    x   <- cbind(data[,x1],.bincode(data[,x1],bk1,include.lowest = T))
  }else{
    bk1 <- as.double(quantile(data[,x1]/data[,x2],probs = seq(0,1,1/bin)))
    x   <- cbind(data[,x1]/data[,x2],.bincode(data[,x1]/data[,x2],bk1,include.lowest = T))
  }
  if(y2 == 0){
    bk2 <- as.double(quantile(data[,y1],probs = seq(0,1,1/bin)))
    y   <- cbind(data[,y1],.bincode(data[,y1],bk2,include.lowest = T))
  }else{
    bk2 <- as.double(quantile(data[,y1]/data[,y2],probs = seq(0,1,1/bin)))
    y   <- cbind(data[,y1]/data[,y2],.bincode(data[,y1]/data[,y2],bk2,include.lowest = T))
  }
  if(z2==0){
    z <- data[,z1]
  }else{
    z <- data[,z1]/data[,z2]
  }
  status <- cbind(x[,2],y[,2],x[,1],y[,1],z)
  num <- matrix(0,nrow=bin,ncol=bin)
  for (i in 1:nrow(status)){
    num[status[i,1],status[i,2]] = num[status[i,1],status[i,2]] + 1
  }
  mat <- matrix(0,nrow=bin,ncol=bin)
  for(i in 1:bin){
    for (j in 1:bin){
      mat[i,j] = median(status[status[,1]==i & status[,2]==j,5])
    }
  }
  library(pheatmap)
  pheatmap(mat,cluster_rows = F,cluster_cols = F)
  status <- subset(status,rowSums(is.finite(status))==5)
  summary(lm(status[,5] ~ status[,3] + status[,4]))
}
```

##########################################################################################
## * Figure 1. Genome-wide R-loop Profiling by Strand-specific R-ChIP
##########################################################################################
## * Fig 1F
The signal intensity profile of R-ChIP data as well as the negative control (WKKD) and input in HEK293T cells.
```{r D210N-WKKD-Input Profile,fig.width=2.46,fig.height=2.49,fig.align = "center",message=F,warning=F}
library(ggplot2)
library(reshape2)
wkkd  <- read.table("Data_MC/Profile_Fig1F/WKKD")
input <- read.table("Data_MC/Profile_Fig1F/Input")
d210n <- read.table("Data_MC/Profile_Fig1F/D210N")
data  <- data.frame(apply(wkkd,2,mean),apply(input,2,mean),apply(d210n,2,mean))
names(data)  = c("WKKD","Input","D210N")
data$Position  = 1:500
data <- melt(data,id.vars = "Position")
ggplot(data,aes(x=Position,y=value,col=variable)) + 
  geom_line(size=1) +
  xlab("Genomic Region (5'-3')") + ylab("RPM") +
  my_clear_theme +
  scale_x_continuous(breaks=seq(0,500,by=100),
                     labels=c("-2K", "", "5'", "3'","","2K")) +
  theme(legend.position=c(.8,.8)) +
  labs(col="")
rm(wkkd,input,d210n,data)
```
## * Fig 1G
Transcription activities in the sense or anti-sense orientation of R-loops.
```{r Strandness (GRO-seq),fig.width=2.05,fig.height=2.39,fig.align = "center",message=F,warning=F}
library(ggplot2)
data <- read.table("Data_MC/Strand_Fig1G/RPKM.data.frame",sep="\t")
data$V2 = factor(data$V2,levels=c("same","oppo"))
ggplot(data,aes(x=V2,y=V1,col=V2)) + 
  geom_boxplot(outlier.colour = "NA",notch = T) +
  coord_cartesian(ylim=c(0,25)) +
  my_clear_theme +
  xlab("") + ylab("GRO-seq (RPKM)") + guides(col=F) +
  scale_x_discrete(labels=c("same" = "Sense", "oppo" = "Antisense"))
rm(data)
```
## * Fig 1H
The strand specificity of R-ChIP signals.
```{r Strandness (R-ChIP),fig.width=2.46,fig.height=2.49,fig.align = "center",message=F,warning=F}
library(ggplot2)
library(reshape2)
sense <- read.table("Data_MC/Strand_Fig1H/D210N")
antisense <- read.table("Data_MC/Strand_Fig1H/D210N.Rev")
data  <- data.frame(apply(sense,2,mean),apply(antisense,2,mean))
names(data) = c("Sense","Antisense")
data$Position = 1:500
data <- melt(data,id.vars = "Position")
ggplot(data,aes(x=Position,y=value,col=variable)) + 
  geom_line(size=1) +
  xlab("Genomic Region (5'-3')") + ylab("RPM") +
  my_clear_theme +
  scale_x_continuous(breaks=seq(0,500,by=100),
                     labels=c("-2K", "", "5'", "3'","","2K")) +
  theme(legend.position=c(.8,.8)) +
  labs(col="")
rm(sense,antisense,data)
```
## * Fig S1C and E
Pairwise comparison of three R-ChIP replicates performed on D210N expressing HEK293T cells.
Pairwise comparison of R-ChIP replicates from ΔHC and D210N expressing HEK293T cells.
```{r R-ChIP Rep HEK293,fig.width=6.46,fig.height=2.53,fig.align = "center",message=F,warning=F}
rep1 <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.HEK293.D210N_rep1.ReadsCoverage",sep="\t")
rep2 <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.HEK293.D210N_rep2.ReadsCoverage",sep="\t")
rep3 <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.HEK293.D210N_rep3.ReadsCoverage",sep="\t")
dele <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.HEK293.135AA.ReadsCoverage",sep="\t")

par(mfrow=c(1,3))
smoothScatter(log2(rep1[,4]),log2(rep2[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="Rep 1 (Log2 Tag Counts)",ylab="Rep 2 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep1[,4],rep2[,4]),2),sep=""))
smoothScatter(log2(rep1[,4]),log2(rep3[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="Rep 1 (Log2 Tag Counts)",ylab="Rep 3 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep1[,4],rep3[,4]),2),sep=""))
smoothScatter(log2(rep2[,4]),log2(rep3[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="Rep 2 (Log2 Tag Counts)",ylab="Rep 3 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep2[,4],rep3[,4]),2),sep=""))

par(mfrow=c(1,3))
smoothScatter(log2(rep1[,4]),log2(dele[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="D210N Rep 1 (Log2 Tag Counts)",ylab="135AA (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep1[,4],dele[,4]),2),sep=""))
smoothScatter(log2(rep2[,4]),log2(dele[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="D210N Rep 2 (Log2 Tag Counts)",ylab="135AA (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep2[,4],dele[,4]),2),sep=""))
smoothScatter(log2(rep3[,4]),log2(dele[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="D210N Rep3 (Log2 Tag Counts)",ylab="135AA (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep3[,4],dele[,4]),2),sep=""))
rm(dele,rep1,rep2,rep3)
```
##
##########################################################################################
## Figure 2. Sequence Features and Genomic Distribution of R-ChIP Signals
##########################################################################################
## * Fig 2A
Distribution of sizes of R-ChIP narrow and broad peakas by the MACS2 program.
```{r Size Cmp (boxplot),fig.width=2.46,fig.height=2.19,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/NarrowVSBroad_Figs2AS2/peakSize",sep="\t")
data <- subset(data,data$V2=="R-ChIP" & data$V3=="HEK293")
data$V1 <- log10(data$V1)
data$V4 <- factor(data$V4,levels=c("Narrow","Broad"))
ggplot(data,aes(x=V4,y=V1,col=V4)) + geom_boxplot(notch = T,outlier.size = .5) +
  scale_y_continuous(breaks=c(2,log10(150),log10(300),log10(500),3,log10(2000)),
                     labels = c(100,150,300,500,1000,2000)) +
  coord_cartesian(ylim=c(2,3.5)) + labs(col="Peak Size") +
  my_clear_theme + xlab("") + ylab("Peak Size")
```
## * Fig 2B
Base composition, G/C content and G/C skew associated with a composite R-loop map.
```{r Sequence (HEK293),fig.width=2.80,fig.height=2.51,fig.align = "center",message=F,warning=F}
library(reshape2)
library(ggplot2)
allseq <- read.table("Data_MC/Sequence/tts.all.nucleotides",sep="\t",header = T)
allseq <- melt(allseq,id.vars = "Pos")
ggplot(data = allseq,aes(x=Pos-400,y=value,col=variable)) + 
  geom_line() + 
#  stat_smooth(method = "loess", span=1/3) + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  labs(col="") + ylim(-0.05,0.7) + xlab("Distance to R-loop summit (nt)")
```
```{r Rloop Sig (HEK293),fig.width=2.52,fig.height=3.67,fig.align = "center",message=F,warning=F}
signal <- read.table("Data_MC/Signal/summit400.signal")
input  <- read.table("Data_MC/Signal/summit400.input")
plot(seq(-400,400,by=1),apply(signal,2,mean) - apply(input,2,mean),type='l',
     xlab="Distance to R-loop summit",ylab="RPM")
rm(signal,input)
```
## * Fig 2C
Percentages of the total R-loop number according to associated consecutive G numbers (G-clusters) in the ±50 bp flanking region of the G/C skew summit in comparison with background.
```{r Consective Gs, fig.width=2.50,fig.height=2.40,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/Sequence/HEK293.consecutiveGs")
data$num = c(2:6,2:6)
ggplot(data,aes(x=num,y=V2/12906,fill=V3)) + 
  geom_bar(stat="identity",position = position_dodge()) +
  my_clear_theme +
  ylab("% R-loop") + xlab("# Consecutive Gs") +
  theme(legend.position = c(.8,.8)) + labs(fill="")
```
## * Fig 2D
Coincidence between potential G-quadruplex (G4) forming regions and R-ChIP mapped R-loops, emphasizing predominant overlap on the non-template DNA strand.
```{r G4 Dis,fig.width=4.45,fig.height=4.81,fig.align = "center",message=F,warning=F}
G4=7724; noG4=12906-G4; Total=12906
plus=7014;Other=Total-plus
minus=1466;overlap=756;plusonly=plus-overlap;minusonly=minus-overlap;Other2 = noG4
iniR=.2 
colors=list(Other='white',G4='#e5f5e0',noG4='#a1d99b',plus='#3182bd',minus='#fec44f')
library('plotrix')
# from outer circle to inner circle
pie(1, radius=iniR, init.angle=90, col=c('white'), border = NA, labels='')
floating.pie(0,0,c(plus,Other),
             radius=4*iniR, startpos=pi/2, col=as.character(colors[c('plus','Other')]),border=NA)
floating.pie(0,0,c(plusonly,overlap,minusonly,Other2),radius=3*iniR,
             startpos=pi/2, col=as.character(colors[c('Other',"minus","minus","Other")]),border=NA)
floating.pie(0,0, c(G4,noG4), radius=2*iniR, startpos=pi/2,
             col=as.character(colors[c('G4',"noG4")]), border = NA)
legend(0.6, 3*iniR, gsub("_"," ",names(colors)[-1]),
       col=as.character(colors[-1]), pch=19,bty='n', ncol=1)
```
## * Fig 2E
R-loop profile relative to sequences that have the potential to form G-quadruplex.
```{r G4 profile,fig.width=2.60,fig.height=2.60,fig.align = "center",message=F,warning=F}
library(ggplot2)
data <- read.table("Data_MC/Sequence/HEK293.all.G4",sep="\t")
ggplot(data,aes(x=V1-401,y=V2/12906)) + 
  geom_line() +
  xlab("Distance to R-Loop Summit") +
  ylab("Proportion of G Quandraplex") + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA))
```
## alternative of Fig 2F (overlapping annotation)
The genomic distribution of R-ChIP mapped R-loops by RChIP in HEK293.
```{r Genomic Dis (HEK293),fig.width=3.14,fig.height=1.04,fig.align = "center",message=F,warning=F}
library(reshape2)
library(ggplot2)
# TSS, TSS/GB, TSS/GB/TTS, TSS/TTS, TTS, TTS/GB, GB, Inter 
RChIP <- c(4512,1800,690,1075,170,297,2177,2185)
End   <- cumsum(RChIP)

Row4 <- c(rep(1,End[4]),rep(0,End[8]-End[4]))
Row3 <- c(rep(0,End[2]),rep(2,End[6]-End[2]),rep(0,End[8]-End[6]))
Row2 <- c(rep(0,End[1]),rep(3,End[3]-End[1]),rep(0,End[5]-End[3]),
          rep(3,End[7]-End[5]),rep(0,End[8]-End[7]))
Row1 <- c(rep(0,End[7]),rep(4,End[8]-End[7]))

data <- as.data.frame(cbind(Row1,Row2,Row3,Row4))
data$Pos <- 1:nrow(data)
data <- melt(data,id.vars = "Pos")
ggplot(data,aes(y=variable,x=Pos/max(Pos),fill=as.factor(value))) + 
  geom_tile() + my_clear_theme +
  scale_fill_manual(values=c("000000","#1f78b4","#ff7f00","#add8e6","#33a02c")) +
  guides(fill=F) + xlab("") + ylab("") +
  scale_y_discrete(labels=c("Intergenic","GeneBody","TTS","TSS")) +
  scale_x_continuous(labels=paste(seq(0,100,by=25),"%",sep="")) +
  geom_vline(xintercept = End/max(End),linetype=1,col="grey",alpha=0.5)
  
```
## * Fig 2G
R-loop signal intensity at different genomic regions
```{r Regional Signal (HEK293),fig.width=2.46,fig.height=2.61,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/Signal/HEK293.RChIP.anno")
data$V13 <- factor(data$V13,levels=c("TSS","Genebody","TTS","Inter"))
ggplot(data,aes(x=V13,y=V5,col=V13)) + 
  geom_boxplot(outlier.colour = "NA",notch = T) +
  coord_cartesian(ylim=c(0,1.5)) +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  guides(col=F) + xlab("") + ylab("R-loop Signal")

```
## * Fig S2B
Overall R-loop signal intensity associated with each of the four R-ChIP groups.
```{bash eval = F, echo = T}
# Peaks
broad=../../data/R-ChIP/HEK293/peaks/broad/all.peaks.bed6
narrow=../../data/R-ChIP/HEK293/peaks/narrow/all.peaks.bed6
cell=HEK293
# For Venn
echo "Common (Broad)"
bedtools intersect -a $broad  -b $narrow -wa -s | sort -u | wc -l
echo "Common (Narrow)"
bedtools intersect -a $narrow -b $broad  -wa -s | sort -u | wc -l
echo "Broad Specific"
bedtools intersect -a $broad  -b $narrow -wa -s -v | sort -u | wc -l
echo "Narrow Specific"
bedtools intersect -a $narrow -b $broad  -wa -s -v | sort -u | wc -l
```
```{r Signal Cmp (Narrow vs Broad),fig.width=3.56,fig.height=1.70,fig.align = "center",message=F,warning=F}
library(ggplot2)
data <- read.table("Data_MC/NarrowVSBroad_Figs2AS2/HEK293.PeakSignal",sep="\t")
ggplot(data,aes(x=V2,y=V1,col=V2)) + 
  geom_boxplot(outlier.colour = "NA") +
  coord_cartesian(ylim=c(0,1.8)) +
    theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  theme(axis.text.x = element_blank()) +
  xlab("") + ylab("R-loop Signal") + labs(col="")
wilcox.test(data[data$V2=="Broad_Common",1],data[data$V2=="Broad_Specific",1])
wilcox.test(data[data$V2=="Narrow_Common",1],data[data$V2=="Narrow_Specific",1])
```
##
##########################################################################################
## Fig 3. Systematic Comparison of R-loops Captured by the Catalytically Dead RNASEH1 versus S9.6
##########################################################################################
## * Fig 3A
Peak Size (RChIP vs DRIP)
```{r Peak Size (vs DRIP),fig.width=2.46,fig.height=2.30,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/peakSize.3tech")
data$V2=factor(data$V2,levels=c("R-ChIP","RDIP","DRIP","DRIPc"))
data <- subset(data,data$V3=="K562")
ggplot(data,aes(x=V3,y=log10(V1),fill=V2)) + 
  geom_boxplot(notch = T,outlier.size = 0.4) + 
  my_clear_theme +
  scale_y_continuous(breaks=c(2,log10(150),log10(300),log10(500),3,log10(2000),4,5),
                     labels = c(100,150,300,500,1000,2000,10000,100000)) +
  ylab("Peak Size") + xlab("") + labs(fill="")
```
## alternative of Fig 3D (overlapping annotation)
The overall genomic distribution of R-loops mapped by R-ChIP and DRIP-seq.
```{r,Genomic Dis (K562),fig.width=3.14,fig.height=1.04,fig.align = "center",message=F,warning=F}
library(reshape2)
library(ggplot2)
# TSS, TSS/GB, TSS/GB/TTS, TSS/TTS, TTS, TTS/GB, GB, Inter 
RChIP <- c(2794,1226,567,937,106,166,1721,2177)
# or
DRIP <- c(229,1741,4039,184,185,804,1125,377)
End   <- cumsum(DRIP)

Row4 <- c(rep(1,End[4]),rep(0,End[8]-End[4]))
Row3 <- c(rep(0,End[2]),rep(2,End[6]-End[2]),rep(0,End[8]-End[6]))
Row2 <- c(rep(0,End[1]),rep(3,End[3]-End[1]),rep(0,End[5]-End[3]),
          rep(3,End[7]-End[5]),rep(0,End[8]-End[7]))
Row1 <- c(rep(0,End[7]),rep(4,End[8]-End[7]))

data <- as.data.frame(cbind(Row1,Row2,Row3,Row4))
data$Pos <- 1:nrow(data)
data <- melt(data,id.vars = "Pos")
ggplot(data,aes(y=variable,x=Pos/max(Pos),fill=as.factor(value))) + 
  geom_tile() + my_clear_theme +
  scale_fill_manual(values=c("000000","#1f78b4","#ff7f00","#add8e6","#33a02c")) +
  guides(fill=F) + xlab("") + ylab("") +
  scale_y_discrete(labels=c("Intergenic","GeneBody","TTS","TSS")) +
  scale_x_continuous(labels=paste(seq(0,100,by=25),"%",sep="")) +
  geom_vline(xintercept = c(0,End/max(End)),linetype=1,col="grey",alpha=0.5)
```
## Fig 3E, S3B (similar to Fig 2B)
## Fig 3F, S3C (similar to Fig 1F)
## * Fig S3A
Pairwise comparison of three R-ChIP replicates performed on D210N expressing HEK293T cells.
```{r R-ChIP Rep K562,fig.width=3.39,fig.height=3.99,message=F,warning=F}
rep1 <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.K562.D210N_rep1.ReadsCoverage",sep="\t")
rep2 <- read.table("Data_MC/RChIPRepli_FigsS1CE3A/chipseq.K562.D210N_rep2.ReadsCoverage",sep="\t")
smoothScatter(log2(rep1[,4]),log2(rep2[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="Rep 1 (Log2 Tag Counts)",ylab="Rep 2 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(rep1[,4],rep2[,4]),2),sep=""))
rm(rep1,rep2)
```
## Fig S4A & 4B
(RDIP) Base composition, G/C content and G/C skew associated with a composite R-loop map.
```{r Nuc Comp RDIP,fig.width=2.94,fig.height=2.28,message=F,warning=F}
library(reshape2)
library(ggplot2)
allseq <- read.table("Data_MC/Sequence/HEK293.RDIP.all.nucleotides",sep="\t",header = T)
allseq <- melt(allseq,id.vars = "Pos")
ggplot(data = allseq,aes(x=Pos-800,y=value,col=variable)) + 
#  geom_line() + 
  stat_smooth(method = "loess", span=1/3) + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  labs(col="") + ylim(-0.05,0.7)
```
```{r RDIP signalfig.width=2.52,fig.height=3.67,fig.align = "center",message=F,warning=F}
signal <- read.table("Data_MC/Signal/RDIP.summit.800.signal")
input  <- read.table("Data_MC/Signal/RDIP.summit.800.input")
plot(seq(-400,400,by=1),apply(signal,2,mean) - apply(input,2,mean),type='l',
     xlab="Distance to R-loop summit",ylab="RPM")
rm(signal,input)
```
(DRIPc) Base composition, G/C content and G/C skew associated with a composite R-loop map.
```{r Nuc Comp DRIPc,fig.width=2.94,fig.height=2.28,message=F,warning=F}
library(reshape2)
library(ggplot2)
allseq <- read.table("Data_MC/Sequence/DRIPc.all.nucleotides",sep="\t",header = T)
allseq <- melt(allseq,id.vars = "Pos")
ggplot(data = allseq,aes(x=Pos-4000,y=value,col=variable)) + 
#  geom_line() + 
  stat_smooth(method = "loess", span=1/3) + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  labs(col="") + ylim(-0.05,0.7)
rm(allseq)
```
```{r DRIPc signal,fig.width=2.52,fig.height=3.67,fig.align = "center",message=F,warning=F}
signal <- read.table("Data_MC/Signal/shuf")
plot(seq(-4000,4000,by=1),apply(signal,2,mean),type='l',
     xlab="Distance to R-loop summit",ylab="RPM")
rm(signal)
```

##
##########################################################################################
## Fig 4. Other R-loop Hotspots in the Human Genome
##########################################################################################
## * Fig 4A, S5A & S5B
R-loop associated tRNAs by R-ChIP in HepG2 and K562 cell lines, and by DRIP-seq
```{r pie venn plot,fig.width=6.34,fig.height=4.21,fig.align = "center",message=F,warning=F}
# HEK293 R-ChIP
Total=360+81+73+92;Lowlyexpressed=81;Nonexpressed=73;Nonalignable=92;Expressed=360
Wrloop=280;WOrloop=Expressed-Wrloop;Other=Total-Wrloop-WOrloop
Samerloop=257;Other2 = Total-Samerloop
Oppositerloop=239;Overlap=216;SingleSame=Samerloop-Overlap;Singleoppo=Oppositerloop-Overlap;Other3=Total-Wrloop
# or K562 R-ChIP
Total=92+83+96+335;Nonalignable=92;Nonexpressed=83;Lowlyexpressed=96;Expressed=335
Wrloop=299;WOrloop=Expressed-Wrloop;Other=Total-Wrloop-WOrloop
Samerloop=280;Other2 = Total-Samerloop
Oppositerloop=269;Overlap=250;SingleSame=Samerloop-Overlap;Singleoppo=Oppositerloop-Overlap;Other3=Total-Wrloop
# plot
iniR=0.4 
colors=list(Other='white',Nonalignable='#e5f5e0',Nonexpressed='#a1d99b',
            Lowlyexpressed='#3182bd',Expressed='#fec44f',Wrloop='#fc9272',
            WOrloop='#9ecae1',Samerloop='#ffeda0',Oppositerloop='#fee0d2')
library('plotrix')
pie(1, radius=iniR, init.angle=90, col=c('white'), border = NA, labels='')
floating.pie(0,0,c(SingleSame,Overlap,Singleoppo,Other3),radius=5*iniR, startpos=pi/2,
             col=as.character(colors[c('Other','Oppositerloop','Oppositerloop','Other')]),border=NA)
floating.pie(0,0,c(Samerloop,Other2),radius=4*iniR,
             startpos=pi/2, col=as.character(colors[c('Samerloop','Other')]),border=NA)
floating.pie(0,0,c(Wrloop,WOrloop,Other),
             radius=3*iniR, startpos=pi/2, col=as.character(colors[c('Wrloop','WOrloop','Other')]),border=NA)
floating.pie(0,0, c(Expressed,Lowlyexpressed,Nonexpressed,Nonalignable), radius=2*iniR, startpos=pi/2,
             col=as.character(colors[c('Expressed','Lowlyexpressed','Nonexpressed','Nonalignable')]), 
             border = NA)
legend(0.6, 3*iniR, gsub("_"," ",names(colors)[-1]),
       col=as.character(colors[-1]), pch=19,bty='n', ncol=1)
# K562 DRIP
Total=92+83+96+335;Nonalignable=92;Nonexpressed=83;Lowlyexpressed=96;Expressed=335
Wrloop=39;WOrloop=Expressed-Wrloop;Other=Total-Wrloop-WOrloop
Unique=8;nonUnique=31;Other2 = Total-Unique-nonUnique
iniR=0.2
colors=list(Other='white',Nonalignable='#e5f5e0',Nonexpressed='#a1d99b',
            Lowlyexpressed='#3182bd',Expressed='#fec44f',Wrloop='#fc9272',
            WOrloop='#9ecae1',Unique='#ffeda0',nonUnique='#fee0d2')
library('plotrix')
pie(1, radius=iniR, init.angle=90, col=c('white'), border = NA, labels='')
floating.pie(0,0,c(Unique,nonUnique,Other2),radius=4*iniR,
             startpos=pi/2, col=as.character(colors[c('Unique',"nonUnique",'Other')]),border=NA)
floating.pie(0,0,c(Wrloop,WOrloop,Other),radius=3*iniR, startpos=pi/2,
             col=as.character(colors[c('Wrloop','WOrloop','Other')]),border=NA)
floating.pie(0,0, c(Expressed,Lowlyexpressed,Nonexpressed,Nonalignable), radius=2*iniR, startpos=pi/2,
             col=as.character(colors[c('Expressed','Lowlyexpressed','Nonexpressed','Nonalignable')]), 
             border = NA)
legend(0.6, 3*iniR, gsub("_"," ",names(colors)[-1]),
       col=as.character(colors[-1]), pch=19,bty='n', ncol=1)
```
##
##########################################################################################
## Fig 5. Elevated RNAPII pausing at TSS allows for increased R-loop formation
##########################################################################################
## Fig 5B & S6B
Pol-II decay and R-loop decay by ChIP-PCR and R-ChIP PCR, respectively
```{r R-ChIP PCR}
library("ggplot2")
library("reshape2")
## read data
polII  <- read.table("Data/ChIP-PCR/PolII.rep",sep="\t",row.names = 1,na.strings = "NA")
polII0 <- read.table("Data/ChIP-PCR/PolII.basal.rep",sep="\t",row.names = 1,na.strings = "NA")
rloop  <- read.table("Data/ChIP-PCR/Rloop.rep",sep="\t",row.names = 1,na.strings = "NA")
rloop0 <- read.table("Data/ChIP-PCR/Rloop.basal.rep",sep="\t",row.names = 1,na.strings = "NA")
## formulate data to data.frame format
# polII
colnames(polII) = seq(0,72,6)
polII      <- melt(t(polII))
polII$gene <- gsub("_[123]","",polII$Var2,perl=T)
genes      <- unique(polII$gene) # will keep the original order
polII      <- polII[order(polII[,"gene"],polII[,"Var1"]),]
# basal polII
basalpolII <- apply(matrix(polII0$V2,nrow=3),2,mean)
names(basalpolII) = genes
# rloop
colnames(rloop) = seq(0,72,6)
rloop      <- melt(t(rloop))
rloop$gene <- gsub("_[123]","",rloop$Var2,perl=T)
rloop      <- rloop[order(rloop[,"gene"],rloop[,"Var1"]),]
# basal rloop
basalrloop <- apply(matrix(rloop0$V2,nrow=3),2,mean)
names(basalrloop) = genes

## (optional) only focus on a subset
selected <- c("ATP5C1","ATP5G1","ATRAID","CLSPN","GADD45A","LENG8",
              "PIGH","PMPCA","SMG5","SNRPD2","THOP1")
polII <- polII[polII$gene%in%selected,]
rloop <- rloop[rloop$gene%in%selected,]
polII <- polII[order(polII[,"gene"],polII[,"Var1"]),]
rloop <- rloop[order(rloop[,"gene"],rloop[,"Var1"]),]
basalrloop <- basalrloop[selected]
basalpolII <- basalpolII[selected]

genes <- unique(rloop$gene) # ascending order
## plot
# polII - plot for each gene (1)
ggplot(polII,aes(x=Var1,y=value)) + 
  geom_point(alpha=0.3,size=0.3) +
  stat_smooth(method = "loess",span=1/3) +
  xlab("Time") + ylab("Input%") +
  facet_wrap(~gene,scales = "free")
# polII - summary plot for all gene (2)
ggplot(polII,aes(x=Var1,y=value)) + 
  xlab("Time After DRB Removal (min)") + ylab("Input %") +
  stat_smooth(data = polII,aes(x=Var1,y=value,col=gene),
              method = "loess",span=1/3,se=F,lwd=0.4,lty=2) +
  stat_smooth(method = "loess",span=1/3) +
  labs(col="Gene") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  scale_color_brewer(palette="Paired")
# polII - take noDRB treatment into consideration (3)
npolII <- polII
for(i in 1:nrow(npolII)){
  npolII[i,"value"] = polII[i,"value"]-basalpolII[polII[i,"gene"]]
}
ggplot(npolII,aes(x=Var1,y=value)) + 
  xlab("Time After DRD Removal (min)") + ylab("ΔInput %") +
  stat_smooth(data = npolII,aes(x=Var1,y=value,col=gene),
              method = "loess",span=1/3,se=F,lwd=0.4,lty=2) +
  stat_smooth(method = "loess",span=1/3) +
  labs(col="Gene") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  scale_color_brewer(palette="Paired")

p <- loess(npolII$value~npolII$Var1,span=1/3)
xl <- seq(0,72,72/1000)
out=predict(p,xl)
xl[which(out==min(abs(out)))]
xl[which(out==-min(abs(out)))]
# rloop - summary plot for all genes (2)
ggplot(rloop,aes(x=Var1,y=value)) + 
  xlab("Time After DRB Removal (min)") + ylab("Input%") +
  stat_smooth(data = rloop,aes(x=Var1,y=value,col=gene),
              method = "loess",span=1/3,se=F,lwd=0.4,lty=2) +
  stat_smooth(method = "loess",span=1/3) +
  labs(col="Gene") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  scale_color_brewer(palette="Paired")
# rloop - take noDRB treatment into consideration (3)
nrloop <- rloop
for(i in 1:nrow(nrloop)){nrloop[i,"value"] = rloop[i,"value"]-basalrloop[rloop[i,"gene"]]}
ggplot(nrloop,aes(x=Var1,y=value)) + 
  xlab("Time After DRB Removal (min)") + ylab("delta Input%") +
  stat_smooth(data = nrloop,aes(x=Var1,y=value,col=gene),
              method = "loess",span=1/3,se=F,lwd=0.4,lty=2) +
  stat_smooth(method = "loess",span=1/3) +
  labs(col="Gene") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  scale_color_brewer(palette="Paired")

r <- loess(nrloop$value~nrloop$Var1,span=1/3)
xl <- seq(0,72,72/1000)
out=predict(r,xl)
xl[which(out==min(abs(out)))]
xl[which(out==-min(abs(out)))]

# polII - in-depth analysis for each gene
# *** cycle ***
# decay: time of the firsr vallay
# half / halfbasal: halflife time w/ or w/o consideration of basal level
# cycle: period
# summit: number of summit
# time- or density- peak first/second/third: time or relative density of the 1st, 2nd or 3rd summits
cycle <- data.frame(genes=genes,decay=0,half=0,halfbasal=0,cycle=0,summit=0,
                    timepeakfirst=0,densitypeakfirst=0,
                    timepeaksecond=0,densitypeaksecond=0,
                    timepeakthird=0,densitypeakthird=0)
row.names(cycle) <- genes
for (i in 1:length(genes)){
  # first vallay & cycle
  gene = genes[i]
  lo <- loess(polII[polII$gene==gene,"value"]~polII[polII$gene==gene,"Var1"],span = 1/3)
  xl <- seq(0,72,72/1000)
  out = predict(lo,xl)
  cycle[gene,"decay"] <- xl[which(diff(diff(out)>0) == 1) + 1][1]
  
  half <- out - (max(out)-(max(out)-out[which(xl==cycle[gene,"decay"])])/2)
  for(k in 1:(length(half)-1)){
    if(half[k]>0 && half[k+1]<0){
      cycle[gene,"half"] = xl[k]
      break
    }
  }
  
  halfbasal <- out-(max(out)-(max(out)-basalpolII[gene])/2)
  for(j in 1:(length(halfbasal)-1)){
    if(halfbasal[j]>0 && halfbasal[j+1]<0){
      cycle[gene,"halfbasal"] = xl[j]
      break
    }
  }
  
  summit <- c(FALSE, diff(diff(out)>0)==-1)
  summit <- xl[summit]
  # include the peaks after first valley
  summit <- subset(summit,summit>cycle[gene,"decay"])
  cycle[gene,"summit"] = length(summit)
  cycle[gene,"cycle"] <- mean(summit[2:length(summit)] - summit[1:(length(summit)-1)])
    
  # the first peak after complete decay
  for(l in summit){
    if(l>cycle[gene,"decay"]){
      cycle[gene,"timepeakfirst"] <- l
      cycle[gene,"densitypeakfirst"] <- out[which(xl==l)] - basalpolII[gene]
#      cycle[gene,"densitypeakfirst"] <- (out[which(xl==l)] - basalpolII[gene]) / basalpolII[gene]
#      cycle[gene,"densitypeakfirst"] <- out[which(xl==l)] / basalpolII[gene]
      cycle[gene,"timepeaksecond"] <- subset(summit,summit>l)[1]
      cycle[gene,"densitypeaksecond"] <- out[which(xl==cycle[gene,"timepeaksecond"])] - basalpolII[gene]
      if ( length (subset(summit,summit>l)) >= 2){
        cycle[gene,"timepeakthird"] <- subset(summit,summit>l)[2]
        cycle[gene,"densitypeakthird"] <- out[which(xl==cycle[gene,"timepeakthird"])] - basalpolII[gene]
      }else{
        cycle[gene,"timepeakthird"] <- cycle[gene,"timepeaksecond"] + cycle[gene,"cycle"]
        cycle[gene,"densitypeakthird"] <- out[which(xl==cycle[gene,"timepeakthird"])] - basalpolII[gene]
      }
      break
    }
  }
  p<- ggplot(polII[polII$gene==gene,],aes(x=Var1,y=value)) + 
    geom_point(alpha =0.3) +
    stat_smooth(method = "loess",span=1/3) +
    xlab("Time") + ylab("Pol II") +
    ggtitle(gene) +
    geom_hline(yintercept=basalpolII[gene]) +
    annotate("text", x=55, y=max(polII[polII$gene==gene,"value"])-0.2,
             label=paste("Period = ",cycle[gene,"cycle"],sep="")) +
    annotate("text", x=55, y=max(polII[polII$gene==gene,"value"])-0.1,
             label=paste("First = ",cycle[gene,"decay"],sep=""))
  print(p)
}

# rloop - in-depth analysis for each gene
# *** decay ***
# a: parameter estimate
# half / halfbasal: halflife time w/ or w/o consideration of basal level
# a- or b- pvalue: significance
# first- / second- or third polIIpeak: R-loop density at time point of PolII summit
decay <- data.frame(genes=genes,a=0,half=0,halfbasal=0,apvalue=0,bpvalue=0,
                    firstpolIIpeak=0,secondpolIIpeak=0,thirdpolIIpeak=0)
row.names(decay) <- genes
for (i in 1:length(genes)){
  gene = genes[i]
  expr = rloop[rloop$gene==gene,"value"]
  time = rloop[rloop$gene==gene,"Var1"]
  
  fit <- nls(expr ~ exp(a * time) + b, start = list(a=0 ,b=0 ))
  summary(fit)
  para <- summary(fit); a <- para$parameters[1,1]; b = para$parameters[2,1]
  point <- predict(fit,list(time=time))
  
  plot(time,expr,main=gene,xlab="Time",ylab="R-Loop Signal",ylim=c(max(0,basalrloop[gene]-0.05),max(expr)))
  lines(time,point,col="red")
  abline(basalrloop[gene],0,lty=2,col="red")
  text(45,max(point)-0.1,paste("y=e^(",round(a,4),"x) + ",round(b,4),sep=""))
  text(45,max(point)-0.3,paste("t1/2 = ",round(log((1+b)-0.5-b)/a,1),sep=""))
  
  decay[gene,"a"] = a
  decay[gene,"half"] = round(log((1+b)-0.5-b)/a,1)
  decay[gene,"halfbasal"] = round(log(1+b-(1+b-basalrloop[gene])/2-b)/a,1)
  decay[gene,"apvalue"] =  para$parameters[1,4]
  decay[gene,"bpvalue"] =  para$parameters[2,4]
  decay[gene,"firstpolIIpeak"] = predict(fit,list(time=cycle[gene,"timepeakfirst"])) - basalrloop[gene]
  decay[gene,"secondpolIIpeak"] = predict(fit,list(time=cycle[gene,"timepeaksecond"])) - basalrloop[gene]
  decay[gene,"thirdpolIIpeak"] = predict(fit,list(time=cycle[gene,"timepeakthird"])) - basalrloop[gene]
}

# Comparason of the half basal
halftime <- as.data.frame(cbind(c(cycle$halfbasal,decay$halfbasal),
                                rep(row.names(cycle),2),
                                rep(c("PolII","R Loop"),each=11)))
names(halftime) <- c("Half","Gene","Class")
ggplot(halftime,aes(x=Class,y=as.numeric(as.vector(Half)),col=Class)) + 
  geom_violin() +
  geom_boxplot(width=0.2) +
  geom_point(position=position_dodge(0.8)) +
  geom_line(aes(x=Class,y=as.numeric(as.vector(Half)),group=Gene),lwd=0.3,linetype=2,col="grey") + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  ylab("Half Life (min)") + xlab("") + 
  theme(axis.text.x = element_blank())
```
## Fig 5E
Signal intensity distribution of overall R-ChIP in response to DRB treatment [DRB(+)] and removal (Post-DRB). 
```{r DRB Rloop Signal,fig.width=2.18,fig.height=2.36,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/DRBrelated/ActiveTSS-Rloop")
names(data)[19:21] = c("DRB(-)","DRB(+)","Post-DRB")
data <- melt(data[,c(19:21)])
ggplot(data,aes(x=variable,y=value,fill=variable)) + 
  geom_boxplot(notch = T,outlier.colour = NA) +
  my_clear_theme +
  coord_cartesian(ylim=c(0,3)) +
  xlab("") + ylab("R-loop Signal") + guides(fill=F)
```
## Fig 5F
Signal intensity distribution of overall GRO-seq in response to DRB treatment [DRB(+)] and removal (Post-DRB). 
```{r GRO-seq Signal,fig.width=2.18,fig.height=2.36,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/DRBrelated/ActiveTSS-Rloop")
names(data)[c(7,9,11)] = c("DRB(-)","DRB(+)","Post-DRB")
data <- melt(data[,c(7,9,11)])
ggplot(data,aes(x=variable,y=value,fill=variable)) + 
  geom_boxplot(notch = T,outlier.colour = NA) +
  my_clear_theme +
  coord_cartesian(ylim=c(0,100)) +
  xlab("") + ylab("GRO-seq Signal") + guides(fill=F)
```
## Fig S6A
Reproducibility of R-ChIP and GRO-seq
```{r R-ChIP Rep DRB HEK293,fig.width=6.46,fig.height=2.53,fig.align = "center",message=F,warning=F}
noDRB1<- read.table("Data/RChIP_Rep/chipseq.D210N_noDRB_V5_rep1.ReadsCoverage",sep="\t")
noDRB2<- read.table("Data/RChIP_Rep/chipseq.D210N_noDRB_V5_rep2.ReadsCoverage",sep="\t")
DRB2h1<- read.table("Data/RChIP_Rep/chipseq.D210N_DRB2h_V5_rep1.ReadsCoverage",sep="\t")
DRB2h2<- read.table("Data/RChIP_Rep/chipseq.D210N_DRB2h_V5_rep2.ReadsCoverage",sep="\t")
post1 <- read.table("Data/RChIP_Rep/chipseq.D210N_postDRB30min_V5_rep1.ReadsCoverage",sep="\t")
post2 <- read.table("Data/RChIP_Rep/chipseq.D210N_postDRB30min_V5_rep2.ReadsCoverage",sep="\t")

par(mfrow=c(1,3))
smoothScatter(log2(noDRB1[,4]),log2(noDRB2[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="noDRB Rep 1 (Log2 Tag Counts)",ylab="noDRB Rep 2 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(noDRB1[,4],noDRB2[,4]),2),sep=""))
smoothScatter(log2(DRB2h1[,4]),log2(DRB2h2[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="DRB 2h Rep 1 (Log2 Tag Counts)",ylab="DRB2h Rep 2 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(DRB2h1[,4],DRB2h2[,4]),2),sep=""))
smoothScatter(log2(post1[,4]),log2(post2[,4]),xlim=c(0,12),ylim=c(0,12),
              xlab="postDRB Rep 2 (Log2 Tag Counts)",ylab="postDRB Rep 2 (Log2 Tag Counts)")
text(3,10,paste("R = ",round(cor(post1[,4],post2[,4]),2),sep=""))
```
##
##########################################################################################
## Fig 6. Correlation of R-loop Levels with G/C Content and RNAPII Pausing at TSSs
##########################################################################################
## Fig 6A
R-loops were divided into three groups according to the signal intensity (green) and the sequence features associated with each group are shown.
Four datasets were tested: HEK293 (3 replicates, MC paper), HEK293.TSS (3 replicates), DRBrelated, DRBrelated.TSS
```{r Rloop Sig vs Sequence,fig.width=2.71,fig.height=2.49,fig.align = "center",message=F,warning=F}
library(reshape2)
library(ggplot2)
# hig vs med vs low
allseq <- read.table("Data_MC/Sequence/HEK293.hig.nucleotides",sep="\t",header = T)
allseq <- melt(allseq,id.vars = "Pos")
ggplot(data = allseq,aes(x=Pos-400,y=value,col=variable)) + 
  geom_line() + 
#  stat_smooth(method = "loess", span=1/3) + 
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  labs(col="") + ylim(-0.05,0.7) + xlab("Distance to R-loop summit (nt)") +
  ylab("Base composition or GC skew")
```
```{r Rloop Sig (H-M-L),fig.width=2.52,fig.height=3.67,fig.align = "center",message=F,warning=F}
signal <- read.table("Data_MC/Signal/High.signal")
input  <- read.table("Data_MC/Signal/High.input.signal")
plot(seq(-400,400,by=1),apply(signal,2,mean) - apply(input,2,mean),type='l',
     xlab="Distance to R-loop summit",ylab="RPM",ylim=c(0,2))
rm(signal,input)
```
## Fig 6B-E
R-loop ~ GC + local expression, local expression ~ R-loop + GC
```{r colorgrid,fig.width=2.94,fig.height=2.27,fig.align = "center",message=F,warning=F}
data <- read.table("Data_MC/DRBrelated/ActiveTSS-Rloop")
# Rloop ~ GC + Expression
BinnedHeatmapPlot(data,z1=19,x1=23,y1=8,bin=5)
# Rloop ~ GC + PP (promoter proximal region expression)
BinnedHeatmapPlot(data,z1=19,x1=23,y1=7,bin=5)
# FC Rloop ~ GC + FC PP (induction)
BinnedHeatmapPlot(data,z1=20,x1=23,y1=9,z2=19,y2=7,bin=5)
# FC Rloop ~ GC + FC PP (decay)
data[,24] = data[,9] - data[,11]; data[data$V24<0,24] = 0
data[,25] = data[,20] - data[,21]; data[data$V25<0,25] = 0
data <- subset(data,data$V20!=0)
BinnedHeatmapPlot(data,z1=25,x1=23,y1=24,y2=9,z2=20,bin=5)

# R-loop retard the following PolII
BinnedErrorBarPlot(data,x1=19,y1=9,y2=7,xlab="Steady State R-loop\nSignal",
                   ylab="Increased GRO-seq\nSignal", bin=12)
BinnedErrorBarPlot(data,x1=20,y1=24,y2=9,xlab="R-loop Signal\n(DRB2h)",
                   ylab="Decreased GRO-seq\nSignal", bin=12)
```
##
##########################################################################################
## Fig 7. Requirement of a Free RNA End for Promoting R-loop Formation
## and Proposed Model for R-loop Initiation and Elongation
##########################################################################################
## Fig 7A
The distance distribution of potential free RNA ends (TSS, TTS, TSS and TTS of tRNA/eRNAs) relative to the summits of R-ChIP mapped R-loops
```{r,Dis2FreeEnd,fig.width=4.63,fig.height=3.82,fig.align = "center",message=F,warning=F}
# summit dis to free end / peak dis to free end
data <- read.table("Data_MC/Freeend/summitdis2freeend")
data <- subset(data,data$V1 < 3000)
hist(data$V1,xlim=c(0,3000),breaks = c(seq(-1,2999,by=200),max(data$V1)),
     freq = T, main="" , col="lightblue", xlab="Distance to Free End")
```
##
##########################################################################################
## Others (not shown in MC paper)
##########################################################################################
## *** tRNA R-loops (related to Fig 4A, B & C)
Differences between Rloop(+) and Rloop(-) tRNAs
```{r cmpr w and w/o group,fig.width=2.94,fig.height=2.27,fig.align = "center",message=F,warning=F}
library("ggplot2")
data <- read.table("Data/tRNA/HEK293.RChIP.uniqtRNAs.bed6",sep="\t",na.strings = "NA")
data <- subset(data,data$V10>=10)
data$V15 = "w"
wo <- which(is.na(data$V12) & is.na(data$V14))
data[wo,"V15"] = "wo"
data$V15 <- factor(data$V15,levels=c("wo","w"))

pvalue = format(wilcox.test(data[data$V15=="wo",8],data[data$V15=="w",8])$p.value,scientific = T)
ggplot(data,aes(x=V15,y=V8,fill=V15)) + geom_violin() + ylim(0,1.1) +
  geom_boxplot(width=0.1,outlier.colour = "NA",fill="white") +
  ylab("Alignable Proportion") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA))

pvalue = format(wilcox.test(data[data$V15=="wo",10],data[data$V15=="w",10])$p.value,scientific = T)
ggplot(data,aes(x=V15,y=V10,fill=V15)) + geom_violin() +
  geom_boxplot(width=0.1,outlier.colour = "NA",fill="white") +
  ylab("GRO-Seq (RPKM)") + ylim(0,170) +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA))

```
Correlation of Rloop signal intensity between sense- and antisense- Rloops
```{r cor between sense & antisense R-loop,fig.width=2.49,fig.height=2.27,fig.align = "center",message=F,warning=F}
data <- read.table("Data/tRNA/HEK293.RChIP.uniqtRNAs.bed6",sep="\t",na.strings = "NA")
data <- subset(data,data$V10>=10)
SA <- data[!is.na(data$V11) & !is.na(data$V13),c(12,14)]
names(SA) <- c("V1","V2")
bk <- as.double(quantile(SA$V1,probs=seq(0,1,0.1)))
SA$V3 <- .bincode(SA$V1,breaks = bk,include.lowest = T)

me_x <- c(); q5_x <- c(); q95_x <- c();
me_y <- c(); q5_y <- c(); q95_y <- c();
for(i in 1:10){
  me_x[i]  <- median(SA[SA$V3==i,"V1"])
  q5_x[i]  <- quantile(SA[SA$V3==i,"V1"],probs = 0.25)
  q95_x[i] <- quantile(SA[SA$V3==i,"V1"],probs = 0.75)
  me_y[i]  <- median(SA[SA$V3==i,"V2"])
  q5_y[i]  <- quantile(SA[SA$V3==i,"V2"],probs = 0.25)
  q95_y[i] <- quantile(SA[SA$V3==i,"V2"],probs = 0.75)
}

mat <- as.data.frame(cbind(me_x,q5_x,q95_x,me_y,q5_y,q95_y))
library(ggplot2)
ggplot(mat,aes(x=me_x,y=me_y)) +
  geom_errorbar(aes(ymin = q5_y,ymax = q95_y),col="blue") + 
  geom_errorbarh(aes(xmin = q5_x,xmax = q95_x),col="blue") +
  geom_point(col="red") +
  xlab("Sense R-Loop Signal") + 
  ylab("Antisense R-Loop Signal") +
  theme(panel.grid.major = element_line(colour="NA"),
        panel.grid.minor = element_line(colour="NA"),
        panel.background = element_rect(fill="NA"),
        panel.border = element_rect(colour="black", fill=NA)) +
  theme(legend.position=c(0,1), legend.justification=c(0,1))

cor(SA$V1,SA$V2)
```