CNS文章代码学习（一）Immunity 三级淋巴结构

发布于：2022-11-08 ⋅ 阅读:(821) ⋅ 点赞:(0)

 

目录

一：样本情况

二：GEO数据处理

三：代码学习

①00_visium_preprocessing.R

②01_MCP_counter.R

第一步 把组织学的注释（用10X的Loup Browser)添加到metadata中

第二步 绘图

第三步 是一个整合 批量绘制所有的图？

③TLS_imprint_signature

原理：Leave-one-out

原理:ROC & AUC

原理：将数据输入ROCit包

 文献代码

一：样本情况

 主要学习对空间转录组的样本处理

二：GEO数据处理

Visium数据的读取若想使用标准读取

格式必须与下面完全一致才可以用标准读取

Seurat::Load10X_Spatial("GSM5924042_frozen_a_1_fuvben/outs")

outs文件夹下应该有以下内容：

spatial文件下(spatial文件夹不可以省略）

注意所有的名字必须一样，不可以带有样本编号等信息，否则难以读取；GEO下载的数据多不标准，需要自己整理成spatial Ranger处理后的标准结果

三：代码学习

①00_visium_preprocessing.R

整理所有样本到指定位置

files <- list.files() #设置路径，然后展示所有文件
[1] "code"                                                  "GSE175540_batch_corrected_bulk_RNA_seq_TPM_25_02_22.csv"
 [3] "GSM5924030_ffpe_c_2"                                     "GSM5924031_ffpe_c_3"                                    
 [5] "GSM5924032_ffpe_c_4"                                     "GSM5924033_ffpe_c_7"                                    
 [7] "GSM5924034_ffpe_c_10"                                    "GSM5924035_ffpe_c_20"                                   
 [9] "GSM5924036_ffpe_c_21"                                    "GSM5924037_ffpe_c_34"                                   
[11] "GSM5924038_ffpe_c_36"                                    "GSM5924039_ffpe_c_39"                                   
[13] "GSM5924040_ffpe_c_45"                                    "GSM5924041_ffpe_c_51"                                   
[15] "GSM5924042_frozen_a_1"                                   "GSM5924043_frozen_a_3"                                  
[17] "GSM5924044_frozen_a_15"                                  "GSM5924045_frozen_a_17"                                 
[19] "GSM5924046_frozen_b_1"                                   "GSM5924047_frozen_b_7"                                  
[21] "GSM5924048_frozen_b_13"                                  "GSM5924049_frozen_b_18"                                 
[23] "GSM5924050_frozen_c_2"                                   "GSM5924051_frozen_c_5"                                  
[25] "GSM5924052_frozen_c_23"                                  "GSM5924053_frozen_c_57"                                 
[27] "RCC_immunity_Tertiary_lymphoid.Rproj" 

slide_list <-files[3:26] #全部空间样本 24个
slide_list_test <- slide_list[1:3] #取随意3个测试代码

正式代码 附带拆解学习

system.time({  #确定运行时间
spatial_list <- sapply(slide_list_test,function(slide){  #一个sapply函数，对slide_list_test中的每一个元素应用自创函数slide
  print(slide)
  raw_data_directory <- paste0("/data/ybk/GEO/STvisium/RCC_immunity_Tertiary_lymphoid/",slide,"/outs/")
  spatial_object <- Seurat::Load10X_Spatial(raw_data_directory) #标准读取
  
  # Collect all genes coded on the mitochondrial genome 找到线粒体基因再计算其比例
  mt.genes <- grep(pattern = "^MT-", x = rownames(spatial_object), value = TRUE) 
  spatial_object$percent.mito <- (Matrix::colSums(spatial_object@assays$Spatial@counts[mt.genes, ])/Matrix::colSums(spatial_object@assays$Spatial@counts))*100
  
  #remove mt genes 移除线粒体基因
  genes_to_keep <- setdiff(names(which(Matrix::rowSums(spatial_object@assays$Spatial@counts )>5)),mt.genes)
  
  #QC去掉些spots和gene
  spatial_object_subset <- subset(spatial_object,features =genes_to_keep, subset = nFeature_Spatial > 300 & percent.mito < 30)
  cat("Spots removed: ", ncol(spatial_object) - ncol(spatial_object_subset), "\n")
  cat("Genes kept: ", length(genes_to_keep),"from",nrow(spatial_object), "\n") 
  
  #使用SCT的标准化（详见Seurat)
  spatial_object_subset <- SCTransform(spatial_object_subset, 
                                       assay = "Spatial", 
                                       verbose = T)
  
  #compute mcp scores计算MCP的score,使用SCT后的数据
  mcp_scores <- MCPcounter.estimate(spatial_object_subset[["SCT"]]@data,featuresType = "HUGO_symbols")
  rownames(mcp_scores) <- make.names(rownames(mcp_scores))#把MCPcounter的type名中间的空格去掉
  cell_types <- rownames(mcp_scores) #就是各种T Bcell
  for(cell_type in cell_types){   #for循环写入meta信息
    spatial_object_subset <- AddMetaData(object= spatial_object_subset,metadata = mcp_scores[cell_type,],col.name = cell_type)
  }
  # Unsupervised analysis 
  spatial_object_subset <- RunICA(spatial_object_subset, assay = "SCT", verbose = FALSE) #这里用ica降维；Seurat官方用的是pca
  spatial_object_subset <- FindNeighbors(spatial_object_subset, reduction = "ica")
  spatial_object_subset <- FindClusters(spatial_object_subset, verbose = FALSE,resolution = 0.4)
  spatial_object_subset <- RunUMAP(spatial_object_subset, reduction = "ica", dims = 1:30)
  return(spatial_object_subset)
  })    #“{”是function的；“）”是sapply的
}) #system.time

# user  system elapsed 可以看到3个样本也运行了十分钟
# 341.571  15.858 625.920

查看一下结果

spatial_list
# $GSM5924030_ffpe_c_2
# An object of class Seurat 
# 32170 features across 4482 samples within 2 assays 
# Active assay: SCT (16085 features, 3000 variable features)
# 1 other assay present: Spatial
# 2 dimensional reductions calculated: ica, umap
# 1 image present: slice1
# 
# $GSM5924031_ffpe_c_3
# An object of class Seurat 
# 31116 features across 4754 samples within 2 assays 
# Active assay: SCT (15558 features, 3000 variable features)
# 1 other assay present: Spatial
# 2 dimensional reductions calculated: ica, umap
# 1 image present: slice1
# 
# $GSM5924032_ffpe_c_4
# An object of class Seurat 
# 31516 features across 3829 samples within 2 assays 
# Active assay: SCT (15758 features, 3000 variable features)
# 1 other assay present: Spatial
# 2 dimensional reductions calculated: ica, umap
# 1 image present: slice1

 meta信息

 Idents信息

 原code中给的保存代码，蛮有意思的，第一次见到save.image()

save.image(paste0(res_folder,Sys.Date(),"_","seurat_processed.RData"))

这里其实就是保存一下这些spatial对象啦

②01_MCP_counter.R

上文已经计算出来了MCPcounter的得分，这里把它们出图

CNS级别的出图

第一步 把组织学的注释（用10X的Loup Browser)添加到metadata中

我的每一个样本的文件夹下储存相应的名为TLS_annotation.csv的注释文件

#ADD annotation 
#再次使用slide_list_test替代slide_list
slide_list_test
# [1] "GSM5924030_ffpe_c_2" "GSM5924031_ffpe_c_3"
# [3] "GSM5924032_ffpe_c_4"

for(slide in slide_list_test){
  #correct ident
  spatial_list[[slide]]$orig.ident <- slide #把orig.ident改过来
  #adding all annotation available
  annot <- paste0(paste0("/data/ybk/GEO/STvisium/RCC_immunity_Tertiary_lymphoid/",slide),"/TLS_annotation.csv")
  for(x in annot){
    annot_table <- read.csv(x)
    rownames(annot_table) <- annot_table[,1]
    annot_table$Barcode <- NULL
    spatial_list[[slide]] <- AddMetaData(object= spatial_list[[slide]],
                                         metadata = annot_table,
                                         col.name = paste0(colnames(annot_table),"_annot"))
  }
}
# 经过与对应文件中的csv验证，这些样本的ano已经被准确添加到对应样本的metadata中
spatial_list$GSM5924030_ffpe_c_2@meta.data["AAACGCTGGGCACGAC-1",]
spatial_list$GSM5924031_ffpe_c_3@meta.data["AAAGTGTGATTTATCT-1",]

                            orig.ident nCount_Spatial
AAAGTGTGATTTATCT-1 GSM5924031_ffpe_c_3           7535
                   nFeature_Spatial percent.mito
AAAGTGTGATTTATCT-1             3589            0
                   nCount_SCT nFeature_SCT    T.cells
AAAGTGTGATTTATCT-1      15618         4105 0.04951051
                   CD8.T.cells Cytotoxic.lymphocytes
AAAGTGTGATTTATCT-1           0             0.3549867
                   B.lineage NK.cells Monocytic.lineage
AAAGTGTGATTTATCT-1 0.6940119        0         0.6839274
                   Myeloid.dendritic.cells Neutrophils
AAAGTGTGATTTATCT-1                       0   0.2408643
                   Endothelial.cells Fibroblasts
AAAGTGTGATTTATCT-1         0.1194506    2.314476
                   SCT_snn_res.0.4 seurat_clusters
AAAGTGTGATTTATCT-1               0               0
                   TLS_2_cat_annot
AAAGTGTGATTTATCT-1             TLS

第二步 绘图

取消坐标轴

#hide spatial axes for seurat plots 
hide_axis <- theme(axis.title.x=element_blank(),
                   axis.text.x=element_blank(),
                   axis.ticks.x=element_blank(),
                   axis.title.y=element_blank(),
                   axis.text.y=element_blank(),
                   axis.ticks.y=element_blank())

绘制，这里挑选了一个样本来做

#figure 1 main 
GSM5924030_ffpe_c_2 <- spatial_list$GSM5924030_ffpe_c_2
Idents(GSM5924030_ffpe_c_2)  <-as.numeric(Idents(GSM5924030_ffpe_c_2)) #如果不用as.numeric，结果是一个因子
#H&E
he_GSM5924030_ffpe_c_2 <- SpatialPlot(GSM5924030_ffpe_c_2,repel = F,label = F,image.alpha=1,alpha = c(0,0), pt.size.factor =0.000001) + geom_point(alpha=0)+
  NoLegend() +
  ggtitle("TLS positive")+  
  theme(text=element_text(size=14))+ 
  theme(text=element_text(face = "bold"));he_GSM5924030_ffpe_c_2  
# MCP scores
cell_types <- make.names(rownames(MCPcounter.estimate(MCPcounterExampleData)))
plot_list_mcp <- lapply(cell_types[c(5,1,2,4,6,8,9,10)],function(x){
  plot_list_mcp <- SpatialPlot(GSM5924030_ffpe_c_2,
                               features = x,
                               image.alpha=0,
                               pt.size.factor = 1.5) +
    scale_fill_viridis_c(option="C") + 
    DarkTheme() +
    hide_axis +    
    theme(text=element_text(size=14))+ 
    theme(text=element_text(face = "bold"))+
    theme(legend.text=element_text(size=7))
  
})
plot_list_mcp

#CA9
ca9 <- SpatialPlot(GSM5924030_ffpe_c_2,
                   features = "CA9",
                   image.alpha=0,
                   pt.size.factor = 1.5) +
  scale_fill_viridis_c(option="C") + 
  DarkTheme() +
  hide_axis +    
  theme(text=element_text(size=14))+ 
  theme(text=element_text(face = "bold"))+
  theme(legend.text=element_text(size=7))

plot_list_mcp[[9]] <- ca9
plot_list_mcp[[10]] <- he_GSM5924030_ffpe_c_2
lay <- rbind(c(10,9,1,2,3,4),
             c(10,9,5,6,7,8))
png(paste0("results_ybk/",Sys.Date(),"_","FIGURE1_.png"),width = 450,height=250,res = 300,units = "mm",bg="transparent")
grid.arrange(grobs = plot_list_mcp, layout_matrix = lay)
dev.off()

出图的效果是所有的图集中在一个list也就是plot_list_mcp当中，grid.arrange函数排图

好看的

第三步 是一个整合 批量绘制所有的图？

我的代码应该还是要把spatial_list换为spatial_list_test

下面为作者源代码 我就没有再测试啦 先用到自己的数据上看看

###这一步应该是对上面函数的整合
sapply(names(spatial_list),function(id){
  print(id)
  #figure 1 main 
  spatial_obj <- spatial_list[[id]]
  Idents(spatial_obj)  <-as.numeric(Idents(spatial_obj)) 
  #H&E
  he <- SpatialPlot(spatial_obj,repel = F,label = F,image.alpha=1,alpha = c(0,0), pt.size.factor =0.00000) +
    NoLegend() +
    ggtitle(id)+  
    theme(text=element_text(size=10))+ 
    theme(text=element_text(face = "bold"))
  #MCP scores
  mcp_keep <- intersect(cell_types[c(5,1,2,4,6,8,9,10)],colnames(spatial_obj@meta.data))
  plot_list_mcp <- lapply(mcp_keep,function(x){
    plot_list_mcp <- SpatialPlot(spatial_obj,
                                 features = x,
                                 image.alpha=0,
                                 pt.size.factor = 1.8) +
      scale_fill_viridis_c(option="C") + 
      DarkTheme() +
      hide_axis +    
      theme(text=element_text(size=10))+ 
      theme(text=element_text(face = "bold"))+ 
      theme(legend.text=element_text(size=8))
  })
  #CA9
  if("CA9" %in% rownames(spatial_obj)){
    ca9 <- SpatialPlot(spatial_obj,
                       features = "CA9",
                       image.alpha=0,
                       pt.size.factor = 1.8) +
      scale_fill_viridis_c(option="C") + 
      DarkTheme() +
      hide_axis +    
      theme(text=element_text(size=10))+ 
      theme(text=element_text(face = "bold"))+
      theme(legend.text=element_text(size=8))
  }else{
    ca9 <- SpatialPlot(spatial_obj,repel = F,label = F,image.alpha=0,alpha = c(0,0), pt.size.factor =0.00000) + NoLegend() + ggtitle("CA9 not expressed")
  }
  plot_list_mcp[[9]] <- ca9
  plot_list_mcp[[10]] <- he
  
  test_null <- sapply(plot_list_mcp,is.null)
  if(any(test_null)){
    check_null <- which(test_null)
    plot_list_mcp[[check_null]] <- SpatialPlot(spatial_obj,repel = F,label = F,image.alpha=0,alpha = c(0,0), pt.size.factor =0.00000) + NoLegend() 
  }
  lay <- rbind(c(10,9,1,2,3,4),
               c(10,9,5,6,7,8))
  grid.arrange(grobs = plot_list_mcp, layout_matrix = lay)
})
dev.off()

③TLS_imprint_signature

原理：Leave-one-out

文中探索这个signature的方法叫做 Leave-one-out，应该说是交叉验证（a method that performs sequential permutations by identifying genes from 3 tumors and validating them on the fourth
one.）。挑选了4个有TLS的冰冻数据样本，计算他们TLS区（手动注释）和非TLS区spot的DEG(使用Seurat的findmarkers函数）。

 

 

随后拿三个样本的DEG，和第四个样本的信息做ROC曲线，并计算AUC（使用ROCit包），这样总共可计算四个样本的箱状图。只有那些至少在三个样本中，AUC都大于0.7的DEG才能被保留作为signature，最后剩下29个。随后还把这29个gene放在剩余有TLS的样本中进行验证，发现除了个别样本（质量差），所有的gene AUC都保持在0.7以上。

 

原理:ROC & AUC

1.ROC曲线的横坐标是伪阳性率FPR（也叫假正类率，False Positive Rate），纵坐标是真阳性率TPR（真正类率，True Positive Rate）。还有伪阴性率（FNR）和伪阴性率（FNR），暂不讨论它俩。

（1）伪阳性率（FPR）：判定为正例却不是真正例的概率，即真负例中判为正例的概率（比如一个人本来没有乙肝，结果试剂盒错把这个人判断为乙肝）

（2）真阳性率（TPR）：判定为正例也是真正例的概率，即真正例中判为正例的概率（也即正例召回率）（比如一个人本来有乙肝，也被试剂盒检测出来了乙肝）

2.列表大概就是这么个情况，想把TP变成TPR要经过计算. roc曲线就是根据TPR和FPR的坐标点绘制出来的。

 

3.ROC曲线用来评价一种分类器的好坏。如果它的AUC超过0.5（一半），就证明分类器比随机猜要好。曲线上的点代表一个个阈值，随着阈值的不同，假阳性和真阳性的概率也会改变。比如乙肝试剂盒的检测阈值越小，很多低携带量的患者就可以被成功检出，真阳性率升高，但也会有正常人被误伤为乙肝阳性，那么假阳性率也会变高。

 

原理：将数据输入ROCit包

文中的ROC包是ROCit，了解一下这个包的输入

 可以看到重要的就是score和class两个参数

文章code中可以看到score是一个样本的（差异）基因*（所有）细胞的表达矩阵；class是在提取样本TLS手动注释的metadata信息，其中中组织学注释为 “TLS” 的细胞会返回True

rocit(
score = spatial_list[[x]][["SCT"]]@data[gene,],    
#取出来一个行位基因列为barcode的SCT表达矩阵
class =  spatial_list[[x]]$TLS_2_cat_annot=="TLS", 
#组织学注释为TLS
negref = F
)

 这里对class举个例子，大概就是这样一个东西。是不是就是所谓的class of the obervations？

结果是一个逻辑值构成的向量，向量的名字是细胞名

 

这里参考一下SPSS在计算AUC时的数据输入：需要一个class列，一个value列。这样在设定不同的阈值时，就会根据value来分配到合适的class。比如，当以Value=0.6作为阈值，那么被试1,2,3,4被分类成阳性样本，其他被试被分类成阴性样本，据此，我们可以计算得到TPR=3/4，FPR=1/4。这样，我们就可以得到一组（TPR,FPR）值，制定多个阈值，我们就可以得到多组（TPR,FPR）值，把这些值（TPR,FPR）绘制出来得到的曲线就是ROC曲线。

SPSS需要这样整理数据：

回到ROCit的输入

假如score是这样一个矩阵

cell1 cell2 cell3
DEG1 12 6 4
DEG2 10 8 2
DEG3 6 5 1

那么class就是:(注意这是个向量，上面的cell是向量的name)

cell1 cell2 cell3
True False False

对于DEG1：

阈值取12，cell1判断为TLS，cell2 cell3非TLS，和class结果一致，那么TPR就是1，FPR是0；

阈值取6，cell1 cell2判断为TLS，cell3非TLS，有一个假阳性，那么TPR就是1/(1+0)，FPR是1/(1+1)；

阈值取4，cell1 cell2 cell3判断为TLS，有两个假阳性，那么TPR就是1/(1+0)，FPR是2/2；

例子有点极端，把AUC搞成1了，但也就说明通过DEG1的值来判断是否为TLS是非常好的；因为高表达的cell1是True，而低表达的cell2 cell3都是False。

以此类推，可以做出每一个差异基因的AUC。

如下图每个肿瘤有多个DEG，每个DEG都有自己的AUC，就形成一个小黑点。

那么leave one out的验证是如何实现的呢？

我猜想应该是用三个样本的DEG高表达基因取交集，得到一个基因list（猜的）。然后用这些gene在第四个样本中的表达量组成的score和第四个样本的TLS注释情况做AUC分析，这样就可以绘出第四个样本的gene-TLS AUC得分箱状图（如下）。

 文献代码

接下来来看看作者的代码究竟是怎么搞的

本文含有隐藏内容，请 开通VIP 后查看

网站公告

今日签到

点亮在社区的每一天
去签到

热门文章

最新发布