单细胞转录组复现（1）
前言二、操作1.看文章2.读入数据，过滤2.Seurat流程3.计算细胞比例4.基因差异4.1.气泡图4.2.热图总结前言
说明：之前做的复现有一点问题，这一版已经是更改过的。放心食用。期刊：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.标准化、高变基因、Scalesc <- 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.6641745sum(bl[c('0','2','3','5','7','8','9')])##文章中图左下角1，6群比例为24.5%，我的也是左下角1，6群##我的结果是0.2451713sum(bl[c('1','6')])##CAFs的比例,此处我的类是"3",文中比例8%，我的结果：0.08302181bl['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的专门分析，假时序分析，我会继续更新，今天就先复现到这儿吧。
目前正在学这方面的知识，欢迎大家和我一起进步。
水平有限，如有错误，请批评指正。
