source("/Users/daga/Documents/Projects/Universal_scripts/Diffbind_visualization_functions.R")
source(file = "/Users/daga/Documents/Projects/Autoimmune/SLE/Scripts/sva_correction_de_analysis_fn.R")
setwd("/Users/daga/Documents/Projects/Autoimmune_all/Scripts/")
library("apeglm")
################# Functions
######### Running differential Functions
to_run_save_de = function (raw_count_data_matrix,metatable,experiment_design,outfn_maplot,outqc,outfn_de_result) {
dds = DESeqDataSetFromMatrix(countData = raw_count_data_matrix,
colData = metatable,
design = experiment_design)
############### For visualisation
dds = DESeq(dds)
res <- results(dds)
##### Correction for p-value
de_res_correct = to_get_corrected_pval_result_table(de_result = res,
plot_hist = FALSE,out_fn = "")
###### SHrinkage
resLFC <- lfcShrink(dds, coef="Condition_Patient_vs_Healthy", type="apeglm")
resNorm <- lfcShrink(dds, coef="Condition_Patient_vs_Healthy", type="normal")
pdf (outfn_maplot,width = 12,height = 6)
par(mfrow=c(1,3))
maplot(res = res,labelsig = FALSE,ylim = c(-3,3),main = "Deseq2")
maplot(res = resLFC ,labelsig = FALSE,ylim = c(-3,3), main = "apelgm")
maplot(res = resNorm,labelsig = FALSE,ylim = c(-3,3), main = "normal")
dev.off()
##### QC plots
pdf(outqc,width = 12,height = 6)
hist(res$pvalue, col="grey", main = "res")
maplot(res,main="res MA Plot",labelsig = FALSE)
volcanoplot(res, lfcthresh=1, sigthresh=0.05, textcx=.8, main = "res")
hist(de_res_correct$pvalue, col="grey", main = "res_correct")
maplot(de_res_correct, main="res_correct MA Plot",labelsig = FALSE)
volcanoplot(de_res_correct, lfcthresh=1, sigthresh=0.05, textcx=.8, main = "res_correct")
dev.off()
## Saving
save(list = c("res","resNorm","de_res_correct"),file = outfn_de_result)
}
###### SLE_only (rread consensus meta creted for diffbind)
disease_meta = read.table("../Consensus_meta_table/Psoriasis_SLE_Graves_JIASF_ex_consensus_metadata_diffbind.txt",
header = T,sep = "\t")
disease_meta[,c("Peaks","bamReads","PeakCaller")] = NULL
names(disease_meta) = c("Sample_id","Age","Condition","Batch","Disease")
disease_meta$Age_group = cut(disease_meta$Age,breaks = seq(0,100,20))
disease_meta$Age_group = factor(disease_meta$Age_group)
disease_meta$Disease_info = sapply(X = disease_meta$Disease,FUN = modify_string,delimtter = "_",position = 2,USE.NAMES = F)
common_Psoriasis_SLE_controls = c("EK19","EK2","EK3")
####### JIA Condition
disease_meta[grepl(pattern = "^JIA.*PB",disease_meta$Sample_id),"Condition"] = "Healthy"
############### Design_matrix
################## DDS rld object
raw_count_file = "../Consensus_peaks_raw_counts/Psoriais_SLE_Graves_JIASF_ex_ovelap3_with_nosummit_with_nosummit_std_chr_only.txt"
raw_count_data = read.table(file = raw_count_file,header = T,sep = "\t")
######### Modifying_table
raw_count_data$Peak_identity = paste(raw_count_data$CHR,":",raw_count_data$START,
"-",raw_count_data$END,sep = "")
rownames(raw_count_data) = raw_count_data$Peak_identity
raw_count_data[,c("Peak_identity","CHR","START","END")] = NULL
########### Subsetting for the each disease
SLE_samples = c(common_Psoriasis_SLE_controls,as.character(disease_meta[disease_meta$Disease_info == "SLE","Sample_id"]))
Psoriasis_samples = as.character(disease_meta[disease_meta$Disease_info == "Psoriasis","Sample_id"])
Graves_samples = as.character(disease_meta[disease_meta$Disease_info == "Graves","Sample_id"])
JIA_samples = as.character(disease_meta[disease_meta$Disease_info == "JIA","Sample_id"])
######## Mattrix
#raw_count_data_mat = as.matrix(raw_count_data)
#summary(as.numeric(raw_count_data_mat))
###### Subsetting df into each disease and performing disease specific DE analysis based on consensus peakset
SLE_raw_count_df = subset(x = raw_count_data,select = SLE_samples)
Psoriasis_raw_count_df = subset(x = raw_count_data,select = Psoriasis_samples)
Graves_raw_count_df = subset(x = raw_count_data,select = Graves_samples)
JIA_raw_count_df = subset(x = raw_count_data,select = JIA_samples)
#####Subsetting df into disease meta df
##To match metadata and raw count data cols for deseq2
SLE_meta_df = disease_meta[match(colnames(SLE_raw_count_df),disease_meta$Sample_id),]
all(SLE_meta_df$Sample_id == colnames(SLE_raw_count_df))
Psoriasis_meta_df = disease_meta[match(colnames(Psoriasis_raw_count_df),disease_meta$Sample_id),]
all(Psoriasis_meta_df$Sample_id == colnames(Psoriasis_raw_count_df))
Graves_meta_df = disease_meta[match(colnames(Graves_raw_count_df),disease_meta$Sample_id),]
all(Graves_meta_df$Sample_id == colnames(Graves_raw_count_df))
JIA_meta_df = disease_meta[match(colnames(JIA_raw_count_df),disease_meta$Sample_id),]
all(JIA_meta_df$Sample_id == colnames(JIA_raw_count_df))
############### SLE
SLE_design_matrix = ~ Age_group +Batch + Condition
to_run_save_de(raw_count_data_matrix = SLE_raw_count_df,
metatable = SLE_meta_df,
experiment_design =SLE_design_matrix,
outfn_maplot = "../Analysis_plot/maplots/SLE_consensus_peak_overlap3_withJIASF_de_maplots.pdf",
outfn_de_result = "../Diff_expression_results/SLE_consensus_peak_overlap3_withJIASF_de_result.RData",
outqc = "../Analysis_plot/maplots/SLE_consensus_peak_withJIASF_overlap3_de_qcplots.pdf")
###### Psoriasis
Psoriasis_design_matrix = ~ Batch + Age_group + Condition
to_run_save_de(raw_count_data_matrix = Psoriasis_raw_count_df,
metatable = Psoriasis_meta_df,
experiment_design =Psoriasis_design_matrix,
outfn_maplot = "../Analysis_plot/maplots/Psoriasis_consensus_peak_overlap3_withJIASF_de_maplots.pdf",
outfn_de_result = "../Diff_expression_results/Psoriasis_consensus_peak_overlap3_withJIASF_de_result.RData",
outqc = "../Analysis_plot/maplots/Psoriasis_consensus_peak_overlap3_withJIASF_de_qcplots.pdf")
####### JIA
JIA_design_matrix = ~ Batch + Condition
to_run_save_de(raw_count_data_matrix = JIA_raw_count_df,
metatable = JIA_meta_df,
experiment_design =JIA_design_matrix,
outfn_maplot = "../Analysis_plot/maplots/JIA_consensus_peak_overlap3_withJIASF_de_maplots.pdf",
outfn_de_result = "../Diff_expression_results/JIA_consensus_peak_overlap3_withJIASF_de_result.RData",
outqc = "../Analysis_plot/maplots/JIA_consensus_peak_overlap3_withJIASF_de_qcplots.pdf")
###### Graves
Graves_design_matrix = ~ Batch + Condition
to_run_save_de(raw_count_data_matrix = Graves_raw_count_df,
metatable = Graves_meta_df,
experiment_design =Graves_design_matrix,
outfn_maplot = "../Analysis_plot/maplots/Graves_consensus_peak_overlap3_withJIASF_de_maplots.pdf",
outfn_de_result = "../Diff_expression_results/Graves_consensus_peak_overlap3_withJIASF_de_result.RData",
outqc = "../Analysis_plot/maplots/Graves_consensus_peak_overlap3_withJIASF_de_qcplots.pdf")
############### Saving LFC shrinkage normal for all diseases in one df
load("../Diff_expression_results/SLE_consensus_peak_withJIASF_de_result.RData")
SLE_resNorm = resNorm
SLE_resNorm_df = as.data.frame(SLE_resNorm)
rm(list = c("res","resNorm","de_res_correct"))
load("../Diff_expression_results/Psoriasis_consensus_peak_withJIASF_de_result.RData")
Psoriasis_resNorm = resNorm
Psoriasis_resNorm_df = as.data.frame(Psoriasis_resNorm)
rm(list = c("res","resNorm","de_res_correct"))
load("../Diff_expression_results/Graves_consensus_peak_withJIASF_de_result.RData")
Graves_resNorm = resNorm
Graves_resNorm_df = as.data.frame(Graves_resNorm)
rm(list = c("res","resNorm","de_res_correct"))
load("../Diff_expression_results/JIA_consensus_peak_withJIASF_de_result.RData")
JIA_resNorm = resNorm
JIA_resNorm_df = as.data.frame(JIA_resNorm)
rm(list = c("res","resNorm","de_res_correct"))
######### Matching the rownames to make common df
all(rownames(SLE_resNorm_df) == rownames(Psoriasis_resNorm_df))
all(rownames(Graves_resNorm_df) == rownames(JIA_resNorm_df))
all(rownames(SLE_resNorm_df) == rownames(Graves_resNorm_df))
resNorm_l2FC_all_disease_df = data.frame("SLE_l2FC_normal" = SLE_resNorm_df$log2FoldChange,
"Psoriasis_l2FC_normal" = Psoriasis_resNorm_df$log2FoldChange,
"Graves_l2FC_normal" = Graves_resNorm_df$log2FoldChange,
"JIA_l2FC_normal" = JIA_resNorm_df$log2FoldChange)
rownames(resNorm_l2FC_all_disease_df) = rownames(SLE_resNorm_df)
save(list = c("resNorm_l2FC_all_disease_df"),file = "../Diff_expression_results/resNorm_l2FC_all_disease_withJIASF.Rdata")