万字长文带你手把手完成一个单细胞转录组文献 图表复现（GSE123837）

做完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与转移相关的基因作为基因集来判断生存分析

GSVA

需要自己制作 brca.gmt 文件哦，就是前面得到的基因集。

生存确实具有显著性，不同的原因可能是

1.作者用的RFS，metabric只有OS，但是都具有显著性2.基因signature有所差异3.作者应该不是用的GSVA来算的，因为没有引用，code里面也没有提=-=，所以我也不知道作者具体用的什么方法来算基因集

References

[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