三阴性乳腺癌提取和分析

##########################################################################################
## 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\\F22\\brca')


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


query <- GDCquery(project = "TCGA-BRCA", 
                  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())

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 <-  count_matrix %>% 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\\F22\\brca')

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


##########################################################################################
## 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\\F22\\brca')


query <- GDCquery(project = "TCGA-BRCA", 
                  data.category = "Clinical", 
                  file.type = "xml")
GDCdownload(query)

clinical <- GDCprepare_clinic(query, clinical.info = "patient")


#合并两类患者，得到肾透明细胞癌的临床信息
survival_data <- as.data.frame(clinical)

write.csv(survival_data , file = 'survival.csv')




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

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


options(stringsAsFactors = F)

library(stringr)

rm(list=ls())

setwd('D:\\SCIwork\\F22\\brca')

load( "mRNA_exprSet.Rda")


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


for (i in 1:length(metadata[,1])) {
  num <- as.numeric(as.character(substring(metadata[i,1],14,15)))
  if (num == 1 ) {metadata[i,2] <- "T"}
  if (num != 1) {metadata[i,2] <- "N"}
}
names(metadata) <- c("id","group")
metadata$group <- as.factor(metadata$group)



metadata_T <- subset(metadata,metadata$group == "T")

metadata_T

metadata_N <- metadata[which(!(metadata$id %in% metadata_T$id)),]

metadata_N$group <- 'Control'

metadata_T$id1 <- substr(x= metadata_T$id, start = 1, stop = 12)



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

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


p <- read.csv('survival.csv',  header = T, sep = ',')

colnames(p)[grep("receptor_status", colnames(p))]
## [1] "breast_carcinoma_estrogen_receptor_status"               
## [2] "breast_carcinoma_progesterone_receptor_status"           
## [3] "lab_proc_her2_neu_immunohistochemistry_receptor_status"  
## [4] "metastatic_breast_carcinoma_estrogen_receptor_status"    
## [5] "metastatic_breast_carcinoma_progesterone_receptor_status"

# examining how many triple-negative receptors samples 
table(p$breast_carcinoma_estrogen_receptor_status == 'Negative' &
        p$breast_carcinoma_progesterone_receptor_status == 'Negative' &
        p$lab_proc_her2_neu_immunohistochemistry_receptor_status == 'Negative')

# extracting tnbc samples 
tnbc_samples <- p[p$breast_carcinoma_estrogen_receptor_status == 'Negative' &
                    p$breast_carcinoma_progesterone_receptor_status == 'Negative' &
                    p$lab_proc_her2_neu_immunohistochemistry_receptor_status == 'Negative', ]

tnbc_samples <- tnbc_samples$bcr_patient_barcode

metadata_TN <- metadata_T[which(metadata_T$id1 %in% tnbc_samples),]

metadata_TN$id1 <- NULL

metadata_TN$group <- 'TNBC'

metadata_TT <- metadata_T[which(!(metadata_T$id1 %in% tnbc_samples)),]

metadata_TT$id1 <- NULL

metadata_TT$group <- 'non-TNBC'