这里是佳奥，我们进入到单细胞转录组的下游分析。
1 vCAF和mCAF细胞亚群特有基因——结合TGCA
作者拿到了四组分类，其中重点关注vCAF和mCAF组，而组内的差异基因便进一步放到TGCA数据库绘制热图。
下载表达矩阵，选出文章基因，计算相关性
library(stringr)
##导入分组
vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
vCAF=unlist(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,', ') )
##UCSC_xena数据库
https://xenabrowser.net/datapages/
##选择选择GDC TCGA Breast Cancer (BRCA)
##选择gene expression RNAseq第一个
##点击download后面的网址下载(设置好路径)
if(F){
library(data.table)
# 文件BRCA.htseq_counts.tsv.gz从UCSC的XENA数据库下载，大于100M所以不提供在这里。
# a=fread('TCGA_BRCA/UCSC_xena/TCGA-BRCA.htseq_counts.tsv.gz',data.table=F)
a[1:4,1:4]
# Ensembl_ID TCGA-3C-AAAU-01A TCGA-3C-AALI-01A TCGA-3C-AALJ-01A
# 1 ENSG00000000003.13 9.348728 8.714246 10.356452
# 2 ENSG00000000005.5 1.584963 1.584963 5.727920
# 3 ENSG00000000419.11 10.874981 10.834471 10.329796
# 4 ENSG00000000457.12 10.121534 11.512247 8.867279
library(org.Hs.eg.db)
library(stringr)
esid=str_split(a$Ensembl_ID,
'[.]',simplify = T)[,1]
columns(org.Hs.eg.db)
head(esid)
rownames(a)=esid
e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" , "SYMBOL" ),keytype = 'ENSEMBL')
##修改大小写
vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL]
mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL]
##取表达矩阵需要的基因
ng=e2s[match(c(vCAF,mCAF),e2s$SYMBOL),1]
mat=a[ng,]
mat=mat[,-1]
save(mat,file = 'TCGA-BRCA-vCAF-and-mCAF-genes-expression.Rdata')
}
load(file = 'TCGA-BRCA-vCAF-and-mCAF-genes-expression.Rdata')
mat[1:4,1:4]
M=cor(t(mat))
colnames(M)=c(vCAF,mCAF)
rownames(M)=c(vCAF,mCAF)
pheatmap::pheatmap(M)
pheatmap::pheatmap(M,cluster_rows = F,cluster_cols = F)
n=t(scale(t( M )))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap::pheatmap(n,cluster_rows = F,cluster_cols = F)
QQ截图20220904151708.png
2 任意基因集在癌症表现——结合TGCA
library(stringr)
vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
vCAF=unlist(str_split(vCAF,', ')[[1]])
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,', ')[[1]])
ECM=c('COL1A1', 'COL1A2','COL3A1')
# https://www.genecards.org/cgi-bin/carddisp.pl?gene=DIPK2B
# CXorf36 --> DIPK2B
endothelial=c('CDH5', 'DIPK2B','TIE1')
if(F){
library(data.table)
# 文件BRCA.htseq_counts.tsv.gz从UCSC的XENA数据库下载，大于100M所以不提供在这里。
# a=fread('../../../TCGA_BRCA/UCSC_xena/TCGA-BRCA.htseq_counts.tsv.gz',data.table=F)
a[1:4,1:4]
# Ensembl_ID TCGA-3C-AAAU-01A TCGA-3C-AALI-01A TCGA-3C-AALJ-01A
# 1 ENSG00000000003.13 9.348728 8.714246 10.356452
# 2 ENSG00000000005.5 1.584963 1.584963 5.727920
# 3 ENSG00000000419.11 10.874981 10.834471 10.329796
# 4 ENSG00000000457.12 10.121534 11.512247 8.867279
library(org.Hs.eg.db)
library(stringr)
esid=str_split(a$Ensembl_ID,
'[.]',simplify = T)[,1]
columns(org.Hs.eg.db)
rownames(a)=esid
a=a[,-1]
e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" , "SYMBOL" ),keytype = 'ENSEMBL')
vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL,]
mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL,]
ng=e2s[match(vCAF,e2s$SYMBOL),1]
vCAF_value=colMeans(a[ng,])
ng=e2s[match(mCAF,e2s$SYMBOL),1]
mCAF_value=colMeans(a[ng,])
##多拿的表达矩阵
ng=e2s[match(ECM,e2s$SYMBOL),1]
ECM_value=colMeans(a[ng,])
ng=e2s[match(endothelial,e2s$SYMBOL),1]
endothelial_value=colMeans(a[ng,])
dat=data.frame(vCAF_value=vCAF_value,
mCAF_value=mCAF_value,
ECM_value=ECM_value,
endothelial_value=endothelial_value )
save(dat,file = 'step2.Rdata')
}
load(file = 'step2.Rdata')
## 两种方法批量绘制散点图（多图组合）
## https://mp.weixin.qq.com/s/SG1VvfL8lVCniGCjW2DqfA
## 画boxplot
library(ggpubr)
colnames(dat)
ggscatter(dat, x = "vCAF_value", y = "endothelial_value",
color = 'black', shape = 21, size = 0.5, # Points color, shape and size
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE,
cor.coeff.args = list(method = "pearson", label.sep = "\n")
)
QQ截图20220904152838.png
QQ截图20220904152901.png
3 多个基因集相关性热图
library(stringr)
vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
vCAF=unlist(str_split(vCAF,', ')[[1]])
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,', ')[[1]])
ECM=c('COL1A1', 'COL1A2','COL3A1')
# https://www.genecards.org/cgi-bin/carddisp.pl?gene=DIPK2B
# CXorf36 --> DIPK2B
endothelial=c('CDH5', 'DIPK2B','TIE1')
# for 参考文献ref32：proliferation signature
proliferation <- c("BAG1","ESR1","FOXA1","GPR160","NAT1","NAT1","MAPT","MLPH","PGR")
# for ref31: stroma-related signatures
stroma <- 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')
# for ref29: endothelial/microvasculature signatures
microvasculature <- 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')
if(F){
library(data.table)
# 文件BRCA.htseq_counts.tsv.gz从UCSC的XENA数据库下载，大于100M所以不提供在这里。
# a=fread('../../../TCGA_BRCA/UCSC_xena/TCGA-BRCA.htseq_counts.tsv.gz',data.table=F)
a[1:4,1:4]
# Ensembl_ID TCGA-3C-AAAU-01A TCGA-3C-AALI-01A TCGA-3C-AALJ-01A
# 1 ENSG00000000003.13 9.348728 8.714246 10.356452
# 2 ENSG00000000005.5 1.584963 1.584963 5.727920
# 3 ENSG00000000419.11 10.874981 10.834471 10.329796
# 4 ENSG00000000457.12 10.121534 11.512247 8.867279
library(org.Hs.eg.db)
library(stringr)
esid=str_split(a$Ensembl_ID,
'[.]',simplify = T)[,1]
columns(org.Hs.eg.db)
rownames(a)=esid
a=a[,-1]
e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" , "SYMBOL" ),keytype = 'ENSEMBL')
vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL]
mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL]
ng=e2s[match(vCAF,e2s$SYMBOL),1]
vCAF_value=colMeans(a[ng,])
ng=e2s[match(mCAF,e2s$SYMBOL),1]
mCAF_value=colMeans(a[ng,])
ng=e2s[match(ECM,e2s$SYMBOL),1]
ECM_value=colMeans(a[ng,])
ng=e2s[match(endothelial,e2s$SYMBOL),1];ng
endothelial_value=colMeans(a[ng,])
ng=e2s[match(proliferation,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
proliferation_value=colMeans(a[ng,])
ng=e2s[match(stroma,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
stroma_value=colMeans(a[ng,])
ng=e2s[match(microvasculature,e2s$SYMBOL),1];ng=ng[!is.na(ng)];ng
microvasculature_value=colMeans(a[ng,])
dat=data.frame(vCAF=vCAF_value,
mCAF=mCAF_value,
ECM=ECM_value,
endothelial=endothelial_value,
microvasculature=microvasculature_value,
stroma=stroma_value,
proliferation=proliferation_value)
save(dat,file = 'step3.Rdata')
}
load(file = 'step3.Rdata')
dat[1:4,1:4]
M=cor( dat)
pheatmap::pheatmap(M)
QQ截图20220904153244.png
4 数据库探索
这是课程的额外部分，这里只提供相关代码的下载。
探索TGCA——BRCA
https://github.com/jmzeng1314/TCGA_BRCA
探索MEMTABRIC数据集
https://github.com/jmzeng1314/METABRIC
探索GTEX数据库
https://github.com/jmzeng1314/gtex_BRCA
5 展望
后续关于10X数据资源：
https://www.10xgenomics.com/cn
或者关注生信技能树的单细胞天地。
这里是佳奥，我们后续的学习再见！
