多个基因集相关性热图

Original 刘小泽单细胞天地 2022-06-06

收录于合集 #第一期单细胞视频笔记汇总 21个

课程笔记

粉丝：有单细胞线上课程吗？

小编：什么？我们的单细胞转录组分析线上课程已经上线好久了，你们竟然都不知道吗，每篇推文后面的课程推荐没人看的吗，小编已哭晕在厕所

好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接）

全网第二个单细胞视频课程预售

这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获！

前言

第四单元第三讲：多个基因集相关性热图

课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
基于前面的两节，这一节变得更容易理解

文章正文还是以乳腺癌为例，绘制了6个数据集和两个亚型基因集之间的相关性

这次还是要对基因集进行操作，因此也是需要进行上一次的colMeans()操作，只不过上次总共得到4个，这次是7个（因为下文将troma1和stroma2合并为一个stroma）

获取各个基因集的基因

一共5篇参考文献，需要将其中各个基因集拿出来：

 1library(stringr)
 2# vCAF基因集
 3vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
 4vCAF=unlist(str_split(vCAF,', '))
 5# mCAF基因集
 6mCAF='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'
 7mCAF=unlist(str_split(mCAF,', '))
 8# ECM基因集
 9ECM=c('COL1A1', 'COL1A2','COL3A1')
10# Endothelial基因集
11endothelial=c('CDH5', 'DIPK2B','TIE1')
12# for 参考文献ref32：proliferation 基因集
13proliferation <- c("BAG1","ESR1","FOXA1","GPR160","NAT1","NAT1","MAPT","MLPH","PGR")
14# for ref31: stroma-related 基因集
15stroma <- c('DCN', 'CSPG2', 'CDH11', 'COL3A1', 'FAP', 'PEDF', 'FBN1', 'PDGFRL', 'CTSK', 'HTRA1', 'ASPN', 'SPARC', 'COL5A2', 'LOXL1', 'MMP2', 'SPON1', 'SFRP4', 'ITGBL1', 'CALD1', 'COPZ2', 'MFAP2', 'ANGPTL2', 'PLAU', 'COL1A2', 'LRRC17', 'C1QTNF3', 'SNAI2', 'PCOLCE', 'POSTN', 'ECM2', 'FBLN1', 'ADAM12', 'MMP11', 'AEBP1', 'PDGFRB', 'GAS1', 'COL6A3', 'RARRES2', 'COL6A1', 'TGFB3', 'NDN', 'C1R', 'LRP1', 'COL10A1', 'DPYSL3', 'OLFML2B', 'MMP14', 'DACT1', 'MGC3047', 'THBS2')
16# for ref29: endothelial/microvasculature 基因集
17microvasculature <- c('ARAP3','ADCY4','ESAM','ERG','SLC43A3','SOX7','PTPRB','PTPRM','AFAP1L1','MMRN2','TENC1','STARD9','COL4A3','LRRK1','PALD1','NPR3','ROBO4','NOTCH4','TIE1','RASIP1','ACVRL1','RAMP2','FAM110D','EGFL7','SMAD6','FGD5','ENG','CASKIN2','ACKR2','SLC9A3R2','CALCRL','HSPA12B','EPAS1','EHD4','LATS2','ICAM1','HBEGF','PLTP','C1orf54','CTTNBP2NL','MYO1B','SLCO2A1','KIFC1','EPHB4','SOX13','DRAM1','PECAM1','ENTPD1','ICAM2','CLDN5','SDPR','CDH5','GPR116','ELTD1','KDR','HILPDA','NPNT')

然后得到7个基因集的表达量

 1# 读入数据TCGA-BRCA.htseq_counts.tsv.gz
 2library(data.table)
 3filepath <- file.choose()
 4a=fread(filepath ,data.table=F)
 5
 6# Ensembl ID切割
 7library(stringr)
 8esid=str_split(a$Ensembl_ID,
 9                 '[.]',simplify = T)[,1]
10# ID转换
11e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" ,  "SYMBOL" ),keytype = 'ENSEMBL')
12vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL,]
13mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL,]
14
15# 获得表达量
16rownames(a)=esid
17a=a[,-1]
18# vCAF
19ng=e2s[match(vCAF,e2s$SYMBOL),1]
20vCAF_value=colMeans(a[ng,])
21# mCAF
22ng=e2s[match(mCAF,e2s$SYMBOL),1]
23mCAF_value=colMeans(a[ng,])
24# ECM
25ng=e2s[match(ECM,e2s$SYMBOL),1]
26ECM_value=colMeans(a[ng,])
27# Endothelial
28ng=e2s[match(endothelial,e2s$SYMBOL),1]
29endothelial_value=colMeans(a[ng,])
30# proliferation(注意这里加上了检测基因名是否存在的代码)
31ng=e2s[match(proliferation,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
32proliferation_value=colMeans(a[ng,])
33# stroma
34ng=e2s[match(stroma,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
35stroma_value=colMeans(a[ng,])
36# microvasculature
37ng=e2s[match(microvasculature,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
38microvasculature_value=colMeans(a[ng,])
39
40dat=data.frame(vCAF=vCAF_value,
41                 mCAF=mCAF_value,
42                 ECM=ECM_value,
43                 endothelial=endothelial_value,
44                 microvasculature=microvasculature_value,
45                 stroma=stroma_value,
46                 proliferation=proliferation_value)
47
48> dat[1:4,1:4] 
49                     vCAF     mCAF      ECM endothelial
50TCGA-A1-A0SP-01A 8.794181 13.07456 17.13670    9.990037
51TCGA-BH-A201-01A 9.549604 14.86951 19.40557   10.895579
52TCGA-E2-A14T-01A 9.600401 13.33230 17.61368   10.475556
53TCGA-AC-A8OS-01A 9.998759 14.07297 19.24798   11.932741