一、写在前面

突发奇想，现在生信数据的存储主要阵地——GEO数据库所收录的第一篇SCI论文做了怎样的工作？于是我们在GEO数据库中进行了检索(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE000001)并找到了这篇文章。这是一篇2000 年 8 月 1 日爱荷华大学N. Hayward& J. Trent团队在《Nature》上在线发表了题为“Molecular classifcation of cutaneous malignant melanoma by gene expression profling”的letter，主要包含了38个样本的探针数据并进行了黑色素瘤的分型。站在上帝视角的我们显然觉得探针的数据分析起来非常的简单(相关性分析、层次聚类谁不会啊)，这类探针数据在如今单细胞、空转等多组学数据满天飞的现在显得有些小儿科(时代的一粒沙，个人的一座山)。但原本计划两天完成的复现工作，实际分析过程我们来回拉扯了一个月。其中许多工具、可视化方法现如今已经逐渐的退出了历史舞台，但大牛们当年能够迈出GEO数据库的第一步，初步提出了表达矩阵的肿瘤分型思想，不可谓不是跨时代*(永远可以相信《Nature》的含金量)。

1.png

由于显示问题，本文内容略有缺失点击查看原文

简介：

人类最常见的一种皮肤相关的癌症是皮肤恶性肿瘤。皮肤黑色素瘤的发病率正在急剧上升，但在晚期疾病的非手术治疗方面几乎没有进展。尽管付出了大量努力试图识别黑色素瘤预后的独立预测因子，但目前(当时，2000年以前)没有公认的组织病理学、分子或免疫组织化学标志物能够定义这种肿瘤的亚型。因此，文章报告通过对一系列样本的基因表达进行分析发现了一种黑色素瘤的亚型。令人惊讶的是，许多用于定义这一亚型的基因在体外形成原始管状网络的侵袭性黑色素瘤中表现出差异性调控，而这种特征出现在某些高度侵袭性转移性黑色素瘤中。转录组分析可以识别出未被认识的皮肤黑色素瘤亚型，并预测可能对疾病进展重要的、可通过实验验证的表型特征。

往期复现工作：

《NC》代码复现| scRNA-seq揭示人纤维化皮肤病中的成纤维细胞异质性和间充质成纤维细胞增加

单细胞复写|3.急性心肌梗塞(AMI)外周血scRNA-Seq分析实战(链接重置)

二、主要内容

（1）皮肤黑色素瘤样本聚类分析

文章收集了38例皮肤样本的转录组数据，包括31例黑色素瘤样本和7例对照样本。fig1展示了几种分析方法的整合，以可视化皮肤黑色素瘤肿瘤样本之间的整体表达模式关系。通过使用所有肿瘤样本的皮尔逊相关系数矩阵，将31个黑色素瘤样本以层次聚类树状图的形式展示出来（fig1a），同时通过多维尺度分析（MDS）展示31例样本的三维分布（fig1b）。分层树状图表明，19个样本紧密聚类在树状图的底部，处于相似度最高的区域。目前还没有一种既定的方法来评估通过聚类预测技术获得的分群关系的可靠性。因此，文章使用两种独立的方法进一步进行测试。第一种方法通过检查强分类基因的频率与此类基因的预期频率（如果表达是随机变化的）以及与将相同样本随机划分为19个和12个新分组时强分类器的频率进行比较，来检查单个基因区分19个主要聚类与其余样本的能力。聚类结果的非随机性是显而易见的。具体而言，许多基因的表达模式在初始样本聚类之间有很大差异，因此可以作为良好的分类器（fig1c，红色三角形）。然而，当样本被分组为相同大小的随机分区时，很难找到对样本进行分类的表达模式（fig1c，蓝线）。因此，在随机形成的簇中，表达行为与基因相对于这些簇的真正随机行为基本上没有区别（fig1c，比较蓝线与空心圆）。第二种方法是基于向数据集引入随机扰动后评估聚类成员资格。对于每个样本，通过引入均值为0、标准偏差为0.15的随机高斯噪声来扰动每个基因的对数比（这是通过计算所有31个样本中单个基因的对数比的中位标准偏差得出的变化估计值）。然后对扰动后的数据集进行层次聚类，并对原始树（fig1a）和扰动树进行比较。比较包括原始树和扰动树分割成k个聚类，然后计算原始树中成对样本聚类在一起但在扰动树中不聚类在一起的配对样本的比例（将此度量称为差异对的比例加权，因为它对较大的聚类给予更大的权重）。对于每个可能的分割，在多个扰动数据集上重复进行比较。给定k（k=2：30），随后对扰动数据集中差异对的有权比例进行平均，从而确定了加权平均差异对（WADPk）。将原始树分割成9组或更少组所得到的聚类具有很强的可重复性（fig1d）。值得注意的是，WADPk的增加几乎与主要包含19个元素的聚类分裂成较小的子聚类完全重合。这些结果有力地支持了这样一种观点，即在本研究中识别出的黑色素瘤样本的主要聚类代表了一个真实且高度可重复的分组。

2.png

（2）两组患者肿瘤倾袭性比较

文中还纳入了两对来自同一患者的标本。它们是M92-001和M93-007（来自同一患者的两个不同样本，分别在手术切除后一年获得），以及TD-1376-3和TC1376-3（同一肿瘤的活检样本和细胞培养物在体外传代三次）。尽管细胞传代次数与聚类组之间没有显著关联，但纳入了TD-1376-3/TC-1376-3这对样本作为细胞培养效果的另一个对照。基于fig1中全部基因表达的线性相关性，fig2和fig3展示了用于指导“基因簇”解释的经验方法。fig2a描绘了用于提取单个基因“加权列表”的统计方法，这些基因在所有实验中的变化方差正确地定义了给定样本簇的边界。fig2b沿垂直轴按降序排列了在19个样本（fig1a和b）中定义主要黑色素瘤簇最有力的基因列表。样本沿水平轴按簇包含情况排序，颜色饱和度与测量到的基因表达比率的幅度成正比。加权基因列表也可用于指导对更大规模的基因表达数据集的分析。fig3a展示了本研究中皮肤黑色素瘤样本的所有数据，以彩色图像的形式呈现，基因按照表达模式的相似性沿垂直轴排列。fig3b-e展示了使用fig2b中“加权列表”中的四个基因（MART-1、CD63、 tropomyosin和WNT5A）。

3.png

最后，在对皮肤黑色素瘤的微阵列分析的同时，研究了一系列具有转移相关特征的葡萄膜黑色素瘤标本，包括侵袭能力和体外血管生成模拟能力。令人惊讶的是，高度侵袭性葡萄膜黑色素瘤细胞与皮肤黑色素瘤样本主要簇列中的相同基因呈强烈负相关（fig2b）。这一观察结果，结合加权列表中基因的已知生物学功能，表明被分配到主要皮肤黑色素瘤簇列中的样本具有较低的运动性和侵袭能力，因为它们下调了与细胞扩散或迁移相关的基因。在主要聚类之外的样本中基因表达趋势则是相反的，这些与其他文献报道一致。进一步使用一系列细胞检测来表征皮肤黑色素瘤迁移性和侵袭性区别。正如对它们基因表达模式的分析所预测的那样，与主要簇组外的黑色素瘤相比，主要簇组内的黑色素瘤的迁移性、侵袭能力和血管生成模拟能力都有所降低。

4.png

三、小结

本研究中的患者群体预后普遍较差，无论是典型的临床因素（如年龄、性别、活检部位）还是体外特征（如传代次数）都与临床结果没有强相关性。相比之下，基于基因表达对这些肿瘤进行分子分类（fig1）可以识别出一种此前未被检测到的癌症亚型。总共有15名患者的生存信息可用，尽管结果在统计学上不显著，但具有一定意义。在紧密关联的19名患者中，有生存信息的10名患者中有3人死亡，而在其余组中，5名患者中有4人死亡。在目前无法确定谁最有可能出现疾病快速进展和死亡风险的情况下，基于基因表达模式对黑色素瘤进行分类是有可能的，可以进一步用于识别对转移过程的各个方面至关重要的基因。目前黑色素瘤样本在临床上通过表达模式进行细分的能力程度仍有待阐明，但是本文为黑色素瘤的其他临床相关亚群研究提供了坚实的分子基础。

四、文章复现

1.复现集锦

5.jpg

2.数据下载和预处理

示例数据集从GEO数据库下载即可。

set.seed(123)
setwd("/home/cwj/project/01_GEO")
suppressMessages({  library(GEOquery)  
library(scatterplot3d)  library(pheatmap)  
library(ggplot2)
  # 安装并加载 ClassComparison 包
  # install.packages("ClassComparison")  
library(ClassComparison)})
# 从GEO数据库获取基因表达矩阵
data <- getGEO("GSE1", destdir="/home/cwj/project/01_GEO",GSEMatrix =TRUE, getGPL=TRUE, AnnotGPL=TRUE)
exp <- exprs(data[[1]])#表达矩阵信息
cli <- pData(data[[1]]) #临床信息
GPL <- fData(data[[1]]) #平台信息
gpl <- GPL[,c(1,3,11)] #提取需要的ID列和gene symbol列write.csv(exp,"GSE1_expeSet.csv")write.csv(cli,"GSE1_metadata.csv")
gpl$'Gene symbol' <- data.frame(sapply(gpl$'Gene symbol',function(x)unlist(strsplit(x,"///"))[1]),stringsAsFactors=F)[,1]
exp <- as.data.frame(exp)
# 样本名替换成结果图中使用的样本名
colnames(exp)<-cli$title 
#增加新列存储基因ID信息
exp$ID <- rownames(exp) 
#为表达矩阵添加Gene symbol名称
exp <-merge(exp,gpl,by="ID") 
dim(exp)

## [1] 8192   41

# 相同基因存在不同转录本，行名重新命名
rownames(exp) <- paste0(rownames(exp),"_",exp[,"Gene symbol"])
#删除无Gene symbol的行exp <- na.omit(exp) 
dim(exp)

## [1] 6275   41

sample_all <- cli$title

3.31例黑色素瘤样本聚类

# fig1a

6.png

# 方法2：基于所有基因表达矩阵先计算样本之间的相关性矩阵，基于相关性，进行聚类划分（与原文一致性最高）

7.png

# 方法3：unique使用剩下的4902个基因表达矩阵先计算样本之间的相关性矩阵，基于相关性，进行聚类划分（与原文一致性次之）

8.png

4.31例黑色素瘤样本三维空间分布

# fig1b：

9.png

5.基因与TNOM评分函数

# fig1c

## Class "TNoM" [package "ClassComparison"]## ## Slots:##                                                                         ## Name:        data   tnomData       nCol       nRow classifier       call## Class:     matrix    numeric    numeric    numeric     factor       call

showClass("fullTNoM")

## Class "fullTNoM" [package "ClassComparison"]## ## Slots:##                                                         ## Name:        dex     fakir       obs       scr      name## Class:   numeric   numeric   numeric   numeric character

bogus <- as.matrix(exp_mali)

## [1] "ordering..."## [1] "matrifying..."## [1] "lapplying..."## [1] "cumsuming..."## [1] "more matrifying..."## [1] "another apply..."

summary(tn)

## TNoM object with 6275 rows and 31 columns## Call: TNoM(data = bogus, classes = splitter) ## ## Number of genes with maximum number of misclassifications:##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 ##    0    0    0    1    1    5   47  177  479 1043 1643 1998  881    0    0    0 ##   16 ##    0

tnf <- update(tn)

## Simulation number: 1 . 2 . 3 . 4 . 5 . 6 . 7 . 8 . 9 . 10 .

plot(tnf)

10.png

6.WADP与K值

# fig1d

111.png

7.top22权重基因筛选

# fig2b

12.png

pheatmap(exp_mali[top22,order], cluster_row = T, cluster_col = F, scale="row", color = colorRampPalette(c("green", "black", "red"))(100), border_color = NA)

13.png

8.其他基因表达模式展示

#fig3 样本*所有基因矩阵的热图

14.png

pheatmap(t(exp_mali[,order]), cluster_row = F, cluster_col = T, scale="column", show_colnames=F, show_rownames=F, color = colorRampPalette(c("green", "red"))(100), border_color = NA)

15.png

# 计算基因和基因之间的相关性，原文是按照每个基因所在簇进行定位,但是尝试将基因分成不同簇后，没有得到关注基因合适的聚类分群，因此调整策略计算关注基因相关性高的基因用于后续展示

16.png

pheatmap(exp_mali[cor_gene,order], cluster_row = F, cluster_col = F,scale="row", annotation_row = anno_row, color = colorRampPalette(c("green", "black", "red"))(100), border_color = NA)

17.png

五、环境版本信息

sessionInfo()

## R version 4.2.2 Patched (2022-11-10 r83330)## Platform: x86_64-pc-linux-gnu (64-bit)## Running under: Ubuntu 20.04.6 LTS## ## Matrix products: default## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0## ## locale:##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              ##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    ##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   ##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       ## ## attached base packages:## [1] stats     graphics  grDevices utils     datasets  methods   base     ## ## other attached packages:## [1] ClassComparison_3.3.2 oompaBase_3.2.9       ggplot2_3.5.1        ## [4] pheatmap_1.0.12       scatterplot3d_0.3-42  GEOquery_2.66.0      ## [7] Biobase_2.58.0        BiocGenerics_0.44.0  ## ## loaded via a namespace (and not attached):##  [1] tidyselect_1.2.1    xfun_0.35           bslib_0.4.1        ##  [4] purrr_1.0.2         splines_4.2.2       colorspace_2.1-1   ##  [7] openintro_2.4.0     vctrs_0.6.5         generics_0.1.3     ## [10] htmltools_0.5.3     yaml_2.3.6          utf8_1.2.4         ## [13] rlang_1.1.4         R.oo_1.25.0         jquerylib_0.1.4    ## [16] pillar_1.9.0        R.utils_2.12.2      glue_1.8.0         ## [19] withr_3.0.1         DBI_1.1.3           RColorBrewer_1.1-3 ## [22] lifecycle_1.0.4     stringr_1.5.1       airports_0.1.0     ## [25] munsell_0.5.1       gtable_0.3.5        R.methodsS3_1.8.2  ## [28] cherryblossom_0.1.0 evaluate_0.18       labeling_0.4.3     ## [31] knitr_1.41          tzdb_0.4.0          fastmap_1.1.0      ## [34] curl_4.3.3          fansi_1.0.6         highr_0.9          ## [37] readr_2.1.5         scales_1.3.0        cachem_1.0.6       ## [40] limma_3.54.0        jsonlite_1.8.3      farver_2.1.2       ## [43] hms_1.1.3           digest_0.6.30       stringi_1.8.4      ## [46] dplyr_1.1.4         grid_4.2.2          cli_3.6.3          ## [49] tools_4.2.2         magrittr_2.0.3      sass_0.4.3         ## [52] tibble_3.2.1        cluster_2.1.4       usdata_0.2.0       ## [55] tidyr_1.3.1         pkgconfig_2.0.3     data.table_1.14.6  ## [58] xml2_1.3.3          rmarkdown_2.18      rstudioapi_0.14    ## [61] R6_2.5.1            compiler_4.2.2

六、参考

https://www.nature.com/articles/35020115

https://cran.rproject.org/web/packages/ClassComparison/ClassComparison.pdf