在前面的教程：居然有如此多种癌症（是时候开启pan-cancer数据挖掘模式），我们把全部的TCGA的33种癌症的表达量矩阵区拆分成为蛋白编码基因和非编码基因这两个不同的表达量矩阵，并且保存成为了rdata文件。
在不同癌症里面，蛋白质编码相关基因数量一直在一万八附近，而非编码基因数量跨度比较大，从一万二到两万七不等。因为非编码基因的组织病人特异性都很强，恰好最近看到了一个很有意思图可以说明这件事，如下所示：
这个图来源于文献：《Transcriptional landscape and clinical utility of enhancer RNAs for eRNA-targeted therapy in cance》，作者展现的是Enhancer RNA (eRNA) 的癌症特异性很强，其实eRNA就是非编码基因的一种。根据特异性可以分成成为3种：
- 652 ubiquitous eRNAs expressed in ≥10 cancer types,
- 3124 intermediately specific eRNAs that are expressed in 2–9 cancer types,
- 5332 cancer- type-specific eRNAs that are expressed in only one cancer type
接下来我们就针对格式化好的非编码基因的表达量矩阵，进行如下所示的图表绘制！
step1: 批量读取前面的Rdata
代码语言：javascript复制rm(list=ls())
options(stringsAsFactors = F)
library(stringr)
library(data.table)
###### step1: 批量读取前面的Rdata ######
fs=list.files('Rdata/',
pattern = 'htseq_counts')
fs
# 首先加载每个癌症的表达量矩阵（提取非蛋白编码基因部分）
# 成为一个 list
dat_list = lapply(fs, function(x){
# x=fs[1]
pro=gsub('.htseq_counts..Rdata','',x)
print(pro)
load(file = file.path('Rdata/',x))
dat=log2(edgeR::cpm(non_mat) 1)
dat[1:4,1:4]
return(dat)
})
step2: 区分基因
代码语言：javascript复制
###### step2: 区分基因 ######
tmp=table(unlist(lapply(dat_list, rownames)))
all_genes=names(tmp) #[tmp==33]
length(all_genes)
ubiquitous_genes=names(tmp)[tmp==33]
length(ubiquitous_genes)
specific_genes =names(tmp)[tmp==1]
length(specific_genes)
intermediately_genes=setdiff(setdiff(all_genes,specific_genes),ubiquitous_genes)
# 在33个癌症都存在的非编码基因不到九千个
# 但是33个癌症总共涉及到31455个非编码基因
# 独特存在于33种癌症的仅仅是一个里面的是 2236个基因。
在33个癌症都存在的非编码基因不到九千个，但是33个癌症总共涉及到31455个非编码基因。独特存在于33种癌症的仅仅是一个里面的是 2236个基因。
step3: 计算各个基因在各个癌症的均值
代码语言：javascript复制
# 但是突然间发现，下面的代码难度更高了：
sm = do.call(cbind,
lapply(dat_list, function(x){
m=rowMeans(x)
ubiquitous_m=rep(0,length(ubiquitous_genes))
specific_m=rep(0,length(specific_genes))
intermediately_m=rep(0,length(intermediately_genes))
cg=intersect(names(m),ubiquitous_genes)
ubiquitous_m[match(cg,ubiquitous_genes)]=m[match(cg,names(m))]
cg=intersect(names(m),specific_genes)
specific_m[match(cg,specific_genes)]=m[match(cg,names(m))]
cg=intersect(names(m),intermediately_genes)
intermediately_m[match(cg,intermediately_genes)]=m[match(cg,names(m))]
return(c(ubiquitous_m,intermediately_m,specific_m))
}))
colnames(sm)=gsub('TCGA-','',gsub('.htseq_counts..Rdata','',fs))
rownames(sm)=c(ubiquitous_genes,intermediately_genes,specific_genes)
pheatmap::pheatmap(sm)
出丑图如下所示：
勉强可以看到 (ubiquitous_genes,intermediately_genes,specific_genes ，下面我们进行初步美化看看：
step4: 美化
代码语言：javascript复制
sm_bak=sm
sm=sm_bak
sm[sm>6]=6
df1=sm[ubiquitous_genes,]
df1=df1[order(rowMeans(df1),decreasing = T),]
p1=pheatmap(t(df1),
cluster_rows = F,
cluster_cols = F,
show_rownames = T,
show_colnames = F)
p1
df2=sm[intermediately_genes,]
n_genes = apply(df2, 1, function(x) sum(x>0))
tmp_list= split(names(n_genes),n_genes)
lapply(tmp_list, length)
ordd=do.call(c,
lapply(tmp_list, function(cg){
# cg=tmp_list[[1]]
ord=do.call(c,lapply(dat_list, function(x){
cg[cg %in% rownames(x)]
}))
length(ord)
return(unique(ord))
}) )
length(ordd)
df2=df2[order( n_genes ,ordd,rowMeans(df2),decreasing = T),]
p2=pheatmap(t(df2),
cluster_rows = F,
cluster_cols = F,
show_rownames = T,
show_colnames = F)
p2
df3=sm[specific_genes,]
ord=do.call(c,lapply(dat_list, function(x){
specific_genes[specific_genes %in% rownames(x)]
}))
p3=pheatmap(t(df3[ord,]),
cluster_rows = F,
cluster_cols = F,
show_rownames = T,
show_colnames = F)
p3
library(cowplot)
plot_grid(p1$gtable,
p2$gtable,
p3$gtable,nrow = 1,
labels=c('ubiquitous_genes',
'intermediately_genes',
'specific_genes') )
我也说不清楚这个是美化还是丑化了：
买家秀和卖家秀差距不是一点点啊！
同理可以应用于蛋白编码基因表达量矩阵
可以看到：
代码语言：javascript复制> length(all_genes)
[1] 19302
> ubiquitous_genes=names(tmp)[tmp==33]
> length(ubiquitous_genes)
[1] 16370
> specific_genes =names(tmp)[tmp==1]
> length(specific_genes)
[1] 123
> intermediately_genes=setdiff(setdiff(all_genes,specific_genes),ubiquitous_genes)
> length(intermediately_genes)
[1] 2809
几乎所有的基因都在全部的癌症内部有表达量，仅仅是123个基因具有癌症特异性，而且它的表达量肯定是超级低！
你能相信它们是蛋白编码基因吗？
这123个基因很有意思，或许是一个课题？
