TCGAbiolinks全解析
2020-04-14 本文已影响0人
养猪场小老板
部分一
第一步、数据下载阶段
####(1.1)、加载包
rm(list=ls())
library(TCGAbiolinks)
library(plyr)
library(limma)
library(biomaRt)
library(SummarizedExperiment)
####(1.2)
####GDCquery()筛选需要的数据，TCGAbiolinks包下载TCGA数据进行表达差异分析
query <- GDCquery(project = "TCGA-STAD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
####(1.3)
####getResults(query, rows, cols)根据指定行名或列名从query中获取结果,此处用来获得样本的barcode
#### 此处共检索出407个barcodes
samplesDown <- getResults(query,cols=c("cases"))
# 从samplesDown中筛选出TP（实体肿瘤）样本的barcodes
# TCGAquery_SampleTypes(barcode, typesample)
# TP代表PRIMARY SOLID TUMOR；
# NT-代表Solid Tissue Normal(其他组织样本可参考学习文档)
##此处共检索出375个TP样本barcodes
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,
typesample = "TP")
####从samplesDown中筛选出NT(正常组织)样本的barcode
####32个NT样本barcodes
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,
typesample = "NT")
####(1.4)
###设置barcodes参数，筛选符合要求的375个肿瘤样本数据和32正常组织数据
queryDown <- GDCquery(project = "TCGA-STAD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts",
barcode = c(dataSmTP, dataSmNT))
#barcode参数：根据传入barcodes进行数据过滤
第二步：GDCdownload()下载GDCquery()得到的结果
####2.1 下载数据，默认存放位置为当前工作目录下的GDCdata文件夹中。
GDCdownload(queryDown,
method = "api",
directory = "GDCdata",
files.per.chunk = 6)
#method ；"API"或者"client"。"API"速度更快，但是容易下载中断。
#directory：下载文件的保存地址。Default: GDCdata。
#files.per.chunk = NULL:使用API下载大文件的时候，可以把文件分成几个小文件来下载，可以解决下载容易中断的问题。
第三步、数据处理
####3.1：GDCprepare()将前面GDCquery()的结果准备成R语言可处理的SE（SummarizedExperiment）文件
####读取下载的数据并将其准备到R对象中，在工作目录生成（save=TRUE）STAD_case.rda文件
#### GDCprepare():Prepare GDC data,准备GDC数据，使其可用于R语言中进行分析
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename =
"STAD_case.rda")
第四步：TCGAanalyze_Preprocessing()对数据进行预处理：
####使用spearman相关系数去除数据中的异常值
# 去除dataPrep1中的异常值，dataPrep1数据中含有肿瘤组织和正常组织的数据
# 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，写入新文件“STAD_dataPrep.csv”
write.csv(dataPrep2,file = "STAD_dataPrep.csv",quote = FALSE)
第五步：TCGAtumor_purity()筛选肿瘤纯度大于60%的肿瘤barcodes
# 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)
####这里筛选的60%的样本数据并没有用成功，这里直接赋值
# filtered 为被过滤的数据， pure_barcodes是我们要的肿瘤数据
####过滤的数据太多，重新分配
purityDATA$pure_barcodes <- dataSmTP
purityDATA$filtered <- dataSmNT
####
Purity.STAD <- purityDATA$pure_barcodes
normal.STAD <- purityDATA$filtered
第六步：将肿瘤表达矩阵与正常组织表达矩阵合并，进行基因注释
#获取375个肿瘤组织样本+32个正常组织样本,共计407个样本
puried_data <-dataPrep2[,c(Purity.STAD,normal.STAD)]
第七步：进行表达矩阵基因注释
#基因注释,需要加载“SummarizedExperiment”包，
#“SummarizedExperiment container”每个由数字或其他模式的类似矩阵的对象表示。
# 行通常表示感兴趣的基因组范围和列代表样品。
library("SummarizedExperiment")
#将结果写入文件“puried.LIHC.cancer.csv”
rownames(puried_data)<-rowData(dataPrep1)$external_gene_name
write.csv(puried_data,file = "puried.STAD.csv",quote = FALSE)
第八步：进行表达矩阵标准化和过滤，得到用于差异分析的表达矩阵
####1、`TCGAanalyze_Normalization()`使用EDASeq软件包标准化mRNA转录本和miRNA。
dataNorm <- TCGAanalyze_Normalization(tabDF = puried_data,
geneInfo = geneInfo,
method = "gcContent")
####2、将标准化后的数据再过滤，去除掉表达量较低(count较低)的基因，得到最终的数据
dataFilt_STAD_final <- TCGAanalyze_Filtering(tabDF = dataNorm,
method = "quantile",
qnt.cut = 0.25)
第九步、差异基因分析
##定义样本的分组，前375个是肿瘤样本，后32个是正常样本
# 定义肿瘤样本分组
mat1 <- dataFilt_STAD_final[,dataSmTP]
mat1 <- log(mat1+1)
# 定义正常组织样本分组
mat2 <- dataFilt_STAD_final[,dataSmNT]
mat2 <- log(mat2+1)
#### 差异分析
Data_DEGs <- TCGAanalyze_DEA(mat1 = mat1,
mat2 = mat2,
Cond1type = "Tumor",
Cond2type = "Normal",
pipeline="limma",
batch.factors = c("TSS"),
voom = TRUE,
contrast.formula = "Mycontrast=Tumor-Normal")
#View(Data_DEGs)
第十步(1)，富集分析
####设置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)
####第十步(2)，富集分析可视化
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")
第十一，聚类分析
res.hc <- TCGAanalyze_Clustering(Data_DEGs,
method = "hclust",
methodHC = "ward.D2")
#########################################################################################
#########################################################################################
#########################################################################################
部分二
1、 下载临床数据
clin.STAD <- GDCquery_clinic("TCGA-STAD", "clinical",save.csv = FALSE )
2、 进行生存分析
###TCGAbiolinks包中自带的进行生存分析的函数是TCGAanalyze_survival，但是我们也有其他的包可以进行生存分析，在第二部分中会进行对比说明
####TCGAanalyze_survival 的用法：
####利用得到的临床数据探索性别对患者生存的影响
TCGAanalyze_survival(clin.STAD,
clusterCol="gender",
risk.table = F,
xlim = c(100,1000),
ylim = c(0.4,1),
conf.int = T,
color = c("Dark2"))
3 、 探索基因表达对生存的影响
####随机选取PGA5基因为例探索单个基因表达的情况对患者生存的影响
# 1、取特定基因（这里以PGA5基因为例）在癌症样本中的表达
# 2、先读取基因表达矩阵
# 3、提取特定基因在癌症样本中的表达，选出的样本都是肿瘤样本
samplesTP <- TCGAquery_SampleTypes(colnames(dataFilt_STAD_final), typesample = c("TP"))
PGA5 <- dataFilt_STAD_final[c("PGA5"),samplesTP]
View(PGA5)
# 4、修改样本名称（原是TCGA-DD-AAD5-01A-11R-A41C-07），
# 但是clin.STAD数据中的样本名称是12位，
# 因为后期要把 clin.STAD 数据和表达数据结合起来，需要将名字简化成TCGA-DD-AAD5
names(PGA5) <- sapply(strsplit(names(PGA5),'-'),function(x) paste0(x[1:3],collapse="-"))
PGA5 <-as.data.frame(cbind(names(PGA5),PGA5))
colnames(PGA5) <- c("submitter_id","PGA5")
# 5、合并PGA5基因与临床数据
clin.STAD <- merge(clin.STAD,PGA5,by="submitter_id")
View(clin.STAD)
# 6、从 clin.STAD 中选取进行生存分析需要的数据组成新的
##数据框（有barcode、生存状态、死亡时间、随访时间、ABCB1基因的表达数据）
df<-subset(clin.STAD,select =c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,PGA5))
View(df)
# 7、去除PGA5基因表达数据缺失的样本
#去掉 NA
df <- df[!is.na(df$PGA5),]
# 8、根据ABCB1基因的表达情况进行分组，取基因的平均值，大于平均值的为H，小于为L
df$PGA5 <- as.numeric(df$PGA5)
df$exp <- ''
df[df$PGA5 >= mean(df$PGA5),]$exp <- "H"
df[df$PGA5 < mean(df$PGA5),]$exp <- "L"
# 9、用 TCGAanalyze_survival 函数进行生存分析
TCGAanalyze_survival(df,
clusterCol="exp",
risk.table = FALSE,
conf.int = FALSE,
color = c("Dark2"))
# 10、除了使用TCGAanalyze_survival 函数，我们还可以使用 survival、survminer这两个包来进行生存分析。我们可以来对比一下他们的区别。
# 10.1 用status表示患者结局，1表示删失，2表示死亡
df2 <- df
# 10.2 将status表示患者结局，1表示删失，2表示死亡
df2[df2$vital_status=='Dead',]$vital_status <- 2
df2[df2$vital_status=='Alive',]$vital_status <- 1
df2$vital_status <- as.numeric(df2$vital_status)
df2$time <- df2$days_to_death
df2$time[which(is.na(df2$time))] <- df2$days_to_last_follow_up[which(is.na(df2$time))]
View(df2)
library(survival)
library(survminer)
# 建模
fit <- survfit(Surv(time, vital_status)~exp, data=df2) # 根据表达建模
# 显示P value
surv_pvalue(fit)$pval.txt
# 画图
ggsurvplot(fit,pval=TRUE)
