基因在任意癌症表達量相關性

2020 年 3 月 30 日
筆記

序言

第四單元第一講：計算基因在任意癌症表達量相關性課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

從題目可以看到，這次的主角有兩個：基因和癌症中的表達量

針對第一個：我們要知道有哪些基因

從這個表中複製基因名，然後放到R中，但要注意它們中間是,分隔，因此要使用str_split 拆分成單獨的字元串：

library(stringr)  vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'  vCAF=unlist(str_split(vCAF,', ')) # 或者直接使用 as.character(str_split(vCAF, ', '))  mCAF='Dcn, Col12a1, Mmp2, Lum, Mrc2, Bicc1, Lrrc15, Mfap5, Col3A1, Mmp14, Spon1, Pdgfrl, Serpinf1, Lrp1, Gfpt2, Ctsk, Cdh11, Itgbl1, Col6a2, Postn, Ccdc80, Lox, Vcan, Col1a1, Fbn1, Col1a2, Pdpn, Col6a1, Fstl1, Col5a2, Aebp1'  mCAF=unlist(str_split(mCAF,', '))    > vCAF  [1] "Esam"   "Gng11"  "Higd1b" "Cox4i2" "Cygb"   "Gja4"   "Eng"  > head(mCAF)  [1] "Dcn"     "Col12a1" "Mmp2"    "Lum"     "Mrc2"    "Bicc1"

看到基因名的開頭大寫，其餘小寫，就說明是小鼠的基因名

針對第二個：如何獲取癌症基因表達量資訊

文章對四種癌症進行了討論：breast cancer, pancreatic ductal adenocarcinoma, lung adenocarcinoma, and renal clear cell carcinoma

目的就是分別畫這樣一張圖：

分析這張圖片：這是一個相關性圖，如果要做相關性的圖，就要有數值型的數據，那麼就是基因表達量了。我們現在有了基因名，缺的就是一個表達矩陣。因此如何獲取表達矩陣就是最大的一個問題了

首先下載乳腺癌的表達矩陣

網址：https://xenabrowser.net/datapages/

但是如果從這裡直接搜索BRCA的話，會有兩個結果（這裡選擇GDC的表達矩陣）：

[GDC TCGA Breast Cancer (BRCA)](https://xenabrowser.net/datapages/?cohort=GDC TCGA Breast Cancer (BRCA)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (20 datasets) 包含了6萬多個基因的表達矩陣（60,489 identifiers X 1217 samples），因為它直接使用的GenCode注釋文件，不管編碼與否都算作基因；基因名用的Ensembl；使用了log2(count+1)

[TCGA Breast Cancer (BRCA)](https://xenabrowser.net/datapages/?cohort=TCGA Breast Cancer (BRCA)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (30 datasets) 其中只有2萬多個基因(20,531 identifiers X 1218 samples)，並且得到的是RSEM標準化表達量（使用了log2(norm_count+1)）；使用SYMBOL基因名

點擊下圖鏈接開始下載：文件大小133M

然後讀入乳腺癌的表達矩陣

使用fread函數

library(data.table)  filepath <- file.choose()# 然後會彈出來一個對話框，找到自己下載的TCGA-BRCA.htseq_counts.tsv.gz，點OK，然後這個文件的路徑就保存在了filepath  a=fread(filepath ,data.table=F)  dim(a)  # [1] 60488  1218  a[1:4,1:4]  # Ensembl_ID TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A201-01A  # 1 ENSG00000000003.13         8.787903        12.064743        11.801304  # 2  ENSG00000000005.5         0.000000         2.807355         4.954196  # 3 ENSG00000000419.11        11.054604        11.292897        11.314017  # 4 ENSG00000000457.12        10.246741         9.905387        11.117643

接著進行ID轉換，Ensembl =》 Symbol ID

需要用到人類的物種注釋包：org.Hs.eg.db

library(org.Hs.eg.db)  # 先看看包的簡介  > org.Hs.eg.db  OrgDb object:  | DBSCHEMAVERSION: 2.1  | Db type: OrgDb  | Supporting package: AnnotationDbi  | DBSCHEMA: HUMAN_DB  | ORGANISM: Homo sapiens  | SPECIES: Human  | EGSOURCEDATE: 2019-Apr26  | EGSOURCENAME: Entrez Gene  | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA  ...  # 再看看這個注釋包里有什麼資訊  > head(ls("package:org.Hs.eg.db"))  [1] "org.Hs.eg"          "org.Hs.eg.db"       "org.Hs.egACCNUM"    "org.Hs.egACCNUM2EG" "org.Hs.egALIAS2EG"  [6] "org.Hs.egCHR"  # 然後看看其中Ensembl的基因是什麼樣子  > head(toTable(org.Hs.egENSEMBL))    gene_id      ensembl_id  1       1 ENSG00000121410  2       2 ENSG00000175899  3       3 ENSG00000256069  4       9 ENSG00000171428  5      10 ENSG00000156006  6      12 ENSG00000196136

發現相對於我們得到TCGA的Emsemble ID，它沒有小數點後面的部分，因此我們也需要切割Ensembl ID =>str_split()

library(stringr)  esid=str_split(a$Ensembl_ID,                   '[.]',simplify = T)[,1]  > head(esid)  [1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" "ENSG00000000938"  rownames(a)=esid

開始進行ID轉換 => select()或bitr()

這裡二者結果一樣

# 第一種方式：官方函數  e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" ,  "SYMBOL" ),keytype = 'ENSEMBL')  dim(e2s)  # [1] 60686     2  #其中很大一部分的Ensemble ID是沒有Symbol對應的。如果出去symbol為NA的值：最後也就剩下25591個基因  nrow(e2s)-sum(is.na(e2s$SYMBOL))  # [1] 25591    # 第二種方式：R包函數  library(clusterProfiler)  gene_tr <- bitr(esid, fromType = "ENSEMBL",                                    toType = "SYMBOL",                                    OrgDb = org.Hs.eg.db)   nrow(gene_tr)  # 25591    identical(e2s$SYMBOL[!is.na(e2s$SYMBOL)],gene_tr$SYMBOL)  # [1] TRUE

這樣我們就同時擁有了Ensembl ID和Symbol ID：在TCGA矩陣中獲取表達量用Emsembl ID，可視化用Symbol ID

# 小鼠基因變大寫，然後挑出來存在於e2s的基因  vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL]  mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL]  # 得到匹配基因的Ensembl ID(總共38個基因)，準備去獲取表達量  ng=e2s[match(c(vCAF,mCAF),e2s$SYMBOL),1]  mat=a[ng,]  mat=mat[,-1]    dim(mat)  # [1]   38 1217  > mat[1:4,1:4]                  TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A201-01A TCGA-E2-A14T-01A  ENSG00000149564        10.279611        10.059344        10.907642        10.458407  ENSG00000127920         9.776433         9.726218        10.948367        10.496854  ENSG00000131097         5.614710         4.857981         5.930737         6.658211  ENSG00000131055         6.022368         6.129283         6.629357         6.475733

最後就是計算相關性，準備繪製熱圖

計算相關性就是利用cor()函數，但是有個問題，它是對行處理還是對列處理？

# 都不用去搜索，自己隨便測試一下就能出來結果(新建一個矩陣，然後對它進行默認的相關性分析)  > matrix(1:10,nrow = 2)       [,1] [,2] [,3] [,4] [,5]  [1,]    1    3    5    7    9  [2,]    2    4    6    8   10  > cor(matrix(1:10,nrow = 2))       [,1] [,2] [,3] [,4] [,5]  [1,]    1    1    1    1    1  [2,]    1    1    1    1    1  [3,]    1    1    1    1    1  [4,]    1    1    1    1    1  [5,]    1    1    1    1    1  # 很明顯，這是對列進行處理

因此，我們如果想看基因之間的相關性，就將上面的mat矩陣轉置一下就可以：

M=cor(t(mat))  colnames(M)=c(vCAF,mCAF)  rownames(M)=c(vCAF,mCAF)  # 然後為了避免高表達量對許多低表達量的遮蓋，我們進行一個標準化處理  n=t(scale(t( M )))  n[n>2]=2  n[n< -2]= -2  pheatmap::pheatmap(n,cluster_rows = F,cluster_cols = F)

典型的相關性熱圖

補充

在Xena資料庫搜索pancreatic 會有5個數據集：

[GDC TCGA Pancreatic Cancer (PAAD)](https://xenabrowser.net/datapages/?cohort=GDC TCGA Pancreatic Cancer (PAAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (14 datasets) 60,489 identifiers X 182 samples；log2(count+1) 下載地址：https://gdc.xenahubs.net/download/TCGA-PAAD.htseq_counts.tsv.gz （20M）
[Pancreatic Cancer (Balagurunathan 2008)](https://xenabrowser.net/datapages/?cohort=Pancreatic Cancer (Balagurunathan 2008)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (2 datasets)
[Pancreatic Cancer (Harada 2008)](https://xenabrowser.net/datapages/?cohort=Pancreatic Cancer (Harada 2008)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (2 datasets)
[Pancreatic Cancer (Jones 2008)](https://xenabrowser.net/datapages/?cohort=Pancreatic Cancer (Jones 2008)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (2 datasets)
[TCGA Pancreatic Cancer (PAAD)](https://xenabrowser.net/datapages/?cohort=TCGA Pancreatic Cancer (PAAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (25 datasets) 20,531 identifiers X 183 samples；log2(RSEM_norm_count+1) 下載地址：https://tcga.xenahubs.net/download/TCGA.PAAD.sampleMap/HiSeqV2.gz

在Xena資料庫搜索lung adenocarcinoma 會有3個數據集：

[GDC TCGA Lung Adenocarcinoma (LUAD)](https://xenabrowser.net/datapages/?cohort=GDC TCGA Lung Adenocarcinoma (LUAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (15 datasets) 60,489 identifiers X 585 samples 下載地址：https://gdc.xenahubs.net/download/TCGA-LUAD.htseq_counts.tsv.gz （63.8M）
[Lung Adenocarcinoma (Ding 2008)](https://xenabrowser.net/datapages/?cohort=Lung Adenocarcinoma (Ding 2008)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (2 datasets)
[TCGA Lung Adenocarcinoma (LUAD)](https://xenabrowser.net/datapages/?cohort=TCGA Lung Adenocarcinoma (LUAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443) (27 datasets) 20,531 identifiers X 576 samples 下載地址：https://tcga.xenahubs.net/download/TCGA.LUAD.sampleMap/HiSeqV2.gz