【原】TCGA的28篇教程- 指定癌癥查看感興趣基因的表達(dá)量

健明 2021-07-14

展開全文

m長期更新列表：
使用R語言的cgdsr包獲取TCGA數(shù)據(jù)（cBioPortal）TCGA的28篇教程- 使用R語言的RTCGA包獲取TCGA數(shù)據(jù) （離線打包版本）TCGA的28篇教程- 使用R語言的RTCGAToolbox包獲取TCGA數(shù)據(jù) （FireBrowse portal）TCGA的28篇教程- 批量下載TCGA所有數(shù)據(jù) （ UCSC的 XENA）TCGA的28篇教程- 數(shù)據(jù)下載就到此為止吧

本教程目錄:

- 文章來源

- 那么首先需要得到上圖的基因列表

- 使用下載這些基因在GBM的芯片數(shù)據(jù)的表達(dá)量。

- 下載感興趣基因的表達(dá)矩陣

- 畫散點圖

- 相關(guān)系數(shù)的總結(jié)

文章來源

這里我重現(xiàn)一篇CELL文章的某個分析要點,，具體見我博客的解讀：http://www./3072.html

在TCGA的GBM的芯片表達(dá)數(shù)據(jù)里面可以看到以下正相關(guān)關(guān)系：

NCR2 expression
upregulation of PDGF-DD-induced NK cell cytokine genes (芯片數(shù)據(jù))
downregulation of tumor cell-cycle genes
greater survival

TCGA GBM mRNA gene expression data obtained using the Affymetrix HT Human Genome U133a microarray platform (n = 539 patients)was downloaded through the UCSC data portal (https://) and matched to the gene expression data (Figures 4C and 4D). 7 of 9 cytokine genes (Figure 4C) and 27 of 34 cell cycle genes (Figure 4D) were matched with the GBM cohort.

分別如下圖

NCR2基因表達(dá)與細(xì)胞因子相關(guān)基因表達(dá)正相關(guān)

NCR2表達(dá)與細(xì)胞周期相關(guān)基因表達(dá)負(fù)相關(guān)

NCR2高表達(dá)的保護(hù)作用

那么首先需要得到上圖的基因列表

在KEGG官網(wǎng)可以查看細(xì)胞周期基因列表：http://www./kegg-bin/show_pathway?hsa04110 , 而細(xì)胞因子是：http://www./kegg-bin/show_pathway?hsa04060

這里可以用R語言本身的包KEGG.db來獲取.

library(KEGG.db)
ls("package:KEGG.db")

##  [1] "KEGG" "KEGGENZYMEID2GO"  "KEGGEXTID2PATHID"
##  [4] "KEGGGO2ENZYMEID"  "KEGGMAPCOUNTS" "KEGGPATHID2EXTID"
##  [7] "KEGGPATHID2NAME"  "KEGGPATHNAME2ID"  "KEGG_dbInfo"  
## [10] "KEGG_dbconn"   "KEGG_dbfile"   "KEGG_dbschema"

cellcycle_genes=KEGGPATHID2EXTID[['hsa04110']]
cytokine_genes=KEGGPATHID2EXTID[['hsa04060']]

使用下載這些基因在GBM的芯片數(shù)據(jù)的表達(dá)量,。

library(cgdsr)
library(DT)
# Get list of cancer studies at server
## 獲取有哪些數(shù)據(jù)集
# save(all_TCGA_studies,all_dataset,all_tables,cellcycle_expr,cytokine_expr,NCR2_expr,file = 'GBM_microarray_TCGA.Rdata')
load(file = 'GBM_microarray_TCGA.Rdata')
mycgds <- CGDS("http://www./public-portal/")
#all_TCGA_studies <- getCancerStudies(mycgds)
all_TCGA_studies[grepl('gbm',all_TCGA_studies$cancer_study_id),]

##  cancer_study_id name
## 55 gbm_tcga_pub2013  Glioblastoma (TCGA, Cell 2013)
## 56  gbm_tcga_pub   Glioblastoma (TCGA, Nature 2008)
## 57   gbm_tcga Glioblastoma Multiforme (TCGA, Provisional)
## 97  lgggbm_tcga_pub Merged Cohort of LGG and GBM (TCGA, Cell 2016)
##  description
## 55 <a href="http://cancergenome./">The Cancer Genome Atlas (TCGA)</a> Glioblastoma project.<br> <a >Raw data via the TCGA Data Portal</a>.
## 56 <a href="http://cancergenome./">The Cancer Genome Atlas (TCGA)</a> Glioblastoma project. 206 primary glioblastoma samples.<br> <a >Raw data via the TCGA Data Portal</a>.
## 57  TCGA Glioblastoma Multiforme; raw data at the <A >NCI</A>.
## 97  Whole-exome sequencing from TCGA LGG and GBM cases

帶有GBM的項目有4個。

下載感興趣基因的表達(dá)矩陣

很明顯那篇文章作者提到了Affymetrix HT Human Genome U133a microarray platform (n = 539 patients),，所以用的是gbm_tcga那個數(shù)據(jù)集,，用下面的代碼下載：

## 獲取在 "gbm_tcga" 數(shù)據(jù)集中有哪些表格（每個表格都是一個樣本列表）
#all_tables <- getCaseLists(mycgds, "gbm_tcga")
DT::datatable(all_tables[,1:4],
  extensions = 'FixedColumns',
  options = list(
 #dom = 't',
 scrollX = TRUE,
 fixedColumns = TRUE
  ))

## 而后獲取可以下載哪幾種數(shù)據(jù)，一般是mutation,，CNV和表達(dá)量數(shù)據(jù)
#all_dataset <- getGeneticProfiles(mycgds, "gbm_tcga")
DT::datatable(all_dataset,
  extensions = 'FixedColumns',
  options = list(
 #dom = 't',
 scrollX = TRUE,
 fixedColumns = TRUE
  ))

all_dataset[4,]

##   genetic_profile_id genetic_profile_name
## 4 gbm_tcga_mrna_U133 mRNA expression (U133 microarray only)
##   genetic_profile_description
## 4 mRNA expression data from the Affymetrix U133 microarray.
##   cancer_study_id genetic_alteration_type show_profile_in_analysis_tab
## 4 740   MRNA_EXPRESSION   false

all_tables[9,1:4]

##   case_list_id  case_list_name
## 9 gbm_tcga_mrna_U133 Tumor Samples with mRNA data (U133 microarray only)
##   case_list_description cancer_study_id
## 9 All samples with mRNA expression data (528 samples) 740

my_dataset <- 'gbm_tcga_mrna_U133'
my_table <- "gbm_tcga_mrna_U133" 
#cellcycle_expr <- getProfileData(mycgds, cellcycle_genes, my_dataset, my_table)
dim(cellcycle_expr)

## [1] 528 124

#cytokine_expr <- getProfileData(mycgds, cytokine_genes, my_dataset, my_table)
dim(cytokine_expr)

## [1] 528 265

#NCR2_expr <- getProfileData(mycgds, 'NCR2', my_dataset, my_table)

畫散點圖

cor_plot <- function(x,y){
  #x=NCR2_expr$NCR2
  #y=cellcycle_expr$CCNB1
  plot(x,y, xlab = 'NCR2', ylab = 'gene')
  model = lm(y ~ x)
  summary(model) 
  int =  model$coefficient["(Intercept)"]
  slope =model$coefficient["x"]
  abline(int, slope,
   lty=1, lwd=2, col="red")   
  r= format(cor(x,y),digits = 4)
  p= format(cor.test(x,y)$p.value,digits = 4)
  title(main = paste0('p value=',p), 
  sub  = paste0('correlation=',r))
}

cor_plot(x=NCR2_expr$NCR2,y=cellcycle_expr$CCNB1)

cor_plot(x=NCR2_expr$NCR2,y=cellcycle_expr$BUB1B)

cor_plot(x=NCR2_expr$NCR2,y=cytokine_expr$CCL1)

cor_plot(x=NCR2_expr$NCR2,y=cytokine_expr$CCL4)

這些圖跟作者的差不多了,，但是很明顯作者的展現(xiàn)方式并不是最好的。因為這樣的散點圖太多了,。其實只需要p值以及相關(guān)性系數(shù)即可,。