學徒帶你一步步從CCLE資料庫裡面根據指定基因在指定細胞系裡面提取表達矩陣進行熱圖可視化

2020 年 3 月 30 日
筆記

昨天生信技能樹發布了學徒作業：學徒作業-在CCLE資料庫裡面根據指定基因在指定細胞系裡面提取表達矩陣很有意思，任務簡單的說就是重複這個圖

下面看一個優秀的學徒作業

首先我們要知道CCLE 是個資料庫。這個資料庫在我心目中和TCGA 還有 GTEx 並稱三大資料庫。TCGA 大家已經很熟悉了，GTEx裡面存儲的主要是正常人體各個組織的轉錄數據。

CCLE 是有1457種癌細胞的表達的資料庫。

下載數據

網站點這裡 CCLE。（https://portals.broadinstitute.org/ccle）

打開後大家需要註冊一下哈然後login

我們可以看到裡面有很多數據，我們挑選RSEM gene和Gene definitions (GENCODE19, GTEx7)(這個用來獲取基因ID) 點擊下載就可以哈。

分割線—————————————————————————————————————————————

數據處理 (獲取表達矩陣)

我們打開R語言，開始導入數據

library(rio)  x1<- import("CCLE_RNAseq_rsem_genes_tpm_20180929.txt")  head(x1)[1:10,1:10]  dim(x1) 數據並不小

boxplot(head(x1[1:10,3:10]))

畫個圖看看數據，看來有點問題，這個稍後我們處理

現在我們導入基因ID數據，就是那個gtf文件哦，一般人肯定是只想到下載表達矩陣文件。

library(data.table)  id <- fread("gencode.v19.genes.v7_model.patched_contigs.gtf.gz",data.table = F)#這一步至關重要，大家可以嘗試變成T 然後看看數據格式    class(id)

dim(id)

id$V9

這裡我們發現第九列內容很多，然後畫紅線的就是我們需要獲取的數據。

看一下第三列

table(id$V3)

第三列有3種，我們這裡需要的是GENE

獲取只有gene的所有行

id1<- id[grepl("gene",id$V3),]

然後我們處理第九列

library(stringr)  options(stringsAsFactors = F)  id1$gene_id<- lapply(id1[,9], function(x){    y=strsplit(x,';')[[1]][2]    strsplit(y,' ')[[1]][3]  }) #把gtf的第9列拆一下獲得EnsembolID  id1$id<- lapply(id1[,9], function(x){    y=strsplit(x,';')[[1]][5]    strsplit(y,' ')[[1]][3]  }) #把gtf的第9列拆一下獲得gene name  idf <- id1[,c(10:11)]  d<- as.data.frame(sapply(idf, function(x) gsub(""", "", x))) #這一步可以幫我們去掉雙引號  x2<- merge(d,x1,"gene_id") #原始矩陣匹配 通過ensembleID  x2 <- x2[,-3]  到這裡我們就獲得了Ensembol ID 對應的 基因ID

準備矩陣行和列

w <- c("CD55", "CD46", "CRIg", "CR1", "Factor H", "Factor I", "FHL1", "C4BP", "Properdin" ,"C1INH")

以上的基因是文章內給的，但是我在矩陣裡面搜了一下

w %in% x2$id

#有部分基因是沒有的。這個時候我懷疑的是我的矩陣有問題，於是我上官網看了一下，官網也沒有。

然後我就把那些沒有的基因去網上分別查了一下。因為很多時候一個基因對應很多名字

詳細可以看生信菜鳥團的這篇文章

w[w %in% x2$id] #這個幾個是有的

w[!w %in% x2$id] #這幾個是沒有的，其實因為名字不一樣。

#CFP=properdin #SERPING1=C1INH #C4BP=C4BPA #VSIG4=CRIg

我手動替換了新的i基因名字，然後我喜歡處理數據框，把它變成了數據框。

w1 <- c("CD55","CD46","VSIG4","CR1","CFH","FHL1" ,"C4BPA", "CFP", "SERPING1")  w2 <- data.frame(id=c("CD55","CD46","VSIG4","CR1","CFH","FHL1" ,"C4BPA", "CFP", "SERPING1"),                   b=rep(1,9))

此處是對細胞名字做準備處理(行名)，先建立向量，然後變成數據框

xb <- c("BT474","BT549",'MDA-MB-231','HCC1937','MDA-MB-361','MDA-MB-436',"AU565","SK-BR-3","MCF-7",'MDA-MB-453')

這裡多了一步是因為我發現所有細胞名字之間沒有 – 作為分隔號

xbb <-  c("BT474","BT549",'MDAMB231','HCC1937','MDAMB361','MDAMB436',"AU565","SKBR3","MCF7",'MDAMB453')  xb1 <- data.frame(n=xb,                    b=rep(1,10))  xb2 <- data.frame(n=xbb,                    b=rep(1,10))

head(x2)

我們發現所有的細胞名字其實還包含了細胞的種類資訊等，但是他們之間都是用 _ 這個符號來分隔。

x3<- merge(w2,x2,"id")

把之前的基因名字與原始矩陣匹配一下

dim(x3)

只剩9行，對應我們需要的9個基因

再將行給基因名字賦值給行名

rownames(x3)=x3$id    x3 <-x3[,-2] #刪除多餘的    x4 <- x3[,-1] #在刪除一列 賦值給x4      sum(table(colnames(x3)))

將細胞的名字全部取出來，變成數據框因為我喜歡處理數據框

w3<- data.frame(n=colnames(x3),             n2=rep(1,1021)) #建立相匹配的列    w3 <- w3[-1,] #刪除多餘  library(tidyr)  w4 <- w3 %>%    separate(n,into = c("n"),sep="\_")

這一步至關重要因為這一步將下劃線分隔符後的所有資訊全部刪除了

將剛剛的做好的細胞名字賦值給表達矩陣

colnames(x4) <- w4$n  x4<- x4[,-1] #刪除多餘

注意這時候，相對應的細胞名字我們還么有篩選

先畫個圖看看

library(pheatmap)    boxplot(x4[,1:10]) 和之前一樣

pheatmap(x4)

這裡我們注意到表達矩陣的值範圍非常大

xf<- log(x4+1)

老規矩我們取個log

boxplot(xf[,1:10])

pheatmap(xf)

好像好一點了，但是還需要進一步標準化

n=scale(xf)  class(n)  n<- as.data.frame(n)  n<- n[,xb2$n]  ##################這裡我們將獲取了需要的細胞種類  boxplot(n)

pheatmap(n)

這次我們就可以明顯看到賦值範圍縮小了。

最後調試一下

color.3 <- colorRampPalette(rev(c("#ff0000", "#ffffff",  "#0000ff")))(100)  pheatmap(n,scale = "row",show_rownames = T,display_numbers = F,           treeheight_row = F,treeheight_col = F,           fontsize_number=50,angle_col = "45",color = color.3,legend = TRUE,title )

由於pheatmap 似乎無法將輕易的將lab位置進行轉換。大家可以用ggplot 試一下。

原圖

是不是很類似？

生信技能樹目前已經公開了三個生信知識庫，記得來關注哦~

每周文獻分享

https://www.yuque.com/biotrainee/weeklypaper

腫瘤外顯子分析指南

https://www.yuque.com/biotrainee/wes

生物統計從理論到實踐

https://www.yuque.com/biotrainee/biostat