cgdsr
library(cgdsr)
library(DT)
mycgds = CGDS("http://www.cbioportal.org/")
test(mycgds)
### all cancer studies
all_tcga_studies <- getCancerStudies(mycgds)
datatable(all_tcga_studies[,c(1,3)])
### get the lggGBM genetic profiles
lgg_genetic <- getGeneticProfiles(mycgds,'lgg_tcga_pan_can_atlas_2018')[,c(1:3)]
datatable(lgg_genetic)
mygenetic <- "lgg_tcga_pan_can_atlas_2018_rna_seq_v2_mrna"
### get the case of lgg
lgg_case <- getCaseLists(mycgds,'lgg_tcga')[,c(1:3)]
datatable(lgg_case)
mycase <- "lgg_tcga_all"
df = getProfileData(mycgds,"TP53",mygenetic,mycase)
df <- na.omit(df)
cldt <- getClinicalData(mycgds,"lgg_tcga_all")
group <- vector(length = 419L)
df_all_data$PERSON_NEOPLASM_CANCER_STATUS <- revalue(factor(df_all_data$PERSON_NEOPLASM_CANCER_STATUS),c(With Tumor = 1,Tumor Free = 0))
gene_name <- "ASCL1"
time_status <- select(df_all_data,c(DAYS_LAST_FOLLOWUP,PERSON_NEOPLASM_CANCER_STATUS,gene_name))
time_status$group <- sapply(time_status[gene_name], function(x){ifelse(x-median(x) <0,"low","high")})
time_status <- time_status[time_status$PERSON_NEOPLASM_CANCER_STATUS !="",]
time_status$PERSON_NEOPLASM_CANCER_STATUS <-revalue(time_status$PERSON_NEOPLASM_CANCER_STATUS,c("Tumor Free" =0,"With Tumor" =1))
time_status$PERSON_NEOPLASM_CANCER_STATUS <- time_status$PERSON_NEOPLASM_CANCER_STATUS %>% droplevels()
time_status$PERSON_NEOPLASM_CANCER_STATUS <- as.integer(time_status$PERSON_NEOPLASM_CANCER_STATUS)
mysurvfit <- survfit(Surv(time_status$DAYS_LAST_FOLLOWUP, time_status$PERSON_NEOPLASM_CANCER_STATUS) ~time_status$group,data = time_status)
mysurvfit2 <- survfit(Surv(time_status$DAYS_LAST_FOLLOWUP, as.integer(time_status$PERSON_NEOPLASM_CANCER_STATUS)) ~time_status$group,data = time_status)
plot(mysurvfit2)
ggsurvplot(mysurvfit2,pval = T)
###### CONSTRUCT A LOOP
everyfit <- function(gene_name = "ASCL1"){
time_status <- select(df_all_data,c(DAYS_LAST_FOLLOWUP,PERSON_NEOPLASM_CANCER_STATUS,gene_name))
time_status$group <- sapply(time_status[gene_name], function(x){ifelse(x-median(x) <0,"low","high")})
time_status <- time_status[time_status$PERSON_NEOPLASM_CANCER_STATUS !="",]
time_status$PERSON_NEOPLASM_CANCER_STATUS <-revalue(time_status$PERSON_NEOPLASM_CANCER_STATUS,c("Tumor Free" =0,"With Tumor" =1))
time_status$PERSON_NEOPLASM_CANCER_STATUS <- time_status$PERSON_NEOPLASM_CANCER_STATUS %>% droplevels()
time_status$PERSON_NEOPLASM_CANCER_STATUS <- as.numeric(time_status$PERSON_NEOPLASM_CANCER_STATUS)
time_status$PERSON_NEOPLASM_CANCER_STATUS[time_status$PERSON_NEOPLASM_CANCER_STATUS ==1] <- 0
time_status$PERSON_NEOPLASM_CANCER_STATUS[time_status$PERSON_NEOPLASM_CANCER_STATUS ==2] <-1
thisfit <- survfit(Surv(time_status$DAYS_LAST_FOLLOWUP, time_status$PERSON_NEOPLASM_CANCER_STATUS) ~time_status$group,data = time_status)
invisible(thisfit)
}
everyfit("NOTCH1") %>% ggsurvplot(pval = T,conf.int=F)
everyfit("MAP2") %>% ggsurvplot(pval = T)
RTCGA
### load packages
library(RTCGA)
library(RTCGA.clinical)
library(RTCGA.rnaseq)
library(dplyr)
library(DT)
#????????????
infoTCGA <- infoTCGA() #??????????????????????????????????????????????????????????????????????????????
# Create the clinical data
#library(RTCGA.clinical)
clin <- survivalTCGA(GBMLGG.clinical)
clin[1:5,] ## status 0 = alive,1 = dead
##library(RTCGA.mRNA) #???????????????
class(GBMLGG.rnaseq) #???????????????????????????????????????
dim(GBMLGG.rnaseq) #??????????????????????????????696????????????20532?????????
GBMLGG.rnaseq[1:5,1:5] ## row for sample col for genes
###read my gene_List
list.files()
library(xlsx)
genelist <- read.xlsx("gene_List.xlsx",1,header = F)
genelist
myvector <- vector(mode = "list",length = length(genelist$X1))
for (i in seq_along(genelist$X2)) {
myvector[[i]] <- grep(genelist$X2[i],names(GBMLGG.rnaseq),value = T)
}
myvector[[2]] <- myvector[[2]][1]
myvector[[3]] <- myvector[[3]][2]
myvector[[5]] <- myvector[[5]][4]
myvector[[7]] <- myvector[[7]][2]
myvector[[8]] <- myvector[[8]][3]
myvector <- unlist(myvector)
#### select
exprSet <- GBMLGG.rnaseq %>% as_tibble() %>%
select("bcr_patient_barcode",myvector) %>%
mutate(bcr_patient_barcode = substr(bcr_patient_barcode,1,12)) %>%
inner_join(clin,by="bcr_patient_barcode")
colnames(exprSet)[2:9] <- sub("\\|[0-9]+","",colnames(exprSet)[2:9])
## survial
library(survival)
library(survminer)
my.surv <- Surv(exprSet$times, exprSet$patient.vital_status)## fist get the surv object
log_rank_p <- apply(exprSet[,2:9], 2, function(value1){
group <- ifelse(value1 > median(value1),"high","low")
kmfit <- survfit(my.surv ~ group,data = exprSet)
data.survdiff <- survdiff(my.surv ~ group)
p.value <- 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
})
log_rank_p < 0.05
names(exprSet)
### select one of the genes in the genelist and plot survival curve
names(exprSet[2:9])
gene = "NOTCH1"
plot_gene <- function(genename = gene){
group <- sapply(exprSet[,match(genename,names(exprSet))][1],function(x){ifelse(x > median(exprSet[[genename]]),"high","low")})
exprSet$group <<- group
kmfit <- survfit(my.surv ~group,data = exprSet)
invisible(kmfit)
}
plot_gene() %>% ggsurvplot(conf.int=FALSE, pval=TRUE ) + ggtitle(gene)
#### save files
tmp <- log_rank_p %>% data.frame(pvalue = .)
write.csv(tmp,"survival_p_value.csv",quote = T,sep = ",")
