查看原文
其他

评估任意基因集在癌症的表现

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


课程笔记




粉丝:有单细胞线上课程吗?

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

好了,戏演完了,下面郑重介绍下我们的单细胞线上课程:(详情戳下方链接) 

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


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

希望大家能有所收获!


目录



前言

第四单元第二讲:评估任意基因集在癌症的表现
课程链接在:http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

上一篇是探索两个细胞亚群(vCAF、mCAF)特有的基因在TCGA中的表现,发现两个亚群的基因都是和TCGA
相关的基因在内部相关,说明了分群的效果不错

目的就是做下面这个图的相关性分析:

可以看到,横坐标的vCAF就是我们前一篇得到的vCAF基因集在TCGA数据集中的表达量,那么纵坐标,就需要去文章里找,作者是拿到了5篇不同参考文献的6个数据集

文章正文放了四张相关性的图,是vCAF和mCAF与第27篇参考文献中的两个乳腺癌数据集进行的比较

然后再来看看第27篇文献的图,其中列出了乳腺癌的ECM和Endothelial的基因集

然后这篇参考文献的作者定义基因集的方法就是:在大部分癌症中都存在的基因就是基因集,因此我们看到,本文使用的基因集中就3个基因

看到其中有一个奇怪的基因名CXorf36,按说人类的基因名都应该是大写,所以拿到GenCard查询一下:https://www.genecards.org/cgi-bin/carddisp.pl?gene=DIPK2B&keywords=CXorf36。发现这个基因目前叫DIPK2B

有了基因集,就去获取4个基因集的表达量

各个基因如下:

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')

还是使用上一篇的GDC乳腺癌TCGA的表达矩阵(60,489 identifiers X 1217 samples

需要注意的是:上一篇我们单纯比较多个基因相关性,所以得到多个基因的表达量然后做个热图就好;但是这次要比较的是两个基因集(每一组内都有不同数量的基因)。作者用散点图来展现,其中的每一个点实际上就是一个样本,但是同一个样本在两个基因集中对应的基因数量不同,不能简单拿任何一个基因进行比较。作者给的方法是:

但是这里我们只是简单对每个样本取个均值

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
19ng=e2s[match(vCAF,e2s$SYMBOL),1]
20vCAF_value=colMeans(a[ng,])
21ng=e2s[match(mCAF,e2s$SYMBOL),1]
22mCAF_value=colMeans(a[ng,])
23ng=e2s[match(ECM,e2s$SYMBOL),1]
24ECM_value=colMeans(a[ng,])
25ng=e2s[match(endothelial,e2s$SYMBOL),1]
26endothelial_value=colMeans(a[ng,])
27
28dat=data.frame(vCAF_value=vCAF_value,
29                 mCAF_value=mCAF_value,
30                 ECM_value=ECM_value,
31                 endothelial_value=endothelial_value )
最后绘制散点图
1library(ggpubr)
2colnames(dat)
3ggscatter(dat, x = "vCAF_value", y = "endothelial_value",
4          color = 'black', shape = 21, size = 0.5# Points color, shape and size
5          add = "reg.line",  # Add regressin line
6          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
7          conf.int = TRUE# Add confidence interval
8          cor.coef = TRUE,  
9          cor.coeff.args = list(method = "pearson",  label.sep = "\n")
10)





seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比

使用seurat3的merge功能整合8个10X单细胞转录组样本

单细胞转录组逐渐沦为科研标配技术

单细胞核测序在人类肾脏上的应用

单细胞数据该怎么利用?

单细胞去除聚类的离群点

如何去学一个R包(下)

如何去学一个R包(上)

是否可以根据10X转录组数据来推断基因组CNV信息呢?

你以为细胞聚在一起就是一类细胞吗

如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程

生信技能树(爆款入门培训课)全国巡讲约你

(南京、南宁见!)全国巡讲第17-18站(生信入门课加量不加价)



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存