library(TCGAbiolinks)
library(SummarizedExperiment)

query_GBM <- GDCquery(project = "TCGA-GBM",
                      data.category = "Transcriptome Profiling",
                      data.type = "Gene Expression Quantification",
                      workflow.type = "HTSeq - FPKM-UQ",
                      legacy = FALSE,
                      sample.type = c("Primary Tumor", "Recurrent Tumor"),
                      access = "open")


GDCdownload(query_GBM)
##### file was created on 2/13/2021  #####

e_gbm <- load(file = "exp.gbm.rda")
e_gbm
colData(data)
assays(data)


exp_Prep <- TCGAanalyze_Preprocessing(object = data, 
                                      cor.cut = 0.6,
                                      datatype = "HTSeq - FPKM-UQ")                      

exp_Filt <- TCGAanalyze_Filtering(tabDF = exp_Prep,
                                  method = "quantile", 
                                  qnt.cut =  0.25)   

t_exp_matrix <- t(exp_Filt)
dim(t_exp_matrix)

REST_position <- which(colnames(t_exp_matrix)  == "ENSG00000084093")

out_cor <- t(cor(t_exp_matrix[, REST_position], t_exp_matrix[, -REST_position]))

out_cor <- as.data.frame(out_cor)

library(data.table)
setDT(out_cor, keep.rownames = "Genes")

out_cor_ranked <- out_cor[order(out_cor$V1),]
head(out_cor_ranked, 100)
colnames(out_cor_ranked)[2] <- "Pearson r"

library(dplyr)

out_cor_ranked %>%
  mutate(p_value = NA) %>%
  mutate(q_value = NA)

for (i in 1:42335) {
  gene_correl <- which(colnames(t_exp_matrix)  == out_cor_ranked$Genes[i])
  out_cor_ranked$p_value[i] = cor.test(t_exp_matrix[, REST_position], t_exp_matrix[, gene_correl])$p.value
}


library(WGCNA)
qqs <- qvalue(out_cor_ranked$p_value)
out_cor_ranked$q_value <- qqs$qvalues
write.csv(as.data.frame(out_cor_ranked), "Output_TCGA_ALL correlations with REST.csv")