主要想介紹的是下載下來數(shù)據(jù)的數(shù)據(jù)處理這一塊际看，不是特別建議RTCGA的下載方法，列出了是因為怕沒有源文件后續(xù)數(shù)據(jù)處理看不懂
一下載數(shù)據(jù)
RTCGA下載是把所有數(shù)據(jù)一起下載下來咬崔，存在不是最新的問題，此為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()#查看所有腫瘤類型每種數(shù)據(jù)分別有多少
DT::datatable(all_TCGA_cancers)
library(RTCGA.clinical)
library(RTCGA.mRNA)
## ?mRNA
## ?clinical
二提取數(shù)據(jù)中的表達矩陣
#提取表達矩陣，已經(jīng)寫好的函數(shù)伏嗜，直接可以得到表達矩陣
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
三處理下載下來的數(shù)據(jù)
expr#通過輸入名字查看表達矩陣的情況
nb_samples <- table(expr$dataset)#看表達矩陣中的每個類型的數(shù)據(jù)都有多少
nb_samples
expr$dataset <- gsub(pattern = ".mRNA", replacement = "", expr$dataset)
#把expr的dataset這一列的數(shù)字中的,mRNA變成空的，更換表達形式從BRCA.mRNA變成BRCA
expr$dataset
expr$bcr_patient_barcode <- paste0(expr$dataset, c(1:590, 1:561, 1:154))
#把病人樣本名稱的這列簡化變成對應(yīng)的腫瘤名稱加序號伐厌，這個數(shù)字是根據(jù)原來的腫瘤數(shù)量算出來的
expr
四利用下載下來的數(shù)據(jù)畫圖阅仔，看組間差異
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
小總結(jié)：處理數(shù)據(jù)最主要的是通過查看數(shù)據(jù)（輸它的名字），了解數(shù)據(jù)的構(gòu)成弧械，然后根據(jù)你的需要個性化的處理數(shù)據(jù)八酒，想實現(xiàn)什么都可以自行百度
方法二：
rm(list=ls())
options(stringsAsFactors = F)
# 注意，并不是說使用 RTCGA.miRNASeq包的數(shù)據(jù)是最佳選擇刃唐，只是因為這個演示起來最方便羞迷。
# 因為GDC官網(wǎng)下載數(shù)據(jù)具有一定門檻，也不是每個人都必須學會的画饥。
getwd()
Rdata_dir='Rdata/'#設(shè)定獲取路徑
Figure_dir='figures/'#設(shè)定獲取路徑
如果開啟下面代碼衔瓮，就會從RTCGA.miRNASeq包里面提取miRNA表達矩陣和對應(yīng)的樣本臨床信息。
library(RTCGA.miRNASeq) #這個RTCGA.miRNASeq包里面包含TCGA數(shù)據(jù)庫里面所有癌癥的miRNAseq的信息
s=rownames(KIRC.miRNASeq)[seq(1,nrow(KIRC.miRNASeq),by=3)]
#rownames(KIRC.miRNASeq)：查看KIRC（腎透明細胞癌）的miRNASeq的行名抖甘，看下圖可以看出有三種數(shù)據(jù)
#seq(1,nrow(KIRC.miRNASeq),by=3)：從第一行起每三行取一個热鞍，取出其中一種數(shù)據(jù)的所有行名（也就是樣本名）
微信截圖_20190714201632.png
此處知識點：TCGA的轉(zhuǎn)錄組數(shù)據(jù)，目前支持3中類型下載，分別是Count, FPKM和FPKM-UQ
轉(zhuǎn)錄組數(shù)據(jù)薇宠，都是采用HTseq-count 進行定量的偷办，之后再轉(zhuǎn)換成FPKM 和 FPKM-UQ
一般推薦采用Count 值，因為：
- 下游進行差異分析的軟件澄港，比如DESeq2, edgeR都是采用Count值椒涯，2. Count值，也可以轉(zhuǎn)換成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)從第一行起每三行取一個回梧，取出其中一種數(shù)據(jù)的所有表達矩陣（也就是樣本名）
##查看樣本可以發(fā)現(xiàn)废岂，只有3到最后才是表達矩陣信息，所以列只取了3:ncol(expr)
得到表達矩陣expr
mi=colnames(expr)#獲取表達矩陣的探針名為mi
expr=apply(expr,1,as.numeric)
#1表示行狱意，2表示列湖苞，c(1,2)表示行和列，此處針對行進行操作详囤，把每行變成一列數(shù)字數(shù)據(jù)
變化前
變化后袒啼，可以看到只取了數(shù)據(jù)，且行列變化了
colnames(expr)=s#把樣本名對應(yīng)列名
rownames(expr)=mi#把探針名對應(yīng)行名
expr[1:4,1:4]
expr=na.omit(expr)#刪除NA值
得到表達矩陣.jpg
expr=expr[apply(expr, 1,function(x){sum(x>1)>10}),]
#對每一行進行function運算纬纪，每行中大于1的數(shù)字之和大于10蚓再，則保留此行（去除一些只測了少數(shù)幾個樣本的行）
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)))]
#得到生存狀態(tài)
## 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')])
#去除這些需要的數(shù)據(jù)變成新的數(shù)據(jù)框
#meta[(grepl('patient.stage_event.pathologic_stage',colnames(meta)))]
## 每次運行代碼，就會重新生成文件包各。
save(expr,meta,
file = file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata') )#按照之前設(shè)定的路徑保存起來
提取出來的臨床信息
## 我們已經(jīng)運行了上面被關(guān)閉的代碼摘仅，而且保存了miRNA表達矩陣和對應(yīng)的樣本臨床信息
# 現(xiàn)在直接加載即可。
load( file =
file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537個病人问畅，但是有593個樣本娃属，每個樣本有 552個miRNA信息。
# 當然护姆，這個數(shù)據(jù)集可以下載原始測序數(shù)據(jù)進行重新比對矾端，可以拿到更多的miRNA信息
