学徒带你一步步从CCLE数据库里面根据指定基因在指定细胞系里面提取表达矩阵进行热图可视化
昨天生信技能树发布了学徒作业：学徒作业-在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 试一下。
原图
是不是很类似？
生信技能树目前已经公开了三个生信知识库，记得来关注哦~
友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
•生信技能树的2019年终总结，你的生物信息学成长宝藏•2020学习主旋律，B站74小时免费教学视频为你领路•全国巡讲全球听（买一得五），你的生物信息学入门课