TCGA數據挖掘（四）：表達差異分析（4）

2019 年 10 月 4 日
筆記

在之前我們的文章：TCGA數據挖掘（三）：表達差異分析中，我們利用的是TCGAbiolinks包中的TCGAanalyze_DEA函數進行差異表達分析，我們也提到可以選擇基於limma或edgeR包進行分析，TCGA數據挖掘（三）：表達差異分析這一講中我們利用的是edgeR包，之後我們在文章：TCGA數據挖掘（四）：表達差異分析（2）和TCGA數據挖掘（四）：表達差異分析（3）中分別也介紹了其他方法的差異分析，包括edgeR和DESeq包，今天這一講，我們就利用TCGAbiolinks包中的TCGAanalyze_DEA函數基於limma包進行差異分析。

數據下載

基因表達數據的下載

數據下載程式碼和之前的一樣，這裡再提供一次。避免出錯不知道原因。

setwd("E:\BioInfoCloud\TCGABiolinks\case_study1")  # 載入相應的包，可能會需要其他包，提示錯誤就安裝缺少的包。  # 因為每個人已經安裝的包不一樣。  library(TCGAbiolinks)  library(plyr)  library(limma)  library(biomaRt)  library(SummarizedExperiment)  # 請求數據。  query <- GDCquery(project = "TCGA-BRCA",                    data.category = "Transcriptome Profiling",                    data.type = "Gene Expression Quantification",                    workflow.type = "HTSeq - Counts")  # 從query中獲取結果表，它可以選擇帶有cols參數的列，並使用rows參數返回若干行。  # 1222個barcode  samplesDown <- getResults(query,cols=c("cases"))  # samplesDown中篩選出TP(主要實體瘤)樣本的barcode  # 1102個barcode  dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown,                                    typesample = "TP")  # 113個barcode  dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,                                    typesample = "NT")    # 選擇100個正常組織和100個腫瘤組織樣本作為研究對象  dataSmTP_short <- dataSmTP[1:100]  dataSmNT_short <- dataSmNT[1:100]  # 根據前面的篩選，再次請求數據  queryDown <- GDCquery(project = "TCGA-BRCA",                        data.category = "Transcriptome Profiling",                        data.type = "Gene Expression Quantification",                        workflow.type = "HTSeq - Counts",                        barcode = c(dataSmTP_short, dataSmNT_short))  #讀取下載的數據並將其準備到R對象中，在工作目錄生成brca_case1.rda文件，  #同時還產生Human_genes__GRCh38_p12_.rda文件  GDCdownload(query = queryDown)  dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename = "brca_case1.rda")

數據處理

# 去除dataPrep1中的異常值，dataPrep數據中含有腫瘤組織和正常組織的數據  dataPrep <- TCGAanalyze_Preprocessing(object = dataPrep1,                                        cor.cut = 0.6,                                        datatype = "HTSeq - Counts")  #使用EDASeq進行標準化  dataNorm <- TCGAanalyze_Normalization(tabDF = dataPrep,                                        geneInfo = geneInfoHT,                                        method = "gcContent")  #根據分位數一處count較低的基因。  dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,                                    method = "quantile",                                      qnt.cut =  0.25)  # 將計數數據轉換為log2計數/百萬(LogCPM)，估計均值-方差關係並使用此關係計算適當的觀測級別權重。  # 然後，數據即可用於線性建模。  v.dataFilt<-voom(dataFilt)  #in order to have 2 batches with multiple samples and avoid batches with one sample  #the user need to call first the get_IDs function on the top of this script if this has not been done already  c1<-which(get_IDs(dataPrep)$tss=="E9")  c2<-which(get_IDs(dataPrep)$tss=="E2")  #taking log transformed data for exploration of batch effects  TCGAbatch_Correction(tabDF = v.dataFilt$E[,c(c1,c2)], batch.factor="TSS", adjustment=NULL)  #if you get error in plot.new():figure margins too large - please use:  #par(mar=c(1,1,1,1))  #and then rerun the command above

#### After exploration of the batches we can continue to work on the original gene matrix for DEA ###  ###############Tumor purity filtering###########  ###vector containing all TCGA barcodes that hhave 60% tumor purity or more  Purity.BRCA<-TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)$pure_barcodes  ################DEA with Molecular subtypes############  ####Molecular subtypes:  ####Filtering data so all samples have a pam50 subtype for BRCA:    #diff contains TCGA samples that have an available molecular subtype  ###Also Applying Tumor purity filtering by making a set difference  #documentation for TCGA_MolecularSubtype is available  ##########Important########  ####in this example no change happens because all the TCGA barcodes  ####have available subtype info  ###setdiff() is here used to remove samples with no subtype info from the TCGA BRCA samples that satisfy  ###tumor purity criteria  diff<-setdiff(Purity.BRCA, TCGA_MolecularSubtype(colnames(dataPrep[,dataSmTP_short]))$filtered)  dataFilt.brca.cancer<-dataPrep[,diff]  dataFilt.brca.normal<-dataPrep[,dataSmNT_short]  dataFilt.brca<-cbind(dataFilt.brca.normal, dataFilt.brca.cancer)  mol_subtypes<-c(rep("Normal", 100), TCGA_MolecularSubtype(colnames(dataFilt.brca.cancer))$subtypes$subtype)  # All barcodes have available molecular subtype info  mol_subtypes<-make.names(mol_subtypes)  # to convert ENSEMBL id to gene name  rownames(dataFilt.brca)<-rowData(dataPrep1)$external_gene_name  # we redo here the normalization and filtering steps  dataNorm.brca <- TCGAanalyze_Normalization(tabDF = dataFilt.brca,                                             geneInfo = geneInfo,                                             method = "gcContent")  # 將標準化後的數據再過濾，得到最終的數據  dataFilt.brca.final <- TCGAanalyze_Filtering(tabDF = dataNorm.brca,                                         method = "quantile",                                         qnt.cut =  0.25)

表達差異分析

這裡利用的是TCGAbiolinks包中的TCGAanalyze_DEA函數，是基於limma包的差異分析。

DEG.brca.TSS <- TCGAanalyze_DEA(MAT=dataFilt.brca.final,                              pipeline="limma",                              batch.factors = c("TSS"),                              Cond1type = "Normal",                              Cond2type = "Tumor",                              fdr.cut = 0.01 ,                              logFC.cut = 1,                              voom = TRUE,                              Condtypes = mol_subtypes,                              contrast.formula = "Mycontrast = (BRCA.LumA+BRCA.LumB)/2 -Normal")  #in this DEA we use Normal as a reference, thus genes with LogFC > 1 are up regulated in the subtypes  #with respect to the normal samples and viceversa for the LogFC < -1.    write.csv(DEG.brca.TSS, "DEGsBRCA_limma_091018.csv", quote = FALSE)  #3241 genes identified    #to check how many with logFC > 5 (for plotting purposes)  DEG.brca.filt.TSS<-DEG.brca.TSS$Mycontrast[which(abs(DEG.brca.TSS$Mycontrast$logFC) >= 6), ]  #47 genes identified

可視化利用火山圖進行可視化。 TCGAVisualize_volcano(DEG.brca.TSS$Mycontrast$logFC, DEG.brca.TSS$Mycontrast$adj.P.Val, filename = "LuminalABvsNormal_FC6.TSS.pdf", xlab = "logFC", names = rownames(DEG.brca.TSS$Mycontrast), show.names = "highlighted", x.cut = 1, y.cut = 0.01, highlight = rownames(DEG.brca.TSS$Mycontrast)[which(abs(DEG.brca.TSS$Mycontrast$logFC) >= 6)], highlight.color = "orange")

在之前利用edgeR包獲得的結果圖如下，我們簡單比較一下2中方法的差異。

很明顯看到2中方法得到的結果是有差異的，這有時候我們會糾結那種方法好，這個就看你自己的研究啦，什麼結果符合自己的用什麼，當然這有點投機取巧的感覺，最好是2中方法得到後取交集。不過下面我們對這2種方法得到的結果進行一個相關性評估。 limma和edgeR結果的相關性比較library(ggplot2) DEGs_limma <- read.csv("DEGsBRCA_limma_091018.csv", row.names = 1) DEGs_edgeR <- read.csv("DEGsBRCA_edgeR_091018.csv", row.names = 1) genes_toTest <- rownames(DEGs_limma)[1:1000] genes_common <- intersect(genes_toTest, rownames(DEGs_edgeR)) dataset <- subset(DEGs_limma, rownames(DEGs_limma) %in% genes_common, select = "Mycontrast.logFC") dataset <- cbind(DEGs_edgeR$Mycontrast.logFC[match(rownames(dataset),rownames(DEGs_edgeR))], dataset) colnames(dataset) <- c("logFC_edgeR", "logFC_limma") pdf("scatterplot_logFC_limma_edgeR_top1000.pdf", width=5, height=3) ggplot(dataset, aes(x=logFC_limma,y=logFC_edgeR))+geom_point()+ xlab("logFC_limma")+ylab("logFC_edgeR")+ theme(legend.title=element_blank(),axis.title=element_text(size=10),axis.text=element_text(size=10), legend.text=element_text(size=20)) dev.off() cor(dataset$logFC_edgeR, dataset$logFC_limma) #0.9936 這2種方法得到的結果雖然有差異，但很小得到的大多數差異基因是一樣的。

Tags: