##########################################################################################
## step1 load package and change Working Directory
###########################################################################################

library(TCGAbiolinks)
library(dplyr)
library(tidyr)
library(tibble)
library(edgeR)
library(limma)
rm(list=ls())
setwd('D:\\SCIwork\\F20ELFN1\\COAD')


##########################################################################################
## step2 download the expresssion data of lncRNA and mRNA
###########################################################################################


query <- GDCquery(project = "TCGA-COAD", 
                  data.category = "Transcriptome Profiling", 
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")



GDCdownload(query, method = "api", files.per.chunk = 50)

library(SummarizedExperiment)

expdat <- GDCprepare(query = query,   save = TRUE, save.filename = "exp.rda")


count_matrix = as.data.frame(assay(expdat))

##########################################################################################

###########################################################################################



rm(list=ls())

setwd('D:\\SCIwork\\F20ELFN1\\COAD')

load('exp.rda')

count_matrix = as.data.frame(assay(data))

count_matrix[1:4,1:4]



fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}


expr <- as.data.frame (apply(count_matrix  , 2, fpkmToTpm))


expr <-  expr %>% rownames_to_column("gene_id")




##########################################################################################

###########################################################################################





setwd("D:\\Originaldata\\GRCH\\Homo_sapiens.GRCh38.90")

load("gtf_df.Rda")

test <- gtf_df[1:5,]

View(test)




mRNA_exprSet <- gtf_df %>%
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
  dplyr::inner_join(expr,by ="gene_id") %>%
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")


save(mRNA_exprSet,file = "mRNA_exprSet.Rda")


mRNA_exprSet <- mRNA_exprSet %>%
  tidyr::separate(gene_id, c("gene_name","gene_id","gene_biotype"),
                  sep = " \\| ")


mRNA_exprSet <- mRNA_exprSet[,-(2:3)]


index <- duplicated(mRNA_exprSet$gene_name)
mRNA_exprSet <- mRNA_exprSet[!index,]
row.names(mRNA_exprSet) <- mRNA_exprSet$gene_name
mRNA_exprSet$gene_name <- NULL


setwd('D:\\SCIwork\\F20ELFN1\\COAD')
save(mRNA_exprSet, file = "mRNA_exprSet.Rda")

##########################################################################################




###########################################################################################

rm(list=ls())


load( "mRNA_exprSet.Rda")



metadata <- data.frame(colnames(mRNA_exprSet ))

colnames(metadata) <- 'barcode'


for (i in 1:length(metadata[,1])) {
  num <- as.numeric(as.character(substring(metadata[i,'barcode'],14,15)))
  if (num == 1 ) {metadata[i,2] <- "T"}
  if (num != 1) {metadata[i,2] <- "N"}
}



names(metadata)[2] <- 'Barcode'

table(metadata$Barcode)

metadata <- subset(metadata, metadata$Barcode == 'T')


mRNA_exprSet <- mRNA_exprSet[,which(colnames(mRNA_exprSet) %in% metadata$barcode)]


setwd('D:\\SCIwork\\F20ELFN1\\COAD')

save(mRNA_exprSet, file = "mRNA_exprSet.Rda")