#使用easyTCGA獲取數據
#清空
rm(list=ls())
gc()
# 安裝bioconductor上面的R包
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
if(!require("BiocManager")) install.packages("BiocManager")
if(!require("TCGAbiolinks")) BiocManager::install("TCGAbiolinks")
if(!require("SummarizedExperiment")) BiocManager::install("SummarizedExperiment")
if(!require("DESeq2")) BiocManager::install("DESeq2")
if(!require("edgeR")) BiocManager::install("edgeR")
if(!require("limma")) BiocManager::install("limma")
# 安裝cran上面的R包
if(!require("survival")) install.packages("survival")
if(!require("broom")) install.packages("broom")
if(!require("devtools")) install.packages("devtools")
if(!require("cli")) install.packages("cli")
#devtools::install_github("ayueme/easyTCGA")
library(easyTCGA)
help(package="easyTCGA")
setwd("F:\\TCGA\\TCGA-COAD")
#下載mRNA、lncRNA和臨床信息
COAD<-getmrnaexpr("TCGA-COAD")#原始下載的count, TPM, FPKM 均沒有經過log2轉化
#下載miRNA
COAD_miRNA<-getmirnaexpr("TCGA-COAD")
#下載copy number variation data
COAD_cnv<-getcnv("TCGA-COAD")
#下載masked somatic mutation 體細胞突變
COAD_snv<-getsnvmaf("TCGA-COAD")
#下載DNA methylation beta value 甲基化數據
getmethybeta("TCGA-COAD")
?
#從下載目錄中打開數據
#差異分析
diff<-diff_analysis(exprset=mrna_expr_counts,#沒有經過log2轉化project="TCGA-COAD",save=F)#批量生存分析
surv<-batch_survival(exprset=mrna_expr_counts,clin=clin_info,is_count = T,optimal_cut = TRUE,project="TCGA-COAD",save_data = FALSE,min_sample_size = 5,print_index = TRUE
)?
#突變分析：瀑布圖
#BiocManager::install("maftools")
library(maftools)
maf<-read.maf(snv,clinicalData=clin_snv)
plotmafSummary(maf)
colnames(clin_snv)
oncoplot(maf=maf,clinicalFeatures=c("ajcc_pathologic_stage","vital_status"),top=10,sortByAnnotation=T
)
?
?
#繪制KM曲線
dim(mrna_expr_counts)
set.seed(123)
colnames(clin_info)
clin<-data.frame(time=clin_info$days_to_last_follow_up,event=clin_info$vital_status)
clin$event<-ifelse(clin$event=="Alive",0,1)
plot_KM(exprset=mrna_expr_counts, marker="CHPF", #基因clin=clin, optimal_cut = TRUE, return_data = TRUE)
?
#正常和癌癥組織基因表達對比箱線圖
rownames(mrna_expr_counts)
plot_gene_paired(exprset=mrna_expr_counts, marker="CHPF", #基因return_data = TRUE)
?
#比較組間基因表達差異
set.seed(123)
group=sample(c(0,1),524,replace = T)
plot_gene(exprset=mrna_expr_counts, marker=c("CHPF","MAOA"), group=group, return_data = TRUE)
?
