做完jimmy老师最优秀的学徒,我居然好久没有做复现任务了,好多函数都忘了。。 趁此机会(2021-07暑期时间)好好复习一下,fighting!
数据来源[1]于细胞学万能老二的期刊《nature cell biology》的一篇文章《Transcriptional diversity and bioenergetic shift in human breast cancer metastasis revealed by single-cell RNA sequencing》[2],提到了乳腺癌转移过程中的mitochondrial oxidative phosphorylation (OXPHOS) 通路激活问题。
同时作者上传了自己的code到github https://github.com/lawsonlab/Single_Cell_Metastasis/blob/master/Single_Cell_Classifier.R,文章总共有1707个细胞,经过筛选后剩下1119个细胞:
作者这里用的是seurat 2.1,但是现在更新到4了=-=
数据下载预处理
作者上传的数据已经是FPKM处理好的,可以看到是3个样品,这里我直接下的raw文件进行合并操作
因为有1707个文件,实在太多=--=,先去吃饭了
结果报错
有2个问题:
1.发现不应该用rbind,应该用cbind=-= 浪费了好多时间2.有一个样本最后一行没有2行,作者传失误了?删去就好啦
文中是1119,有可能是上一步id转换已经粗暴的去掉了=-=
所以需要查看gencode的gtf,查看核糖体基因
结果还是不对,原文是3个数据集分开处理的。明天我再试试分开,结果是否一致=-=
这里不写循环的原因是因为分组的编号不一样,并且作者已经做好基因注释
H1
H1还剩下246个cells
H2
结果h2的读取,遇到了日期基因=-=,作者传的文件还是挺多坑的,需要去除掉这些日期基因
可以看到,与h1相比,少了28个基因,因为这28个因为作者的疏忽变成日期基因了
h2分组的时候可以看到有12个cells命名跟之前的明明不一致,这个只能回到gse的网页去查看所属分组
发现是M,所以全部算作M
H2还剩下396个细胞
H10
同样需要去查看命名
剩下470个cells
小结
总共三个样本最后所剩cell数为1112,与原文1119也是有所差异,与自己直接进行三个样本放在一起id转换也是有所差异。所以个人觉得原因如下
1.作者自己上传H2数据混入了日期基因,导致H2少了28个基因的表达,所以自己分样本结果比原文少2.不能粗暴的基因id转换,因为会去除很多基因,例如线粒体基因=-= 作者把60155个基因全部进行了注释,所以以后做单细胞,最好是通过gtf基因注释
下面是合并3个sc,这里与文章有所不同,因为文章用的是seurat2,现在3可以直接多数据通过anchor合并了,所以结果也许也会不同=--=
这里作者用cell cycle作为batch effects,所以需要去除影响
通过作者给的阈值 logfc>0.25|p<0.05来进行筛选top15基因
火山图
热图
可以看到与作者的结果有所差异
top15基因
作者说直接到网站下载的生存数据but可能付费用户才能下吧,我实在是找不到,所以我还是到metabric下载的数据(这里偷个懒,用Jimmy老师github的数据=-=)
最后有394个basal-like样本,与原文的360个样本接近
后面将top15与转移相关的基因作为基因集来判断生存分析
需要自己制作 brca.gmt 文件哦,就是前面得到的基因集。
生存确实具有显著性,不同的原因可能是
1.作者用的RFS,metabric只有OS,但是都具有显著性2.基因signature有所差异3.作者应该不是用的GSVA来算的,因为没有引用,code里面也没有提=-=,所以我也不知道作者具体用的什么方法来算基因集
[1]
数据来源: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123837
[2]
《Transcriptional diversity and bioenergetic shift in human breast cancer metastasis revealed by single-cell RNA sequencing》: https://doi.org/10.1038/S41556-020-0477-0