单细胞转录组学习笔记-21-基因在任意癌症表达量相关性

刘小泽写于19.9.6-第四单元第一讲：计算基因在任意癌症表达量相关性
笔记目的：根据生信技能树的单细胞转录组课程探索smart-seq2技术相关的分析技术
课程链接在：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)(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个数据集：

在Xena数据库搜索lung adenocarcinoma 会有3个数据集：

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,242评论 5赞 459
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 81,769评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 141,484评论 0赞 319
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,133评论 1赞 263
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,007评论 4赞 355
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,080评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,496评论 3赞 381
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,190评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,464评论 1赞 290
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,549评论 2赞 309
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,330评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,205评论 3赞 312
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,567评论 3赞 298
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,889评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,160评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,475评论 2赞 341
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,650评论 2赞 335