# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: tcgaAnalysisPythonR
#     language: python
#     name: tcgaanalysispythonr
# ---

import os, numpy, warnings
import pandas as pd

os.environ['R_HOME'] = '/home/gdpoore/anaconda3/envs/tcgaAnalysisPythonR/lib/R'
warnings.filterwarnings('ignore')
# %config InlineBackend.figure_format = 'retina'

# %reload_ext rpy2.ipython

# + language="R"
#
# require(ggplot2)
# require(snm)
# require(limma)
# require(edgeR)
# require(dplyr)
# require(edgeR)
# require(pvca)
# require(lme4)
# require(doMC)
#
# numCores <- detectCores()
# registerDoMC(cores=numCores)

# + language="R"
# load("validationSplitVbDataAndMetadata.RData")

# + language="R"
# print(dim(split1MetadataQC))
# print(dim(split2MetadataQC))
#
# print(dim(split1VbDataQC))
# print(dim(split2VbDataQC))

# + language="R"
# qcMetadataS1 <- split1MetadataQC # metadataSamplesAllQCAML
# qcDataS1 <- split1VbDataQC # vbDataBarnDFReconciledQCAML
#
# # Set up design matrix
# covDesignNormS1 <- model.matrix(~0 + sample_type +
#                                   data_submitting_center_label +
#                                   platform +
#                                   experimental_strategy +
#                                   tissue_source_site_label +
#                                   portion_is_ffpe,
#                               data = qcMetadataS1)
# print(colnames(covDesignNormS1))
# colnames(covDesignNormS1) <- gsub('([[:punct:]])|\\s+','',colnames(covDesignNormS1))
# print(colnames(covDesignNormS1))
#
# # Set up counts matrix
# countsS1 <- t(qcDataS1) # DGEList object from a table of counts (rows=features, columns=samples)
#
# # Normalize using edgeR and then plug into voom
# dgeS1 <- DGEList(counts = countsS1)
# keep <- filterByExpr(dgeS1, covDesignNormS1)
# dgeS1 <- dgeS1[keep,,keep.lib.sizes=FALSE]
# print("Now normalizing data...")
# dgeS1 <- calcNormFactors(dgeS1, method = "TMM")
# print("Now applying voom on normalized data...")
# vdgeS1 <- voom(dgeS1, design = covDesignNormS1, plot = TRUE, save.plot = TRUE, normalize.method="none")

# + language="R"
# qcMetadataS2 <- split2MetadataQC # metadataSamplesAllQCAML
# qcDataS2 <- split2VbDataQC # vbDataBarnDFReconciledQCAML
#
# # Set up design matrix
# covDesignNormS2 <- model.matrix(~0 + sample_type +
#                                   data_submitting_center_label +
#                                   platform +
#                                   experimental_strategy +
#                                   tissue_source_site_label +
#                                   portion_is_ffpe,
#                               data = qcMetadataS2)
# print(colnames(covDesignNormS2))
# colnames(covDesignNormS2) <- gsub('([[:punct:]])|\\s+','',colnames(covDesignNormS2))
# print(colnames(covDesignNormS2))
#
# # Set up counts matrix
# countsS2 <- t(qcDataS2) # DGEList object from a table of counts (rows=features, columns=samples)
#
# # Normalize using edgeR and then plug into voom
# dgeS2 <- DGEList(counts = countsS2)
# keep <- filterByExpr(dgeS2, covDesignNormS2)
# dgeS2 <- dgeS2[keep,,keep.lib.sizes=FALSE]
# print("Now normalizing data...")
# dgeS2 <- calcNormFactors(dgeS2, method = "TMM")
# print("Now applying voom on normalized data...")
# vdgeS2 <- voom(dgeS2, design = covDesignNormS2, plot = TRUE, save.plot = TRUE, normalize.method="none")

# + language="R"
#
# print(dim(t(vdgeS1$E)))
# print(dim(t(vdgeS2$E)))

# + language="R"
# taxaIntersect <- intersect(colnames(t(vdgeS1$E)), colnames(t(vdgeS2$E)))
# print(length(taxaIntersect))
#
# vdgeS1Filt <- t(vdgeS1$E)[,taxaIntersect]
# vdgeS2Filt <- t(vdgeS2$E)[,taxaIntersect]
# print(dim(vdgeS1Filt))
# print(dim(vdgeS2Filt))

# + language="R"
#
# # Apply
# bio.var.sample.typeS1 <- model.matrix(~sample_type, #sample_type, # histological_diagnosis_label and disease_type tried but cause function to fail
#                         data=qcMetadataS1)
# adj.varS1 <- model.matrix(~data_submitting_center_label +
#                           platform +
#                           experimental_strategy +
#                           tissue_source_site_label +
#                           portion_is_ffpe,
#                         data=qcMetadataS1)
# colnames(bio.var.sample.typeS1) <- gsub('([[:punct:]])|\\s+','',colnames(bio.var.sample.typeS1))
# colnames(adj.varS1) <- gsub('([[:punct:]])|\\s+','',colnames(adj.varS1))
# print(dim(adj.varS1))
# print(dim(bio.var.sample.typeS1))
# print(dim(vdgeS1Filt))
# print(dim(covDesignNormS1))

# + language="R"
#
# # Apply
# bio.var.sample.typeS2 <- model.matrix(~sample_type, #sample_type, # histological_diagnosis_label and disease_type tried but cause function to fail
#                         data=qcMetadataS2)
# adj.varS2 <- model.matrix(~data_submitting_center_label +
#                           platform +
#                           experimental_strategy +
#                           tissue_source_site_label +
#                           portion_is_ffpe,
#                         data=qcMetadataS2)
# colnames(bio.var.sample.typeS2) <- gsub('([[:punct:]])|\\s+','',colnames(bio.var.sample.typeS2))
# colnames(adj.varS2) <- gsub('([[:punct:]])|\\s+','',colnames(adj.varS2))
# print(dim(adj.varS2))
# print(dim(bio.var.sample.typeS2))
# print(dim(vdgeS2Filt))
# print(dim(covDesignNormS2))

# + language="R"
# snmDataObjSampleTypeWithExpStrategyFAS1 <- snm(raw.dat = t(vdgeS1Filt), 
#                                             bio.var = bio.var.sample.typeS1, 
#                                             adj.var = adj.varS1, 
#                                             rm.adj=TRUE,
#                                             verbose = TRUE,
#                                             diagnose = TRUE)
# snmDataSampleTypeWithExpStrategyFAS1 <- t(snmDataObjSampleTypeWithExpStrategyFAS1$norm.dat)
#
# print(dim(snmDataSampleTypeWithExpStrategyFAS1))

# + language="R"
# snmDataObjSampleTypeWithExpStrategyFAS2 <- snm(raw.dat = t(vdgeS2Filt), 
#                                             bio.var = bio.var.sample.typeS2, 
#                                             adj.var = adj.varS2, 
#                                             rm.adj=TRUE,
#                                             verbose = TRUE,
#                                             diagnose = TRUE)
# snmDataSampleTypeWithExpStrategyFAS2 <- t(snmDataObjSampleTypeWithExpStrategyFAS2$norm.dat)
#
# print(dim(snmDataSampleTypeWithExpStrategyFAS2))

# + language="R"
# save(snmDataSampleTypeWithExpStrategyFAS1, snmDataSampleTypeWithExpStrategyFAS2,
#      split1MetadataQC, split2MetadataQC, 
#      file = "validationSplits_snmDataSampleTypeWithExpStrategyFAS1S2.RData")
# -

# # PCA plotting to visually examine batch effects and batch correction

# + language="R"
# pcaPlotting <- function(pcaObject,pcChoices, dataLabels, factorString, titleString){
#   require(ggbiplot)
#   theme_update(plot.title = element_text(hjust = 0.5))
#   g <- ggbiplot(pcaObject,pcChoices, obs.scale = 1, var.scale = 1,
#                 groups = dataLabels, ellipse = TRUE,
#                 alpha = 0.2,
#                 circle = TRUE,var.axes=FALSE) + 
#     scale_color_discrete(name = factorString) +
#     theme_bw() + 
#     #theme(legend.direction = "horizontal", legend.position = "top") +
#     ggtitle(titleString) + theme(plot.title = element_text(hjust = 0.5))
#   
#   print(g)
# }

# + language="R"
# unnormalizedPCAPlotFA <- pcaPlotting(pcaObject = prcomp(t(vdge$E)),
#                                    pcChoices = c(1,2),
#                                    dataLabels = qcMetadata$data_submitting_center_label,
#                                    factorString = "Batch",
#                                    titleString = "PCA w/o Batch Correction")

# + language="R"
# snmPCAPlotSampleTypeFA <- pcaPlotting(pcaObject = prcomp(snmDataSampleTypeWithExpStrategyFA),
#                           pcChoices = c(1,2),
#                           dataLabels = qcMetadata$data_submitting_center_label,
#                           factorString = "Sequencing Center",
#                           titleString = "PCA w/ SNM Correction\n(Target: Sample Type)")

# +
# # %%R 
# snmPCAPlotGender <- pcaPlotting(pcaObject = prcomp(snmDataGenderWithAML),
#                           pcChoices = c(1,2),
#                           dataLabels = qcMetadata$data_submitting_center_label,
#                           factorString = "Sequencing Center",
#                           titleString = "PCA w/ SNM Correction\n(Target: Gender)")

# + language="R"
# ggsave(plot = unnormalizedPCAPlotFA, 
#        filename = "unnormalizedPCAPlotFA_DecreasedOpacity.png",
#        width = 16.2,
#        height = 5.29,
#        units = "in",
#        dpi = "retina")
#
# ggsave(plot = snmPCAPlotSampleTypeFA, 
#        filename = "snmPCAPlotSampleTypeFA_DecreasedOpacity.png",
#        width = 16.2,
#        height = 5.29,
#        units = "in",
#        dpi = "retina")
#
# # save(snmDataGenderWithAML, metadataSamplesAllQCAML, 
# #      vbDataBarnDFReconciledQCAML, 
# #      file = "amlVbDataAndMetadataAndSNMByGender.RData")

# +
# # %%R
# snmDataObjGenderWithAML <- snm(raw.dat = vdge$E, 
#                                 bio.var = bio.var.gender, 
#                                 adj.var = adj.var, 
#                                 rm.adj=TRUE,
#                                 verbose = TRUE,
#                                 diagnose = TRUE)
# snmDataGenderWithAML <- t(snmDataObjGenderWithAML$norm.dat)

# print(dim(snmDataGenderWithAML))
# -

# # PVCA using key filtered metadata features (i.e. narrowing down the extended version of this)

# + language="R"
# # Implement PVCA
# # From extended model, remove variables that contribute very little if at all:
# # ethnicity, gender, reference_genome
# pct_threshold <- 0.8
# metaPVCAExtendedFiltered <- metadataSamplesAllQC[,c("sample_type",
#                                             "disease_type",
#                                             "data_submitting_center_label",
#                                             "platform",
#                                             "experimental_strategy", 
#                                             "tissue_source_site_label",
#                                             "portion_is_ffpe")]
# print(dim(metaPVCAExtendedFiltered))
# print(dim(snmDataSampleTypeWithExpStrategy))
# print(dim(vbDataBarnDFReconciledQC))

# + language="R"
# pvcaVbRawNoVoomNoSNM_ExtendedFiltered_FA <- PVCA(counts = t(vbDataBarnDFReconciledQC), 
#                                                meta = metaPVCAExtendedFiltered, 
#                                                threshold = pct_threshold,
#                                                inter = FALSE)
# save(pvcaVbRawNoVoomNoSNM_ExtendedFiltered_FA, file = "pvcaVbRawNoVoomNoSNM_ExtendedFiltered_FA.RData")
# PlotPVCA(pvcaVbRawNoVoomNoSNM_ExtendedFiltered_FA, "Raw count data")

# + language="R"
# pvcaVoomNoSNM_ExtendedFiltered_FA <- PVCA(counts = vdge$E,
#                                            meta = metaPVCAExtendedFiltered,
#                                            threshold = pct_threshold,
#                                            inter = FALSE)
# save(pvcaVoomNoSNM_ExtendedFiltered_FA, file = "pvcaVoomNoSNM_ExtendedFiltered_FA.RData")
# PlotPVCA(pvcaVoomNoSNM_ExtendedFiltered_FA, "Voom Normalized")

# + language="R"
# pvcaSampleWithExpStrategySNM_ExtendedFiltered_FA <- PVCA(counts = t(snmDataSampleTypeWithExpStrategyFA), 
#                                                        meta = metaPVCAExtendedFiltered,
#                                                        threshold = pct_threshold,
#                                                        inter = FALSE)
# save(pvcaSampleWithExpStrategySNM_ExtendedFiltered_FA, 
#      file = "pvcnoaSampleWithExpStrategySNM_ExtendedFiltered_FA.RData")
# PlotPVCA(pvcaSampleWithExpStrategySNM_ExtendedFiltered_FA, 
#          "Voom Normalized & SNM Corrected Plus Exp Strategy (Target is Sample Type)")

# + language="R"
# 1+2
# -

# # Examining sample and taxa ratio changes due to batch correction

# + language="R"
# require(ggplot2)
# require(matrixStats)
# divSNMDataSampleType <- snmDataSampleType / t(snmDataObjSampleType$raw.dat)
# taxaMedians <- data.frame(Medians = colMedians(divSNMDataSampleType), 
#                           Taxa = colnames(divSNMDataSampleType),
#                           pval = factor(ifelse(snmDataObjSampleType$pval <=0.05, 
#                                                yes = "P-value <= 0.05", no = "P-value > 0.05")))
# sampleMedians <- data.frame(Medians = rowMedians(divSNMDataSampleType), 
#                             Samples = rownames(divSNMDataSampleType),
#                             SeqCenter = metadataSamplesAllQC$data_submitting_center_label,
#                             SampleType = metadataSamplesAllQC$sample_type,
#                             CancerType = metadataSamplesAllQC$disease_type)
# gt <- ggplot(taxaMedians, aes(x = reorder(Taxa, -Medians), y = Medians, fill = pval)) + 
# geom_bar(stat = "identity") +
# theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
# labs(y = "Median of Normalizing Ratios Per Taxa", x = "Samples", fill = "ANOVA Result Per Taxa") 
#
# gs <- ggplot(sampleMedians, aes(x = reorder(Samples, -Medians), y = Medians, fill = CancerType)) + 
#   geom_bar(stat = "identity") + coord_flip() +
#   theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
#   scale_y_log10() + labs(y = "Median of Normalizing Ratios Per Sample", x = "Samples", fill='Cancer Type') 

# + language="R"
# gt

# + language="R"
# ggsave(plot = gt, 
#        filename = "snmNormMedianPerTaxaPval.png",
#        width = 8.5,
#        height = 6,
#        units = "in",
#        dpi = "retina")

# + language="R"
# require(pheatmap)
# pheatmap(snmDataSampleTypeLMFit$coefficients,
#          clustering_distance_rows = "correlation",
#          clustering_distance_cols = "correlation",
#          show_rownames = FALSE,
#          show_colnames = FALSE,
#          filename = "snmLMFitCoefCorr.png")

# +
# # %%R
# save(snmDataObjPathStage, snmDataPathStage, metadataSamplesAllQCPath, file = "snmResultsPathBinned.RData")
# -



StarcoderdataJupyterScriptsDedupFiltered