目录
第一步 把组织学的注释（用10X的Loup Browser)添加到metadata中
④展示TLS_imprint_signature和其它marker关系
一：样本情况
主要学习对空间转录组的样本处理
二：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表达矩阵（该函数结合了sapply函数，gene为一个单独的基因
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 |
那么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在第四个样本中的基因取交集，获得最后用来做ROC的DEG。对每一个DEG取出score，和第四个样本的TLS注释情况class一起做AUC分析。这样就可以绘出第四个样本的gene-TLS AUC得分箱状图（如下）。一个点代表一个基因，拥有一个AUC值
文献代码
接下来来看看作者的代码究竟是怎么搞的
作者的源代码是"b_1","a_15","a_4","b_17"四个样本，我这里就接上文的三个样本
最后的效果就是两个样本的DEG与第三个样本的表达矩阵gene取交集获得score，并用第三个样本的TLS注释(class)来做ROC曲线。
Differential expression in each samples
# Differential expression in each samples 用举例的三个样本
TLS_pos_ids <- c("GSM5924030_ffpe_c_2",
"GSM5924031_ffpe_c_3",
"GSM5924032_ffpe_c_4")
# 差异基因分析
TLS_list_all_samples <- lapply(TLS_pos_ids,function(y){
x <- spatial_list[[y]]
to_use <- "TLS"
if(length(to_use)!=0){
tls_MAST <- FindMarkers(x,group.by = "TLS_2_cat_annot",ident.1 = to_use,test.use = "MAST")
}else{
return(NA)
}
})
TLS_list_all_samples #三个差异基因的List，对应上方的三个样本Perform discovery on 2 tumors and validate on 1（原作为3 tumors and validate on 1)
hide_ig <- F #这是给一个逻辑值
x <- TLS_pos_ids[c(1)] #给下文X的值举例 正文不需要此代码 跑循环的时候不运行
cross_validation <- sapply(TLS_pos_ids[c(3,2,1)],function(x){
print(x) #x就是对应的3个样本之一
sub <- TLS_pos_ids[TLS_pos_ids!=x] #sub是除了x的另外两个样本
gene_list <- lapply(TLS_list_all_samples[TLS_pos_ids!=x],function(y){ #y是除了x的另外两个样本的差异基因list，也就是sub对应的样本差异基因
rownames(y)[which(y$avg_log2FC > 1 & y$p_val_adj < 0.05)] #按照阈值筛选差异基因
}) ##lapply end 返回的结果是一个list，list中是两个样本的差异基因
tls_sign <- unique(unlist(gene_list)) #把上面筛选的两个样本的所有的差异基因去重，得到筛选后的DEG
intersect_tls_sign <- intersect(tls_sign,rownames(spatial_list[[x]])) #将DEG与x的所有基因取交集（保证上文推导的差异基因都在X中出现）
if(hide_ig){ #默认不运行本段
intersect_tls_sign <- grep(x =intersect_tls_sign, pattern = "IG",value = T,invert = T) #invert=T以后返回的是未匹配到的值
intersect_tls_sign <- grep(x =intersect_tls_sign, pattern = "JCHAIN",value = T,invert = T) #value=T后返回的是匹配到的值，F则返回位置
}
#ROC
#intersect_tls_sign是通过除X的另外几个样本得到的DEG与样本X的所有基因取交集后得到的一串基因名（signature）
rocs <- sapply(intersect_tls_sign,function(gene){ #d对signature中的所有基因挨个应用函数function(gene)
roc_empirical <- rocit(score = spatial_list[[x]][["SCT"]]@data[gene,], #取出来一个：行为一个基因，列为barcode的SCT表达矩阵
class = spatial_list[[x]]$TLS_2_cat_annot=="TLS", #组织学注释为TLS
negref = F)
roc_empirical$AUC
}) #输出的结果rocs是一串ROC值的向量，向量的name是基因
return(list(
"tls_sign"=intersect_tls_sign,
"auc"=sort(rocs)))
}) ###supply end
运行结果解读：
直接来看是个List，class的结果却不是
class(cross_validation)
# [1] "matrix" "array" 
tls_sign是利用该样本进行验证时的TLS signature，auc展示对应基因的值
如果把上文的cross_validation <- sapply(TLS_pos_ids[c(3,2,1)],function(x){改为lapply,则返回结果是list。展示为标准的List结构
class(cross_validation)
# [1] "list"
如果注释掉一些代码
则尽管使用的是sapply，返回的结果仍然是list
class(cross_validation)
# [1] "list" 
既注释掉return，又使用lapply
class(cross_validation)
# [1] "list"
上述的例子涉及Return和sapply和lapply的区分
在return中，如果你自己编了一个函数，试图使这个函数返回多个值，例如返回均值，方差，如果在自定义的函数里，分两行写了均值，方差，将只返回方差。如果要实现这个目的，请使用list。
返回值return:
返回值就是函数给出的结果。打个比方，编写一个函数，就像自己攒一个机器，例如现在攒好 一台豆浆机，该豆浆机要求输入大豆，输入的大豆就是参数， 返回的结果，就是豆浆。如果该豆浆机需要不停地输入大豆， 而不能产出豆浆，这样的机器就一定会被扔掉。函数也是一样的， 需要给出返回值。 R中默认的情况是将最后一句作为返回值。但是为了函数的可读性起见，应该尽量指名返回值。返回值用return()函数给出。 函数在内部处理过程中，一旦遇到return()，就会终止运行， 将return()内的数据作为函数处理的结果给出。
结合上述四种情况可以看出，在没有设置返回值时，只能默认返回最后一个运行值（rocs)。
想要同时返回两个值，请用list连接
不想返回一个列表，请使用sapply
Get list with threshold 按照0.7的阈值筛选
tls_sign <- sort(table(unlist(sapply(cross_validation["auc",],function(auc) names(which(auc > 0.7))))))
tls_sign <- names(which(tls_sign >=3))
tls_sign
# [1] "IGHA1" "IGHG1" "IGHG2" "IGHG3" "IGHM" "IGKC" "IGKV4-1" "IGLC1" "JCHAIN" "POU2AF1"Venn diagramm 图有点丑，不过可以出图
tls_sign_over_07 <- sapply(cross_validation["auc",],function(auc) auc[which(auc > 0.7)])
library(VennDiagram)
toto <- sapply(tls_sign_over_07,function(x) names(x)) #结果是所有经过0.7筛选的基因名
pl <- venn.diagram(
x = toto,
category.names = names(toto),
filename = NULL,cex=4,
fill=c(viridis(3)),
output=T
)
pl
png(paste0(res_folder,Sys.Date(),"_FIGURE_TLS_signature_supp_venn.png"),width = 7,height=5,units = "in",res = 500,bg = 'transparent')
grid.draw(pl)
dev.off()
注意这里
y1=sapply(tls_sign_over_07,function(x) names(x))
y2=sapply(tls_sign_over_07,function(x){
names(x)}) #两种写法等价对上述代码进行可视化
#visualise on all tumours 在所有的肿瘤上对TLS——signature可视化
opt<-"C" #提供调色板选项
#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())
val <-lapply(names(spatial_list)[1:3],function(x){ #我的list中只有三个样本，源代码1:12
print(x)
intersect_tls_sign <- intersect(tls_sign,rownames(spatial_list[[x]]))
spatial_list[[x]] <- AddMetaData(spatial_list[[x]],apply(as.matrix(spatial_list[[x]][["SCT"]]@data[intersect_tls_sign,]),2,mean),col.name = "TLS_signature")
#ROC
if(any(spatial_list[[x]]$TLS_2_cat_annot=="TLS")){
roc_empirical <- rocit(score = spatial_list[[x]]$TLS_signature,
class = spatial_list[[x]]$TLS_2_cat_annot=="TLS",
negref = F)
auc <- roc_empirical$AUC
}else{
auc <- NA
}
SpatialPlot(spatial_list[[x]],image.alpha = 0, features = "TLS_signature") +
scale_fill_viridis_c(option=opt,limits = c(0,4.1)) +
DarkTheme() +
hide_axis +
ggtitle(label =paste0("AUC = ", round(auc,2)))
})
grid.arrange(grobs=val,ncol=3)
png(paste0(res_folder,Sys.Date(),"_FIGURE_TLS_signature_supp_heatmap.png"),width = 24,height=10,units = "in",res = 500,bg = 'transparent')
grid.arrange(grobs=val,ncol=6)
dev.off()
绘制差异基因火山图
（还是只有三个样本
#figure supp TLS signature 差异基因的火山图
#plot volcano
names(TLS_list_all_samples) <- TLS_pos_ids
all_volcanoes <- lapply(TLS_pos_ids[c(3,2,1)], function(ids){
diff_exp <- TLS_list_all_samples[[ids]]
diff_exp <- diff_exp[,c("avg_log2FC","p_val_adj")]
colnames(diff_exp) <- c("log2FC","pval")
plot_volcano(dataset = diff_exp,
pval_FC=c(0.05,2),
#label_to_id = tls_sign,
hide_label = F,
ylabel = "-log10 (Adjusted P value)") + xlim(lims=c(-5,5)) +
ggtitle(label=ids) +
theme(text=element_text(size=15))+
theme(text=element_text(face = "bold"))
})
png(paste0(res_folder,Sys.Date(),"_FIGURE_TLS_signature_supp.png"),width = 10,height=10,units = "in",res = 500)
grid.arrange(grobs=all_volcanoes,n=2)
dev.off()
上门的代码有一段plot_volcano可能是作者source的自建函数。我这里替换一个火山图函数
#figure supp TLS signature 差异基因的火山图
#plot volcano
names(TLS_list_all_samples) <- TLS_pos_ids #给差异基因的list命名
all_volcanoes <- lapply(TLS_pos_ids[c(3,2,1)], function(ids){
diff_exp <- TLS_list_all_samples[[ids]]
diff_exp <- diff_exp[,c("avg_log2FC","p_val_adj")]
colnames(diff_exp) <- c("log2FC","pval")
EnhancedVolcano(diff_exp,
lab = rownames(diff_exp),
x="log2FC",
y="pval",
pCutoff=0.05,
FCcutoff=2)+
# plot_volcano(dataset = diff_exp,
# pval_FC=c(0.05,2),
# #label_to_id = tls_sign,
# hide_label = F,
# ylabel = "-log10 (Adjusted P value)") +
xlim(lims=c(-5,5)) +
ggtitle(label=ids) +
theme(text=element_text(size=15))+
theme(text=element_text(face = "bold"))
})
all_vocanno
grid.arrange(grobs=all_volcanoes,n=1)
BOXPLOT AUC的箱装图
构建数据
library(stringr)
toto <- data.frame(unlist(cross_validation["auc",]))
toto$id <- str_split(rownames(toto),"\\.",simplify = T)[,1]
colnames(toto)[1] <- "values"
绘图
对比原文看一下
④展示TLS_imprint_signature和其它marker关系
使用scale.data
#figure supp 绘制共表达图的代码？？？？好像是的 下面画了Sup 4 figure
opt <- "C"
plot_list <- list()
for(spat in names(spatial_list)[c(1:3)]){ #原文是c(1:2,4:12)，这里进行替换1:3
spat <- names(spatial_list[1]) ###举例，不运行循环
print(spat)
library("patchwork")
pat <- list.files(path = "~/projects/visium/data/ihc",pattern =spat,full.names =T)
ihc <- readPNG(pat,native=T)
# 把IGKC和JCHAIN 和 intersect_tls_sign的gene计算平均值（SCT data）计入两个新的metadata
spatial_list[[spat]] <- AddMetaData(spatial_list[[spat]],apply(as.matrix(spatial_list[[spat]][["SCT"]]@data[c("IGKC","JCHAIN"),]),2,mean),col.name = "Pan_IG")
intersect_tls_sign <- intersect(tls_sign,rownames(spatial_list[[spat]]))
spatial_list[[spat]] <- AddMetaData(spatial_list[[spat]],apply(as.matrix(spatial_list[[spat]][["SCT"]]@data[intersect_tls_sign,]),2,mean),col.name = "TLS_Imprint")
#ROC
if(any(spatial_list[[spat]]$TLS_2_cat_annot=="TLS")){
print("compute_roc")
roc_empirical <- rocit(score = spatial_list[[spat]]$TLS_Imprint,
class = spatial_list[[spat]]$TLS_2_cat_annot=="TLS",
negref = F)
auc <- roc_empirical$AUC
}else{
auc <- NA
}
#compute spatial cor ????
TLS_Imprint <- as.vector(scale(spatial_list[[spat]]$TLS_Imprint)) #提取metadata，
Pan_IG <- as.vector(scale(spatial_list[[spat]]$Pan_IG))
spatial_list[[spat]]@assays$SCT@scale.data <- rbind(spatial_list[[spat]]@assays$SCT@scale.data,TLS_Imprint)#把结果合并到scale.data中(3000个基因+1)
spatial_list[[spat]]@assays$SCT@scale.data <- rbind(spatial_list[[spat]]@assays$SCT@scale.data,Pan_IG)#为什么是scale.data呢
sign_list <- c("TLS_Imprint",
"MZB1",
"Pan_IG",
"IGHM",
"IGHG1",
"IGHA1",
"COL1A1",
"CXCL12")
spat_cor_int <- intersect(sign_list,rownames(spatial_list[[spat]]@assays$SCT@scale.data)) #取交集
spat_cor_ext <- setdiff(sign_list,rownames(spatial_list[[spat]]@assays$SCT@scale.data)) #取差集
# spatial_cor <- RunMoransI(data = spatial_list[[spat]]@assays$SCT@scale.data[spat_cor_int,],
# pos = GetTissueCoordinates(object = spatial_list[[spat]]))
# spatial_cor[spat_cor_ext,] <-c(NA,NA)
he <- SpatialPlot(spatial_list[[spat]],repel = F,label = F,image.alpha=1,alpha = c(0,0), pt.size.factor =0) + NoLegend()
blank <- SpatialPlot(spatial_list[[spat]],repel = F,label = F,image.alpha=0,alpha = c(0,0), pt.size.factor =0) + NoLegend()+
inset_element(p = ihc,
left = 0,
bottom = 0,
right = 1,
top = 1)
all_plots <-lapply(sign_list,function(sign){
p <- SpatialPlot(spatial_list[[spat]],image.alpha = 0, features = sign)
p <- p+ scale_fill_viridis_c(option="C")
p <- p+ DarkTheme()
# p <- p+ annotate("text",
# x=max(GetTissueCoordinates(object = spatial_list[[spat]])[,2]) - 10,
# y = max(GetTissueCoordinates(object = spatial_list[[spat]])[,1]) +30,
# size=4,
# vjust=1,
# hjust=1,
# col="white"),
# label = paste0("Spatial cor: ", round(spatial_cor[sign,"observed"],digits = 2)))
if(sign=="TLS_Imprint" & spat %in% names(spatial_list)[c(1,2,4,6:9)]){
p <- p+ annotate("text",
x=min(GetTissueCoordinates(object = spatial_list[[spat]])[,2]) + 100,
y=min(GetTissueCoordinates(object = spatial_list[[spat]])[,1]) +30,
size=4,
vjust=1,
hjust=1,
col="white",
label = paste0("AUC = ", round(auc,2)))
}
p <- p + hide_axis +theme(legend.title=element_text(size=15))
})
plot_list[[spat]] <- (he | blank | all_plots[[1]] | all_plots[[2]] |all_plots[[3]] ) / ( all_plots[[4]] | all_plots[[5]] | all_plots[[6]] |all_plots[[7]]| all_plots[[8]])
}
pdf(paste0(res_folder,Sys.Date(),"_FIGURE_4_supp.pdf"),width = 15,height=15)
plot_list
dev.off()
绘图用scale.data
最后出的一些图没有太明白和上文的有什么区别
扫一扫
473
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
