查看原文
其他

CNS图表复现07—原来这篇文章有两个单细胞表达矩阵

生信技能树 单细胞天地 2022-06-06

分享是一种态度



回顾

我们的CNS图表复现之旅已经开始,前面6讲是;

如果你也想加入交流群,自己去:你要的rmarkdown文献图表复现全套代码来了(单细胞)找到我们的拉群小助手哈。

正文

文章的第一次分群按照 :

  • immune (CD45+,PTPRC),
  • epithelial/cancer (EpCAM+,EPCAM),
  • stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

的表达量分布,可以拿到如下所示的:

第一次分群

文章提到的各大亚群细胞数量是:(epithelial cells [n = 5,581], immune cells [n = 13,431], stromal cells [n = 4,249]).

我发现自己怎么都重复不出来,因为唯一的质量控制步骤,细胞数量就开始不吻合了:

  
  # 这里没有绝对的过滤标准
  # Filter cells so that remaining cells have nGenes >= 500 and nReads >= 50000
  # 26485 features across 27489 samples
  main_tiss_filtered <- subset(x=main_tiss, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
  main_tiss_filtered
  # 26485 features across 21444 samples

也就是说27489个细胞,被过滤成为了 21444细胞,而文章的细胞数量是23261,但是我明明是跟文章作者一模一样的处理过程。

真的是百思不得其解。

后来认真看了看它提供的文件列表:

 1.4G Aug 30 12:04 S01_datafinal.csv
   11M Aug 30 11:34 S01_metacells.csv
  3.5K Aug 30 11:35 by_sample_ratios_driver_muts_3.30.20.csv
  3.3K Aug 30 11:35 neo-osi_metadata.csv
  184M Aug 30 11:35 neo-osi_rawdata.csv
   13K Aug 30 11:34 samples_x_genes_3.30.20.csv

原来是有两个表达矩阵文件,neo-osi_rawdata.csv 和 S01_datafinal.csv 文件,都需要走流程~

需要注意的是这个文章作者并没有上传这些文件到GEO上面,仅仅是上传了原始数据在:https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA591860,这个原始数据接近4T的文件。

这些csv文件呢,作者上传到了谷歌云盘,我们人工搬运到了百度云,在交流群可以拿到。如果你也想加入交流群,自己去:你要的rmarkdown文献图表复现全套代码来了(单细胞)找到我们的拉群小助手哈。

我漏掉的这个 neo-osi_rawdata.csv 走流程拿到的是  2070 个单细胞

# 26485 features across 3507 samples within 1 assay 
  osi_object_filtered <- subset(x=osi_object, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
  osi_object_filtered 
  # 26485 features across 2070 samples within 1 assay 

然后把两个单细胞表达矩阵构建好的Seurat 对象合并起来,代码如下:

main_tiss_filtered
osi_object_filtered
main_tiss_filtered1 <- merge(x = main_tiss_filtered, y = osi_object_filtered)
main_tiss_filtered1 

得到的就差不多是文章里面的 23514 个单细胞啦,文章的细胞数量是23261。

> main_tiss_filtered1
An object of class Seurat 
26485 features across 23514 samples within 1 assay 
Active assay: RNA (26485 features, 0 variable features)

全部代码可以复制粘贴的

总共是  270 行,基本上涵盖了我们单细胞天地的方方面面知识点。

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
pro='first'

# 单细胞项目:来自于30个病人的49个组织样品,跨越3个治疗阶段

# Single-cell RNA sequencing (scRNA-seq) of
# metastatic lung cancer was performed using 49 clinical biopsies obtained from 30 patients before and during
# targeted therapy. 

# before initiating systemic targeted therapy (TKI naive [TN]), 
# at the residual disease (RD) state,
# acquired drug resistance (progression [PD]).

# 首先读取第一个单细胞转录组表达矩阵文件:
if(T){
  
  # Load data 
  dir='./'
  Sys.time()
  raw.data <- read.csv(paste(dir,"Data_input/csv_files/S01_datafinal.csv", sep=""), header=T, row.names = 1)
  Sys.time()
  dim(raw.data) 
  raw.data[1:4,1:4]
  head(colnames(raw.data)) 
  # Load metadata 
  metadata <- read.csv(paste(dir,"Data_input/csv_files/S01_metacells.csv", sep=""), row.names=1, header=T)
  head(metadata) 
  
  # Find ERCC's, compute the percent ERCC, and drop them from the raw data.
  
  erccs <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = TRUE)
  percent.ercc <- Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
  fivenum(percent.ercc)
  ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE)
  raw.data <- raw.data[-ercc.index,]
  dim(raw.data) 
  # Create the Seurat object with all the data (unfiltered)
  main_tiss <- CreateSeuratObject(counts = raw.data)
  # add rownames to metadta 
  row.names(metadata) <- metadata$cell_id
  # add metadata to Seurat object 
  main_tiss <- AddMetaData(object = main_tiss, metadata = metadata)
  main_tiss <- AddMetaData(object = main_tiss, percent.ercc, col.name = "percent.ercc")
  # Head to check
  head(main_tiss@meta.data)
  
  # Calculate percent ribosomal genes and add to metadata
  ribo.genes <- grep(pattern = "^RP[SL][[:digit:]]", x = rownames(x = main_tiss@assays$RNA@data), value = TRUE)
  percent.ribo <- Matrix::colSums(main_tiss@assays$RNA@counts[ribo.genes, ])/Matrix::colSums(main_tiss@assays$RNA@data)
  fivenum(percent.ribo)
  main_tiss <- AddMetaData(object = main_tiss, metadata = percent.ribo, col.name = "percent.ribo")
  main_tiss
  
  # 唯一的质量控制步骤:
  # 这里没有绝对的过滤标准
  # Filter cells so that remaining cells have nGenes >= 500 and nReads >= 50000
  # 26485 features across 27489 samples
  main_tiss_filtered <- subset(x=main_tiss, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
  main_tiss_filtered
  # 26485 features across 21444 samples
}
main_tiss_filtered
# 26485 features across 21444 samples

# 然后读取第二个单细胞转录组表达矩阵文件:
if(T){
  # load and clean rawdata
  # Load data 
  dir='./'
  Sys.time()
  osi.raw.data <- read.csv(paste(dir,"Data_input/csv_files/neo-osi_rawdata.csv", sep=""), row.names = 1)
  
  Sys.time()
 
  colnames(osi.raw.data) <- gsub("_S.*.homo""", colnames(osi.raw.data))
  osi.raw.data[1:4,1:4]
  tail( osi.raw.data[ 1:4])
  dim(osi.raw.data)
  
  # drop sequencing details from gene count table
  # 因为是HTseq这个软件跑出来的表达矩阵
  row.names(osi.raw.data)[grep("__", row.names(osi.raw.data))]
  osi.raw.data <- osi.raw.data[-grep("__", row.names(osi.raw.data)),]
  
  #Make osi.metadata by cell from osi.metadata by plate
  
  osi.metadata <- read.csv(paste(dir, "Data_input/csv_files/neo-osi_metadata.csv", sep = ""))
  osi.meta.cell <- as.data.frame(colnames(osi.raw.data))
  osi.meta.cell <- data.frame(do.call('rbind', strsplit(as.character(osi.meta.cell$`colnames(osi.raw.data)`),'_',fixed=TRUE)))
  rownames(osi.meta.cell) <- paste(osi.meta.cell$X1, osi.meta.cell$X2, sep = "_")
  colnames(osi.meta.cell) <- c("well""plate")
  osi.meta.cell$cell_id <- rownames(osi.meta.cell)
  library(dplyr)
  osi.metadata <- left_join(osi.meta.cell, osi.metadata, by = "plate")
  rownames(osi.metadata) <- osi.metadata$cell_id
  head(osi.metadata)
  
  unique(osi.metadata$plate)
  
  
  # Find ERCC's, compute the percent ERCC, and drop them from the raw data.
  
  erccs <- grep(pattern = "^ERCC-", x = rownames(x = osi.raw.data), value = TRUE)
  percent.ercc <- Matrix::colSums(osi.raw.data[erccs, ])/Matrix::colSums(osi.raw.data)
  ercc.index <- grep(pattern = "^ERCC-", x = rownames(x = osi.raw.data), value = FALSE)
  osi.raw.data <- osi.raw.data[-ercc.index,]
  dim(osi.raw.data)
  
  
  # Create the Seurat object with all the data (unfiltered)
  
  osi_object <- CreateSeuratObject(counts = osi.raw.data)
  osi_object <- AddMetaData(object = osi_object, metadata = osi.metadata)
  osi_object <- AddMetaData(object = osi_object, percent.ercc, col.name = "percent.ercc")
  # Changing nUMI column name to nReads
  colnames(osi_object@meta.data)[colnames(osi_object@meta.data) == 'nUMI'] <- 'nReads'
  head(osi_object@meta.data)
  
  
  # Calculate percent ribosomal genes and add to osi.metadata
  
  ribo.genes <- grep(pattern = "^RP[SL][[:digit:]]", x = rownames(x = osi_object@assays$RNA@data), value = TRUE)
  percent.ribo <- Matrix::colSums(osi_object@assays$RNA@data[ribo.genes, ])/Matrix::colSums(osi_object@assays$RNA@data)
  osi_object <- AddMetaData(object = osi_object, metadata = percent.ribo, col.name = "percent.ribo")
  osi_object
  # 26485 features across 3507 samples within 1 assay 
  osi_object_filtered <- subset(x=osi_object, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
  osi_object_filtered 
  # 26485 features across 2070 samples within 1 assay 
  VlnPlot(osi_object_filtered, features = "nFeature_RNA")
  VlnPlot(osi_object_filtered, features = "nCount_RNA", log = TRUE)

}
main_tiss_filtered
osi_object_filtered
main_tiss_filtered1 <- merge(x = main_tiss_filtered, y = osi_object_filtered)
main_tiss_filtered1 
raw_sce <- main_tiss_filtered1
if(F){
  # Drop any samples with 10 or less cells
  
  main_tiss_filtered@meta.data$sample_name <- as.character(main_tiss_filtered@meta.data$sample_name)
  sample_name <- as.character(main_tiss_filtered@meta.data$sample_name)
  # Make table 
  tab.1 <- table(main_tiss_filtered@meta.data$sample_name) 
  # Which samples have less than 10 cells 
  samples.keep <- names(which(tab.1 > 10))
  metadata_keep <- filter(main_tiss_filtered@meta.data, sample_name %in% samples.keep)
  # Subset Seurat object 
  tiss_subset <- subset(main_tiss_filtered, cells=as.character(metadata_keep$cell_id))
  tiss_subset
  raw_sce <- tiss_subset
}

# Gene-expression profiles of 23,261 cells were retained after quality control filtering.



raw_sce 
rownames(raw_sce)[grepl('^mt-',rownames(raw_sce),ignore.case = T)]
rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]

raw_sce[["percent.mt"]] <- PercentageFeatureSet(raw_sce, pattern = "^MT-")
fivenum(raw_sce[["percent.mt"]][,1])
rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce),ignore.case = T)]
C<-GetAssayData(object = raw_sce, slot = "counts")
percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
fivenum(percent.ribo)
raw_sce <- AddMetaData(raw_sce, percent.ribo, col.name = "percent.ribo")

plot1 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "percent.ercc")
plot2 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

VlnPlot(raw_sce, features = c("percent.ribo""percent.ercc"), ncol = 2)
VlnPlot(raw_sce, features = c("nFeature_RNA""nCount_RNA"), ncol = 2)
VlnPlot(raw_sce, features = c("percent.ribo""nCount_RNA"), ncol = 2)
raw_sce
 
sce=raw_sce
sce
#sce <- NormalizeData(sce, normalization.method =  "LogNormalize",  scale.factor = 1e6)
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",  scale.factor = 1e5)
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)
# 跟文章保持一致,第一次分群,resolution 采取 0.5
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
table(sce@meta.data$RNA_snn_res.0.2) 
sce <- FindClusters(sce, resolution = 0.8)
table(sce@meta.data$RNA_snn_res.0.8) 
sce <- FindClusters(sce, resolution = 0.5)
table(sce@meta.data$RNA_snn_res.0.5) 

# 不同的 resolution 分群的结果不一样
library(gplots)
tab.1=table(sce@meta.data$RNA_snn_res.0.2,sce@meta.data$RNA_snn_res.0.8) 
balloonplot(tab.1)

# Check clustering stability at given resolution  
# Set different resolutions 
res.used <- seq(0.1,1,by=0.2)
res.used
# Loop over and perform clustering of different resolutions 
for(i in res.used){
  sce <- FindClusters(object = sce, verbose = T, resolution = res.used)
}
# Make plot 
library(clustree)
clus.tree.out <- clustree(sce) +
  theme(legend.position = "bottom") + 
  scale_color_brewer(palette = "Set1") +
  scale_edge_color_continuous(low = "grey80", high = "red")

clus.tree.out 

res.used <- 0.5
sce <- FindClusters(object = sce, verbose = T, resolution = res.used)


set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
phe=data.frame(cell=rownames(sce@meta.data),
               cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne'
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.pdf'))

sce <- RunUMAP(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "umap",label=T)
ggplot2::ggsave(filename = paste0(pro,'_umap.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'))
sce.first=sce
save(sce.first,sce.markers,file = 'first_sce.Rdata')

如果你看不懂代码,请自行回顾:单细胞基础10讲

更多个性化汇总代码见:单细胞初级8讲和高级分析8讲

初级分析

高级分析



如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存