作者：Squirrelity

2022-07-18 補充說明
最近R更新了，很多包都用不了，如果遇到報錯或者是運行不了有可能是因為版本問題。

一、加載對應的R包
這里用到十三個包（距離上次更新之后又新增了不少方法/包）：

library("TCGAbiolinks") library("plyr") library("limma") library("biomaRt") library("SummarizedExperiment") library("stringr") library("ggplotify") library("patchwork") library("cowplot") library('DESeq2') library('edgeR') library("dplyr") library("rtracklayer")

光下載都費了不少功夫www，下面把install代碼放出來(｡･?･)ﾉﾞ直接install.packages()日常失敗我就不說了
大部分包都可以在bioconductor中找到，有遺漏的可以去官網找下載代碼
https://bioconductor.org/

#plyr下載 唯一可以直接install的包啊哈哈 install.packages("plyr") #if(!require(stringr)) #install.packages（）下載法 if(!require(stringr))install.packages('stringr') if(!require(ggplotify))install.packages("ggplotify") if(!require(patchwork))install.packages("patchwork") if(!require(cowplot))install.packages("cowplot")#更新BiocManager3.15，R版本為4.2.0 if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install(version = "3.15")#tccgbiolinks包穩定版本安裝 if (!requireNamespace("BiocManager", quietly=TRUE))install.packages("BiocManager") BiocManager::install("TCGAbiolinks")#limma包下載 if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install("limma")#biomaRt包下載 if (!require("BiocManager", quietly = TRUE))install.packages("biomaRt") BiocManager::install("biomaRt")#最麻煩的SummarizedExperiment包：force = TRUE是根據warning后來加的if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install("SummarizedExperiment",force = TRUE) #edgeR包下載 if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install("edgeR") #DESeq2包下載 if (!requireNamespace("BiocManager", quietly=TRUE))install.packages("BiocManager") BiocManager::install("DESeq2")

在使用之前需要加載BIocManager，代碼參考這個：https://bioconductor.org/install/

if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager") BiocManager::install(version = "3.15")

若提示warning：A version of this package for your version of R might be available elsewhere,see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

可以的解決方法有：

1.乖乖在官網下載cran（最不推薦）
2.在Rstudio上方的tools-global options處找到packages，修改默認的global cran，選擇apply-ok即可
3.找到包的官網，用官網提供的代碼下載（首推???）

二、下載數據
在下載之前設置工作路徑：
dir.create（）創建目錄，getwd（）獲取工作路徑，setwd（）設置工作路徑，由于TCGA下載下來的數據包都挺大的，建議還是選一個比較富裕的盤來作為工作路徑。

dir.create("D:\\BioInfoCloud\\TCGABiolinks\\case_study") setwd("D:\\BioInfoCloud\\TCGABiolinks\\case_study")

這里用到的是R包TCGAbiolinks：
可以參照R包TCGAbiolinks官網使用http://www.bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/casestudy.html#Case_study_n_1:_Pan_Cancer_downstream_analysis_BRCA
示例：

query <- GDCquery(project = "TCGA-BRCA",data.category = "Transcriptome Profiling",data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts")

project選出的是腫瘤項目，而里面用到的都是縮寫，詳見https://blog.csdn.net/Squirrelity/article/details/124259330?spm=1001.2014.3001.5501

建議去TCGA官網repository一邊對照著選所需要的樣本
https://portal.gdc.cancer.gov/repository?facetTab=cases

#下載到本地 GDCdownload(query = query, method = "api") #查看下載的數據 View(query) BRCA.Rnaseq.SE <- GDCprepare( query = query, save = TRUE, save.filename = "brca.rda") BRCAMatrix <- assay(BRCA.Rnaseq.SE,"unstranded") #記住這個文件BRCA.RNAseq_CorOutliers BRCA.RNAseq_CorOutliers <- TCGAanalyze_Preprocessing(BRCA.Rnaseq.SE)

三、ID轉換
下載下來的BRCA.RNAseq_CorOutliers為ENTREZID，而我們肯定是需要圖片顯示基因名而不是ENTREZID，因此進行數據轉換，這里用到包括但不限于的dplyr包和rtracklayer包

library("dplyr") library("rtracklayer")

ID轉換分為四步：
1.導入數據：BRCA.RNAseq_CorOutliers和人類基因組注釋文件；

data=read.table(BRCA.RNAseq_CorOutliers,header=T,sep='\t') #把行名改為gene_id，與gtf保持一致 colnames(data)[1] <- "gene_id"

對照人類基因組注釋文件，對BRCA.RNAseq_CorOutliers進行ID轉換
其中，人類基因組注釋文件參考http://www.360doc.com/content/21/1028/10/77506210_1001626502.shtml

#處理人類基因組注釋文件的數據 gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz') gtf <- as.data.frame(gtf) save(gtf,file="人類基因組注釋文件.Rda") gtf <-load(file="人類基因組注釋文件.Rda") #根據條件篩選基因（大篩選） a = dplyr::filter(gtf,type=="gene") dim（a） #只要gene_name,gene_id,gene_biotype這三行 b = dplyr::select(a,c(gene_name,gene_id,gene_biotype))

2.數據預處理

#ENTREZID帶有，這里去除小數點及后邊的數字（我用excel處理的，ctrl+F無字符替換.*） data1 <- separate(data,gene_id,into = c("gene_id"),sep="\\.")

3.數據處理

#根據gene id 合并文件 c = dplyr::inner_join(b,data,by="gene_id") #去掉2,3列，基因名再去重 d=select(c,-gene_id,-gene_biotype) data1=distinct(d,gene_name,.keep_all = T) #把行名由數字換成基因 rownames(data1)<- data1[,1] data1<-data1[,-1]

4.數據后處理

#如下載的數據取了log2（count-1）這里再返回count data2 <- 2^data1 -1 write.csv(data2,file="data2.csv") data2 <- read.csv("data2.csv") #重新用read打開整行的-會變成.因此需要恢復原來的行名 colnames(data2) <- colnames(BRCA.RNAseq_CorOutliers)

四、差異表達分析
1.參考網址
代碼示例參考這個包的文檔http://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html
TCGA可視化教程
https://www.jianshu.com/p/d3e481f0187a
https://cloud.tencent.com/developer/article/1778874
http://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/casestudy.html

2.數據預處理

#處理后可得熱圖判斷樣本相似性 dataPrep <- TCGAanalyze_Preprocessing(object = BRCA.RNAseq_CorOutliers, cor.cut = 0.6 )

*上面代碼生成的圖如下所示。組內的樣本相似性都很高，符合預期。

3.對數據進行標準化處理+質控+差異化分析
**TCGAanalyze_LevelTab（）**將差異表達基因在正常和腫瘤組織中的表達量數據添加到差異表達分析結果中的主要用法：

TCGAanalyze_LevelTab(FC_FDR_table_mRNA, typeCond1, typeCond2, TableCond1, TableCond2, typeOrder = TRUE)

#數據標準化 dataNorm <- TCGAanalyze_Normalization(tabDF = data2,geneInfo = geneInfoHT,method = "gcContent" ) #數據質控 dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,method = "quantile", qnt.cut = 0.25 ) #分組：NT正常組織組 TP癌癥組織組*# selection of normal samples "NT" samplesNT <- TCGAquery_SampleTypes(barcode = colnames(dataFilt),typesample = c("NT") )) # selection of tumor samples "TP" samplesTP <- TCGAquery_SampleTypes(barcode = colnames(dataFilt), typesample = c("TP") )# 差異表達分析 dataDEGs <- TCGAanalyze_DEA(mat1 = dataFilt[,samples.solid.tissue.normal],mat2 = dataFilt[,samples.primary.tumour],Cond1type = "Normal",Cond2type = "Tumor",fdr.cut = 0.01 ,logFC.cut = 2,method = "glmLRT",pipeline = "edgeR" ) #在正常和腫瘤樣本中差異基因的表達值 dataDEGsFiltLevel <- TCGAanalyze_LevelTab(FC_FDR_table_mRNA = dataDEGs,typeCond1 = "Tumor",typeCond2 = "Normal",TableCond1 = dataFilt[,samplesTP1],TableCond2 = dataFilt[,samplesNT1] )

得到的dataDEGsFiltLevel文件按logFC絕對值排序可得最顯著的top差異表達基因（excel表處理）

五、可視化
ps：得到的圖片有的可以直接看，有的保存在工作路徑上了。
1.PCA主成分分析

TCGAvisualize_PCA()實現主成分分析的主要用法:

TCGAvisualize_PCA(dataFilt, dataDEGsFiltLevel, ntopgenes, group1, group2)

#標準化 dataNorm <- TCGAbiolinks::TCGAanalyze_Normalization(data2, geneInfo)#質量控制 dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,method = "quantile", qnt.cut = 0.25)#選擇正常樣本 group1 <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("NT")) #選擇癌癥樣本 group2 <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("TP"))#Principal Component Analysis plot for ntop selected DEGs pca.top200 <- TCGAvisualize_PCA(dataFilt,dataDEGsFiltLevel, ntopgenes = 200,group1, group2)

上面代碼生成的圖如下所示。

2.火山圖

#為了做圖的需要，突出顯示logFC≥8的gene名稱 DEG.BRCA.filt<-dataDEGs[which(abs(dataDEGs$logFC) >= 8), ] str(DEG.BRCA.filt) #'data.frame': 29 obs. of 5 variables: #說明共有29個基因滿足|logFC|≥8TCGAVisualize_volcano(dataDEGs$logFC, dataDEGs$FDR,filename = "TumorvsNormal_FC8.edgeR.pdf", xlab = "logFC",names = rownames(dataDEGs), show.names = "highlighted",x.cut = 1, y.cut = 0.01, highlight = rownames(dataDEGs)[which(abs(DEG.LIHC.edgeR$logFC) >= 8)],highlight.color = "orange",title = "volcano plot by edgeR")

上面代碼生成的圖如下所示。突出顯示了logFC≥8的gene名稱

3.GO功能分析條形圖
TCGAbiolinks 輸出條形圖，其中包含三個本體的主要類別（分別為GO：生物過程、GO：細胞成分和GO：分子功能）的基因數量。

ansEA <- TCGAanalyze_EAcomplete(TFname = "DEA genes Normal Vs Tumor",RegulonList = dataDEGs$gene_name ) TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),GOBPTab = ansEA$ResBP,GOCCTab = ansEA$ResCC,GOMFTab = ansEA$ResMF,PathTab = ansEA$ResPat,nRGTab = dataDEGs$gene_name,nBar = 10 )

上面代碼生成的圖如下所示。

總結

以上是生活随笔為你收集整理的TCGA差异表达分析|2022.5.1更新的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。