学徒带你一步步从CCLE数据库里面根据指定基因在指定细胞系里面提取表达矩阵进行热图可视化

生信菜鸟团

发布于 2020-03-30 14:27:10

4.1K0

发布于 2020-03-30 14:27:10

文章被收录于专栏：生信菜鸟团生信菜鸟团

昨天生信技能树发布了学徒作业：学徒作业-在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 试一下。