Skip to content

FLEXIs.md

Latest commit

 

History

History
970 lines (798 loc) · 42 KB

FLEXIs.md

File metadata and controls

970 lines (798 loc) · 42 KB
Raw

RNAs in Extracellular Vesicles

Wilfried Guiblet 2025-04-11

I. FLEXIs

I.1 Global GCcontent and Minimum Free Energy

# Create negative control - random introns shorter than 200 bp
awk '{if ($3-$2 <=200) print $0}' references/mm10.GENCODE_VM23.UCSC.introns.bed > references/mm10.GENCODE_VM23.UCSC.introns200bp.bed
shuf -n 200 references/mm10.GENCODE_VM23.UCSC.introns200bp.bed > inputs/mm10.GENCODE_VM23.UCSC.introns200bp.subset.bed


# Get sequence of introns
python scripts/FormatFLEXIs.py inputs/FLEXI_up_EVs.tsv inputs/FLEXI_up_EVs.bed
python scripts/FormatFLEXIs.py inputs/FLEXI_up_cells.tsv inputs/FLEXI_up_cells.bed
bedtools getfasta -s -fi references/mm10.fa -bed inputs/FLEXI_up_EVs.bed > inputs/FLEXI_up_EVs.fa
bedtools getfasta -s -fi references/mm10.fa -bed inputs/FLEXI_up_cells.bed > inputs/FLEXI_up_cells.fa
bedtools getfasta -s -fi references/mm10.fa -bed inputs/mm10.GENCODE_VM23.UCSC.introns200bp.subset.bed > inputs/introns200bp.fa


# Scrambled fasta sequence. This step could be removed: it is now done automatically for ZScore computations
python scripts/Fasta_Scrambler.py inputs/FLEXI_up_EVs.fa inputs/FLEXI_up_EVs.scrambled.fa
python scripts/Fasta_Scrambler.py inputs/FLEXI_up_cells.fa inputs/FLEXI_up_cells.scrambled.fa
python scripts/Fasta_Scrambler.py inputs/introns200bp.fa inputs/introns200bp.scrambled.fa


# Run Global GCcontent and MFE
python scripts/RNAfoldMeasures.py --FastaFile inputs/FLEXI_up_EVs.fa --FastaFileScrambled inputs/FLEXI_up_EVs.scrambled.fa --OutFile results/FLEXI_up_EVs.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/FLEXI_up_cells.fa --FastaFileScrambled inputs/FLEXI_up_cells.scrambled.fa --OutFile results/FLEXI_up_cells.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/introns200bp.fa --FastaFileScrambled inputs/introns200bp.scrambled.fa --OutFile results/introns200bp.folding.csv
FLEXI_up_EVs <- read.csv('output_files/FLEXI_up_EVs.folding.csv', row.names=1)
FLEXI_up_cells <- read.csv('output_files/FLEXI_up_cells.folding.csv', row.names=1)
introns200b <- read.csv('output_files/introns200bp.folding.csv', row.names=1)

head(FLEXI_up_EVs)
##                               Gcontent GCcontent   MFE MFE_G4 MFE_scrambled
## chr1:75185846-75185950(-)    0.3557692 0.6057692 -40.8  -40.8         -38.8
## chr1:88277286-88277373(-)    0.4252874 0.7701149 -40.4  -40.4         -44.0
## chr10:128276624-128276747(+) 0.3414634 0.5934959 -46.9  -46.9         -50.4
## chr10:128911870-128911967(+) 0.3402062 0.5876289 -40.4  -46.9         -28.8
## chr10:39533711-39533808(+)   0.3092784 0.5773196 -33.3  -33.3         -28.8
## chr10:78405480-78405565(-)   0.3176471 0.5764706 -26.1  -26.1         -35.1
##                              MFE_G4_scrambled     Zscore  Zscore_G4
## chr1:75185846-75185950(-)               -38.8 -1.2795071 -1.3361742
## chr1:88277286-88277373(-)               -44.0  0.9140650  0.7433490
## chr10:128276624-128276747(+)            -50.4 -1.0418883 -0.7408230
## chr10:128911870-128911967(+)            -28.8 -2.4940312 -4.6810875
## chr10:39533711-39533808(+)              -28.8 -1.3226812 -0.9315546
## chr10:78405480-78405565(-)              -35.1  0.1481971  0.1972206

ggplot() +
  geom_boxplot(data = FLEXI_up_EVs, aes(x='EV', y=Gcontent)) +
  geom_boxplot(data = FLEXI_up_cells, aes(x='Cells', y=Gcontent)) +
  geom_boxplot(data = introns200b, aes(x='Introns200bp', y=Gcontent)) +
  labs(title='G content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'Cells', 'Introns200bp'))

ggplot() +
  geom_boxplot(data = FLEXI_up_EVs, aes(x='EV', y=GCcontent)) +
  geom_boxplot(data = FLEXI_up_cells, aes(x='Cells', y=GCcontent)) +
  geom_boxplot(data = introns200b, aes(x='Introns200bp', y=GCcontent)) +
  labs(title='GC content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'Cells', 'Introns200bp'))

custom_labels <- c('EV', 'EV_G4', 'Cells', 'Cells_G4', 'Introns200bp', 'Introns200bp_G4')
ggplot() +
  geom_boxplot(data = FLEXI_up_EVs, aes(x=1, y=MFE)) +
  geom_boxplot(data = FLEXI_up_EVs, aes(x=2, y=MFE_G4)) +
  geom_boxplot(data = FLEXI_up_cells, aes(x=3, y=MFE)) +
  geom_boxplot(data = FLEXI_up_cells, aes(x=4, y=MFE_G4)) +
  geom_boxplot(data = introns200b, aes(x=5, y=MFE_G4)) +
  geom_boxplot(data = introns200b, aes(x=6, y=MFE_G4)) +
  labs(title='MFE', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'introns200bp')
ggplot() +
  geom_boxplot(data = FLEXI_up_EVs, aes(x=1, y=Zscore)) +
  #geom_boxplot(data = FLEXI_up_cells, aes(x=2, y=Zscore)) +
  geom_boxplot(data = introns200b, aes(x=2, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'introns200bp')
ggplot() +
  geom_boxplot(data = FLEXI_up_EVs, aes(x=1, y=Zscore_G4)) +
  #geom_boxplot(data = FLEXI_up_cells, aes(x=2, y=Zscore)) +
  geom_boxplot(data = introns200b, aes(x=2, y=Zscore_G4)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore G4', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

FLEXI_up_EVs <- read.csv('output_files/FLEXI_up_EVs.folding.csv', row.names=1)
FLEXI_up_cells <- read.csv('output_files/FLEXI_up_cells.folding.csv', row.names=1)
introns200b <- read.csv('output_files/introns200bp.folding.csv', row.names=1)

# Add sample identifiers and concatenate
concatenated_df <- bind_rows(
    FLEXI_up_EVs %>% mutate(sample = "EVs"),
    #FLEXI_up_cells %>% mutate(sample = "Cells"),
    introns200b %>% mutate(sample = "Introns200bp")
)

pwc <- concatenated_df %>%
  wilcox_test(GCcontent ~ sample, p.adjust.method = "holm")

pwc <- pwc %>% add_xy_position(x = "sample")
ggboxplot(concatenated_df, x = "sample", y = "GCcontent") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_pwc_label(pwc),
    caption = '',
    x = ''
  ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme_bw()

II. 5’ UTRs peaks

II.1 Sliding GCcontent and Minimum Free Energy

WC1 <- read.csv('output_files/dedNA_WC1.5UTR.MFE.sliding.csv')
WC2 <- read.csv('output_files/dedNA_WC2.5UTR.MFE.sliding.csv')
WC3 <- read.csv('output_files/dedNA_WC3.5UTR.MFE.sliding.csv')
WC4 <- read.csv('output_files/dedNA_WC4.5UTR.MFE.sliding.csv')
WC5 <- read.csv('output_files/dedNA_WC5.5UTR.MFE.sliding.csv')
EV1 <- read.csv('output_files/dedNA_EV1.5UTR.MFE.sliding.csv')
EV2 <- read.csv('output_files/dedNA_EV2.5UTR.MFE.sliding.csv')
EV3 <- read.csv('output_files/dedNA_EV3.5UTR.MFE.sliding.csv')
EV4 <- read.csv('output_files/dedNA_EV4.5UTR.MFE.sliding.csv')
EV5 <- read.csv('output_files/dedNA_EV5.5UTR.MFE.sliding.csv')

group.colors <- c('WC' = 'darkgrey',
                  'EV' = 'lightblue')

head(WC1)
##    X Total_Gcontent Total_GCcontent Total_MFE Total_MFE_G4 Total_ZScore
## 1 11      0.2887538       0.6337386 -2.788298    -2.979787  -0.92819211
## 2 12      0.2862209       0.6301925 -2.800000    -2.991489  -0.05066709
## 3 13      0.2877406       0.6301925 -2.904255    -3.095745  -0.16686855
## 4 14      0.2902736       0.6337386 -3.064894    -3.256383  -0.23113506
## 5 15      0.2857143       0.6312057 -3.086170    -3.297872  -0.38241061
## 6 16      0.2882472       0.6261398 -3.151064    -3.368085  -0.28442673
##   Total_ZScore_G4
## 1      -0.1052432
## 2      -0.1112357
## 3      -0.1739058
## 4      -0.6317949
## 5      -0.3401235
## 6      -0.3531138

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_ZScore, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_ZScore, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_ZScore, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_ZScore, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_ZScore, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_ZScore, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_ZScore, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_ZScore, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_ZScore, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_ZScore, color = 'EV')) +
  labs(title='Zscore in sliding window', x="Position in 5'UTRs") +
  ylim(-3, 3) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='') 

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_ZScore_G4, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_ZScore_G4, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_ZScore_G4, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_ZScore_G4, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_ZScore_G4, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_ZScore_G4, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_ZScore_G4, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_ZScore_G4, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_ZScore_G4, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_ZScore_G4, color = 'EV')) +
  labs(title='Zscore G4  in sliding window', x="Position in 5'UTRs") +
  ylim(-3, 3) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='') 

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_MFE, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_MFE, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_MFE, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_MFE, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_MFE, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_MFE, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_MFE, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_MFE, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_MFE, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_MFE, color = 'EV')) +
  labs(title='MFE in sliding window', x="Position in 5'UTRs") +
  #ylim(0.4, 0.7) +
  theme_bw() +
  scale_color_manual(values=group.colors, name='')

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_MFE_G4, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_MFE_G4, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_MFE_G4, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_MFE_G4, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_MFE_G4, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_MFE_G4, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_MFE_G4, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_MFE_G4, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_MFE_G4, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_MFE_G4, color = 'EV')) +
  labs(title='MFE G4 in sliding window', x="Position in 5'UTRs") +
  #ylim(0.4, 0.7) +
  theme_bw() +
  scale_color_manual(values=group.colors, name='') 

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_GCcontent, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_GCcontent, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_GCcontent, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_GCcontent, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_GCcontent, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_GCcontent, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_GCcontent, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_GCcontent, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_GCcontent, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_GCcontent, color = 'EV')) +
  labs(title='GCcontent in sliding window', x="Position in 5'UTRs") +
  ylim(0.4, 0.7) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',) +
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='')

ggplot() +
  geom_line(data = WC1, aes(x = X, y=Total_Gcontent, color = 'WC')) +
  geom_line(data = WC2, aes(x = X, y=Total_Gcontent, color = 'WC')) +
  geom_line(data = WC3, aes(x = X, y=Total_Gcontent, color = 'WC')) +
  geom_line(data = WC4, aes(x = X, y=Total_Gcontent, color = 'WC')) +
  geom_line(data = WC5, aes(x = X, y=Total_Gcontent, color = 'WC')) +
  geom_line(data = EV1, aes(x = X, y=Total_Gcontent, color = 'EV')) +
  geom_line(data = EV2, aes(x = X, y=Total_Gcontent, color = 'EV')) +
  geom_line(data = EV3, aes(x = X, y=Total_Gcontent, color = 'EV')) +
  geom_line(data = EV4, aes(x = X, y=Total_Gcontent, color = 'EV')) +
  geom_line(data = EV5, aes(x = X, y=Total_Gcontent, color = 'EV')) +
  labs(title='Gcontent in sliding window', x="Position in 5'UTRs") +
  ylim(0.1, 0.5) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',) +
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='')

II.2 Peaks Length

WC1 <- read.csv('output_files/dedNA_WC1.5UTR.lengths', sep = '\t', header = FALSE)
WC2 <- read.csv('output_files/dedNA_WC2.5UTR.lengths', sep = '\t', header = FALSE)
WC3 <- read.csv('output_files/dedNA_WC3.5UTR.lengths', sep = '\t', header = FALSE)
WC4 <- read.csv('output_files/dedNA_WC4.5UTR.lengths', sep = '\t', header = FALSE)
WC5 <- read.csv('output_files/dedNA_WC5.5UTR.lengths', sep = '\t', header = FALSE)
EV1 <- read.csv('output_files/dedNA_EV1.5UTR.lengths', sep = '\t', header = FALSE)
EV2 <- read.csv('output_files/dedNA_EV2.5UTR.lengths', sep = '\t', header = FALSE)
EV3 <- read.csv('output_files/dedNA_EV3.5UTR.lengths', sep = '\t', header = FALSE)
EV4 <- read.csv('output_files/dedNA_EV4.5UTR.lengths', sep = '\t', header = FALSE)
EV5 <- read.csv('output_files/dedNA_EV5.5UTR.lengths', sep = '\t', header = FALSE)

group.colors <- c('WC' = 'darkgrey',
                  'EV' = 'lightblue')

ggplot() +
  geom_density(data = WC1, aes(x = V1, color = 'WC'), size = 1) +
  geom_density(data = WC2, aes(x = V1, color = 'WC'), size = 1) +
  geom_density(data = WC3, aes(x = V1, color = 'WC'), size = 1) +
  geom_density(data = WC4, aes(x = V1, color = 'WC'), size = 1) +
  geom_density(data = WC5, aes(x = V1, color = 'WC'), size = 1) +
  geom_density(data = EV1, aes(x = V1, color = 'EV'), size = 1) +
  geom_density(data = EV2, aes(x = V1, color = 'EV'), size = 1) +
  geom_density(data = EV3, aes(x = V1, color = 'EV'), size = 1) +
  geom_density(data = EV4, aes(x = V1, color = 'EV'), size = 1) +
  geom_density(data = EV5, aes(x = V1, color = 'EV'), size = 1) +
  xlim(0, 500) + 
  theme_classic(base_size = 20) +
  labs(y = 'Density',
       x = 'Length (bp)',
       title = "Peaks in 5'UTRs") +
  scale_color_manual(values=group.colors, name='')

III. 5’UTR peaks Up- and Downregulated in EVs

III.1 Differential Expression Analysis

DiffExp <- read.table('output_files/merged_peaks.5UTR.DiffExp.tsv')

# Calculate -log10(p-value)
DiffExp$logP <- -log10(DiffExp$padj)

# Define thresholds for significance (adjust as needed)
logFC_cutoff <- 1
pval_cutoff <- 0.05

# Label significant points
DiffExp$Significance <- "Not Significant"
DiffExp$Significance[DiffExp$log2FoldChange > logFC_cutoff & DiffExp$padj < pval_cutoff] <- "Upregulated"
DiffExp$Significance[DiffExp$log2FoldChange < -logFC_cutoff & DiffExp$padj < pval_cutoff] <- "Downregulated"


# Plot volcano
ggplot(DiffExp, aes(x = log2FoldChange, y = logP)) +
    geom_point(aes(color = Significance), alpha = 0.6) +
    scale_color_manual(values = c("red", "gray", "blue")) +   # Color for each significance category
    labs(title = "5' UTRs Peaks",
         x = "Log2 Fold Change",
         y = "-Log10 P-value") +
    geom_vline(xintercept = c(-logFC_cutoff, logFC_cutoff), linetype = "dashed") +  # Threshold lines for logFC
    geom_hline(yintercept = -log10(pval_cutoff), linetype = "dashed") +             # Threshold line for p-value
    theme_minimal()

III.2 Sliding GCcontent and Minimum Free Energy

upregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.upregulated.200bp.MFE.sliding.csv')
downregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.downregulated.200bp.MFE.sliding.csv')

group.colors <- c('Downregulated' = 'red',
                  'Upregulated' = 'blue')

head(upregulated)
##    X Total_Gcontent Total_GCcontent Total_MFE Total_MFE_G4 Total_ZScore
## 1 26      0.3622573       0.6908107 -16.27910    -17.23373  -0.17352946
## 2 27      0.3626963       0.6920300 -16.27910    -17.25154  -0.17530933
## 3 28      0.3634767       0.6925178 -16.25622    -17.24289  -0.13898513
## 4 29      0.3627939       0.6919325 -16.24254    -17.21147  -0.12510613
## 5 30      0.3633792       0.6920788 -16.20149    -17.18731  -0.08957844
## 6 31      0.3640133       0.6922251 -16.27985    -17.26970  -0.13224817
##   Total_ZScore_G4
## 1      -0.1799672
## 2      -0.1701080
## 3      -0.1434926
## 4      -0.1527102
## 5      -0.1325192
## 6      -0.1522364

ggplot() +
  geom_line(data = upregulated, aes(x = X, y=Total_ZScore, color = 'Upregulated')) +
  geom_line(data = downregulated, aes(x = X, y=Total_ZScore, color = 'Downregulated')) +
  labs(title='ZScore in sliding window', x="Position in peak") +
  ylim(-3, 3) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='')

ggplot() +
  geom_line(data = upregulated, aes(x = X, y=Total_ZScore_G4, color = 'Upregulated')) +
  geom_line(data = downregulated, aes(x = X, y=Total_ZScore_G4, color = 'Downregulated')) +
  labs(title='ZScore G4 in sliding window', x="Position in peak") +
  ylim(-3, 3) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='') 

ggplot() +
  geom_line(data = upregulated, aes(x = X, y=Total_GCcontent, color = 'Upregulated')) +
  geom_line(data = downregulated, aes(x = X, y=Total_GCcontent, color = 'Downregulated')) +
  labs(title='GCcontent in sliding window', x="Position in peak") +
  ylim(0.4, 0.8) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',) +
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='') 

ggplot() +
  geom_line(data = upregulated, aes(x = X, y=Total_Gcontent, color = 'Upregulated')) +
  geom_line(data = downregulated, aes(x = X, y=Total_Gcontent, color = 'Downregulated')) +
  labs(title='Gcontent in sliding window', x="Position in peak") +
  ylim(0.1, 0.6) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',) +
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  theme_bw() +
  scale_color_manual(values=group.colors, name='')

III.3 Peaks Length

upregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.upregulated.bed', sep = '\t', header = FALSE)
downregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.downregulated.bed', sep = '\t', header = FALSE)


group.colors <- c('Downregulated' = 'red',
                  'Upregulated' = 'blue')

ggplot() +
  geom_density(data = upregulated, aes(x = V3-V2, color = 'Upregulated'), size = 1) +
  geom_density(data = downregulated, aes(x = V3-V2, color = 'Downregulated'), size = 1) +

  xlim(0, 500) + 
  theme_classic(base_size = 20) +
  labs(y = 'Density',
       x = 'Length (bp)',
       title = "Peaks in 5'UTRs") +
  scale_color_manual(values=group.colors, name='')

III.4 Global GCcontent and Minimum Free Energy

Upregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.upregulated.folding.csv', row.names=1)
Downregulated <- read.csv('output_files/merged_peaks.5UTR.DiffExp.downregulated.folding.csv', row.names=1)
Random <- read.csv('output_files/merged_peaks.5UTR.DiffExp.shuff.folding.csv', row.names=1)

custom_labels <- c('Upregulated', 'Downregulated', 'Random')
head(Upregulated)
##                           Gcontent GCcontent    MFE  MFE_G4      Zscore
## chr7:100493313-100493476 0.2760736 0.6503067  -52.7  -52.70  0.69032535
## chr11:94211421-94211538  0.4700855 0.7606838  -60.7  -60.70 -2.02682725
## chr3:96238462-96239032   0.3631579 0.6754386 -274.7 -274.70  0.04861162
## chr13:45389726-45389817  0.5054945 0.7362637  -40.6  -52.10 -0.74755219
## chr1:167284965-167285077 0.4642857 0.7767857  -66.5  -66.50 -1.55664327
## chr3:95658738-95658899   0.2857143 0.6708075  -57.1  -57.78  0.56601310
##                           Zscore_G4
## chr7:100493313-100493476  0.6531920
## chr11:94211421-94211538  -0.3233781
## chr3:96238462-96239032    0.4226899
## chr13:45389726-45389817  -1.2877757
## chr1:167284965-167285077 -0.7155077
## chr3:95658738-95658899    0.1682145

ggplot() +
  geom_boxplot(data = Upregulated, aes(x='Upregulated', y=GCcontent)) +
  geom_boxplot(data = Downregulated, aes(x='Downregulated', y=GCcontent)) +
  geom_boxplot(data = Random, aes(x='Random', y=GCcontent)) +
  labs(title='GC content', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

ggplot() +
  geom_boxplot(data = Upregulated, aes(x='Upregulated', y=Gcontent)) +
  geom_boxplot(data = Downregulated, aes(x='Downregulated', y=Gcontent)) +
  geom_boxplot(data = Random, aes(x='Random', y=Gcontent)) +
  labs(title='G content', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

ggplot() +
  geom_boxplot(data = Upregulated, aes(x='Upregulated', y=MFE)) +
  geom_boxplot(data = Downregulated, aes(x='Downregulated', y=MFE)) +
  geom_boxplot(data = Random, aes(x='Random', y=MFE)) +
  labs(title='MFE', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels) +
  ylim(-100, 0) 

ggplot() +
  geom_boxplot(data = Upregulated, aes(x='Upregulated', y=MFE_G4)) +
  geom_boxplot(data = Downregulated, aes(x='Downregulated', y=MFE_G4)) +
  geom_boxplot(data = Random, aes(x='Random', y=MFE_G4)) +
  labs(title='MFE G4', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels) +
  ylim(-100, 0)

ggplot() +
  geom_boxplot(data = Upregulated, aes(x=1, y=Zscore)) +
  geom_boxplot(data = Downregulated, aes(x=2, y=Zscore)) +
  geom_boxplot(data = Random, aes(x=3, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

ggplot() +
  geom_boxplot(data = Upregulated, aes(x=1, y=Zscore_G4)) +
  geom_boxplot(data = Downregulated, aes(x=2, y=Zscore_G4)) +
  geom_boxplot(data = Random, aes(x=3, y=Zscore_G4)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore G4', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

IV. Tailer analysis

module load python
#pip install jla-tailer
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_WC1.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_WC2.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_WC3.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_WC4.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_WC5.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_EV1.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_EV2.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_EV3.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_EV4.bam
Tailer -a gencode.vM23.annotation.gtf rawfiles/dedNA_EV5.bam

mv rawfiles/dedNA_*_tail.csv Tailer/gencode_genes/

### IMPORTANT : MODIFY LINE 171 TO RUN ON NON-GENES:

vim /home/guibletwm/.local/lib/python3.10/site-packages/Tailer/TailerFunctions.py

### 

# Extract entries of interest


awk 'NR==FNR {regions[$1,$2,$3,$6]; next} ($1,$4-1,$5,$7) in regions' references/45S.bed references/rmsk.mm10.gtf| awk '{gsub(/transcript_id/, "gene_name"); print}'  > references/45S.gtf
sed 's/4\.5SRNA/45SRNA/g' references/45S.gtf > temp.txt && mv temp.txt  references/45S.gtf


awk 'NR==FNR {regions[$1,$2,$3,$6]; next} ($1,$4-1,$5,$7) in regions' inputs/BEDfinal/tRNA.bed references/rmsk.mm10.gtf | awk '{gsub(/transcript_id/, "gene_name"); print}' > references/tRNA.gtf
sed -i 's/gene_name "\([^"]*\)"/gene_name "\1"/g; s/-/_/g; s/_/_/g' references/tRNA.gtf


head references/rmsk.mm10.gtf > references/rmsk.mm10.test.gtf
Tailer -a references/rmsk.mm10.test.gtf rawfiles/dedNA_WC1.bam


Tailer -a references/45S.gtf rawfiles/dedNA_WC1.bam
Tailer -a references/45S.gtf rawfiles/dedNA_WC2.bam
Tailer -a references/45S.gtf rawfiles/dedNA_WC3.bam
Tailer -a references/45S.gtf rawfiles/dedNA_WC4.bam
Tailer -a references/45S.gtf rawfiles/dedNA_WC5.bam
Tailer -a references/45S.gtf rawfiles/dedNA_EV1.bam
Tailer -a references/45S.gtf rawfiles/dedNA_EV2.bam
Tailer -a references/45S.gtf rawfiles/dedNA_EV3.bam
Tailer -a references/45S.gtf rawfiles/dedNA_EV4.bam
Tailer -a references/45S.gtf rawfiles/dedNA_EV5.bam

mv rawfiles/dedNA_*_tail.csv Tailer/rmsk.45S/

Tailer -a references/tRNA.gtf rawfiles/dedNA_WC1.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_WC2.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_WC3.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_WC4.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_WC5.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_EV1.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_EV2.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_EV3.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_EV4.bam
Tailer -a references/tRNA.gtf rawfiles/dedNA_EV5.bam


mv rawfiles/dedNA_*_tail.csv Tailer/rmsk.tRNA/

V.Mir

Part I

grep -w Mir1981 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 > Mir.EV.bed
grep -w Mir1982 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6418 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6931 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6958 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6963 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6968 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6977 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6979 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6981 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir6998 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7012 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir702 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7054 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7060 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7067 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7068 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir7070 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed
grep -w Mir877 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.EV.bed


grep -w Mir1231 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 > Mir.WholeCells.bed
grep -w Mir1933 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir3064 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir3103 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir5107 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir6537 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir6983 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir6990 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir6994 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir704 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
grep -w Mir7091 references/gencode.vM23.annotation.bed | awk '{if ($8=="transcript") print $0}' | cut -f 1-6 >> Mir.WholeCells.bed
# Get sequence
bedtools getfasta -s -fi references/mm10.fa -bed Mir.EV.bed > inputs/Mir.EVs.fa
bedtools getfasta -s -fi references/mm10.fa -bed Mir.WholeCells.bed > inputs/Mir.WholeCells.fa


# Scrambled fasta sequence. This step could be removed: it is now done automatically for ZScore computations
python scripts/Fasta_Scrambler.py inputs/Mir.EVs.fa inputs/Mir.EVs.scrambled.fa
python scripts/Fasta_Scrambler.py inputs/Mir.WholeCells.fa inputs/Mir.WholeCells.scrambled.fa


# Run Global GCcontent and MFE
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.EVs.fa --FastaFileScrambled inputs/Mir.EVs.scrambled.fa --OutFile results/Mir.EVs.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.WholeCells.fa --FastaFileScrambled inputs/Mir.WholeCells.scrambled.fa --OutFile results/Mir.WholeCells.folding.csv
Mir_EVs <- read.csv('output_files/Mir.EVs.folding.csv', row.names=1)
Mir_WholeCells <- read.csv('output_files/Mir.WholeCells.folding.csv', row.names=1)

head(Mir_EVs)
##                               Gcontent GCcontent   MFE MFE_G4    Zscore
## chr1:184822406-184822488(-)  0.2926829 0.5365854 -30.3 -30.30 -2.344104
## chr10:80828796-80828870(+)   0.3108108 0.5945946 -26.8 -27.38 -1.348830
## chr5:137529479-137529594(-)  0.2782609 0.5565217 -37.0 -37.00 -1.283028
## chr11:102999847-102999921(-) 0.3108108 0.5675676 -29.4 -29.40 -2.972589
## chr15:89185464-89185537(-)   0.3013699 0.6301370 -28.3 -28.30 -2.048988
## chr15:103350454-103350525(-) 0.2957746 0.5492958 -27.1 -27.10 -2.111003
##                              Zscore_G4
## chr1:184822406-184822488(-)  -2.507128
## chr10:80828796-80828870(+)   -1.299429
## chr5:137529479-137529594(-)  -1.356480
## chr11:102999847-102999921(-) -2.542852
## chr15:89185464-89185537(-)   -1.875746
## chr15:103350454-103350525(-) -2.222638

ggplot() +
  geom_boxplot(data = Mir_EVs, aes(x='EV', y=Gcontent)) +
  geom_boxplot(data = Mir_WholeCells, aes(x='Cells', y=Gcontent)) +
  labs(title='G content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'Cells'))

ggplot() +
  geom_boxplot(data = Mir_EVs, aes(x='EV', y=GCcontent)) +
  geom_boxplot(data = Mir_WholeCells, aes(x='Cells', y=GCcontent)) +
  labs(title='GC content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'Cells'))

custom_labels <- c('EV', 'EV_G4', 'Cells', 'Cells_G4')
ggplot() +
  geom_boxplot(data = Mir_EVs, aes(x=1, y=MFE)) +
  geom_boxplot(data = Mir_EVs, aes(x=2, y=MFE_G4)) +
  geom_boxplot(data = Mir_WholeCells, aes(x=3, y=MFE)) +
  geom_boxplot(data = Mir_WholeCells, aes(x=4, y=MFE_G4)) +
  labs(title='MFE', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'Cells')
ggplot() +
  geom_boxplot(data = Mir_EVs, aes(x=1, y=Zscore)) +
  geom_boxplot(data = Mir_WholeCells, aes(x=2, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'Cells')
ggplot() +
  geom_boxplot(data = Mir_EVs, aes(x=1, y=Zscore_G4)) +
  geom_boxplot(data = Mir_WholeCells, aes(x=2, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore G4', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

Part II

# REMOVED Mirt1 AND Mir100hg > too long


while read entry; do 
    grep -w "$entry" references/gencode.vM23.annotation.bed | awk '$8=="transcript"' | cut -f 1-6 | awk -v OFS='\t' -v e="$entry" '{print $0, e}'; 
done < inputs/Mir.Up_in_EVs.txt > inputs/Mir.Up_in_EVs.bed


while read entry; do 
    grep -w "$entry" references/gencode.vM23.annotation.bed | awk '$8=="transcript"' | cut -f 1-6 | awk -v OFS='\t' -v e="$entry" '{print $0, e}'; 
done < inputs/Mir.Up_in_EVs.No_mirtrons.txt > inputs/Mir.Up_in_EVs.No_mirtrons.bed


while read entry; do 
    grep -w "$entry" references/gencode.vM23.annotation.bed | awk '$8=="transcript"' | cut -f 1-6 | awk -v OFS='\t' -v e="$entry" '{print $0, e}'; 
done < inputs/Mir.Up_in_cells.Top100.txt > inputs/Mir.Up_in_cells.Top100.bed


# Also remove Mir17hg, Mir124-2hg, Mir142hg because length

while read entry; do 
    grep -w "$entry" references/gencode.vM23.annotation.bed | awk '$8=="transcript"' | cut -f 1-6 | awk -v OFS='\t' -v e="$entry" '{print $0, e}'; 
done < inputs/Mir.Up_in_cells.All.txt > inputs/Mir.Up_in_cells.All.bed
# Get sequence
bedtools getfasta -s -fi references/mm10.fa -bed inputs/Mir.Up_in_EVs.bed  | awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' > inputs/Mir.Up_in_EVs.fa
bedtools getfasta -s -fi references/mm10.fa -bed inputs/Mir.Up_in_EVs.No_mirtrons.bed  | awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' > inputs/Mir.Up_in_EVs.No_mirtrons.fa
bedtools getfasta -s -fi references/mm10.fa -bed inputs/Mir.Up_in_cells.Top100.bed  | awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' > inputs/Mir.Up_in_cells.Top100.fa
bedtools getfasta -s -fi references/mm10.fa -bed inputs/Mir.Up_in_cells.All.bed  | awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' > inputs/Mir.Up_in_cells.All.fa


# Run Global GCcontent and MFE
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.Up_in_EVs.fa  --OutFile results/Mir.Up_in_EVs.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.Up_in_EVs.No_mirtrons.fa  --OutFile results/Mir.Up_in_EVs.No_mirtrons.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.Up_in_cells.Top100.fa  --OutFile results/Mir.Up_in_cells.Top100.folding.csv
python scripts/RNAfoldMeasures.py --FastaFile inputs/Mir.Up_in_cells.All.fa  --OutFile results/Mir.Up_in_cells.All.folding.csv
Mir.Up_in_EVs <- read.csv('output_files/Mir.Up_in_EVs.folding.csv', row.names=1)
Mir.Up_in_EVs.No_mirtrons <- read.csv('output_files/Mir.Up_in_EVs.No_mirtrons.folding.csv', row.names=1)
Mir.Up_in_cells.Top100 <- read.csv('output_files/Mir.Up_in_cells.Top100.folding.csv', row.names=1)
Mir.Up_in_cells.All <- read.csv('output_files/Mir.Up_in_cells.All.folding.csv', row.names=1)


head(Mir.Up_in_EVs)
##                              Gcontent GCcontent   MFE MFE_G4    Zscore
## chr11:78073004-78073070(+)  0.2424242 0.3787879 -28.6  -28.6 -6.653518
## chr11:79726399-79726511(+)  0.2589286 0.4553571 -39.7  -39.7 -4.515450
## chr1:153899925-153900036(+) 0.2432432 0.6126126 -38.3  -38.3 -2.803481
## chr8:23142572-23142662(-)   0.3111111 0.6555556 -63.8  -63.8 -9.624782
## chr15:84951525-84951623(-)  0.2959184 0.6122449 -51.6  -51.6 -4.774950
## chr19:6264642-6264728(+)    0.3255814 0.6046512 -47.3  -47.3 -5.646973
##                             Zscore_G4
## chr11:78073004-78073070(+)  -7.388169
## chr11:79726399-79726511(+)  -4.180001
## chr1:153899925-153900036(+) -2.675875
## chr8:23142572-23142662(-)   -8.886640
## chr15:84951525-84951623(-)  -5.367994
## chr19:6264642-6264728(+)    -4.217688

ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x='EV', y=Gcontent)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x='EV no mirtrons', y=Gcontent)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x='Cells Top100', y=Gcontent)) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x='Cells All', y=Gcontent)) +
  labs(title='G content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'EV no mirtrons', 'Cells Top100', 'Cells All'))

ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x='EV', y=GCcontent)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x='EV no mirtrons', y=GCcontent)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x='Cells Top100', y=GCcontent)) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x='Cells All', y=GCcontent)) +
  labs(title='GC content', x='') +
  theme_bw() +
  scale_x_discrete(limits = c('EV', 'EV no mirtrons', 'Cells Top100', 'Cells All'))

custom_labels <- c('EV', 'EV no mirtrons', 'Cells Top100', 'Cells All')
ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x=1, y=MFE)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x=2, y=MFE)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x=3, y=MFE)) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x=4, y=MFE)) +
  labs(title='MFE', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x=1, y=MFE)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x=2, y=MFE)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x=3, y=MFE), outliers = FALSE) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x=4, y=MFE)) +
  labs(title='MFE', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'EV no mirtrons', 'Cells Top100', 'Cells All')
ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x=1, y=Zscore)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x=2, y=Zscore)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x=3, y=Zscore)) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x=4, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)

custom_labels <- c('EV', 'EV no mirtrons', 'Cells Top100', 'Cells All')
ggplot() +
  geom_boxplot(data = Mir.Up_in_EVs, aes(x=1, y=Zscore_G4)) +
  geom_boxplot(data = Mir.Up_in_EVs.No_mirtrons, aes(x=2, y=Zscore)) +
  geom_boxplot(data = Mir.Up_in_cells.Top100, aes(x=3, y=Zscore)) +
  geom_boxplot(data = Mir.Up_in_cells.All, aes(x=4, y=Zscore)) +
  geom_hline(aes(yintercept=1.96), col='black', linetype = 'dashed',)+
  geom_hline(aes(yintercept=-1.96), col='black', linetype = 'dashed',)+
  labs(title='Zscore G4', x='') +
  theme_bw() +
  scale_x_discrete(limits = custom_labels)