###### Selection de genes stables
######         2020.01.14 

# 1. Methode : 
#   a. Recherche de genes stables "graphiquement" séparemment dans chaque réplicat puis on sélectionne les genes en communs stable 
#   dans les trois réplicats
#   b. séparation des genes en trois set de genes selon expression
#   c. calcul du SD de ces gènes selectionnés dans chacun des sets
#   d. sélection des trois gènes de chaque set dont le sd est le plus faible et présent dans les trois réplicats 
#   e. plot de ces gènes pour vérifier. Eventuellement en enlever quelque uns 

# 2. load librairies and set env   
setwd("~/RMI2/Projet_TDD/20200114_Normalisation")
library(ggplot2)
library(reshape2)
library(DESeq2)
library(stringr)
library(dplyr)

all_data=read.delim("data/HTSeq_count_stats_all_libraries.csv", header=TRUE, sep=",", row.names=1, stringsAsFactors = TRUE) 

####### MACROPHAGES #######
# 3. Macro actives --------------------------------------------------------

tdata <- data.frame(all_data[-(1:5),
                             grep("^LPS(no)?_macro_(0h_Triptolide|3h_Triptolide)_(m|i).*exon$",colnames(all_data), perl = TRUE)])
data = apply(tdata, 2, function(x) as.numeric(as.integer(x)))
rownames(data) = rownames(tdata)

# if I want to calculate a sd, i need use RPM values
RPMfactor = apply(data, 2, sum)
RPMfactor = 100000/RPMfactor

dataRPM = data
for (i in 1:ncol(data))
{
  dataRPM[,i] = data[,i]*RPMfactor[i]
}


# replicate 1 
data_rep1 = dataRPM[,c("LPS_macro_0h_Triptolide_m3_4_exon","LPS_macro_3h_Triptolide_m3_6_exon")]
plot(x=log10(data_rep1[,1]), y=log10(data_rep1[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep1", ylab = "3h Triptolide rep1",
     main = "Normalization verification in Activated macrophage")
abline(a=0, b=1, col="dark blue")

rep1_data_stable_act = data_rep1[data_rep1[,1] > 0.2 & 
                                 data_rep1[,1] + 0.3 * data_rep1[,1] < data_rep1[,2] &
                                 data_rep1[,2] - 0.9 * data_rep1[,2] < data_rep1[,1],]

points(x=log10(rep1_data_stable_act[,1]), y=log10(rep1_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 2 
data_rep2 = dataRPM[,c("LPS_macro_0h_Triptolide_m4_4_exon","LPS_macro_3h_Triptolide_m4_6_exon")]
plot(x=log10(data_rep2[,1]), y=log10(data_rep2[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Activated macrophage")
abline(a=0, b=1, col="dark blue")

rep2_data_stable_act = data_rep2[data_rep2[,1] > 0.2 & 
                                 data_rep2[,1] + 0.3 * data_rep2[,1] < data_rep2[,2] &
                                 data_rep2[,2] - 0.9 * data_rep2[,2] < data_rep2[,1],]

points(x=log10(rep2_data_stable_act[,1]), y=log10(rep2_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 3 
data_rep3 = dataRPM[,c("LPS_macro_0h_Triptolide_m5_4_exon","LPS_macro_3h_Triptolide_m5_6_exon")]
plot(x=log10(data_rep3[,1]), y=log10(data_rep3[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Activated macrophage")
abline(a=0, b=1, col="dark blue")

rep3_data_stable_act = data_rep3[data_rep3[,1] > 0.2 &
                                 data_rep3[,1] + 0.3 * data_rep3[,1] < data_rep3[,2] &
                                 data_rep3[,2] - 0.9 * data_rep3[,2] < data_rep3[,1],]

points(x=log10(rep3_data_stable_act[,1]), y=log10(rep3_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

macroAStable <- intersect(intersect(rownames(rep1_data_stable_act),rownames(rep2_data_stable_act)),rownames(rep3_data_stable_act))

# 4. Macro resting --------------------------------------------------------
# replicate 1 
data_rep1 = dataRPM[,c("LPSno_macro_0h_Triptolide_i3_4_exon","LPSno_macro_3h_Triptolide_i3_6_exon")]
plot(x=log10(data_rep1[,1]), y=log10(data_rep1[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep1", ylab = "3h Triptolide rep1",
     main = "Normalization verification in Resting macrophage")
abline(a=0, b=1, col="dark blue")
rep1_data_stable_act = data_rep1[data_rep1[,1] > 0.2 &
                                 data_rep1[,1] + 0.3 * data_rep1[,1] < data_rep1[,2] &
                                 data_rep1[,2] - 0.9 * data_rep1[,2] < data_rep1[,1],]
points(x=log10(rep1_data_stable_act[,1]), y=log10(rep1_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 2 
data_rep2 = dataRPM[,c("LPSno_macro_0h_Triptolide_i4_4_exon","LPSno_macro_3h_Triptolide_i4_6_exon")]
plot(x=log10(data_rep2[,1]), y=log10(data_rep2[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Resting macrophage")
abline(a=0, b=1, col="dark blue")
rep2_data_stable_act = data_rep2[data_rep2[,1] > 0.2 &
                                   data_rep2[,1] + 0.3 * data_rep2[,1] < data_rep2[,2] &
                                   data_rep2[,2] - 0.9 * data_rep2[,2] < data_rep2[,1],]
points(x=log10(rep2_data_stable_act[,1]), y=log10(rep2_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 3 
data_rep3 = dataRPM[,c("LPSno_macro_0h_Triptolide_i5_4_exon","LPSno_macro_3h_Triptolide_i5_6_exon")]
plot(x=log10(data_rep3[,1]), y=log10(data_rep3[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Activated macrophage")
abline(a=0, b=1, col="dark blue")

rep3_data_stable_act = data_rep3[data_rep3[,1] > 0.2 & 
                                   data_rep3[,1] + 0.3 * data_rep3[,1] < data_rep3[,2] &
                                   data_rep3[,2] - 0.9 * data_rep3[,2] < data_rep3[,1],]

points(x=log10(rep3_data_stable_act[,1]), y=log10(rep3_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

macroRStable <- intersect(intersect(rownames(rep1_data_stable_act),rownames(rep2_data_stable_act)),rownames(rep3_data_stable_act))
macroStable <- intersect(macroRStable, macroAStable)

# sd calcul and gene selection
db = as.data.frame(dataRPM[rownames(dataRPM) %in% macroStable,])
db$genes <- rownames(db)
db = melt(db, id.vars = "genes")

set1  = db[db$value<=10&db$value>1,]
set2 = db[db$value<=100&db$value>10,]
set3 = db[db$value<=1000&db$value>100,]

set1 <- group_by(set1, genes)
set1 <- summarise(set1, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set1 = set1[order(-set1$n,set1$sd),]

set2 <- group_by(set2, genes)
set2 <- summarise(set2, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set2 = set2[order(-set2$n,set2$sd),]

set3 <- group_by(set3, genes)
set3 <- summarise(set3, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set3 = set3[order(-set3$n,set3$sd),]

stable = dataRPM[c(as.character(set1$genes[1:13]),
                   as.character(set2$genes[1:4]), 
                   as.character(set3$genes[1:2])
),]
pdf(file = "results/00_selected_stables_genes_macrophages.pdf")
# resting
plot(x = log10(dataRPM[,"LPSno_macro_0h_Triptolide_i3_4_exon"]), y = log10(dataRPM[,"LPSno_macro_3h_Triptolide_i3_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Resting macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPSno_macro_0h_Triptolide_i3_4_exon"]),
       y = log10(stable[,"LPSno_macro_3h_Triptolide_i3_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"LPSno_macro_0h_Triptolide_i4_4_exon"]), y = log10(dataRPM[,"LPSno_macro_3h_Triptolide_i4_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep4", ylab = "3h Triptolide rep4",
     main = "Normalization verification in Resting macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPSno_macro_0h_Triptolide_i4_4_exon"]),
       y = log10(stable[,"LPSno_macro_3h_Triptolide_i4_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"LPSno_macro_0h_Triptolide_i5_4_exon"]), y = log10(dataRPM[,"LPSno_macro_3h_Triptolide_i5_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep5", ylab = "3h Triptolide rep5",
     main = "Normalization verification in Resting macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPSno_macro_0h_Triptolide_i5_4_exon"]),
       y = log10(stable[,"LPSno_macro_3h_Triptolide_i5_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

# activated
plot(x = log10(dataRPM[,"LPS_macro_0h_Triptolide_m3_4_exon"]), y = log10(dataRPM[,"LPS_macro_3h_Triptolide_m3_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in activated macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPS_macro_0h_Triptolide_m3_4_exon"]),
       y = log10(stable[,"LPS_macro_3h_Triptolide_m3_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"LPS_macro_0h_Triptolide_m4_4_exon"]), y = log10(dataRPM[,"LPS_macro_3h_Triptolide_m4_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep4", ylab = "3h Triptolide rep4",
     main = "Normalization verification in activated macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPS_macro_0h_Triptolide_m4_4_exon"]),
       y = log10(stable[,"LPS_macro_3h_Triptolide_m4_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"LPS_macro_0h_Triptolide_m5_4_exon"]), y = log10(dataRPM[,"LPS_macro_3h_Triptolide_m5_6_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep5", ylab = "3h Triptolide rep5",
     main = "Normalization verification in activated macrophages")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"LPS_macro_0h_Triptolide_m5_4_exon"]),
       y = log10(stable[,"LPS_macro_3h_Triptolide_m5_6_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)
dev.off()
# LYMPOCYTES -----------------------------------------------------
# Lympho Activated

tdata <- data.frame(all_data[-(1:5),
                             grep("^(R|A).*(0h_Triptolide|3h_Triptolide)_exon$",colnames(all_data), perl = TRUE)])
data = apply(tdata, 2, function(x) as.numeric(as.integer(x)))
rownames(data) = rownames(tdata)

# if I want to calculate a sd, i need use RPM values
RPMfactor = apply(data, 2, sum)
RPMfactor = 100000/RPMfactor

dataRPM = data
for (i in 1:ncol(data))
{
  dataRPM[,i] = data[,i]*RPMfactor[i]
}


# replicate 1 
data_rep1 = dataRPM[,c("A2_7_Activated_0h_Triptolide_exon","A2_9_Activated_3h_Triptolide_exon")]
plot(x=log10(data_rep1[,1]), y=log10(data_rep1[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep1", ylab = "3h Triptolide rep1",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

rep1_data_stable_act = data_rep1[data_rep1[,1]>0.2 & 
                                   data_rep1[,1]+0.4*data_rep1[,1]<data_rep1[,2] &
                                   data_rep1[,2]-0.90*data_rep1[,2] < data_rep1[,1],]

points(x=log10(rep1_data_stable_act[,1]), y=log10(rep1_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 2 
data_rep2 = dataRPM[,c("A3_7_Activated_0h_Triptolide_exon","A3_9_Activated_3h_Triptolide_exon")]
plot(x=log10(data_rep2[,1]), y=log10(data_rep2[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

rep2_data_stable_act = data_rep2[data_rep2[,1]>0.2 & 
                                   data_rep2[,1]+0.7*data_rep2[,1]<data_rep2[,2] & 
                                   data_rep2[,2]-0.9*data_rep2[,2] < data_rep2[,1],]

points(x=log10(rep2_data_stable_act[,1]), y=log10(rep2_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 3 
data_rep3 = dataRPM[,c("A4_7_Activated_0h_Triptolide_exon","A4_9_Activated_3h_Triptolide_exon")]
plot(x=log10(data_rep3[,1]), y=log10(data_rep3[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

rep3_data_stable_act = data_rep3[data_rep3[,1]>0.2 &
                                   data_rep3[,1]+0.6*data_rep3[,1]<data_rep3[,2] &
                                   data_rep3[,2]-0.9*data_rep3[,2] < data_rep3[,1],]

points(x=log10(rep3_data_stable_act[,1]), y=log10(rep3_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

LymphoAStable <- intersect(intersect(rownames(rep1_data_stable_act),rownames(rep2_data_stable_act)),rownames(rep3_data_stable_act))
stablesLymphoA = rownames(stable)

# 6. Lympho Resting -------------------------------------------------------
# replicate 1 
data_rep1 = dataRPM[,c("R2_4_Resting_0h_Triptolide_exon","R2_6_Resting_3h_Triptolide_exon")]
plot(x=log10(data_rep1[,1]), y=log10(data_rep1[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep1", ylab = "3h Triptolide rep1",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

rep1_data_stable_act = data_rep1[data_rep1[,1]>0.2 & 
                                   data_rep1[,1]+0.5*data_rep1[,1]<data_rep1[,2] & 
                                   data_rep1[,2]-0.9*data_rep1[,2] < data_rep1[,1],]

points(x=log10(rep1_data_stable_act[,1]), y=log10(rep1_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 2 
data_rep2 = dataRPM[,c("R3_4_Resting_0h_Triptolide_exon","R3_6_Resting_3h_Triptolide_exon")]
plot(x=log10(data_rep2[,1]), y=log10(data_rep2[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

rep2_data_stable_act = data_rep2[data_rep2[,1]>0.2 & 
                                   data_rep2[,1]+0.6*data_rep2[,1]<data_rep2[,2] &
                                   data_rep2[,2]-0.9*data_rep2[,2] < data_rep2[,1],]

points(x=log10(rep2_data_stable_act[,1]), y=log10(rep2_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

# replicate 3 
data_rep3 = dataRPM[,c("R4_4_Resting_0h_Triptolide_exon","R4_6_Resting_3h_Triptolide_exon")]
plot(x=log10(data_rep3[,1]), y=log10(data_rep3[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

rep3_data_stable_act = data_rep3[data_rep3[,1]>0.2 & 
                                   data_rep3[,1]+0.6*data_rep3[,1]<data_rep3[,2] & 
                                   data_rep3[,2]-0.9*data_rep3[,2] < data_rep3[,1],]

points(x=log10(rep3_data_stable_act[,1]), y=log10(rep3_data_stable_act[,2]),
       col = alpha("red", 1) , pch = 20, cex = 1)
legend("bottomright", legend = "stables genes", col = c("red"), pch = 20 )

plot(x=log10(data_rep1[,1]), y=log10(data_rep1[,2]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep1", ylab = "3h Triptolide rep1",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

LymphoRStable <- intersect(intersect(rownames(rep1_data_stable_act),rownames(rep2_data_stable_act)),rownames(rep3_data_stable_act))
LymphoStable <- intersect(LymphoRStable, LymphoAStable)

# sd calcul and gene selection
db = as.data.frame(dataRPM[rownames(dataRPM) %in% LymphoStable,])
db$genes <- rownames(db)
db = melt(db, id.vars = "genes")

set1  = db[db$value<=10&db$value>1,]
set2 = db[db$value<=100&db$value>10,]
set3 = db[db$value<=1000&db$value>100,]

set1 <- group_by(set1, genes)
set1 <- summarise(set1, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set1 = set1[order(-set1$n,set1$sd),]

set2 <- group_by(set2, genes)
set2 <- summarise(set2, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set2 = set2[order(-set2$n,set2$sd),]

set3 <- group_by(set3, genes)
set3 <- summarise(set3, mean=mean(value), sd = sd(value)/mean(value), n = length(value))
set3 = set3[order(-set3$n,set3$sd),]

stable = dataRPM[c(as.character(set1$genes[1:10]),
                   as.character(set2$genes[1:10]), 
                   as.character(set3$genes[1:10])
),]


# resting

pdf(file = "results/00_selected_stables_genes_lymphocytes.pdf")
plot(x = log10(dataRPM[,"R2_4_Resting_0h_Triptolide_exon"]), y = log10(dataRPM[,"R2_6_Resting_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"R2_4_Resting_0h_Triptolide_exon"]),
       y = log10(stable[,"R2_6_Resting_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"R3_4_Resting_0h_Triptolide_exon"]), y = log10(dataRPM[,"R3_6_Resting_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"R3_4_Resting_0h_Triptolide_exon"]),
       y = log10(stable[,"R3_6_Resting_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"R4_4_Resting_0h_Triptolide_exon"]), y = log10(dataRPM[,"R4_6_Resting_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep4", ylab = "3h Triptolide rep4",
     main = "Normalization verification in Resting Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"R4_4_Resting_0h_Triptolide_exon"]),
       y = log10(stable[,"R4_6_Resting_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

# activated

plot(x = log10(dataRPM[,"A2_7_Activated_0h_Triptolide_exon"]), y = log10(dataRPM[,"A2_9_Activated_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep2", ylab = "3h Triptolide rep2",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"A2_7_Activated_0h_Triptolide_exon"]),
       y = log10(stable[,"A2_9_Activated_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"A3_7_Activated_0h_Triptolide_exon"]), y = log10(dataRPM[,"A3_9_Activated_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep3", ylab = "3h Triptolide rep3",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"A3_7_Activated_0h_Triptolide_exon"]),
       y = log10(stable[,"A3_9_Activated_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)

plot(x = log10(dataRPM[,"A4_7_Activated_0h_Triptolide_exon"]), y = log10(dataRPM[,"A4_9_Activated_3h_Triptolide_exon"]),
     col = alpha("black", 0.1) , pch = 20, cex = 1,
     xlab = "0h Triptolide rep4", ylab = "3h Triptolide rep4",
     main = "Normalization verification in Activated Lympho")
abline(a=0, b=1, col="dark blue")

points(x = log10(stable[,"A4_7_Activated_0h_Triptolide_exon"]),
       y = log10(stable[,"A4_9_Activated_3h_Triptolide_exon"]),
       col = alpha("red", 1) , pch = 20, cex = 1)
dev.off()

# 7. Save Data ------------------------------------------------------------

stablesGenes = list(LymphoStable,macroStable)
save(stablesGenes, file = "results/00_stableGenesInEachCondition.RData")