Show the code

# library
library(BSgenome.Umaydis.ENSMBL.UM1) # forged ustilago genome
library(ggplot2)
library(BindingSiteFinder)
library(rtracklayer)
library(ComplexHeatmap)
library(GenomicFeatures)
library(forcats)
library(tidyr)
library(dplyr)
library(tidyverse)
library(GenomicRanges)
library(magick)
library(magrittr)
library(gridExtra)
library(IRanges)
library(Biostrings)
library(ggpp)
library(gginnards)
library(ggrepel)
library(ggpubr)
library(ggforce)
library(ggrastr)
library(viridis)
library(reshape2)
library(gprofiler2)
library(ggsci)
library(ggh4x)
library(ggplotify)
library(gridExtra)
library(circlize)
library(EnrichedHeatmap)
library(UpSetR)
library(kableExtra)
library(cowplot)
library(rstatix)
library(beeswarm)
library(clusterProfiler)
library(ggseqlogo)
library(tidyHeatmap)
library(paletteer)
library(ggvenn)
library(colorspace)
library(ggpointdensity)
library(lookup)
library(rstatix)
library(ggrepel)

Show the code

load("C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/Kathis_lab/hyperTRIBE/hyperTRIBE.rds")

Show the code

# reproducible editing sites in Khd4-Ada-Gfp and control-Ada
hyper.rep

GRangesList object of length 2:
$khd4
GRanges object with 377 ranges and 6 metadata columns:
        seqnames    ranges strand |          id     gene.id no.reads.rep1
           <Rle> <IRanges>  <Rle> | <character> <character>     <integer>
    [1]       10    122633      - |   10_122633  UMAG_03680            52
    [2]       10    166332      + |   10_166332  UMAG_03696            43
    [3]       10    195059      + |   10_195059  UMAG_10492            36
    [4]       10    195431      + |   10_195431  UMAG_10492            34
    [5]       10    198031      + |   10_198031  UMAG_10492            47
    ...      ...       ...    ... .         ...         ...           ...
  [373]        9    541920      - |    9_541920  UMAG_03569            38
  [374]        9    542067      - |    9_542067  UMAG_03569            60
  [375]        9    614585      + |    9_614585  UMAG_03597            88
  [376]        9    666949      - |    9_666949  UMAG_03615           117
  [377]        9    698817      - |    9_698817  UMAG_03623            23
        perc.editing.rep1 no.reads.rep2 perc.editing.rep2
                <numeric>     <integer>         <numeric>
    [1]               5.8            49               8.2
    [2]               9.3            41               9.8
    [3]              11.1            21               9.5
    [4]              14.7            41               9.8
    [5]              10.6            48              10.4
    ...               ...           ...               ...
  [373]              15.8            38               7.9
  [374]              10.0            57               8.8
  [375]               8.0            97               6.2
  [376]               5.1           113               9.7
  [377]               8.7            25               8.0
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

$ctrl
GRanges object with 1129 ranges and 6 metadata columns:
         seqnames    ranges strand |          id     gene.id no.reads.rep1
            <Rle> <IRanges>  <Rle> | <character> <character>     <integer>
     [1]       10    109281      - |   10_109281  UMAG_03674            34
     [2]       10    125437      + |   10_125437  UMAG_03681            41
     [3]       10    142066      - |   10_142066  UMAG_10487            58
     [4]       10    169565      + |   10_169565  UMAG_10490            26
     [5]       10    171252      - |   10_171252  UMAG_03698            47
     ...      ...       ...    ... .         ...         ...           ...
  [1125]        9    652510      - |    9_652510  UMAG_03610            33
  [1126]        9    666949      - |    9_666949  UMAG_03615           125
  [1127]        9    711406      - |    9_711406  UMAG_03627            53
  [1128]        9     74604      + |     9_74604  UMAG_03418            49
  [1129]        9     83223      + |     9_83223  UMAG_11071            85
         perc.editing.rep1 no.reads.rep2 perc.editing.rep2
                 <numeric>     <integer>         <numeric>
     [1]               5.9            22              13.6
     [2]               7.3            38               5.3
     [3]              12.1            59              10.2
     [4]              11.5            26              11.5
     [5]              19.1            40               5.0
     ...               ...           ...               ...
  [1125]               9.1            37              10.8
  [1126]              11.2           118               6.8
  [1127]              11.3            44               6.8
  [1128]               8.2            58              12.1
  [1129]               5.9            95               5.3
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

Show the code

# Venn diaram
venn.list = lapply(hyper.rep, function(hyper){
  return(hyper$id)
})

names(venn.list) = c("Khd4-Ada-Gfp", "control-Ada")
# FigS6A and B

ggvenn(venn.list, names(venn.list), 
       fill_color = c("dodgerblue3", ada.col),  
       set_name_size = 6) +
  ggtitle("Overlap between Khd4-Ada-Gfp and control-Ada")+
  theme(plot.title = element_text(size = 20, hjust = 0.5, vjust = 3))

Show the code

# Identifying the RNA sequence motif

################################################################################
#AUACCC
################################################################################

# get gene sequences
gene.seqs = getSeq(Umaydis, genes.GR)

# find AUACCC motifs
motif.auaccc = DNAString("ATACCC")

# find all AUACCC motifs in the genome
auaccc.gr = vmatchPattern(motif.auaccc, Umaydis)

# count motifs in genes
genes.GR$AUACCC = countOverlaps(genes.GR, auaccc.gr)
#table(genes.GR$AUACCC)

# assign motifs to genes
ov = findOverlaps(auaccc.gr, genes.GR)

# check whether motifs overlap with more than one gene
table(duplicated(from(ov)))


FALSE  TRUE 
 1975    18

Show the code

# keep only motifs in genes & duplicate motifs that overlap with 2 genes, so they count for both genes
auaccc.gr = auaccc.gr[from(ov)]
auaccc.gr$gene.id = genes.GR$gene_id[to(ov)]

# count motifs in transcript regions
#findOverlaps(hits, tx.regions) %>% to(.) %>% table(.)
#findOverlaps(auaccc.gr, tx.regions, select="first") %>% names(tx.regions)[.] %>% table(.)


################################################################################
#AGAUCU
################################################################################

# find AGAUCU motifs
motif.ctrl = DNAString("AGATCT")
agaucu.gr = vmatchPattern(motif.ctrl, Umaydis)
genes.GR$AGAUCU = countOverlaps(genes.GR, agaucu.gr)
table(genes.GR$AGAUCU)


   0    1    2    3    4    5    6    7    8   11 
4000 1990  672  189   50   19    7    2    1    1

Show the code

# assign motifs to genes
ov = findOverlaps(agaucu.gr, genes.GR)

# check whether motifs overlap with more than one gene
#table(duplicated(from(ov)))

# keep only motifs in genes & duplicate motifs that overlap with 2 genes, so they count for both genes
agaucu.gr = agaucu.gr[from(ov)]
agaucu.gr$gene.id = genes.GR$gene_id[to(ov)]

# count motifs in transcript regions
#findOverlaps(hits, tx.regions) %>% to(.) %>% table(.)
#findOverlaps(agaucu.gr, tx.regions, select="first") %>% names(tx.regions)[.] %>% table(.)

0.1 De novo motif discovery analysis on sequences carrying unique Khd4–Ada–Gfp editing sites.

XSTREME analysis was done using the random genome sequences of lenght 501-nt as a background and unique-Khd4-Ada-Gfp editing sites (501-nt) as an input.

First, the Khd4-Ada-Gfp editing sites were extended on each side by 250nt to make a 501-nt window.

Show the code

# keep only id column
hyper.spe = lapply(hyper.rep, function (x){
  mcols(x) = x[,c(1)]
  names(x) = x$id
  return(x)
})

# split the reproducible editing sites that are unique to Khd4 and Control , and shared
hyper.specific = list(
  Khd4.specific  = hyper.spe$khd4[!(hyper.spe$khd4$id %in% hyper.spe$ctrl$id)],
  shared = hyper.spe$khd4[hyper.spe$khd4$id %in% hyper.spe$ctrl$id],
  ctrl.specific  = hyper.spe$ctrl[!(hyper.spe$ctrl$id %in% hyper.spe$khd4$id)]
)

# as a GRanges list
hyper.specific = as(hyper.specific, "GRangesList")

# extend the window by 250 bp on both sides
hyper_250 = lapply(hyper.specific, function(x){
  x = x+250 # window 250 bp
})

# get fastafiles
# get sequence
khd4.specific = Biostrings::getSeq(Umaydis, hyper_250$Khd4.specific)
khd4.ctrl.shared = Biostrings::getSeq(Umaydis, hyper_250$shared)
ctrl.specific = Biostrings::getSeq(Umaydis, hyper_250$ctrl.specific)

# export the fasta files
output = "C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/TRIBE_data/Editing_events_all/Khd4/bedgraphs/MEME/fasta/"

writeXStringSet(khd4.specific, filepath = paste0(output,"khd4.specific_250.fasta"))
writeXStringSet(khd4.ctrl.shared, filepath = paste0(output,"khd4.ctrl.shared_250.fasta"))
writeXStringSet(ctrl.specific, filepath = paste0(output,"ctrl.specific_250.fasta"))

Scatter plot comparing the relative enrichment ratio of enriched motifs with the percentage of each motif in the tested sequences.

Show the code

# making a scatter plot
# unique to Khd4-Ada-Gfp
Motif.frame.khd4 = data.frame(Motif = 
c("HMHHKAUACCC","HMKAUACCCVH","UGSAGGCGCURSWGS","UHCCUCUYCHCCMUC","CAGCARCARCADCAG","CCUCKAAACSNCGC","CYAAKSMUGCUYCUR","CYCCCW","CUUGAYUUGUWWYCC","YYUUGUAYUUGUA","UYGCUGCUCUC", "ACMUCUCCAAGVYCA","CUCCUGAUCCUGCCU", "AGCGSUCRUUUGASC", "CCUCUKU", "CAACGCUCCAC","CCAGCAACAGCARSM", "GCAGAUGRUAYCCAK", "RASAAGMVSAWCAAG", "AVUGGUGCUGU"),
absolute = c(60.9, 59.3, 16.0, 65.8, 13.6, 37.4, 11.1, 37.9, 2.5, 6.2, 18.5, 6.2, 5.3, 3.7, 16.0, 50.6, 29.2, 3.7, 38.3, 11.1),
EnR = c(9.40,9.60,8.35,1.46,4.01,1.80,3.30,1.60,37.98,4.82,2.08, 4.34, 5.06, 4.93, 1.81, 1.27, 1.43, 3.88   , 1.10, 1.05),
set=rep("khd4",20))

pl.khd4 = ggplot(Motif.frame.khd4, aes(x=absolute, y=EnR, color=set))+
  geom_point(alpha=1, size=2)+
  scale_color_manual(values = khd4.col)+
  coord_cartesian(ylim=c(0,10))+
  labs(subtitle= "Khd4-Ada-Gfp",
       x="sequences with motif (%)",
       y = "motif enrichment ratio")+
  theme_paper()+
  theme(legend.position = "none")

pl.khd4

Make the sequence logo of the most enriched two motifs

Show the code

################################################################################
# De novo discoevered motif on Khd4-Ada-Gfp edited transcripts
################################################################################

# motif logo with u in it
# Motif-1 TP% = 60.9   ENR-Ratio = 9.40 HMHHKAUACCC
A= c(0.220238,0.380952,0.369048,0.392857,0.041667,1000000,0.000000,1000000,0.000000,0.000000,0.059524
)
C= c(0.440476,0.321429,0.351190,0.261905,0.107143,0.000000,0.000000,0.000000,0.976190,1000000,0.940476
)
G = c(0.107143,0.113095,0.029762,0.119048,0.333333,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000
)
U = c(0.232143,0.184524,0.250000,0.226190,0.517857,0.000000,1000000,0.000000,0.023810,0.000000,0.000000
)
 
pwDF1bU = rbind(A,C,G,U)


# MoUif-2 UP% = 59.3   ENR-RaUio = 9.60 HMKAUACCCVH

A= c(0.366890,0.432382,0.026357,0.999729,0.000082,0.999729,0.019779,0.000082,0.050456,0.234389,0.213404
)
C= c( 0.304485,0.245517,0.122745,0.000103,0.000103,0.000103,0.927513,0.999751,0.949376,0.374600,0.306852
)
G = c(0.032919,0.116138,0.316355,0.000090,0.006656,0.000090,0.006656,0.000090,0.000090,0.211350,0.159894
)
U = c(0.295706,0.205963,0.534542,0.000078,0.993160,0.000078,0.046052,0.000078,0.000078,0.179661,0.319849
)
 
pwDF2bU = rbind(A,C,G,U)

# make ggseqlogo
m1 = ggseqlogo(pwDF1bU, col_scheme="nucleotide", font="helvetica_bold")

m2 = ggseqlogo(pwDF2bU, col_scheme="nucleotide", font="helvetica_bold")

pl = ggarrange(m1,m2,
              labels = LETTERS[1:2],
              nrow = 1)

pl

Next I want to select mRNA targets that are specific to Khd4-Ada-Gfp and control-Ada data

Show the code

# Khd4.spe
Khd4.spe= hyper.rep$khd4[!(hyper.rep$khd4$gene.id %in% hyper.rep$ctrl$gene.id)]

# mRNAs specific to control
ctrl.spe.tx = hyper.rep$ctrl[!(hyper.rep$ctrl$gene.id %in% hyper.rep$khd4$gene.id)]

# make a list with mRNAs specific to each data set
# list
hyper.spe = list(
  Khd4_Ada_Gfp = Khd4.spe,
  control_Ada = ctrl.spe.tx
) %>% as(., "GRangesList")

0.2 Distribution of editing sites from the nearest AUACCC motif

Only editing sites on mRNAs specific to each dataset were considered

Show the code

pf = lapply(names(hyper.spe), function(i){
  #i = "khd4"
  hyper = hyper.spe[[i]]
  head(hyper)
  
  # calculate distance to nearest AUACCC motif
  dists.auacc = distanceToNearest(hyper, auaccc.gr)
  
  # keep only pairs on same genes
  sel = hyper$gene.id[from(dists.auacc)] == auaccc.gr$gene.id[to(dists.auacc)]
  table(sel)
  dists.auacc= subset(dists.auacc, sel)
  mcols(dists.auacc)$motif = "AUACCC"
  
  # calculate distance to nearest AGAUCU motif
  dists.agaucu = distanceToNearest(hyper, agaucu.gr)
  
  # keep only pairs on same genes
  sel = hyper$gene.id[from(dists.agaucu)] == agaucu.gr$gene.id[to(dists.agaucu)]
  table(sel)
  dists.agaucu= subset(dists.agaucu, sel)
  mcols(dists.agaucu)$motif = "AGAUCU"
  
  # get stats of distances
  summary(mcols(dists.auacc)$distance)
  summary(mcols(dists.agaucu)$distance)
 
  pf = rbind(mcols(dists.auacc), mcols(dists.agaucu)) %>% as.data.frame()
  pf$set = i
  
  return(pf)
})

pf = do.call(rbind, pf)

# FigS7B

# Histogram
pl4 = ggplot(pf, aes(x=distance, fill=set)) +
  geom_histogram( binwidth = 150, aes(y = ..density..), position="dodge") +
  coord_cartesian(xlim=c(0, 2500)) +
  scale_fill_manual(values=c(Khd4_Ada_Gfp = "dodgerblue4", control_Ada = "grey60")) +
  labs(x="distance of editing sites to the nearest motif (nt)",
       y="motif density")+
  theme_bw() +
  facet_wrap(~factor(motif, levels=c("AUACCC", "AGAUCU")))+myTheme1

pl4

Show the code

# median distance between editing sites and the motif
#Khd4-Ada-Gfp
Khd4_auaccc= filter(pf, set == "Khd4_Ada_Gfp", motif == "AUACCC") 
med_K_auaccc = median(Khd4_auaccc$distance)


#Control-Ada
Ctrl_auaccc= filter(pf, set == "control_Ada", motif == "AUACCC") 
med_C_auaccc = median(Ctrl_auaccc$distance)


data.frame(Khd4_Ada_Gfp = med_K_auaccc,
           control_Ada = med_C_auaccc) %>% 
  kbl(caption = "median distance between auaccc and editing events",html_font = "Cambria", bold = T, color="black") %>% 
  kable_styling(bootstrap_options = "striped", "hover", full_width = F, html_font = "Cambria") %>%
  column_spec(1, bold = T)

median distance between auaccc and editing events
Khd4_Ada_Gfp	control_Ada
120	863

0.3 Studying the enrichment of motifs in mRNAs specific to each datasets

Here, we consider only the editing events on Khd4-bound mRNAs and on mRNAs specific to control-Ada. Khd4-bound mRNAs- All mRNAs that carry reproducible Khd4-Ada-Gfp editing sites with atleast one AUACCC motif. Control-Ada mRNAs- mRNAs that have reproducible control-Ada editing sites but no Khd4-Ada-Gfp sites

Show the code

# Khd4-bound mRNAs
Khd4_boundRNA= hyper.rep$khd4[(hyper.rep$khd4$gene.id %in% auaccc.gr$gene.id)]

# mRNAs specific to control
ctrl.spe.tx = hyper.rep$ctrl[!(hyper.rep$ctrl$gene.id %in% hyper.rep$khd4$gene.id)]

#check if the numbers are similar to the venn diagram
#n_distinct(khd4.spe.tx$gene.id)

# make a list with mRNAs specific to each data set
Khd4.spe= hyper.rep$khd4[!(hyper.rep$khd4$gene.id %in% hyper.rep$ctrl$gene.id)]


# list
hyper.rep.txspe = list(
  khd4 = Khd4_boundRNA,
  ctrl = ctrl.spe.tx
) %>% as(., "GRangesList")

1 Analyis of the high-confidence targets

Show the code

# only Khd4
khd4.rep.edit=hyper.rep$khd4
motifs=genes.GR
khd4.rep.edit$AUACCC=lookup(khd4.rep.edit$gene.id, motifs$gene_id,motifs$AUACCC)

#granges of high confidence Khd4
high.conf.khd4=subset(khd4.rep.edit,AUACCC>0)

high.conf = high.conf.khd4

1.1 Transcript region enrichment of motif

Since Khd4 binds AUACCC specifically, we first want to check the transcript region that the AUACCC motif is enriched in high-confidence mRNA targets. Besides using the AGAUCU motif as a control, I also want to use the all hyphae expressed AUACCC-containing mRNAs as a putative targets

Furthermore, I will use the already available AUACCC and AGAUCU motif position in U. maydis genome, as these already include the transcript region locations of the motif.

Show the code

# load DGE data to identify hyphae expressed genes
# Use the lookup function to identify geneIDs overlapping with the high confident targets  
Fi_wtvskd_DEG$high.conf_targets=lookup(Fi_wtvskd_DEG$gene.ID,high.conf$gene.id, high.conf$gene.id )

# auaccc motif
Fi_wtvskd_DEG$auaccc=lookup(Fi_wtvskd_DEG$gene.ID,auaccc.gr$gene.id, auaccc.gr$gene.id )

# agaucu motif
Fi_wtvskd_DEG$agaucu=lookup(Fi_wtvskd_DEG$gene.ID,agaucu.gr$gene.id, agaucu.gr$gene.id )

# only expressed genes
Fi_wtvskd_BM = Fi_wtvskd_DEG[Fi_wtvskd_DEG$baseMean>10,]

Show the code

auaccc_DGE = Fi_wtvskd_BM[!is.na(Fi_wtvskd_BM$auaccc),] 
  
agaucu_DGE = Fi_wtvskd_BM[!is.na(Fi_wtvskd_BM$agaucu),]

Show the code

# AUACCC
auaccc.mot=read.csv("C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/TRIBE_data/Editing_events_all/Khd4/bedgraphs/AUACCCmotif.csv")

auaccc.mot$high.conf.targets = lookup::lookup(auaccc.mot$Gene_ID,high.conf$gene.id,high.conf$gene.id)

auaccc.mot$put.targets = lookup::lookup(auaccc.mot$Gene_ID, auaccc_DGE$gene.ID, auaccc_DGE$gene.ID)

# AGAUCU
agaucu.mot=read.csv("C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/TRIBE_data/Editing_events_all/Khd4/bedgraphs/AGAUCUmotif.csv")

agaucu.mot$high.conf.targets = lookup::lookup(agaucu.mot$GeneID,high.conf$gene.id,high.conf$gene.id)

agaucu.mot$put.targets = lookup::lookup(agaucu.mot$GeneID, auaccc_DGE$gene.ID, auaccc_DGE$gene.ID)

Show the code

# highconfident targets (auaccc motif)
df_high_conf_auaccc = auaccc.mot[!is.na(auaccc.mot$high.conf.targets),] %>% .[,"Position"] %>% table(.) %>% as.data.frame %>% set_colnames(.,c("variable", "Freq")) %>% mutate(motif= "AUACCC", set = "Khd4_bound_mRNA")

# highconfident targets (agaucu motif)
df_high_conf_agaucuc = agaucu.mot[!is.na(agaucu.mot$high.conf.targets),] %>% .[,"Position"] %>% table(.) %>% as.data.frame %>% set_colnames(.,c("variable", "Freq")) %>% mutate(motif= "AGAUCU", set = "Khd4_bound_mRNA")

# putative targets (auaccc motif)
df_put_auaccc = auaccc.mot[!is.na(auaccc.mot$put.targets),] %>% .[,"Position"] %>% table(.) %>% as.data.frame %>% set_colnames(.,c("variable", "Freq")) %>% mutate(motif= "AUACCC", set = "putative_targets")

# putative targets (agaucu motif)
df_put_agaucuc = agaucu.mot[!is.na(agaucu.mot$put.targets),] %>% .[,"Position"] %>% table(.) %>% as.data.frame %>% set_colnames(.,c("variable", "Freq")) %>% mutate(motif= "AGAUCU", set = "putative_targets")

df = rbind(df_high_conf_auaccc,df_high_conf_agaucuc,df_put_auaccc,df_put_agaucuc)

ggplot(df, aes(x=factor(set, levels=c("Khd4_bound_mRNA", "putative_targets")), y = Freq, fill = factor(variable, levels=c("5UTR", "ORF", "3UTR"))))+
  geom_bar(position = position_fill(), stat = "identity")+
  scale_fill_brewer(palette = "Greens")+
  facet_wrap(~factor(motif, levels=c("AUACCC", "AGAUCU")))+
  labs(x="",
       y = "% with motif",
       fill= "transcript region")+
  myTheme1

1.2 Positional analysis of the AUACCC motif around the start and stop codon

Here I want to plot the meta profile of the AUACCC motif in Khd4-bound mRNA and all AUACCC-containing mRNAs. The goal is to assess the specific enrichment of this motif around the start and stop codons, with a particular focus on the 5’ and 3’ UTR regions.

Show the code

# stop windows
# stop codon 
stops= flank(genes.GR, -303,start = FALSE) %>% flank(.,-3) 
stops=flank(stops,-1, start= FALSE) 

# start windows
# start codon 
starts=flank(genes.GR,-303, start=TRUE) %>% flank(., -3, start=FALSE)
starts=flank(starts,-1, start= TRUE)

Show the code

#all AUACCC motif high.conf
# granges of all AUACCC motif in the high conf target
sel.highconf=as.vector(auaccc.gr$gene.id  %in%  high.conf$gene.id)# filter auaccc.gr genes that are also found in high confident Khd4 targets
allHC.auaccc.gr=auaccc.gr[sel.highconf] # subset auaccc.gr for high confident Khd4 targets based on gene id

Show the code

# all high confident targets with AGAUCU
sel.agaucu.HC=as.vector(agaucu.gr$gene.id %in% high.conf$gene.id)
HT.agaucu.gr=agaucu.gr[sel.agaucu.HC]

Show the code

# putative targets with AUACCC
sel.auaccc.put=as.vector(auaccc.gr$gene.id %in% auaccc_DGE$gene.ID)
put.auaccc.gr=auaccc.gr[sel.auaccc.put]

Show the code

# putative targets with AGAUCU
sel.agaucu.put=as.vector(agaucu.gr$gene.id %in% auaccc.gr$gene.id)
put.agaucu.gr=agaucu.gr[sel.agaucu.put]

Show the code

get.relative.distance = function(ref.site, motifs, windowsize,name){
  overlap=findOverlaps(ref.site+windowsize, motifs)
  rel.dist=start(motifs)[to(overlap)]-start(ref.site)[from(overlap)]
  sel= as.vector(strand(ref.site)[from(overlap)]=="-")
  rel.dist [sel] = rel.dist [sel]*-1
  pf=data.frame(rel.dist = rel.dist, set=name) 
  return(pf)
}

Show the code

# AUACCC
# stop window
pf1=rbind(
    get.relative.distance(stops, put.auaccc.gr,500, "putative targets"),
    get.relative.distance(stops, allHC.auaccc.gr,500, "Khd4_bound mRNA")
)

dens1_stop_auaccc=ggplot(pf1, aes(rel.dist,col=set))+
      geom_density(adjust =1.5,size=1,alpha=0.5)+
      scale_color_manual(values=c("Khd4_bound mRNA" = "darkred",  "putative targets" = "grey30" ))+
      coord_cartesian(xlim = c(-500,300),ylim = c(0,0.007))+
      labs(x="",
           y="")+
      myTheme1
      

# start window
pf2=rbind(
    get.relative.distance(starts, put.auaccc.gr,500, "putative targets"),
    get.relative.distance(starts, allHC.auaccc.gr,500, "Khd4_bound mRNA")
)

dens1_start_auaccc=ggplot(pf2, aes(rel.dist,col=set))+
      geom_density(adjust =1.5,size=1,alpha=0.5)+
      scale_color_manual(values=c("Khd4_bound mRNA" = "darkred",  "putative targets" = "grey30" ))+
      coord_cartesian(xlim = c(-300,500),ylim = c(0,0.007))+
      labs(x="",
           y="")+
      myTheme1
      

auaccc_dense_pl =ggarrange(dens1_start_auaccc, dens1_stop_auaccc, common.legend = TRUE) 

auaccc_dense_pl = annotate_figure(auaccc_dense_pl, top = text_grob("AUACCC"))             

# AGAUCU
# stop window
pf3=rbind(
    get.relative.distance(stops, put.agaucu.gr,500, "putative targets"),
    get.relative.distance(stops, HT.agaucu.gr,500, "Khd4_bound mRNA")
)

dens1_stop_agaucuc=ggplot(pf3, aes(rel.dist,col=set))+
      geom_density(adjust =1.5,size=1,alpha=0.5)+
      scale_color_manual(values=c("Khd4_bound mRNA" = "darkred",  "putative targets" = "grey30" ))+
      coord_cartesian(xlim = c(-500,300),ylim = c(0,0.007))+
      labs(x="",
           y="")+
      myTheme1
      


# start window

pf4=rbind(
    get.relative.distance(starts, put.agaucu.gr,500, "putative targets"),
    get.relative.distance(starts, HT.agaucu.gr,500, "Khd4_bound mRNA")
)

dens1_start_agaucu=ggplot(pf4, aes(rel.dist,col=set))+
      geom_density(adjust =1.5,size=1,alpha=0.5)+
      scale_color_manual(values=c("Khd4_bound mRNA" = "darkred",  "putative targets" = "grey30" ))+
      coord_cartesian(xlim = c(-300,500),ylim = c(0,0.007))+
      labs(x="",
           y="")+
      myTheme1
      
agaucu_dense_pl = ggarrange(dens1_start_agaucu, dens1_stop_agaucuc, common.legend = TRUE)

agaucu_dense_pl = annotate_figure(agaucu_dense_pl, top=text_grob("AGAUCU"))

pl = ggarrange(auaccc_dense_pl,agaucu_dense_pl, nrow= 2,common.legend = TRUE )

annotate_figure(pl, bottom = text_grob("position relative to start or stop codon (nt)"), left = "motif density")

Show the code

sessionInfo()

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_Germany.utf8  LC_CTYPE=English_Germany.utf8   
[3] LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C                    
[5] LC_TIME=English_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] lookup_1.0                        ggpointdensity_0.1.0             
 [3] colorspace_2.1-0                  ggvenn_0.1.10                    
 [5] paletteer_1.6.0                   tidyHeatmap_1.8.1                
 [7] ggseqlogo_0.2                     clusterProfiler_4.10.1           
 [9] beeswarm_0.4.0                    rstatix_0.7.2                    
[11] cowplot_1.1.3                     kableExtra_1.4.0                 
[13] UpSetR_1.4.0                      EnrichedHeatmap_1.32.0           
[15] circlize_0.4.16                   ggplotify_0.1.2                  
[17] ggh4x_0.2.8                       ggsci_3.2.0                      
[19] gprofiler2_0.2.3                  reshape2_1.4.4                   
[21] viridis_0.6.5                     viridisLite_0.4.2                
[23] ggrastr_1.0.2                     ggforce_0.4.2                    
[25] ggpubr_0.6.0                      ggrepel_0.9.5                    
[27] gginnards_0.2.0                   ggpp_0.5.8-1                     
[29] gridExtra_2.3                     magrittr_2.0.3                   
[31] magick_2.8.3                      lubridate_1.9.3                  
[33] stringr_1.5.1                     purrr_1.0.2                      
[35] readr_2.1.5                       tibble_3.2.1                     
[37] tidyverse_2.0.0                   dplyr_1.1.4                      
[39] tidyr_1.3.1                       forcats_1.0.0                    
[41] GenomicFeatures_1.54.4            AnnotationDbi_1.64.1             
[43] Biobase_2.62.0                    ComplexHeatmap_2.18.0            
[45] BindingSiteFinder_2.0.0           ggplot2_3.5.1                    
[47] BSgenome.Umaydis.ENSMBL.UM1_2.0.0 BSgenome_1.70.2                  
[49] rtracklayer_1.62.0                BiocIO_1.12.0                    
[51] Biostrings_2.70.1                 XVector_0.42.0                   
[53] GenomicRanges_1.54.1              GenomeInfoDb_1.38.8              
[55] IRanges_2.36.0                    S4Vectors_0.40.2                 
[57] BiocGenerics_0.48.1              

loaded via a namespace (and not attached):
  [1] fs_1.6.3                    matrixStats_1.2.0          
  [3] bitops_1.0-7                enrichplot_1.22.0          
  [5] HDO.db_0.99.1               httr_1.4.7                 
  [7] RColorBrewer_1.1-3          doParallel_1.0.17          
  [9] tools_4.3.1                 backports_1.5.0            
 [11] utf8_1.2.4                  R6_2.5.1                   
 [13] lazyeval_0.2.2              ggdist_3.3.2               
 [15] GetoptLong_1.0.5            withr_3.0.1                
 [17] prettyunits_1.2.0           cli_3.6.1                  
 [19] scatterpie_0.2.3            labeling_0.4.3             
 [21] Rsamtools_2.18.0            systemfonts_1.1.0          
 [23] yulab.utils_0.1.5           gson_0.1.0                 
 [25] DOSE_3.28.2                 svglite_2.1.3              
 [27] rstudioapi_0.16.0           RSQLite_2.3.7              
 [29] generics_0.1.3              gridGraphics_0.5-1         
 [31] shape_1.4.6.1               crosstalk_1.2.1            
 [33] car_3.1-2                   distributional_0.4.0       
 [35] dendextend_1.17.1           GO.db_3.18.0               
 [37] Matrix_1.6-4                ggbeeswarm_0.7.2           
 [39] fansi_1.0.6                 abind_1.4-5                
 [41] lifecycle_1.0.4             yaml_2.3.7                 
 [43] carData_3.0-5               SummarizedExperiment_1.32.0
 [45] qvalue_2.34.0               SparseArray_1.2.3          
 [47] BiocFileCache_2.10.2        blob_1.2.4                 
 [49] promises_1.3.0              crayon_1.5.3               
 [51] lattice_0.22-5              KEGGREST_1.42.0            
 [53] pillar_1.9.0                knitr_1.48                 
 [55] fgsea_1.28.0                rjson_0.2.21               
 [57] codetools_0.2-20            fastmatch_1.1-4            
 [59] glue_1.6.2                  ggfun_0.1.5                
 [61] data.table_1.15.4           vctrs_0.6.5                
 [63] png_0.1-8                   treeio_1.26.0              
 [65] gtable_0.3.5                rematch2_2.1.2             
 [67] cachem_1.0.8                xfun_0.45                  
 [69] mime_0.12                   S4Arrays_1.2.0             
 [71] tidygraph_1.3.1             iterators_1.0.14           
 [73] nlme_3.1-165                ggtree_3.10.1              
 [75] bit64_4.0.5                 progress_1.2.3             
 [77] filelock_1.0.3              vipor_0.4.7                
 [79] DBI_1.2.3                   DESeq2_1.42.1              
 [81] tidyselect_1.2.1            bit_4.0.5                  
 [83] compiler_4.3.1              curl_5.2.1                 
 [85] xml2_1.3.6                  DelayedArray_0.28.0        
 [87] plotly_4.10.4               shadowtext_0.1.4           
 [89] scales_1.3.0                rappdirs_0.3.3             
 [91] digest_0.6.33               rmarkdown_2.27             
 [93] htmltools_0.5.8.1           pkgconfig_2.0.3            
 [95] MatrixGenerics_1.14.0       highr_0.11                 
 [97] dbplyr_2.5.0                fastmap_1.1.1              
 [99] rlang_1.1.1                 GlobalOptions_0.1.2        
[101] htmlwidgets_1.6.4           shiny_1.9.1                
[103] farver_2.1.2                jsonlite_1.8.8             
[105] BiocParallel_1.36.0         GOSemSim_2.28.1            
[107] RCurl_1.98-1.13             polynom_1.4-1              
[109] GenomeInfoDbData_1.2.11     patchwork_1.2.0            
[111] munsell_0.5.1               Rcpp_1.0.12                
[113] ape_5.8                     stringi_1.8.3              
[115] ggraph_2.2.1                zlibbioc_1.48.0            
[117] MASS_7.3-60                 plyr_1.8.9                 
[119] parallel_4.3.1              graphlayouts_1.1.1         
[121] splines_4.3.1               hms_1.1.3                  
[123] locfit_1.5-9.10             igraph_2.0.3               
[125] ggsignif_0.6.4              biomaRt_2.58.2             
[127] XML_3.99-0.16.1             evaluate_0.24.0            
[129] httpuv_1.6.15               tzdb_0.4.0                 
[131] foreach_1.5.2               tweenr_2.0.3               
[133] polyclip_1.10-6             clue_0.3-65                
[135] xtable_1.8-4                broom_1.0.6                
[137] restfulr_0.0.15             tidytree_0.4.6             
[139] later_1.3.2                 aplot_0.2.3                
[141] memoise_2.0.1               GenomicAlignments_1.38.0   
[143] cluster_2.1.6               timechange_0.3.0