祖传的单个10x样本的seurat标准代码

其实单细胞领域进展太快，我那些课程内容关于R包相关的代码基本上过时了，因为R语言本身都经历了一个超级大的变革！考虑到不能把粉丝带歪，我早就全部公开了系列视频课程。还创立了《单细胞天地》这个公众号 ：

全网第一个单细胞课程（免费基础课程）

• 免费学习地址在B站：https://www.bilibili.com/video/av38741055 ，欢迎提问弹幕交流！
• 务必听课后完成结业考核20题：https://mp.weixin.qq.com/s/lpoHhZqi-_ASUaIfpnX96w
• 课程配套资料文档在：https://docs.qq.com/doc/DT2NwV0Fab3JBRUx0

技能树出品的第二个单细胞课程（进阶课程，仍然免费）

• 详情请自行阅读介绍 https://mp.weixin.qq.com/s/bLfO-8ri_SNUepGs4UwRQw
• 本课程长期答疑文档，https://docs.qq.com/doc/DT0FxbEpHYU5ZVlpu

因为课程涉及到知识点太多，所以我拆分成为了5个子课程，欢迎B站提问弹幕交流！全部链接是：

• 「生信技能树」单细胞进阶数据处理之文献导读，链接是：https://www.bilibili.com/video/BV17f4y1R7N8
• 「生信技能树」使用10X单细胞转录组数据探索免疫治疗，链接是：https://www.bilibili.com/video/BV1xD4y1S74P
• 「生信技能树」单细胞基因组数据拷贝数变异分析流程，链接是：https://www.bilibili.com/video/BV1Yf4y1R75R
• 「生信技能树」云服务器处理单细胞转录组数据，链接是：https://www.bilibili.com/video/BV154411Z7DU
• 「生信技能树」使用Smart-seq2单细胞转录组数据探索小鼠性腺发育，链接是：https://www.bilibili.com/video/BV1454y1q77Z

也希望可以帮助到你。

这里做一个统一的代码更新

复制粘贴就可以使用的代码哦，单个10x样本的seurat标准代码如下：

### ---------------###### Create: Jianming Zeng### Date: 2020-08-27 15:00:26### Email: [e-mail]### Blog: http://www.bio-info-trainee.com/### Forum: http://www.biotrainee.com/thread-1376-1-1.html### CAFS/SUSTC/Eli Lilly/University of Macau### Update Log: 2020-08-27 First version###### ---------------
rm(list=ls())options(stringsAsFactors = F)library(Seurat)pro='S1'# 搞清楚你的10x单细胞项目的cellranger输出文件夹哦hp_sce <- CreateSeuratObject(Read10X('scRNAseq_10_s1/filtered_feature_bc_matrix/'), pro)
hp_scerownames(hp_sce)[grepl('^mt-',rownames(hp_sce))]rownames(hp_sce)[grepl('^Rp[sl]',rownames(hp_sce))]
hp_sce[["percent.mt"]] <- PercentageFeatureSet(hp_sce, pattern = "^mt-")fivenum(hp_sce[["percent.mt"]][,1])rb.genes <- rownames(hp_sce)[grep("^Rp[sl]",rownames(hp_sce))]C<-GetAssayData(object = hp_sce, slot = "counts")percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100hp_sce <- AddMetaData(hp_sce, percent.ribo, col.name = "percent.ribo")
plot1 <- FeatureScatter(hp_sce, feature1 = "nCount_RNA", feature2 = "percent.mt")plot2 <- FeatureScatter(hp_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")CombinePlots(plots = list(plot1, plot2))
VlnPlot(hp_sce, features = c("percent.ribo", "percent.mt"), ncol = 2)VlnPlot(hp_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)VlnPlot(hp_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)hp_sce
hp_sce1 <- subset(hp_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)hp_sce1

sce=hp_sce1scesce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)GetAssay(sce,assay = "RNA")sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)# 步骤 ScaleData 的耗时取决于电脑系统配置（保守估计大于一分钟）sce <- ScaleData(sce)sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)ElbowPlot(sce)sce <- FindNeighbors(sce, dims = 1:15)sce <- FindClusters(sce, resolution = 0.2)table([e-mail]$RNA_snn_res.0.2)sce <- FindClusters(sce, resolution = 0.8)table([e-mail]$RNA_snn_res.0.8)

set.seed(123)sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)DimPlot(sce,reduction = "tsne",label=T)phe=data.frame(cell=rownames([e-mail]), res2=[e-mail]$RNA_snn_res.0.2, res8=[e-mail]$RNA_snn_res.0.8, cluster =[e-mail]$seurat_clusters)head(phe)table(phe$cluster)tsne_pos=Embeddings(sce,'tsne')DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')
head(phe)table(phe$cluster)head(tsne_pos)dat=cbind(tsne_pos,phe)save(dat,file=paste0(pro,'_for_tSNE.pos.Rdata'))load(file=paste0(pro,'_for_tSNE.pos.Rdata'))library(ggplot2)p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster), geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()print(p)theme= theme(panel.grid =element_blank()) + ## 删去网格 theme(panel.border = element_blank(),panel.background = element_blank()) + ## 删去外层边框 theme(axis.line = element_line(size=1, colour = "black"))p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))print(p)ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2.pdf'))
plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt")plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")CombinePlots(plots = list(plot1, plot2))ggplot2::ggsave(filename = paste0(pro,'_CombinePlots.pdf'))
VlnPlot(sce, features = c("percent.ribo", "percent.mt"), ncol = 2)ggplot2::ggsave(filename = paste0(pro,'_mt-and-ribo.pdf'))VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)ggplot2::ggsave(filename = paste0(pro,'_counts-and-feature.pdf'))VlnPlot(sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)

table([e-mail]$seurat_clusters)
for( i in unique([e-mail]$seurat_clusters) ){ markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25) print(x = head(markers_df)) markers_genes = rownames(head(x = markers_df, n = 5)) VlnPlot(object = sce, features =markers_genes,log =T ) ggsave(filename=paste0(pro,'_VlnPlot_subcluster_',i,'_sce.markers_heatmap.pdf')) FeaturePlot(object = sce, features=markers_genes ) ggsave(filename=paste0(pro,'_FeaturePlot_subcluster_',i,'_sce.markers_heatmap.pdf'))}
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
DT::datatable(sce.markers)write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))library(dplyr)top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)DoHeatmap(sce,top10$gene,size=3)ggsave(filename=paste0(pro,'_sce.markers_heatmap.pdf'))# FeaturePlot( sce, top10$gene )# ggsave(filename=paste0(pro,'_sce.markers_FeaturePlot.pdf'),height = 49)

library(SingleR)sce_for_SingleR <- GetAssayData(sce, slot="data")mouseImmu <- ImmGenData()pred.mouseImmu <- SingleR(test = sce_for_SingleR, ref = mouseImmu, labels = mouseImmu$label.main)
mouseRNA <- MouseRNAseqData()pred.mouseRNA <- SingleR(test = sce_for_SingleR, ref = mouseRNA, labels = mouseRNA$label.fine )
cellType=data.frame(seurat=[e-mail]$seurat_clusters, mouseImmu=pred.mouseImmu$labels, mouseRNA=pred.mouseRNA$labels)sort(table(cellType[,2]))sort(table(cellType[,3]))table(cellType[,2:3])save(sce,cellType, file=paste0(pro,'_output.Rdata') )
rm(list=ls())options(stringsAsFactors = F)library(Seurat)pro='S1'library(SingleR)load(file=paste0(pro,'_output.Rdata') )
cg=names(tail(sort(table(cellType[,2]))))cellname=cellType[,2]cellname[!cellname %in% cg ]='other'table([e-mail]$seurat_clusters,cellname)table(cellname)# Epithelial cells, 0,2,5,7,8,12# Endothelial cells, 10# Fibroblasts , 1,3,6,9,14,15# Macrophages , 4,13,16# NKT , 11# other ,4,11# Stromal cells,cl=[e-mail]$seurat_clustersps=ifelse(cl %in% c(0,2,5,7,8,12),'epi', ifelse(cl %in% c(1,3,6,9,14,15),'fibro', ifelse(cl %in% c(4,13,16),'macro', ifelse(cl %in% c(10),'endo', ifelse(cl %in% c(11),'NKT','other' )))))table(ps)cellname=paste(cl,ps)[e-mail]$cellname=cellnameDimPlot(sce,reduction = "tsne",label=T,group.by = 'cellname')ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_raw.pdf'))dat$cluster=cellnamelibrary(ggplot2)p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster), geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()print(p)theme= theme(panel.grid =element_blank()) + ## 删去网格 theme(panel.border = element_blank(),panel.background = element_blank()) + ## 删去外层边框 theme(axis.line = element_line(size=1, colour = "black"))p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))print(p)ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.2_singleR_pretty.pdf'))


这里面包括了我一直强调的单细胞基础10讲相关内容 ：

• 01. 上游分析流程
• 02.课题多少个样品，测序数据量如何
• 03. 过滤不合格细胞和基因（数据质控很重要）
• 04. 过滤线粒体核糖体基因
• 05. 去除细胞效应和基因效应
• 06.单细胞转录组数据的降维聚类分群
• 07.单细胞转录组数据处理之细胞亚群注释
• 08.把拿到的亚群进行更细致的分群
• 09.单细胞转录组数据处理之细胞亚群比例比较

还有一些个性化汇总。

如果你需要的是CNS文章全部图表，发表级别的代码，那就必须看CNS文章自带的GitHub哈，比如：你要的rmarkdown文献图表复现全套代码来了（单细胞）