部分一
第一步、數(shù)據(jù)下載階段
####(1.1)、加載包
rm(list=ls())
library(TCGAbiolinks)
library(plyr)
library(limma)
library(biomaRt)
library(SummarizedExperiment)
####(1.2)
####GDCquery()篩選需要的數(shù)據(jù)墓造，TCGAbiolinks包下載TCGA數(shù)據(jù)進行表達差異分析
query <- GDCquery(project = "TCGA-STAD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
####(1.3)
####getResults(query, rows, cols)根據(jù)指定行名或列名從query中獲取結(jié)果,此處用來獲得樣本的barcode
#### 此處共檢索出407個barcodes
samplesDown <- getResults(query,cols=c("cases"))
# 從samplesDown中篩選出TP（實體腫瘤）樣本的barcodes
# TCGAquery_SampleTypes(barcode, typesample)
# TP代表PRIMARY SOLID TUMOR；
# NT-代表Solid Tissue Normal(其他組織樣本可參考學(xué)習(xí)文檔)
##此處共檢索出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)
###設(shè)置barcodes參數(shù)棚菊，篩選符合要求的375個腫瘤樣本數(shù)據(jù)和32正常組織數(shù)據(jù)
queryDown <- GDCquery(project = "TCGA-STAD",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts",
barcode = c(dataSmTP, dataSmNT))
#barcode參數(shù)：根據(jù)傳入barcodes進行數(shù)據(jù)過濾
第二步：GDCdownload()下載GDCquery()得到的結(jié)果
####2.1 下載數(shù)據(jù)，默認(rèn)存放位置為當(dāng)前工作目錄下的GDCdata文件夾中。
GDCdownload(queryDown,
method = "api",
directory = "GDCdata",
files.per.chunk = 6)
#method ；"API"或者"client"迫靖。"API"速度更快，但是容易下載中斷兴使。
#directory：下載文件的保存地址系宜。Default: GDCdata。
#files.per.chunk = NULL:使用API下載大文件的時候发魄，可以把文件分成幾個小文件來下載盹牧，可以解決下載容易中斷的問題。
第三步励幼、數(shù)據(jù)處理
####3.1：GDCprepare()將前面GDCquery()的結(jié)果準(zhǔn)備成R語言可處理的SE（SummarizedExperiment）文件
####讀取下載的數(shù)據(jù)并將其準(zhǔn)備到R對象中汰寓，在工作目錄生成（save=TRUE）STAD_case.rda文件
#### GDCprepare():Prepare GDC data,準(zhǔn)備GDC數(shù)據(jù)，使其可用于R語言中進行分析
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename =
"STAD_case.rda")
第四步：TCGAanalyze_Preprocessing()對數(shù)據(jù)進行預(yù)處理：
####使用spearman相關(guān)系數(shù)去除數(shù)據(jù)中的異常值
# 去除dataPrep1中的異常值苹粟，dataPrep1數(shù)據(jù)中含有腫瘤組織和正常組織的數(shù)據(jù)
# TCGAanalyze_Preprocessing(object, cor.cut = 0, filename = NULL,
# width = 1000, height = 1000, datatype = names(assays(object))[1])
# 函數(shù)功能描述：Array Array Intensity correlation (AAIC) and correlation boxplot to define outlier
dataPrep2 <- TCGAanalyze_Preprocessing(object = dataPrep1,
cor.cut = 0.6,
datatype = "HTSeq - Counts")
#將預(yù)處理后的數(shù)據(jù)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毛好，
# 這里設(shè)置cpe=0.6（cpe是派生的共識度量望艺，是將所有方法的標(biāo)準(zhǔn)含量歸一化后的均值純度水平，以使它們具有相等的均值和標(biāo)準(zhǔn)差）
# 篩選腫瘤純度大于等于60%的樣本數(shù)據(jù)
purityDATA <- TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)
####這里篩選的60%的樣本數(shù)據(jù)并沒有用成功肌访，這里直接賦值
# filtered 為被過濾的數(shù)據(jù)找默， pure_barcodes是我們要的腫瘤數(shù)據(jù)
####過濾的數(shù)據(jù)太多，重新分配
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”每個由數(shù)字或其他模式的類似矩陣的對象表示。
# 行通常表示感興趣的基因組范圍和列代表樣品旨剥。
library("SummarizedExperiment")
#將結(jié)果寫入文件“puried.LIHC.cancer.csv”
rownames(puried_data)<-rowData(dataPrep1)$external_gene_name
write.csv(puried_data,file = "puried.STAD.csv",quote = FALSE)
第八步：進行表達矩陣標(biāo)準(zhǔn)化和過濾咧欣，得到用于差異分析的表達矩陣
####1、`TCGAanalyze_Normalization()`使用EDASeq軟件包標(biāo)準(zhǔn)化mRNA轉(zhuǎn)錄本和miRNA轨帜。
dataNorm <- TCGAanalyze_Normalization(tabDF = puried_data,
geneInfo = geneInfo,
method = "gcContent")
####2魄咕、將標(biāo)準(zhǔn)化后的數(shù)據(jù)再過濾，去除掉表達量較低(count較低)的基因蚌父，得到最終的數(shù)據(jù)
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)喝滞，富集分析
####設(shè)置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, #顯示條形圖的數(shù)量
filename = "TCGAvisualize_EAbarplot_Output.pdf")
第十一右遭，聚類分析
res.hc <- TCGAanalyze_Clustering(Data_DEGs,
method = "hclust",
methodHC = "ward.D2")
#########################################################################################
#########################################################################################
#########################################################################################
部分二
1、 下載臨床數(shù)據(jù)
clin.STAD <- GDCquery_clinic("TCGA-STAD", "clinical",save.csv = FALSE )
2缤削、 進行生存分析
###TCGAbiolinks包中自帶的進行生存分析的函數(shù)是TCGAanalyze_survival窘哈，但是我們也有其他的包可以進行生存分析，在第二部分中會進行對比說明
####TCGAanalyze_survival 的用法：
####利用得到的臨床數(shù)據(jù)探索性別對患者生存的影響
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數(shù)據(jù)中的樣本名稱是12位骇窍，
# 因為后期要把 clin.STAD 數(shù)據(jù)和表達數(shù)據(jù)結(jié)合起來，需要將名字簡化成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基因與臨床數(shù)據(jù)
clin.STAD <- merge(clin.STAD,PGA5,by="submitter_id")
View(clin.STAD)
# 6像鸡、從 clin.STAD 中選取進行生存分析需要的數(shù)據(jù)組成新的
##數(shù)據(jù)框（有barcode、生存狀態(tài)、死亡時間只估、隨訪時間、ABCB1基因的表達數(shù)據(jù)）
df<-subset(clin.STAD,select =c(submitter_id,vital_status,days_to_death,days_to_last_follow_up,PGA5))
View(df)
# 7着绷、去除PGA5基因表達數(shù)據(jù)缺失的樣本
#去掉 NA
df <- df[!is.na(df$PGA5),]
# 8蛔钙、根據(jù)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 函數(shù)進行生存分析
TCGAanalyze_survival(df,
clusterCol="exp",
risk.table = FALSE,
conf.int = FALSE,
color = c("Dark2"))
# 10彬向、除了使用TCGAanalyze_survival 函數(shù)兼贡，我們還可以使用 survival、survminer這兩個包來進行生存分析娃胆。我們可以來對比一下他們的區(qū)別遍希。
# 10.1 用status表示患者結(jié)局，1表示刪失里烦，2表示死亡
df2 <- df
# 10.2 將status表示患者結(jié)局凿蒜，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) # 根據(jù)表達建模
# 顯示P value
surv_pvalue(fit)$pval.txt
# 畫圖
ggsurvplot(fit,pval=TRUE)
