前言

说明：之前做的复现有一点问题，这一版已经是更改过的。放心食用。

期刊：Cancer
标题：Single-Cell T tanscriptomic Analysis of T umor-Derived Fibroblasts and Normal Tissue-Resident Fibroblasts Reveals Fibroblast Heterogeneity in Breast Cancer。

文章链接

摘要：主要对小鼠乳腺癌CAFs和正常小鼠成纤维细胞进行单细胞测序，并对CAFs的异质性进行了描述。并未涉及机制研究，算是一篇纯的不行的生信文章。适合入门。
废话不多说，开办。

# 一、数据下载文章的数据一把藏在dataability里面，但是这篇文章的数据在方法中。除了部分supplementary中的图没有数据外，正文数据都可以下载。下载解压后是txt文件。

数据下载

二、操作

1.看文章

当然是先看文章了，然后根据文章的方法设置参数，当然你也可以自己设置参数，用文章的参数是为了检验自己有没有别的地方做错了。

2.读入数据，过滤

这里文献提供的是counts矩阵。（部分文章提供的是稀疏矩阵，需要read10x读取）
代码如下：

rm(list = ls())
options(stringsAsFactors = F)
library(dplyr)
library(Seurat)
library(patchwork)
library(R.utils)
library(textshape)
library(monocle)
##读取CSV文件，如果是三文件的10x数据可直接创建seurat对象
a <- read.table(file = '4T1_viable_raw_counts.txt',header = T,sep = '\t',fill = T)

2.Seurat流程

seurat标准流程包括：（复杂的还应该包括多个数据集的整合）
创建seurat对象，
过滤，
（在此我没有对线粒体基因进行过滤，但是整体效果与文章大致相同。没做线粒体基因过滤的原因呢，主要是基因名字里面没有找到’Mt-‘开头的基因，如果要过滤，可以网上搜索小鼠线粒体基因，本人较懒，所以，哈哈）
normalize，
findvariablefeatures，
scale，
runpca，
findneighbors，
findclusters
可视化
代码如下：

# 一、Seurat流程 -------------------------------------------------------------------
sc <- CreateSeuratObject(counts = a,project = 'mouse_viable',min.cells = 5)#1.过滤 #这里'Mt'开头的基因并非全是线粒体基因
sc[['percent.mt']] <- PercentageFeatureSet(sc,pattern = '^Mt')VlnPlot(sc,features = c('percent.mt','nFeature_RNA','nCount_RNA'),group.by = 'orig.ident')
sc <- subset(sc,subset=nFeature_RNA<=7800&nFeature_RNA>500)#2.标准化、高变基因、Scale
sc <- NormalizeData(sc)
sc <- FindVariableFeatures(sc,selection.method = 'vst')
sc <- ScaleData(sc,features = rownames(sc)) #这里对所有的基因进行scale，是方便画热图，否则部分基因画不出#3.PCA线性降维
sc <- RunPCA(sc,npcs = 30)
ElbowPlot(sc,ndims = 30,reduction = 'pca')#4.聚类
sc <- FindNeighbors(sc,dims = 1:15)
sc <- FindClusters(sc,resolution = 0.3)#5.非线性降维及可视化
sc <- RunUMAP(sc,dims = 1:15)
sc <- RunTSNE(object = sc,dims = 1:15)
DimPlot(sc,reduction = 'umap',label = T)
DimPlot(sc,reduction = 'tsne',label = T)

文章的图

复现的图：可以看见，除了各类的位置不一样，类的形状，大小基本与原图一致

有的人就会问了，你敢不敢再过分一点啊？
敢，那我们来把文中的细胞比例算一下！！

3.计算细胞比例

文中讲到免疫细胞66.4%，上皮细胞24.5%，CAFs 8%。那我不安排一下？

代码如下：

bl <- table(sc@meta.data$seurat_clusters)/sum(table(sc@meta.data$seurat_clusters))
bl <- as.array(bl)
##文章中的图为左上角的7类比例为66.4%（和我的cluster名字有点区别）
#我的图中是0，2，3，5，7，8，9,比例为0.6641745
sum(bl[c('0','2','3','5','7','8','9')])
##文章中图左下角1，6群比例为24.5%，我的也是左下角1，6群
##我的结果是0.2451713
sum(bl[c('1','6')])
##CAFs的比例,此处我的类是"3",文中比例8%，我的结果：0.08302181
bl['4']

结果基本（完全）一摸一样，如果你非要说不一样，那只能说我的比他的要精确

4.基因差异

4.1.气泡图

文章的图我就不放了，放一下自己做的图吧。和文章的图是一模一样啊

代码如下：

gens <- c('Ptprc','Itgam','Cd14','Arg2','Mrc1','Fcgr2b','Irf7','Nusap1','Mki67','S100a8','Ly6g','Cd79a','Cd19','Nkg7','Cd3d','Thy1','Pdpn','Pdgfra','Epcam','Pecam1','Mcam','Pdgfrb','Cspg4','Rgs5')
DotPlot(object = sc,features = gens,group.by = 'seurat_clusters')+RotatedAxis()#顺便把E图的气泡图一做
dotgene <- c('Bmp1','Tgfbr2','Tgfbr3','Il11ra1','Il33','Cxcl14','Cxcl12','Cxcl1','C4b','C3','C1s1','C1ra','Sfrp4','Sfrp2','Sfrp1','Has1','Chl1','Cdh11')
DotPlot(sc,features = dotgene)+RotatedAxis()

4.2.热图

哇，心态小崩，这么多基因，纯手敲的，文章中竟然不能复制图中的文字。。。
文中的图应该是用pheatmap做的，所以颜色及排版显得不同。有兴趣的也可以去尝试。主要是提取scale.data中的矩阵进行。

heatgene <- c('Col1a1','Col1a2','Col3a1','Col4a1','Col4a2','Col5a1','Col5a2','Col5a3','Col6a1','Col6a2','Col6a3','Col8a1','Col12a1','Col14a1','Col15a1','Col16a1','Adamts2','Ctsk','Lox','Loxl1','Loxl2','Loxl3','Mmp2','Mmp3','P3h1','P3h3','P3h4','P4ha1','P4ha2','P4ha3','P4hb','Plod1','Plod2','Plod3','Bgn','Lum','Ogn','Prelp','Prg4','Dcn','Aspn','Vcan','Spon1','Postn','Tnc','Cthrc1','Thbs2','Fbln1','Fbn2','Spon2','Fbn1','Dpt','Fndc1','Cyr61','Nid1','Fbln2')
DoHeatmap(sc,features = heatgene)

有的人又会问了，你敢不敢还过分一点啊？
敢，且听下回分解。

总结

关于文中的注释问题，可以自己根据marker基因进行手动修改。这个需要一定的知识背景。
其实这篇文章的分析到这才只是开始，后面还有对CAFs的专门分析，假时序分析，我会继续更新，今天就先复现到这儿吧。

目前正在学这方面的知识，欢迎大家和我一起进步。
水平有限，如有错误，请批评指正。

单细胞转录组文章复现系列(一)——seurat相关推荐

seurat提取表达矩阵_Hemberg-lab单细胞转录组数据分析
单细胞RNA-seq简介混合RNA-seq2000年末的重大技术突破,取代微阵列表达芯片被广泛使用通过混合大量细胞获取足够RNA用于建库测序,来定量每个基因的平均表达水平用于比较转录组,例如比较 ...
Seurat 单细胞转录组测序数据分析教程（二）——python(scanpy)
Seurat 单细胞转录组测序数据分析教程(二)--python(scanpy) 文章参考至scanpy官网,做了一个更详细的解读. 数据由来自健康捐赠者的 3k PBMC组成,可从 10x Geno ...
代码分析 | 单细胞转录组质控详解
前言 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅综述:三万字 ...
单细胞转录组数据整合分析专题研讨会(2019.11)
2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...
有了易生信，导师再也不用担心我的单细胞转录组整合分析啦
2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...
单细胞转录组专题研讨会第二期
单细胞转录组之前是跟常规转录组一起开课的,但后来因为涉及内容多,也需要更专业的讲解,8月份第一次单飞独自开课,邀请中科院单细胞分析算法开发博士倾力授课,一鸣惊人,取得了很好的效果. 现于2019年11 ...
单细胞转录组单飞第二期开课啦！！
单细胞转录组之前是跟常规转录组一起开课的,但后来因为涉及内容多,也需要更专业的讲解,8月份第一次单飞独自开课,邀请中科院单细胞分析算法开发博士倾力授课,一鸣惊人,取得了很好的效果. 现于2019年11 ...
易生信群体和单细胞转录组专题第6期于5月10日在北京开课了
群体转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,分析思路简洁清晰,是入门生信,解决最常见问题的首选. 单细胞分析是近几年的明星技术,多次被Nature.Science评为 ...
中科院单细胞分析算法开发博士带你做单细胞转录组分析
" 福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组和Python课程的线上直播课.报名参加线上直播课的老师可在1年内选择参加同课程的一次 ...

单细胞转录组文章复现系列(一)——seurat

单细胞转录组复现（1）

前言