主要想介绍的是下载下来数据的数据处理这一块，不是特别建议RTCGA的下载方法，列出了是因为怕没有源文件后续数据处理看不懂
一下载数据
RTCGA下载是把所有数据一起下载下来，存在不是最新的问题，此为2015年的
# Load the bioconductor installer.
source("https://bioconductor.org/biocLite.R")
# Install the main RTCGA package
biocLite("RTCGA")
# Install the clinical and mRNA gene expression data packages
biocLite("RTCGA.clinical") ## 14Mb
biocLite('RTCGA.rnaseq') ## (612.6 MB)
biocLite("RTCGA.mRNA") ## (85.0 MB)
biocLite('RTCGA.mutations') ## (103.8 MB)
library(RTCGA)
## Welcome to the RTCGA (version: 1.8.0).
all_TCGA_cancers=infoTCGA()#查看所有肿瘤类型每种数据分别有多少
DT::datatable(all_TCGA_cancers)
library(RTCGA.clinical)
library(RTCGA.mRNA)
## ?mRNA
## ?clinical
二提取数据中的表达矩阵
#提取表达矩阵，已经写好的函数，直接可以得到表达矩阵
expr <- expressionsTCGA(BRCA.mRNA, OV.mRNA, LUSC.mRNA,
extract.cols = c("GATA3", "PTEN", "XBP1","ESR1", "MUC1"))
## Warning in flatten_bindable(dots_values(...)): '.Random.seed' is not an
三处理下载下来的数据
expr#通过输入名字查看表达矩阵的情况
nb_samples <- table(expr$dataset)#看表达矩阵中的每个类型的数据都有多少
nb_samples
expr$dataset <- gsub(pattern = ".mRNA", replacement = "", expr$dataset)
#把expr的dataset这一列的数字中的,mRNA变成空的，更换表达形式从BRCA.mRNA变成BRCA
expr$dataset
expr$bcr_patient_barcode <- paste0(expr$dataset, c(1:590, 1:561, 1:154))
#把病人样本名称的这列简化变成对应的肿瘤名称加序号，这个数字是根据原来的肿瘤数量算出来的
expr
四利用下载下来的数据画图，看组间差异
library(ggpubr)
## Loading required package: ggplot2
## Loading required package: magrittr
# GATA3
ggboxplot(expr, x = "dataset", y = "GATA3",
title = "GATA3", ylab = "Expression",
color = "dataset", palette = "jco")
Rplotsd.jpeg
加上P值
my_comparisons <- list(c("BRCA", "OV"), c("OV", "LUSC"))
ggboxplot(expr, x = "dataset", y = "GATA3",
title = "GATA3", ylab = "Expression",
color = "dataset", palette = "jco")+
stat_compare_means(comparisons = my_comparisons)
Rplotwe.jpeg
label.select.criteria <- list(criteria = "`y` > 3.9 & `x` %in% c('BRCA', 'OV')")
ggboxplot(expr, x = "dataset",
y = c("GATA3", "PTEN", "XBP1"),
combine = TRUE,
color = "dataset", palette = "jco",
ylab = "Expression",
label = "bcr_patient_barcode", # column containing point labels
label.select = label.select.criteria, # Select some labels to display
font.label = list(size = 9, face = "italic"), # label font
repel = TRUE # Avoid label text overplotting
)
Rplot12.jpeg
小总结：处理数据最主要的是通过查看数据（输它的名字），了解数据的构成，然后根据你的需要个性化的处理数据，想实现什么都可以自行百度
方法二：
rm(list=ls())
options(stringsAsFactors = F)
# 注意，并不是说使用 RTCGA.miRNASeq包的数据是最佳选择，只是因为这个演示起来最方便。
# 因为GDC官网下载数据具有一定门槛，也不是每个人都必须学会的。
getwd()
Rdata_dir='Rdata/'#设定获取路径
Figure_dir='figures/'#设定获取路径
如果开启下面代码，就会从RTCGA.miRNASeq包里面提取miRNA表达矩阵和对应的样本临床信息。
library(RTCGA.miRNASeq) #这个RTCGA.miRNASeq包里面包含TCGA数据库里面所有癌症的miRNAseq的信息
s=rownames(KIRC.miRNASeq)[seq(1,nrow(KIRC.miRNASeq),by=3)]
#rownames(KIRC.miRNASeq)：查看KIRC（肾透明细胞癌）的miRNASeq的行名，看下图可以看出有三种数据
#seq(1,nrow(KIRC.miRNASeq),by=3)：从第一行起每三行取一个，取出其中一种数据的所有行名（也就是样本名）
微信截图_20190714201632.png
此处知识点：TCGA的转录组数据，目前支持3中类型下载，分别是Count, FPKM和FPKM-UQ
转录组数据，都是采用HTseq-count 进行定量的，之后再转换成FPKM 和 FPKM-UQ
一般推荐采用Count 值，因为：
- 下游进行差异分析的软件，比如DESeq2, edgeR都是采用Count值，2. Count值，也可以转换成FPKM和FPKM-UQ
expr <- expressionsTCGA(KIRC.miRNASeq)#提取表达矩阵
dim(expr)
expr[1:40,1:4]
expr
expr=as.data.frame(expr[seq(1,nrow(expr),by=3),3:ncol(expr)])##expr[seq(1,nrow(expr),by=3)从第一行起每三行取一个，取出其中一种数据的所有表达矩阵（也就是样本名）
##查看样本可以发现，只有3到最后才是表达矩阵信息，所以列只取了3:ncol(expr)
得到表达矩阵expr
mi=colnames(expr)#获取表达矩阵的探针名为mi
expr=apply(expr,1,as.numeric)
#1表示行，2表示列，c(1,2)表示行和列，此处针对行进行操作，把每行变成一列数字数据
变化前
变化后，可以看到只取了数据，且行列变化了
colnames(expr)=s#把样本名对应列名
rownames(expr)=mi#把探针名对应行名
expr[1:4,1:4]
expr=na.omit(expr)#删除NA值
得到表达矩阵.jpg
expr=expr[apply(expr, 1,function(x){sum(x>1)>10}),]
#对每一行进行function运算，每行中大于1的数字之和大于10，则保留此行（去除一些只测了少数几个样本的行）
library(RTCGA.clinical)
meta <- KIRC.clinical#得到临床信息
tmp=as.data.frame(colnames(meta))#得到列名的矩阵
meta[(grepl('patient.bcr_patient_barcode',colnames(meta)))]#得到样本名
meta[(grepl('patient.days_to_last_followup',colnames(meta)))]#得到生存时间
meta[(grepl('patient.days_to_death',colnames(meta)))]
#得到生存时间（死亡病人的）
meta[(grepl('patient.vital_status',colnames(meta)))]
#得到生存状态
## patient.race # patient.age_at_initial_pathologic_diagnosis # patient.gender
# patient.stage_event.clinical_stage
meta=as.data.frame(meta[c('patient.bcr_patient_barcode','patient.vital_status',
'patient.days_to_death','patient.days_to_last_followup',
'patient.race',
'patient.age_at_initial_pathologic_diagnosis',
'patient.gender' ,
'patient.stage_event.pathologic_stage')])
#去除这些需要的数据变成新的数据框
#meta[(grepl('patient.stage_event.pathologic_stage',colnames(meta)))]
## 每次运行代码，就会重新生成文件。
save(expr,meta,
file = file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata') )#按照之前设定的路径保存起来
提取出来的临床信息
## 我们已经运行了上面被关闭的代码，而且保存了miRNA表达矩阵和对应的样本临床信息
# 现在直接加载即可。
load( file =
file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537个病人，但是有593个样本，每个样本有 552个miRNA信息。
# 当然，这个数据集可以下载原始测序数据进行重新比对，可以拿到更多的miRNA信息
