TCGAbiolinks 是一个用于TCGA数据综合分析的R/BioConductor软件包，能够通过GDC Application Programming Interface (API)访问 National Cancer Institute (NCI) Genomic Data Commons (GDC) ，来搜索、下载和准备TCGA相关数据，以便在R中进行分析。
0、镜像和R包安装和加载
#设置镜像
rm(list=ls())
options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
#意思是，如果没有安装BiocManager包，那就安装BiocManager
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
BiocManager::install("SummarizedExperiment")
BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
library(SummarizedExperiment)
1、TCGA数据下载
library(TCGAbiolinks)#加载包
query <- GDCquery(project = "TCGA-STAD", #肿瘤类型，##TCGAbiolinks:::getGDCprojects()$project_id#可查看
data.category = "Transcriptome Profiling",#数据范畴#转录组数据
data.type = "Gene Expression Quantification",#选定要下载的数据类型#gene表达数据
workflow.type = "HTSeq - Counts")#选定要下载RNAseq的Counts文件
#######################################################################################################
#下载rna-seq的counts数据
#data.type = "Gene Expression Quantification"
#下载miRNA数据
#data.type = "miRNA Expression Quantification"
#下载Copy Number Variation数据
#data.type = "Copy Number Segment"
#######################################################################################################
samplesDown <- getResults(query,cols=c("cases"))
#getResults(query, rows, cols)根据指定行名或列名从query中获取结果,此处用来获得样本的barcode，即，样本名
# 此处共检索出407个barcodes
head(samplesDown)#407样本名
[1] "TCGA-F1-A448-01A-11R-A24K-31" "TCGA-D7-6522-01A-11R-1802-13"
[3] "TCGA-VQ-A94R-01A-11R-A414-31" "TCGA-RD-A8MV-01A-11R-A36D-31"
[5] "TCGA-BR-4191-01A-02R-1131-13" "TCGA-BR-8291-01A-11R-2343-13"
#######################################################################################################
# TCGAquery_SampleTypes(barcode, typesample) #筛选肿瘤和正常组织的barcode，方便后面分组做差异分析
# TP代表PRIMARY SOLID TUMOR；NT-代表Solid Tissue Normal（其他组织样本可参考学习文档）
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,
typesample = "TP")
#View(dataSmTP )# 从samplesDown中筛选出375个TP样本barcodes
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,
typesample = "NT")
#View(dataSmNT)# 从samplesDown中筛选出32个正常样本barcodes
#######################################################################################################
###设置barcodes参数，筛选符合要求的375个肿瘤样本数据和32正常组织数据
#barcode参数：根据传入barcode进行数据过滤
queryDown <- GDCquery(project = "TCGA-SATD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts",
barcode = c(dataSmTP, dataSmNT))
#######################################################################################################
GDCdownload(queryDown, #上面queryDown的结果
method = "api",#两种方法，api或者gdc-client
files.per.chunk = 6)#分为多个片段下载,可以解决下载容易中断的问题。
##下载完成，前面如果已经下载了，这里不会再次下载
2、数据整理
########################################################################################
library(SummarizedExperiment)
#GDCprepare()将前面GDCquery（）的结果准备成R语言可处理的SE（SummarizedExperiment）文件
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename =
"STAD_case.Rda")
count_matrix=assay(dataPrep1)#数组
write.csv(count_matrix,file = paste("TCGA-STAD","Counts.csv",sep = "-"))#文件保存
###########################################################################################
# 去除dataPrep1中的异常值，dataPrep1数据中含有肿瘤组织和正常组织的数据
##TCGAanalyze_Preprocessing()对数据进行预处理：
#使用spearman相关系数去除数据中的异常值
#TCGAanalyze_Preprocessing(object, cor.cut = 0, filename = NULL,width = 1000, height = 1000, datatype = names(assays(object))[1])
# 函数功能描述：
#Array Array Intensity correlation (AAIC) and correlation boxplot to define outlier
dataPrep2 <- TCGAanalyze_Preprocessing(object = dataPrep1,
cor.cut = 0.6,
datatype = "HTSeq - Counts")
#将预处理后的数据dataPrep2，写入新文件“LIHC_dataPrep.csv”
write.csv(dataPrep2,file = "STAD_dataPrep.csv",quote = FALSE)
###########################################################################################
# TCGAtumor_purity(barcodes, estimate, absolute, lump, ihc, cpe);
#使用来自5种方法的5个估计值作为阈值对TCGA样本进行过滤，这5个值是estimate, absolute, lump, ihc, cpe;
#这里设置cpe=0.6（cpe是派生的共识度量，是将所有方法的标准含量归一化后的均值纯度水平，以使它们具有相等的均值和标准差）
#筛选肿瘤纯度大于等于60%的样本数据
purityDATA <- TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)
# filtered 为被过滤的数据， pure_barcodes是我们要的肿瘤数据
Purity.STAD<-purityDATA$pure_barcodes
normal.STAD<-purityDATA$filtered
#获取肿瘤纯度大于60%的375个肿瘤组织样本+32个正常组织样本,共计407个样本
puried_data <-dataPrep2[,c(Purity.STAD,normal.STAD)]
###########################################################################################
#基因注释,需要加载“SummarizedExperiment”包，“SummarizedExperiment container”每个由数字或其他模式的类似矩阵的对象表示。行通常表示感兴趣的基因组范围和列代表样品。
library("SummarizedExperiment")
rowData(dataPrep1) #传入数据dataPrep1必须为
# SummarizedExpensembl_gene_id external_gene_name original_ensembl_gene_id
# <character> <character>
#ENSG00000000003 ENSG00000000003 TSPAN6 ENSG00000000003.13
#ENSG00000000005 ENSG00000000005 TNMD ENSG00000000005.5
#ENSG00000000419 ENSG00000000419 DPM1 ENSG00000000419.11
#ENSG00000000457 ENSG00000000457 SCYL3 ENSG00000000457.12
#ENSG00000000460 ENSG00000000460 C1orf112 ENSG00000000460.15
#将结果写入文件“puried.STAD.cancer.csv”
rownames(puried_data)<-rowData(dataPrep1)$external_gene_name
write.csv(puried_data,file = "puried.STAD.csv",quote = FALSE)
###########################################################################################
#进行表达矩阵标准化和过滤，得到用于差异分析的表达矩阵
#TCGAanalyze_Normalization（）# 使用EDASeq软件包标准化mRNA转录本和miRNA。
#TCGAanalyze_Normalization()执行EDASeq包中的如下功能：
#1. EDASeq::newSeqExpressionSet
#2. EDASeq::withinLaneNormalization
#3. EDASeq::betweenLaneNormalization
#4. EDASeq::counts
dataNorm <- TCGAanalyze_Normalization(tabDF = puried_data,#RNAseq表达矩阵，行代表基因，列代表样本|
geneInfo = geneInfo,
method = "gcContent")
#|geneInfo|关于geneLength和gcContent的20531个基因的矩阵，“geneInfoHT”和“geneInfo”可选。
#method |选择标准化的方法，基于’gcContent’ 或 ’geneLength’的标准化方法可选
#将标准化后的数据再过滤，去除掉表达量较低（count较低）的基因，得到最终的数据
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,
method = "quantile", #用于过滤较低count数的基因的方法，有’quantile’, ’varFilter’, ’filter1’, ’filter2’
qnt.cut = 0.25)#过滤的阈值
str(dataFilt)
write.csv(dataFilt,file = "TCGA_STAD_final.csv",quote = FALSE)
3、基因差异表达分析
#首先读入表达矩阵文件
dataFilt_LIHC_final <- read.csv("TCGA_STAD_final.csv", header = T,check.names = FALSE）
# 定义行名
rownames(dataFilt_LIHC_final) <- dataFilt_LIHC_final[,1]
dataFilt_LIHC_final <- dataFilt_LIHC_final[,-1]
View(dataFilt_STAD_final)
#定义样本的分组，肿瘤组和正常组
#目前只知道根据barcode14和15位数字0—9为肿瘤，大于就为正常样本
metadata <- data.frame(colnames(dataFilt_STAD_final))#View(metadata)
for (i in 1:length(metadata[,1])) {
num <- as.numeric(substring(metadata[i,1],14,15))
if (num %in% seq(1,9)) {metadata[i,2] <- "Tum"}
if (num %in% seq(10,29)) {metadata[i,2] <- "Nor"}
}
colnames(metadata)[2] <- c("fenzu")
ifelse(metadata$fenzu%in%c("Tum")==TRUE,"Tum","Nor")#发现从376行开始是正常的
# 定义肿瘤样本分组
mat1 <- dataFilt_STAD_final[,1-375]
mat1 <- log(mat1+1)
# 定义正常组织样本分组
mat2 <- dataFilt_STAD_final[,376-407]
mat2 <- log(mat2+1)
#数据准备好了，开始差异表达分析
Data_DEGs <- TCGAanalyze_DEA(mat1 = mat1,#肿瘤组织的表达矩阵
mat2 = mat2,#正常样本的表达矩阵
Cond1type = "Tumor",#mat1中的样品分组信息
Cond2type = "Normal",#mat2中的样品分组信息
pipeline="limma",#用limma包还是edgeR
batch.factors = c("TSS"),#批处理纠正选项
voom = TRUE,
contrast.formula = "Mycontrast=Tumor-Normal")
View(Data_DEGs)#结果展示
######################################################################################
#设置logFC，挑选表达有差异的基因进行富集分析
Data_DEGs_high_expr <- Data_DEGs[Data_DEGs$logFC >=1,]
Genelist <- rownames(Data_DEGs_high_expr)
ansEA <- TCGAanalyze_EAcomplete(TFname="DEA genes Normal Vs Tumor", Genelist)
#富集分析结果可视化##
TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
GOBPTab = ansEA$ResBP,
GOCCTab = ansEA$ResCC,
GOMFTab = ansEA$ResMF,
PathTab = ansEA$ResPat,
nRGTab = Genelist,
nBar = 10, #显示条形图的数量
filename = "TCGAvisualize_EAbarplot_Output.pdf")
#结果没富集出来提示
#Error in mtext(mainLab, side = 3, line = -1, outer = TRUE, font = 2) : plot.new has not been called yet