R version: R Under development (unstable) (2023-10-22 r85388)
Bioconductor version: 3.19
Package version: 1.19.0

1 介绍

微阵列(microarray)技术的发展和RNA测序技术的快速应用为检测生物样品的转录谱(或转录组)提供了平台(Cieślik and Chinnaiyan 2018)。转录组学分析通常侧重于样品不同分组之间基因的“差异表达(differential expression)”,然而随着公开可用RNA数据的快速增长,越来越多的人开始使用相对的方法来量化样品与特定的基因标签(gene signatures)代表的特征之间的一致性(Cieślik and Chinnaiyan 2018)。对于用突变或融合基因进行肿瘤亚型分类与鉴定,以及寻找驱动癌症进展的基因组损伤而言,基因组的测序十分重要,而转录组学分析则可以提供携带这些突变的细胞的形态或表型信息。

癌症是一种异质性疾病,其具有多种临床和病理亚型。以乳腺癌为例,临床上主要依靠激素受体(雌激素受体:ER;孕酮受体:PR)的表达或Erb-B2受体酪氨酸激酶(HER2)的过度表达进行分型,这些特征的癌症分型有可直接靶向治疗的药物。一个常见的例子就是临床上使用PAM50(prediction analysis of microarray 50)标签来区分不同的乳腺癌亚型(Parker et al. 2009, @cieslik18)。许多癌症的亚型分类主要依赖于识别大型患者队列中的复发性突变(recurrent mutations),通过全基因组或全外显子组测序来获得具有临床意义的亚型(Cancer Genome Atlas Research Network 2013, @papaemmanuil16)

流行的“relative approach”有single-sample gene set enrichment analysis (ssGSEA) (Barbie et al. 2009),用户可以通过GenePattern web-tool使用该方法(http://software.broadinstitute.org/cancer/software/genepattern)。另一种常见的方法是gene set variation analysis (GSVA) (Hänzelmann, Castelo, and Guinney 2013),以R/Bioconductor package的形式实现(https://bioconductor.org/packages/release/bioc/html/GSVA.html)。该package中还包括了ssGSEA,PLAGE (Tomfohr, Lu, and Kepler 2005)和z-score (Lee et al. 2008)方法。ssGSEA和GSVA都使用了Kolmogorov-Smirnov like random-walk statistic来将归一化(normalised)的基因排名转换成样本得分。这一“归一化(normalised)”的步骤使得这些方法不适用于单一样本(single-sample),而且样本组成带来的variations(比如样本由不同的肿瘤亚型组成)会对样本得分造成影响。其次,这些方法产生的样本得分具有不同的数值范围(range),为阐释结果带来一些困难。为了克服这些问题,我们开发了一种单样本的基因集打分方法singscore (Foroutan et al. 2018) (http://bioconductor.org/packages/singscore/), 它利用了给定基因集中基因的排名并相对这些基因排名的最大值与最小值进行归一化(normalization)。

The Cancer Genome Atlas(TCGA)提供了上千的病人转录组数据,通常这些数据还有配对的基因组或表观基因组数据(通常为DNA甲基化数据)。这些数据可以帮助我们找到突变所对应的功能上的变化,探索由于表观因素或者转录调控而引起的肿瘤异质性。在这里,我们展示了singscore方法 (Foroutan et al. 2018) 可以利用NPM1c mutation, KMT2A (MLL) 基因融合(gene fusions), 和 PML-RARA 基因融合(gene fusions)的转录组gene signatures对TCGA中的AML样本进行分类。在不需要参数拟合或估计的情况下,用singscore对基因集打分可以区分出携带有这些突变的样本。在这里,我们将介绍基因集打分不仅可以用于识别差异,还可以衡量临床结果不同的AML亚型之间的相对相似性。

2 生物问题简述

与大多数癌症一样,急性髓细胞白血病(acute myeloid leukemia,AML)是一种具有许多亚型的异质性疾病。通过对TCGA AML基因组数据的分析,有人根据特定“驱动突变(driver mutations)”的存在与否对AML进行了分类,总结并扩展了之前已定义的临床上的亚型 (Cancer Genome Atlas Research Network 2013)。最近一项主要针对基因组数据的研究进一步完善了具有临床意义的AML亚型 (Papaemmanuil et al. 2016),其中囊括了许多共同发生以及相互排斥的突变。

值得注意的是,在临床AML样本中最常见的突变之一是NPM1基因第12外显子上的移码突变(frameshift mutation)(Papaemmanuil et al. 2016)。这种突变导致核磷蛋白(nucleophosmin)的异常定位与胞质积累(cytoplasmic accumulation),因此这种突变经常被称作NPM1c突变 (Brunetti et al. 2018)。如 Verhaak et al. (2005) 所述,NPM1c突变与同源框结构域(Hox)转录因子家族的活性失调相关,而该转录因子家族对于发育模式(developmental patterning)是必需的。最近的研究进一步证实了该突变在疾病进展中的作用,NPM1c的缺失会导致AML细胞的分化(Brunetti et al. 2018)

AML中的复发性遗传损伤还包括赖氨酸甲基转移酶2A(KMT2A,以前被称为MLL)融合基因,KMT2A基因(KMT2A -PTD)内的部分串联重复,以及早幼粒细胞白血病蛋白(PML)和视黄酸受体α(RARA)之间的基因融合(PML-RARA)。鉴于NPM1c突变对Hox基因家族的失调作用,具有MLL基因融合的AML样本显示出Hox家族基因的表达失调成为了一个值得探究的问题 (Hess 2004, @ross04)。然而,具有MLL-PTD的样本似乎显示出与MLL-融合样本相对不同的表型 (Ross et al. 2004)。尽管有充分的证据证明NPM1c突变和其它遗传损伤在阻碍AML细胞分化中的作用,但具有PML-RARA融合的样本被诊断为一种叫做急性早幼粒细胞白血病(acute promyelocytic leukemia,APL)的AML亚型。这种临床AML亚型与French-American-British (FAB)分类系统中的FAB-M3相关。由于早幼粒细胞阶段的分化阻滞,该亚型的细胞显示出特有的形态 (Thé et al. 1991)

在本分析流程中,我们证明了singscore方法 (Foroutan et al. 2018)从转录组数据中区分不同肿瘤“驱动突变(driver mutations)”的能力。我们使用了已定义的NPM1c突变的gene signature(Verhaak et al. 2005),以及PML-RARA基因融合和MLL融合的gene signatures。后两个gene signatures来自于儿童AML样本,尽管儿童AML和成人AML之间的突变谱差别较大(Ma et al. 2018),但研究表明该signatures能够很好地区分具有类似基因损伤的成人AML样本(Ross et al. 2004)。这些gene signatures已存在于MSigDB数据库(molecular signatures database)中(Liberzon et al. 2015),通过这些gene signatures,我们证明了singscore的双向打分方法可以根据不同的突变将TCGA AML样本分类,且结果具有良好的精确度(precision)和召回率(recall)。该方法的优点之一在于它是能够根据表型特征(phenotypic signatures)将样本映射到二维或者更高维度的空间中。通过比较NPM1c和KMT2A-/MLL- 融合signatures的得分,我们展示了这种分类可能是由于Hox家族失调产生的共同的下游效应而导致的。我们还将NPM1c突变signature与PML-RARA signature进行比较,这两个亚型的明显分离反映了它们显著不同的表型和相互排斥的突变。

3 数据下载与处理

我们可以通过Genomic Data Commons(GDC)获得多种不同预处理的TCGA数据。转录组数据有count值和FPKM值,有用上四分位数标准化(upper quantile normalisation)前后的值。经过其它方式预处理后的数据可以在 www.cbioportal.orgfirebrowse.org 找到。GDC的数据使用了STAR的“two-pass”模式进行序列比对,使用了STAR进行定量。用户在使用 GDC data transfer tool 下载GDC中的数据时,可以通过GDC portal选择自己感兴趣的特定文件进行下载。下载之后需要读取并整合这些文件之后才可以进行下游分析。以上这些步骤(包括下载)可以通过R包TCGAbiolinks (Colaprico et al. 2015, @R-TCGAbiolinks)完成。这个包支持使用GDC API和GDC client进行数据下载,我们会使用该包完成下载、注释并将数据整合成SummarizedExperiment R对象。

在开始分析之前需要完成一下几步数据处理: 1. 创建查询列表以下载特定的文件; 2. 执行查询以下载数据; 3. 将下载好的数据读入R中; 4. 过滤掉表达量低的基因; 5. 校正样本组成带来的偏差并对基因长度进行标准化 (Foroutan et al. 2018)

3.1 查询GDC数据库

任何数据分析的第一步是先确定好数据的版本和下载数据的方式。getGDCInfo()函数可以返回GDC数据库中数据的版本号和发布时间。

library(SingscoreAMLMutations)
library(TCGAbiolinks)

#get GDC version information
gdc_info = getGDCInfo()
gdc_info
## $commit
## [1] "3d0c5b1fd25804e51093b02b870742d9042e75bd"
## 
## $data_release
## [1] "Data Release 38.0 - August 31, 2023"
## 
## $status
## [1] "OK"
## 
## $tag
## [1] "4.0.1"
## 
## $version
## [1] 1

接下来我们需要创建一个查询命令,从GDC找到并下载特定的文件。这一步和使用GDC portal创建MANIFEST文件很类似。查询命令的第一个参数(project)指定了项目的名称(使用getGDCprojects()函数或者进入https://portal.gdc.cancer.gov/projects可以获得所有项目名称),TCGA急性髓细胞样白血病(acute myeloid leukemia,AML)数据属于TCGA-LAML项目。接下来还需要指定数据分类(data.category)、数据类型(data.type)和工作流程类型(workflow.type)。下面这个查询命令指定了count-level的转录组数据。输入各参数的值可以从“query”vignette文档的“searching arguments”部分找到(使用vignette("query", package = "TCGAbiolinks")命令)。查询命令最终会返回一个包含了文件名和相关注释的数据框。

这里我们选择count-level的数据而不是FPKM是因为我们先要对基因进行过滤。通常,FPKM的计算在基因过滤之后,以确保FPKM计算时使用的library sizes大小正确。如果library sizes足够大,在无法获得count-level数据的情况下使用FPKM值也是合理的。

#form a query for the RNAseq data
query_rna = GDCquery(
  #getGDCprojects()
  project = 'TCGA-LAML',
  #TCGAbiolinks:::getProjectSummary('TCGA-LAML')
  data.category = 'Transcriptome Profiling',
  data.type = 'Gene Expression Quantification',
  workflow.type = 'STAR - Counts'
)

#extract results of the query
rnaseq_res = getResults(query_rna)
dim(rnaseq_res)
## [1] 151  29
colnames(rnaseq_res)
##  [1] "id"                        "data_format"              
##  [3] "cases"                     "access"                   
##  [5] "file_name"                 "submitter_id"             
##  [7] "data_category"             "type"                     
##  [9] "file_size"                 "created_datetime"         
## [11] "md5sum"                    "updated_datetime"         
## [13] "file_id"                   "data_type"                
## [15] "state"                     "experimental_strategy"    
## [17] "version"                   "data_release"             
## [19] "project"                   "analysis_id"              
## [21] "analysis_state"            "analysis_submitter_id"    
## [23] "analysis_workflow_link"    "analysis_workflow_type"   
## [25] "analysis_workflow_version" "sample_type"              
## [27] "is_ffpe"                   "cases.submitter_id"       
## [29] "sample.submitter_id"

3.2 下载TCGA AML RNA-seq数据(read counts)

GDCdownload函数可以执行查询并使用GDC API下载数据。如果需要下载的文件很大(如RNA-seq的read数据或者甲基化数据),应该把GDCdownload函数里的下载方法切换成“client”。这里我们暂时将数据存储在“GDCdata”文件夹中,用户应该指定自己的存储路径以妥善保存好数据。GDCdownload函数通过这种参数设定的方法使我们能够将不同类型的数据存储在同一个文件夹内进行管理与使用。这里,下载后我们可以看到count-level数据被存储在TEMPDIR/GDCdata/TCGA-LAML/harmonized/Transcriptome_Profiling/Gene_Expression_Quantification/路径下。

datapath = file.path(tempdir(), 'GDCdata')
GDCdownload(query_rna, directory = datapath) #(size: 170MB)

3.3 将count-level数据读入R

GDCprepare函数可以将下载好的数据读入R并处理成RangedSummarizedExperiment对象(来自SummarizedExperiment包),该对象可以同时储存read count、基因注释和临床信息。在调用该函数的时候临床信息自动下载并整合入RangedSummarizedExperiment对象中。RangedSummarizedExperiment对象和ExpressionSet很类似,但它还可以使用基因坐标进行索引以及存储具有相同结构的多个数据矩阵。关于数据特征(feature)的注释被存储在一个RDA/RDATA文件中。

aml_se = GDCprepare(query_rna, directory = datapath)

RangedSummarizedExperiment对象中包含了60660个特征(features)和150个样本。用rowData(se)colData(se)函数可以获得特征和样本的注释信息,用assay(se)函数可以获得counts数据。TCGA数据往往包含一些福尔马林固定、石蜡包埋(formalin-fixed paraffin-embedded,FFPE)的样本,这些样本需要被过滤掉以避免protocol不同而引入的误差。这一步在本数据集中不需要,因为这一步骤是针对实体瘤(solid tumors)而非白血病(leukemias)的。

aml_se
## class: RangedSummarizedExperiment 
## dim: 60660 151 
## metadata(1): data_release
## assays(6): unstranded stranded_first ... fpkm_unstrand fpkm_uq_unstrand
## rownames(60660): ENSG00000000003.15 ENSG00000000005.6 ...
##   ENSG00000288674.1 ENSG00000288675.1
## rowData names(10): source type ... hgnc_id havana_gene
## colnames(151): TCGA-AB-2937-03A-01T-0734-13
##   TCGA-AB-2863-03A-01T-0734-13 ... TCGA-AB-2979-03B-01T-0760-13
##   TCGA-AB-2918-03A-01T-0740-13
## colData names(50): barcode patient ... released sample.aux

3.4 过滤掉counts数低的基因

edgeR包提供了过滤所需的数据标准化与转换方法。这些方法要求数据存储在DGEList对象中,因此我们需要将SummarizedExperiment对象转换成DGEList对象。

library(SummarizedExperiment)
library(edgeR)

#remove ENSEMBL ID version numbers
rownames(aml_se) <- gsub('\\.[0-9]*', '', rownames(aml_se))

aml_dge = DGEList(counts = assay(aml_se), genes = rowData(aml_se))

本分析流程中,我们过滤掉了在大多数样本中表达量很低的基因。这一步骤是差异表达分析的标准步骤,因为这些基因导致离差(dispersion)估计发生偏离。在以排序为基础的方法中这一步也很有必要,因为重复的排名(rank duplication)会减弱方法的判别能力。通常我们只选择那些在一定比例样本中CPM值高于一定阈值的基因。之所以用CPM值进行过滤而不是raw counts是因为CPM值对总reads数,也就是library sizes进行了均一化,因此CPM值相对于raw counts而言是无偏的。阈值的选择不是绝对的,假设AML数据的library sizes范围是1860~4970(百万reads),那么CPM = 1意味着read counts的范围为19~50(reads)。在这里我们保留了在50%的样本中CPM值超过1(CPM > 1)的基因。在特定情况下可能其它过滤低表达基因的方法更加适用。 Chen, Lun, and Smyth (2016)Law et al. (2016) 根据实验设计来过滤基因,他们先在各个组内确定表达量过低的基因,然后在所有样本中过滤掉这些基因。这种过滤策略适合那些样本较少的数据集,AML数据有足够多的样本因此我们直接对所有样本进行过滤。在图1中我们可以看到过滤后logCPMs值的分布更加接近理想的log-normal分布。

prop_expressed = rowMeans(edgeR::cpm(aml_dge) > 1)
keep = prop_expressed > 0.5

op = par(no.readonly = TRUE)
par(mfrow = c(1, 2))
hist(edgeR::cpm(aml_dge, log = TRUE), main = 'Unfiltered', xlab = 'logCPM')
abline(v = log(1), lty = 2, col = 2)
hist(edgeR::cpm(aml_dge[keep, ], log = TRUE), main = 'Filtered', xlab = 'logCPM')
abline(v = log(1), lty = 2, col = 2)
过滤前后AML数据logCPM值的直方图. 过滤之后数据中的零值减少了。在大部分样本里CPM值小于1(logCPM < 0)的基因被过滤掉,使得最后得到了近似log-normal分布的logCPM值。

Figure 1: 过滤前后AML数据logCPM值的直方图
过滤之后数据中的零值减少了。在大部分样本里CPM值小于1(logCPM < 0)的基因被过滤掉,使得最后得到了近似log-normal分布的logCPM值。

par(op)
#subset the data
aml_dge = aml_dge[keep, , keep.lib.sizes = FALSE]
aml_se = aml_se[keep, ]

3.5 转换成FPKM值并归一化(normalisation)

Singscore要求同一个样本内的基因表达量是可以互相比较的,因此我们需要对基因长度进行归一化(Oshlack and Wakefield 2009)。数据可以被转换成transcripts per million (TPM)或者reads/fragments per kilobase per million (RPKM/FPKM),它们都对基因长度进行了归一化。因此针对这两种转换方法,只要library size足够大,singscore得到的结果应当是类似的。edgeR包中的calcNormFactors函数提供了三种对基因长度归一化的方法,在不指定的情况下TMM normalisation是默认的方法。 Chen, Lun, and Smyth (2016)Law et al. (2016) 讨论了在下游分析如差异表达分析之前进行数据归一化的意义。通常归一化的目的是使得样本之间的比较是有意义的。Singscore的分析是针对单个样本内的因此无需对整体样本进行归一化。同理,该方法也不受其它转换方式(如log transformation)的影响。这里我们仅出于可视化目的使用TMM normalisation。

将数据转换成TPM或者RPKM/FPKM值需要先计算所有基因的长度。基因长度的计算依赖于序列比对过程和定量参数。TCGA转录组数据使用了STAR进行序列比对并用STAR进行定量(流程细节见https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/)。STAR统计了比对到每个基因外显子区域的read数量,因此有效基因长度应当是基因的所有外显子区域长度之和。GENCODE v36被用于定量过程中的注释,因此在计算基因长度时仍然应该用这个注释文件。

#download v36 of the GENCODE annotation
library(BiocFileCache) 
gencode_file = 'gencode.v36.annotation.gtf.gz'
gencode_link = paste(
  'ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_36',
  gencode_file,
  sep = '/'
  )
bfc <- BiocFileCache() 
gencode_path <- bfcrpath(bfc, gencode_link) 

rtracklayer这个R包提供了解析GTF文件的函数。

library(rtracklayer)
library(plyr)

gtf = import.gff(gencode_path, format = 'gtf', genome = 'GRCm38.71', feature.type = 'exon')
#split records by gene to group exons of the same gene
grl = reduce(split(gtf, elementMetadata(gtf)$gene_id))
gene_lengths = ldply(grl, function(x) {
  #sum up the length of individual exons
    return(c('gene_length' = sum(width(x))))
}, .id = 'ensembl_gene_id')

我们获取了基因的Ensembl IDs和gene type注释信息以便于下游分析。直接从GTF文件里获得的Ensembl IDs末尾带有版本号,我们需要将其去掉以转换为正式的Ensembl IDs。

#extract information on gene biotype
genetype = unique(elementMetadata(gtf)[, c('gene_id', 'gene_type')])
colnames(genetype)[1] = 'ensembl_gene_id'
gene_lengths = merge(genetype, gene_lengths)

#remove ENSEMBL ID version numbers
gene_lengths$ensembl_gene_id = gsub('\\.[0-9]*', '', gene_lengths$ensembl_gene_id)
saveRDS(gene_lengths, file = 'gene_lengths_gencode_v36.rds')
gene_lengths
## DataFrame with 60660 rows and 3 columns
##       ensembl_gene_id      gene_type gene_length
##           <character>    <character>   <integer>
## 1     ENSG00000000003 protein_coding        4536
## 2     ENSG00000000005 protein_coding        1476
## 3     ENSG00000000419 protein_coding        1207
## 4     ENSG00000000457 protein_coding        6883
## 5     ENSG00000000460 protein_coding        5970
## ...               ...            ...         ...
## 60656 ENSG00000288669 protein_coding        3987
## 60657 ENSG00000288670         lncRNA        1807
## 60658 ENSG00000288671 protein_coding         213
## 60659 ENSG00000288674 protein_coding        9620
## 60660 ENSG00000288675 protein_coding        1680

SummarizedExperiment对象可以存储基因的注释信息,因此我们需要把Ensembl IDs、gene length和gene type加进去。类似的,我们也要把这些注释信息加到DGEList对象中。

#allocate rownames for ease of indexing
rownames(gene_lengths) = gene_lengths$ensembl_gene_id
rowData(aml_se)$gene_length = gene_lengths[rownames(aml_se), 'gene_length']
rowData(aml_se)$gene_biotype = gene_lengths[rownames(aml_se), 'gene_type']

#annotate gene lengths for the DGE object
aml_dge$genes$length = gene_lengths[rownames(aml_dge), 'gene_length']

计算normalisation factors之后我们就可以将数据转换成RPKM/FPKM值了。只要特征数量(行)和样本数量(列)相同,SummarizedExperiment对象可以同时存储多层数据。因此我们直接将FPKM值存在已建立的SummarizedExperiment对象aml_se中。

aml_dge_tmm = calcNormFactors(aml_dge, method = 'TMM')

#compute FPKM values and append to assays
assay(aml_se, 'logFPKM_TMM') = rpkm(aml_dge_tmm, log = TRUE)
aml_se
## class: RangedSummarizedExperiment 
## dim: 17146 151 
## metadata(1): data_release
## assays(7): unstranded stranded_first ... fpkm_uq_unstrand logFPKM_TMM
## rownames(17146): ENSG00000000419 ENSG00000000457 ... ENSG00000288663
##   ENSG00000288670
## rowData names(12): source type ... gene_length gene_biotype
## colnames(151): TCGA-AB-2937-03A-01T-0734-13
##   TCGA-AB-2863-03A-01T-0734-13 ... TCGA-AB-2979-03B-01T-0760-13
##   TCGA-AB-2918-03A-01T-0740-13
## colData names(50): barcode patient ... released sample.aux

3.6 对样本注释突变信息

读者们需要注意的是,本分析流程中我们用了从原始TCGA AML文献中(Cancer Genome Atlas Research Network 2013)提取的突变列表(Supplemental Table 01 at https://gdc.cancer.gov/node/876)而不是标准TCGA流程中的突变信息(可以在National Cancer Institute Genomic Data Commons获得),两者存在一定差异。针对我们所关心的遗传病变(NPM1c, KMT2A-MLL, KMT2A-PTD and PML-RARA),病人通过以下几个标签被分类:

  • Patient ID: TCGA Patient ID
  • NPM1c: 当“NPM1”列包含有“p.W287fs”或“p.W288fs”时,返回TRUE
  • KMT2A-fusion: 当“MLL-partner”列包含有“MLL-”或“-MLL”时,返回TRUE(注意MLL的official gene symbol是KMT2A
  • KMT2A-PTD: 当“MLL-PTD”列包含有“exons”时,返回TRUE
  • PML-RARA: 当“PML-RARA”列包含有“PML-RARA”时,返回TRUE
data(AMLPatientMutationsTCGA)
patient_mutations = AMLPatientMutationsTCGA
patient_mutations = patient_mutations[colnames(aml_se), ] # order samples
aml_mutations = colnames(patient_mutations) # get mutation labels
colData(aml_se) = cbind(colData(aml_se), patient_mutations)
colData(aml_se)[, aml_mutations]
## DataFrame with 151 rows and 4 columns
##                              NPM1c.Mut KMT2A.Fusion KMT2A.PTD  PML.RARA
##                              <logical>    <logical> <logical> <logical>
## TCGA-AB-2937-03A-01T-0734-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2863-03A-01T-0734-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2941-03A-01T-0740-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2910-03A-01T-0740-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2856-03A-01T-0736-13     FALSE        FALSE     FALSE     FALSE
## ...                                ...          ...       ...       ...
## TCGA-AB-2994-03A-01T-0735-13     FALSE        FALSE     FALSE      TRUE
## TCGA-AB-2956-03A-01T-0740-13     FALSE         TRUE     FALSE     FALSE
## TCGA-AB-2955-03A-01T-0734-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2979-03B-01T-0760-13     FALSE        FALSE     FALSE     FALSE
## TCGA-AB-2918-03A-01T-0740-13     FALSE        FALSE     FALSE     FALSE

3.7 匹配Entrez IDs

Ensembl注释(Ensembl IDs)在基因组上有更高的覆盖度,因此适合用于识别突变(variant calling)等类似的分析过程(Zhao and Zhang 2015)。而RefSeq注释(Entrez IDs)对RNA-seq分析更适用,通常RNA-seq分析需要稳定的注释信息以便于未来的比较(Wu, Phan, and Wang 2013)。因此我们选择将Entrez IDs和Ensembl IDs匹配后去掉那些没有匹配到Entrez ID的基因。

匹配可以通过biomaRt的bioconductor R包完成,它可以提供最新的注释信息。匹配还可以通过一年更新两次的R包org.Hs.eg.db(???)完成,它提供了更加稳定的注释信息,提高了分析的可重现性。R包AnnotationDbi(Pagès et al. 2023)中的mapIds函数可以实现匹配功能。

library(org.Hs.eg.db)

rowData(aml_se)$entrezgene = mapIds(
  org.Hs.eg.db,
  keys = rownames(aml_se),
  keytype = 'ENSEMBL',
  column = 'ENTREZID',
  multiVals = 'asNA'
  )
gene_annot = rowData(aml_se)

匹配到多个Entrez ID的Ensembl IDs返回值为NAs,接下来我们就可以去掉这些行以保证一对一的匹配。类似的,我们还需要去掉匹配到多个Ensembl IDs的Entrez IDs。最后我们得到了Ensembl ID和Entrez ID一对一匹配的数据。

#select genes with mapped Entrez IDs
keep = !is.na(gene_annot$entrezgene)

#select genes with unique Entrez IDs
dup_entrez = gene_annot$entrezgene[duplicated(gene_annot$entrezgene)]
keep = keep & !gene_annot$entrezgene %in% dup_entrez

#Biotype of discarded genes (due to non-unique mapping)
head(sort(table(gene_annot[!keep, 'gene_biotype']), decreasing = TRUE), n = 10)
## 
##                             lncRNA               processed_pseudogene 
##                               1758                                511 
##                                TEC transcribed_unprocessed_pseudogene 
##                                231                                211 
##                     protein_coding   transcribed_processed_pseudogene 
##                                131                                 57 
##             unprocessed_pseudogene                           misc_RNA 
##                                 57                                 23 
##     transcribed_unitary_pseudogene                            Mt_tRNA 
##                                  8                                  6
#subset the data
aml_se = aml_se[keep, ]

4 转录组标签(signatures)预测突变状态

在这里,来自 Verhaak et al. (2005) 的标签(signature)被用来预测NPM1c突变的状态。我们需要量化标签(signature)内的基因和它们在样本中表达水平的一致性(concordance)。如标签(signature)内的上调基因(up-regulated genes)在样本中高表达,下调基因(down-regulated genes)在样本中低表达就会产生较高的分数。接下来,这个分数就可以被用来预测样本的突变状态。

首先,我们先将标签(signature)从MSigDB中下载下来并使用R包GSEABase(???)读入GeneSet对象中。接着我们用R/Bioconductor包singscore来给Verhaak标签(signature)打分,singscore包内的一些函数可用于后续的可视化与诊断。最后,我们对最终分数使用了逻辑回归模型(logistic regression model)来预测突变状态。

4.1 下载标签并载入

Verhaak et al. (2005) 的标签包含有上调和下调的基因列表。许多标签都是这种形式以增强其对样本的区分力。MSigDB将这样的标签(signatures)分成两部分并分别用后缀“_UP”和“_DN”标记。这里,我们直接将标签名“VERHAAK_AML_WITH_NPM1_MUTATED”和"_UP“或”_DN"连接以形成下载链接。

#create signature names
verhaak_names = paste('VERHAAK_AML_WITH_NPM1_MUTATED', c('UP', 'DN'), sep = '_')
verhaak_names
## [1] "VERHAAK_AML_WITH_NPM1_MUTATED_UP" "VERHAAK_AML_WITH_NPM1_MUTATED_DN"

利用刚才创造的链接下载后,"_UP“和”_DN“分别可以得到一个XML文件,在这里我们指定好输出XML文件的文件名(参数verhaak_files)。mapply函数可以用来循环下载类似的“链接-输出”参数对。

#generate URLs
verhaak_links = paste0(
  'http://software.broadinstitute.org/gsea/msigdb/download_geneset.jsp?geneSetName=',
  verhaak_names,
  '&fileType=xml'
  )

#download files
verhaak_files = paste0(verhaak_names, '.xml')
verhaak_path <- bfcrpath(bfc, verhaak_links) 

GSEABase包中的函数可以用来读取、解析和处理标签(signatures)。getBroadSets函数用来读取MSigDB XML文件中的标签(signatures)并生成一个GeneSet对象。XML文件还提供了来自原始实验(即HG-U133A)的Gene symbols,Entrez IDs和affymetrix chip IDs。我们仅读入了Entrez IDs因为它能直接和我们的数据匹配起来。GSEABase包中的mapIdentifiers函数可以帮助你进行ID转换,之所以用这个函数而不是AnnotationDbi包中的mapIds是因为这个函数在转换ID的同时保留了GeneSet对象。

library(GSEABase)

verhaak_sigs = getBroadSets(verhaak_path, membersId = 'MEMBERS_EZID')
verhaak_sigs
## GeneSetCollection
##   names: VERHAAK_AML_WITH_NPM1_MUTATED_UP, VERHAAK_AML_WITH_NPM1_MUTATED_DN (2 total)
##   unique identifiers: 10051, 10135, ..., 9828 (437 total)
##   types in collection:
##     geneIdType: SymbolIdentifier (1 total)
##     collectionType: BroadCollection (1 total)

为了使得分析过程中的索引更加方便,我们把SummarisedExperiment对象的行名改成Entrez IDs。

rownames(aml_se) = rowData(aml_se)$entrezgene

4.2 使用标签给样本打分

Singscore 是一种基于排名的,对基因集(gene set)在单个样本中富集程度的度量方法。对不同的基因集(gene set)的打分利用的样本基因表达量是相同的,因此我们只要根据表达量对样本进行一次排序就可以对不同的基因集(gene set)进行打分。在Singscore包中我们将这两步分开以节约计算资源。rankGenes函数可以对以numeric matrix、numeric data fame、ExpressionSet对象、DGEList对象或者SummarizedExperiment对象为存储形式的表达谱数据计算排名。用户需要指定排序过程中重复排名的处理方法(tiesMethod参数),默认方法为“min”,如有10个基因根据其表达量得到的排名均为1,则这10个基因的排名记为1,第11个基因的排名为11。我们推荐在RNA-seq数据里使用该方法,因为RNA-seq数据里通常有很多基因的表达量为0。这样可以减少零值对排名带来的影响,不过在数据质量控制时过滤掉低表达基因仍然是必不可少的步骤(见3.4部分)。

library(singscore)

#apply the rankGenes method to each version of the dataset, excluding counts
aml_ranked = rankGenes(assay(aml_se, 'logFPKM_TMM'))

Singscore针对不同的基因标签(gene signature)有三种计算模式。第一种模式适用于基因标签内有上调(up)和下调(down)基因列表。MSigDB中许多标签(signature),包括这里我们用的 Verhaak et al. (2005) 的标签(signature)都是这种形式的。这个模式下,上调和下调基因列表应分别传递给upSetdownSet。如果基因标签中的基因都是上调或者下调的,用户只要把基因列表传递给upSet参数。在这第二种模式下,下调的基因的得分会直接被倒转(如果得分已经被中心化,则直接取相反数,否则取“1-score”)。如果用户不确定标签(signature)内基因的组成,比如里面的基因可能是上调或者下调的,那么可以使用第三种模式。用户把基因列表指定给参数upSet后将knownDirection参数设定成FALSE即可。

scores默认被中心化,这样前两种模式得到的分数(scores)范围分别是\([-1, 1]\)\([-0.5, 0.5]\)。负的得分表明基因存在相反的富集,即预期上调的基因实际上在样本中是下调的,反之同理。第三种模式得到的分数(scores)无法被中心化,它的范围是\([0, 1]\)。在这个模式下,分数越高表明基因的排名离开中位数越远。如果我们对得分(scores)进行了中心化,负的得分并不能代表相反的富集因此并不适用。在这里,中心化只是为了更好地解释结果。

我们使用默认参数给给NPM1c突变标签(Verhaak Signature)打分,由于标签内含有上调和下调基因列表因此我们使用第一种模式。最后函数会返回一个data frame,里面有对上调基因列表、下调基因列表和所有列表的打分(scores)和离差(dispersion)。在这种模式下,所有列表的离差(dispersion)是上调基因列表和下调基因列表离差(dispersion)的平均值。函数会返回一个warning指明出现在标签(signature)内但不包括在表达谱数据内的基因名称/ID。

#apply the scoring function
verhaak_scores = simpleScore(aml_ranked,
                             upSet = verhaak_sigs[[1]],
                             downSet = verhaak_sigs[[2]])
## Warning in checkGenes(upSet, rownames(rankData)): 29 genes missing: 10265, 108,
## 10924, 11025, 11026, 1672, 200315, 2215, 3201, 3215, 3216, 3569, 3627, 3759,
## 50486, 6346, 6348, 6364, 643332, 6648, 6966, 717, 8337, 8351, 861, 8843, 9518,
## 9627, 9997
## Warning in checkGenes(downSet, rownames(rankData)): 28 genes missing: 10232,
## 10267, 2122, 221981, 2258, 23532, 24141, 25907, 2697, 28526, 3047, 3386, 3848,
## 3934, 4070, 445, 4680, 4681, 5457, 5790, 6091, 653067, 653145, 7102, 7441,
## 8277, 862, 8788

应该注意的是,singscores由两部分组成,即代表富集程度的分数(score)和对排名的离散程度的估计。在理想情况下,所有预期上调的基因都具有高表达,因此将基因从低表达到高表达排序后得到的值越大,这样的值应当位于分布的右端。Singscore旨在量化这种排名的分布,因此它计算了标签(signature)中基因排名的平均值和离差(dispersion)。平均值的计算和其它单样本打分方法类似,但是我们认为用平均值和离差(dispersion)两个数据来观察标签(signature)内基因排名的分布更为合理。默认且推荐的离差(dispersion)统计方法是具有非参数统计特性的绝对中位差(median absolute deviation,MAD)。其它方法还有四分差(inter-quartile range,IQR),用户可以将IQR函数作为参数传递给dispersionFun

head(verhaak_scores)
##                               TotalScore TotalDispersion       UpScore
## TCGA-AB-2937-03A-01T-0734-13 -0.23170496        5067.156 -0.1519431299
## TCGA-AB-2863-03A-01T-0734-13 -0.04380560        5736.921 -0.0469140870
## TCGA-AB-2941-03A-01T-0740-13  0.04910721        5651.301  0.0068156552
## TCGA-AB-2910-03A-01T-0740-13  0.12058917        5519.720 -0.0008802139
## TCGA-AB-2856-03A-01T-0736-13  0.19180029        4736.907  0.1337476636
## TCGA-AB-2874-03A-01T-0735-13 -0.06609889        6439.673 -0.0557252863
##                              UpDispersion    DownScore DownDispersion
## TCGA-AB-2937-03A-01T-0734-13     4097.165 -0.079761828       6037.147
## TCGA-AB-2863-03A-01T-0734-13     5132.020  0.003108487       6341.822
## TCGA-AB-2941-03A-01T-0740-13     5659.084  0.042291555       5643.517
## TCGA-AB-2910-03A-01T-0740-13     6682.819  0.121469388       4356.620
## TCGA-AB-2856-03A-01T-0736-13     4505.621  0.058052624       4968.193
## TCGA-AB-2874-03A-01T-0735-13     5992.669 -0.010373608       6886.677

4.3 对Verhaak标签的诊断

singscoreR包提供了一系列可视化工具来探索基因标签(gene signature)。比如,这些工具可以用来观察双向标签中上调基因列表和下调基因列表的重要程度,观察基因标签中单个基因对不同样本类别的区分能力,以及还可以用来探索最终得分和离差(dispersion)之间的关系。注释可以直接被叠加在图上,Singscore支持连续型和离散型注释。连续型注释可以以向量(vector)形式或者用字符串(string)指定注释位于data frame的某一列。

针对上调基因、下调基因或者所有基因,我们开始探索分数(score)和离差(dispersion)之间的关系。plotDispersion函数可以用来生成带有注释的图。注释可以是离散型或者连续型变量,如果注释被合并在最后的得分data frame中则直接指定列名即可。singscore包中所有的画图函数都可以通过指定isInteractiveTRUE而生成可交互的动态图。

#relative size of text in the figure
relSize = 1.2

#create annotation
mutated_gene = rep('Other', ncol(aml_se))
mutated_gene[aml_se$NPM1c.Mut] = 'NPM1c Mut'
mutated_gene[aml_se$KMT2A.Fusion | aml_se$KMT2A.PTD] = 'MLL Fusion/PTD'
p1 = plotDispersion(verhaak_scores, annot = mutated_gene, textSize = relSize)
p1
使用NPM1c signature上调和/或下调基因列表得到的分数可以将携带有NPM1c mutations或*MLL* fusions/PTDs的样本区分出来. 以Scores和基因集内基因排名的median absolute deviation (MAD)为轴作图,图上可以看出scores值接近0的样本点拥有更高的MADs值。

Figure 2: 使用NPM1c signature上调和/或下调基因列表得到的分数可以将携带有NPM1c mutations或MLL fusions/PTDs的样本区分出来
以Scores和基因集内基因排名的median absolute deviation (MAD)为轴作图,图上可以看出scores值接近0的样本点拥有更高的MADs值。

2显示下调的基因更能区分出具有NPM1c突变的样本。而且我们可以看到该标签(signature)还能够将MLLKMT2A)fusions和PTDs从其它样本中区分开来。在下调基因标签的打分中,具有NPM1c突变的样本分数最高,其次是MLLfusions和PTDs样本。在上调基因标签的打分中,分数也有这样的趋势,而且尽管分数的范围变大了,它区分不同类型样本的能力只发生了略微减弱。事实上,上调基因能够更好地区分出那些没有我们感兴趣的突变的样本(标记为“Other”),这些样本大多数具有负的得分。上调基因标签的打分为负分意味着这些基因在样本中的表达量低于中位数,即可能在样本中是下调的基因。观察离差(MAD)的趋势我们发现,分数越接近0值则离差越大。有三种情况会导致得分接近0值:1)标签中基因的表达量在中位数附近;2)标签中基因的表达量平均分布在表达谱两端;3)标签中基因的表达量均匀分布在整个表达谱中(可能性最大)。后两种情况会导致离差(dispersion)值很大。为了进一步弄清这个问题,我们选择了三个样本并对上调和下调基因标签中基因的排名分布进行作图。三个样本分别为总体得分最高的样本、总体得分最低的样本和总体离差(dispersion)值最高的样本。

library(gridExtra)
library(ggplot2)

#select samples with the properties required
select_samples = c(
  'Max Total Score' = which.max(verhaak_scores$TotalScore),
  'Min Total Score' = which.min(verhaak_scores$TotalScore),
  'Max Dispersion' = which.max(verhaak_scores$TotalDispersion)
  )

#plotRankDensity applied to each sample
rank_plots = lapply(names(select_samples), function(x) {
  #get the sample index
  aml_sample = select_samples[x]
  #drop = FALSE is required to ensure the data.frame is intact
  p1 = plotRankDensity(rankData = aml_ranked[, aml_sample, drop = FALSE],
                       upSet = verhaak_sigs[[1]],
                       downSet = verhaak_sigs[[2]],
                       textSize = relSize)
  
  #overwrite title of the plot with the description of the sample
  #this is possible because singscore uses ggplot2
  p1 = p1 + ggtitle(paste(x,  mutated_gene[aml_sample], sep = '\n')) +
    guides(colour = guide_legend(ncol = 1))
  return(p1)
})

#create a multipanel plot
grid.arrange(grobs = rank_plots, nrow = 1)
基因排名的分布显示了应该上调的基因在高分样本中上调,在低分样本中下调;应该下调的基因在高分样本中下调,在低分样本中上调. 从左到右分别是总体得分最高的样本、总体得分最低的样本和总体离差(dispersion)值最高的样本。barcode图显示了标签内每个基因的排名及其密度函数。上调基因的颜色为绿色而下调基因为紫色。在最高分和最低分样本中,排名两端呈现正态分布(Gaussian distributions)。总体离差(dispersion)最高的样本得分接近0,因为下调基因的排名呈现均匀分布而上调基因的排名呈现双峰分布。

Figure 3: 基因排名的分布显示了应该上调的基因在高分样本中上调,在低分样本中下调;应该下调的基因在高分样本中下调,在低分样本中上调
从左到右分别是总体得分最高的样本、总体得分最低的样本和总体离差(dispersion)值最高的样本。barcode图显示了标签内每个基因的排名及其密度函数。上调基因的颜色为绿色而下调基因为紫色。在最高分和最低分样本中,排名两端呈现正态分布(Gaussian distributions)。总体离差(dispersion)最高的样本得分接近0,因为下调基因的排名呈现均匀分布而上调基因的排名呈现双峰分布。

3显示上调和下调基因列表分别对总体得分最高的样本作出了贡献,这些基因分别位于排名分布的两端。类似地,在得分最低的样本中,上调基因位于排名分布的左端而下调基因位于排名分布的右端从而产生了很低的得分。正如前一张图中看到的那样,上调基因列表提高了NPM1c突变样本和other样本之间的区分度,可能相比下调基因列表能更好地指示野生型NPM1c样本。最后,离差(dispersion)值最高的样本中,下调基因的排名呈均匀分布而上调基因的排名呈双峰分布,双峰分别位于谱的两端。

综上,这些图显示了上调基因和下调基因列表在区分NPM1c突变、MLLfusion/PTD和野生型样本时都很重要。这些图能帮助用户在应用之前先确定标签(signature)的重要性,有时候标签(signature)是在特定情境下产生的,本身带有偏差因而只能在特定情况下使用。这些图还能帮助验证这些标签(signature)的适用情境。

4.4 用Verhaak标签预测突变

突变状态可以通过logit函数的逻辑回归模型(logistic regression model)进行预测。这个模型相对于将标签(signature)中每一个基因作为变量的方法而言更加简便, Verhaak et al. (2005) 的标签(signature)中的436个基因在逻辑回归中可以产生436个预测变量(predictor)。由于样本数少于预测变量(predictor)数,一些特征选取(feature selection)方法可能导致信息的丢失。此外,基因水平的模型在选择更多基因标签(gene signatures)上存在限制。而singscore的模型在某种程度上会继承它非参数的特性。比如该模型对于任何保留了基因排序的数据转换都是稳健的。该模型主要的不足在于丢失了一定的准确性,因为它将436个基因的信息限制到两个分类中。

不管怎样,我们的目的不是讨论哪个模型应该被用于预测突变状态,而是展示singscore和转录组标签(signature)区分样本的能力。因此我们在这里用训练模型的数据评价该模型的表现。我们先把分数和突变注释合并到一起。

library(reshape2)

#ggplot theme
rl = 1.2
current_theme = theme_minimal() +
  theme(
    panel.border = element_rect(colour = 'black', fill = NA),
    panel.grid.minor = element_blank(),
    axis.title = element_text(size = rel(rl) * 1.1),
    axis.text = element_text(size = rel(rl)),
    plot.title = element_text(size = rel(rl)),
    strip.background = element_rect(fill = NA, colour = 'black'),
    strip.text = element_text(size = rel(rl)),
    legend.text = element_text(size = rel(rl)),
    legend.title = element_text(size = rel(rl), face = 'italic'),
    legend.position = 'bottom',
    legend.direction = 'horizontal'
  )

#create a dataframe with the data required: scores and mutations
scoredf = as.data.frame(colData(aml_se)[, aml_mutations])
scoredf$Score = verhaak_scores$TotalScore
scoredf$Dispersion = verhaak_scores$TotalDispersion

在训练模型之前,我们可以先可视化一下,看 Verhaak et al. (2005) 的标签(signature)得到的分数区分突变型和野生型样本的能力。图4展示了分数可以区分开NPM1c突变型样本和野生型样本,MLL (KMT2A) fusions和PML-RARA fusions的情况也是如此。

#restructure the data for ploting
plotdf = melt(
  scoredf,
  id.var = c('Score', 'Dispersion'),
  variable.name = 'Mutation',
  value.name = 'Status'
  )
#convert TRUE-FALSE values to Mut-WT
plotdf$Status = factor(plotdf$Status, labels = c('WT', 'Mut'))
p1 = ggplot(plotdf, aes(Mutation, Score, fill = Status)) +
  geom_boxplot(position = 'dodge', alpha = 0.6) +
  scale_fill_brewer(palette = 'Set2') +
  current_theme +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
p1
NPM1c标签的得分可以区分突变样本与野生型样本. NPM1c标签得分的Boxplot(按照不同突变类型分类)。

Figure 4: NPM1c标签的得分可以区分突变样本与野生型样本
NPM1c标签得分的Boxplot(按照不同突变类型分类)。

为了量化以上的结果,我们将每个基因损伤(genomic lesion)作为因变量(response variable),样本的分数作为预测变量(predictor)进行逻辑回归并列出回归系数估计,标准差,z值和p值。

#fit GLMs to each mutation
glms = lapply(aml_mutations, function(x) {
  #generate a formula for the fit
  form = as.formula(paste0(x, ' ~ Score'))
  glm1 = glm(form, data = scoredf, family = binomial(link = 'logit'))
  return(glm1)
})
names(glms) = aml_mutations

#extract coefficients
coefs = lapply(glms, function(x) coef(summary(x)))
ldply(coefs, function(x) x[2, ], .id = 'Mutation')
##       Mutation  Estimate Std. Error    z value     Pr(>|z|)
## 1    NPM1c.Mut 13.068968   2.183133  5.9863351 2.146220e-09
## 2 KMT2A.Fusion  4.472692   2.177279  2.0542582 3.995071e-02
## 3    KMT2A.PTD  1.675708   2.111781  0.7935046 4.274839e-01
## 4     PML.RARA -6.492618   2.423885 -2.6785997 7.393071e-03

NPM1c突变和分数显著相关,MLL gene fusion也和分数显著相关,这显示了它们在Hox基因失调中共有的作用。有趣的是,PML-RARA fusion的系数是负值,这很可能反映了急性早幼粒细胞白血病(acute promyelocytic leukemia, APL)相对于其它AML亚型而言独特的细胞形态或表型特征。

4.5 评估预测模型的表现

以上统计数字只提供了模型本身的信息而非模型的表现好坏。精确率(precision)和召回率(recall)可以衡量模型表现的好坏。dcanr(Bhuva 2023)提供了函数进行这类计算。

library(dcanr)

#assess sensitivity and specificity
prec_rec = ldply(glms, function(glm1) {
  #predict mutations for the data used in training and convert to binary format
  prediction = as.numeric(predict(glm1) > 0)
  observed = glm1$y
  prec = performanceMeasure(prediction, observed, 'precision')
  recall = performanceMeasure(prediction, observed, 'recall')
  f1 = performanceMeasure(prediction, observed)
  return(c('Precision' = prec, 'Recall' = recall, 'F1' = f1))
}, .id = 'Mutation')
prec_rec
##       Mutation Precision    Recall        F1
## 1    NPM1c.Mut       0.7 0.5675676 0.6268657
## 2 KMT2A.Fusion       NaN 0.0000000 0.0000000
## 3    KMT2A.PTD       NaN 0.0000000 0.0000000
## 4     PML.RARA       NaN 0.0000000 0.0000000

这里我们只能计算对NPM1c突变预测的精确率(precision),具有其它突变类型的样本在此情景下的预测结果是野生型(wild-types),因此无法计算它们的精确率。同理,具有其他突变类型的样本的召回率(recall)是0。NPM1c突变的模型的精确率(precision)、召回率(recall)、F1-score和预期的一样高。就像4.2部分提到的,singscores是由两部分组成的分数,因此在考虑排名的离差(dispersion)后,模型表现应该更加优秀(图2)。

以下的模型显示,使用了singscores中的score及离差(dispersion)后,模型的预测能力显著提升。

#include dispersion in the model
glm_npm1c = glm('NPM1c.Mut ~ Score + Dispersion',
               data = scoredf,
               family = binomial(link = 'logit'))

#evaluate performance of the new model
prediction = as.numeric(predict(glm_npm1c) > 0)
observed = glm_npm1c$y
c(
  'Precision' = performanceMeasure(prediction, observed, 'precision'),
  'Recall' = performanceMeasure(prediction, observed, 'recall'),
  'F1' = performanceMeasure(prediction, observed)
  )
## Precision    Recall        F1 
## 0.7777778 0.7567568 0.7671233

5 多个基因标签(signature)的landscapes

通常我们会对两个独立的或者相互依赖的表型之间的关系感兴趣,比如EMT过程中epithelial和mesenchymal表型。大多数标签(signature)的作用是为了便于量化一些难以衡量的分子表型。因此,我们可以用对应的标签(signature)来探索两个表型之间的关系。 Foroutan et al. (2017) 首次介绍了使用molecular signature landscapes来探索标签(signature)之间的关系,这些标签(signature)均和EMT及TGF\(\beta\)导致的EMT相关。之后, Cursons et al. (2018) 计算了EMT表型标签(signature)的singscores得分并用signature landscape证明了将miR-200c转染到mesenchymal细胞系之后发生了向epithelial表型的转变。 Foroutan et al. (2018) 使用signature landscapes在epithelial-mesenchymal尺度上划分乳腺癌亚型并把它包含在了singscore包中。这里,我们展示了signature landscapes的进一步应用,即使用代表不同突变的转录组标签(signature)来划分AML样本。

5.1 Ross MLL fusion signature vs. Verhaak signature landscape

我们现在使用 Ross et al. (2004) 的MLL-fusion标签(signatures)来给TCGA AML样本打分。不像NPM1c标签,这个标签中的基因可以用来区分具有MLL-fusion基因的样本。我们按照4.1中的流程下载并处理标签。

#create signature names
rossmll_name = 'ROSS_AML_WITH_MLL_FUSIONS'
#generate URLs
rossmll_link = paste0(
  'http://software.broadinstitute.org/gsea/msigdb/download_geneset.jsp?geneSetName=',
  rossmll_name,
  '&fileType=xml'
  )

#download files
rossmll_file = paste0(rossmll_name, '.xml')
rossmll_path <- bfcrpath(bfc, rossmll_link) 
rossmll_sig = getBroadSets(rossmll_path, membersId = 'MEMBERS_EZID')
rossmll_sig
## GeneSetCollection
##   names: ROSS_AML_WITH_MLL_FUSIONS (1 total)
##   unique identifiers: 10113, 10479, ..., 9961 (80 total)
##   types in collection:
##     geneIdType: SymbolIdentifier (1 total)
##     collectionType: BroadCollection (1 total)

该标签包含有上调和下调基因以区分突变型和野生型。我们用4.2部分提到过的singscore的第三种模式进行打分,该模式不需要提前知道或确定基因是上调还是下调的。simpleScore函数里的knownDirection参数应当设定成FALSE4.2部分已经计算过的基因排名可以在这里再次使用,以计算新的标签的分数。

rossmll_scores = simpleScore(aml_ranked, rossmll_sig[[1]], knownDirection = FALSE)

我们可以用plotScoreLandscape函数画hexbin图来可视化分数的分布。使用 Verhaak et al. (2005) 标签及 Ross et al. (2004) MLL-fusions标签得到的分数被传递给了该函数。两个分数应该基于相同的样本计算后得到。这两个分数所对应的名字应当传递给scorenames参数。textSize参数可以用来指定图中文字相对于图片的大小,这在应对不同poster、publications和presentations的要求时非常有用。

p_mll_npm1c = plotScoreLandscape(
  verhaak_scores,
  rossmll_scores,
  scorenames = c('VERHAAK_AML_WITH_NPM1_MUTATED', 'ROSS_AML_WITH_MLL_FUSIONS'),
  textSize = relSize
  )
p_mll_npm1c
MLL fusion和NPM1c标签的landscape. AML中两个标签存在正相关。

Figure 5: MLL fusion和NPM1c标签的landscape
AML中两个标签存在正相关。

5显示,尽管两个标签之间只有17个共同基因(标签大小分别为81和436个基因),两个标签的分数之间仍具有很强的正相关(Spearman’s \(\rho\) = 0.593)。在这样的分析中,我们可能想像 Cursons et al. (2018) 文章里那样把新的数据点投影在图上。或者我们可能希望将已有的数据点投影在图上以观察样本的分类情况。这里我们将MLL fusions、MLL PTDs、PML-RARA fusions和NPM1c mutations四类样本的数据点进行投影。首先我们需要先建立注释。

#new annotation - modify previously used annotations
mutated_gene[aml_se$KMT2A.Fusion] = 'MLL Fusion'
mutated_gene[aml_se$KMT2A.PTD] = 'MLL PTD'
mutated_gene[aml_se$PML.RARA] = 'PML-RARA'

projectScoreLandscape函数可以用来将数据点投影在已有的图上。这个函数使用plotScoreLandscape函数产生的p_mll_npm1c变量并将新的数据点投影上去。计算新数据点分数的基因标签必须和已有图中使用的标签相同。在这里我们直接用之前已经计算好的分数来画新的数据点。subSamples参数用于选择想要投影的数据点子集,这里我们选择了感兴趣的四种样本。

select_aml = !mutated_gene %in% 'Other'
#project above mutations onto the landscape
p1 = projectScoreLandscape(p_mll_npm1c,
                           verhaak_scores,
                           rossmll_scores,
                           subSamples = select_aml,
                           annot = mutated_gene[select_aml])
p1 + theme(legend.box = 'vertical')
Landscape展示了MLL fusions和PTDs之间有很大不同. 不同的突变在landscape上占据了不同的位置,提示了它们具有不同的分子特征。

Figure 6: Landscape展示了MLL fusions和PTDs之间有很大不同
不同的突变在landscape上占据了不同的位置,提示了它们具有不同的分子特征。

6显示MLL fusions和MLL PTDs两组样本在两个维度上被分开。在MLL fusion标签的方向上,MLL PTDs相比MLL fusions的得分更低。这些样本在两个维度的平面上没有形成cluster,这和 Ross et al. (2004) 的结论是一致的。在图上进行样本注释的投射有助于帮助我们解释landscape不同部分的意义。

5.2 Ross PML-RARA fusion signature vs. Verhaak signature landscape

之前的分析显示PML-RARA样本和其它的样本有明显不同。这里我们用来自 Ross et al. (2004)PML-RARA标签重复了刚才的分析以验证这种不同。该标签在MSigDB中的名字为“ROSS_AML_WITH_PML_RARA_FUSION”。我们下载该标签并对所有样本进行打分,用该分数和 Verhaak et al. (2005) 标签得到的分数画landscape图。最后我们将样本投影到图上(见5.1部分)。该标签的和MLL fusions标签来自同一个研究结果因此我们仍然使用相同的参数进行打分。

PML-RARA标签和NPM1c标签之间没有关联. L型的landscape提示了这两个突变背后的分子机制是互斥的。

Figure 7: PML-RARA标签和NPM1c标签之间没有关联
L型的landscape提示了这两个突变背后的分子机制是互斥的。

7展现了和MLL fusion标签完全不同的landscape。PML-RARA标签使得PML-RARA和NPM1c样本完全分开,PML-RARA样本是唯一在PML-RARA标签中得到高分的样本,而且两个标签之间没有显著相关。在样本背景介绍部分,PML-RARA fusion是AML中一种已知的亚型——急性早幼粒细胞白血病(acute prolmyelocytic leukemia,APL)的诊断标志,这种亚型独特的细胞表型反映了早幼粒细胞分化的阻滞(Thé et al. 1991)。这种特性反映在了L型的landscape中,以及这些基因组损伤驱动白血病不同的机制。

6 总结

singscore包通过R/Bioconductor环境给使用者提供了一个用户友好的界面来进行基因集打分(gene set scoring)。 TCGABiolinks包使用户相对容易地访问大型的临床相关数据集,如TCGA及其相关注释。 singscore中包含的诊断和绘图功能允许用户调查感兴趣的基因集(gene set)以确定它们区分样本之间差异的能力。然后不同的基因集(gene set)可以组合在一起以探索不同的细胞表型在大型研究中的关系。以上分析展示了当恰当的基因集被使用时,singscore计算得到的结果可以用于样本分类,结果可以进一步被用到缺少基因组数据的大型转录组数据中。

使用到的R包

此工作流程用到了许多Bioconductor项目3.19版本中的包,适用于R Under development (unstable) (2023-10-22 r85388)或更高版本。本工作流程中用到的所有包已列在下表中:

sessionInfo()
## R Under development (unstable) (2023-10-22 r85388)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] BiocFileCache_2.11.1         dbplyr_2.4.0                
##  [3] dcanr_1.19.0                 gridExtra_2.3               
##  [5] reshape2_1.4.4               singscore_1.23.0            
##  [7] GSEABase_1.65.0              graph_1.81.0                
##  [9] annotate_1.81.0              XML_3.99-0.14               
## [11] org.Hs.eg.db_3.18.0          AnnotationDbi_1.65.0        
## [13] plyr_1.8.9                   rtracklayer_1.63.0          
## [15] edgeR_4.1.0                  limma_3.59.0                
## [17] SummarizedExperiment_1.33.0  Biobase_2.63.0              
## [19] GenomicRanges_1.55.0         GenomeInfoDb_1.39.0         
## [21] IRanges_2.37.0               S4Vectors_0.41.1            
## [23] BiocGenerics_0.49.0          MatrixGenerics_1.15.0       
## [25] matrixStats_1.0.0            TCGAbiolinks_2.31.0         
## [27] ggplot2_3.4.4                SingscoreAMLMutations_1.19.0
## [29] BiocStyle_2.31.0            
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3          jsonlite_1.8.7             
##   [3] magrittr_2.0.3              magick_2.8.1               
##   [5] farver_2.1.1                rmarkdown_2.25             
##   [7] BiocIO_1.13.0               zlibbioc_1.49.0            
##   [9] vctrs_0.6.4                 memoise_2.0.1              
##  [11] Rsamtools_2.19.0            RCurl_1.98-1.12            
##  [13] htmltools_0.5.6.1           S4Arrays_1.3.0             
##  [15] progress_1.2.2              curl_5.1.0                 
##  [17] SparseArray_1.3.0           sass_0.4.7                 
##  [19] bslib_0.5.1                 cachem_1.0.8               
##  [21] GenomicAlignments_1.39.0    igraph_1.5.1               
##  [23] lifecycle_1.0.3             iterators_1.0.14           
##  [25] pkgconfig_2.0.3             Matrix_1.6-1.1             
##  [27] R6_2.5.1                    fastmap_1.1.1              
##  [29] GenomeInfoDbData_1.2.11     digest_0.6.33              
##  [31] colorspace_2.1-0            RSQLite_2.3.1              
##  [33] labeling_0.4.3              filelock_1.0.2             
##  [35] fansi_1.0.5                 httr_1.4.7                 
##  [37] abind_1.4-5                 compiler_4.4.0             
##  [39] rngtools_1.5.2              bit64_4.0.5                
##  [41] withr_2.5.1                 downloader_0.4             
##  [43] BiocParallel_1.37.0         DBI_1.1.3                  
##  [45] hexbin_1.28.3               R.utils_2.12.2             
##  [47] biomaRt_2.59.0              rappdirs_0.3.3             
##  [49] DelayedArray_0.29.0         rjson_0.2.21               
##  [51] tools_4.4.0                 R.oo_1.25.0                
##  [53] glue_1.6.2                  restfulr_0.0.15            
##  [55] grid_4.4.0                  generics_0.1.3             
##  [57] BSgenome_1.71.0             gtable_0.3.4               
##  [59] tzdb_0.4.0                  R.methodsS3_1.8.2          
##  [61] tidyr_1.3.0                 data.table_1.14.8          
##  [63] hms_1.1.3                   xml2_1.3.5                 
##  [65] utf8_1.2.4                  XVector_0.43.0             
##  [67] ggrepel_0.9.4               foreach_1.5.2              
##  [69] pillar_1.9.0                stringr_1.5.0              
##  [71] dplyr_1.1.3                 lattice_0.22-5             
##  [73] bit_4.0.5                   tidyselect_1.2.0           
##  [75] locfit_1.5-9.8              Biostrings_2.71.1          
##  [77] knitr_1.44                  bookdown_0.36              
##  [79] xfun_0.40                   statmod_1.5.0              
##  [81] stringi_1.7.12              yaml_2.3.7                 
##  [83] TCGAbiolinksGUI.data_1.23.0 evaluate_0.22              
##  [85] codetools_0.2-19            tibble_3.2.1               
##  [87] BiocManager_1.30.22         cli_3.6.1                  
##  [89] xtable_1.8-4                munsell_0.5.0              
##  [91] jquerylib_0.1.4             Rcpp_1.0.11                
##  [93] png_0.1-8                   parallel_4.4.0             
##  [95] readr_2.1.4                 blob_1.2.4                 
##  [97] prettyunits_1.2.0           mclust_6.0.0               
##  [99] doRNG_1.8.6                 bitops_1.0-7               
## [101] scales_1.2.1                purrr_1.0.2                
## [103] crayon_1.5.2                rlang_1.1.1                
## [105] KEGGREST_1.43.0             rvest_1.0.3

致谢

工作流程中的结果全部或部分基于TCGA Research Network(https://www.cancer.gov/tcga)产生的数据。感谢来自The Walter and Eliza Hall Institute的Christoffer Flensburg在TCGA AML mutation data上的帮助。

参考文献

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2023. Rmarkdown: Dynamic Documents for R. https://github.com/rstudio/rmarkdown.

Barbie, David A, Pablo Tamayo, Jesse S Boehm, So Young Kim, Susan E Moody, Ian F Dunn, Anna C Schinzel, et al. 2009. “Systematic Rna Interference Reveals That Oncogenic Kras-Driven Cancers Require Tbk1.” Nature 462 (7269): 108.

Bhuva, Dharmesh D. 2023. Dcanr: Differential Co-Expression/Association Network Analysis. https://doi.org/10.18129/B9.bioc.dcanr.

Brunetti, Lorenzo, Michael C Gundry, Daniele Sorcini, Anna G Guzman, Yung-Hsin Huang, Raghav Ramabadran, Ilaria Gionfriddo, et al. 2018. “Mutant Npm1 Maintains the Leukemic State Through Hox Expression.” Cancer Cell 34 (3): 499–512.

Cancer Genome Atlas Research Network. 2013. “Genomic and Epigenomic Landscapes of Adult de Novo Acute Myeloid Leukemia.” New England Journal of Medicine 368 (22): 2059–74.

Chen, Yunshun, Aaron TL Lun, and Gordon K Smyth. 2016. “From Reads to Genes to Pathways: Differential Expression Analysis of Rna-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline.” F1000Research 5.

Cieślik, Marcin, and Arul M Chinnaiyan. 2018. “Cancer Transcriptome Profiling at the Juncture of Clinical Translation.” Nature Reviews Genetics 19 (2): 93.

Colaprico, Antonio, Tiago C Silva, Catharina Olsen, Luciano Garofano, Claudia Cava, Davide Garolini, Thais S Sabedot, et al. 2015. “TCGAbiolinks: An R/Bioconductor Package for Integrative Analysis of Tcga Data.” Nucleic Acids Research 44 (8): e71–e71.

Cursons, Joseph, Katherine A Pillman, Kaitlin G Scheer, Philip A Gregory, Momeneh Foroutan, Soroor Hediyeh-Zadeh, John Toubia, et al. 2018. “Combinatorial Targeting by microRNAs Co-Ordinates Post-Transcriptional Control of Emt.” Cell Systems 7 (1): 77–91.

Foroutan, Momeneh, Dharmesh D Bhuva, Ruqian Lyu, Kristy Horan, Joseph Cursons, and Melissa J Davis. 2018. “Single Sample Scoring of Molecular Phenotypes.” BMC Bioinformatics 19 (1): 404.

Foroutan, Momeneh, Joseph Cursons, Soroor Hediyeh-Zadeh, Erik W Thompson, and Melissa J Davis. 2017. “A Transcriptional Program for Detecting Tgf\(\beta\)-Induced Emt in Cancer.” Molecular Cancer Research 15 (5): 619–31.

Hänzelmann, Sonja, Robert Castelo, and Justin Guinney. 2013. “GSVA: Gene Set Variation Analysis for Microarray and Rna-Seq Data.” BMC Bioinformatics 14 (1): 7.

Hess, Jay L. 2004. “MLL, Hox Genes, and Leukemia: The Plot Thickens.” Blood 103 (8): 2870–1.

Law, Charity W, Monther Alhamdoosh, Shian Su, Gordon K Smyth, and Matthew E Ritchie. 2016. “RNA-Seq Analysis Is Easy as 1-2-3 with Limma, Glimma and edgeR.” F1000Research 5.

Lee, Eunjung, Han-Yu Chuang, Jong-Won Kim, Trey Ideker, and Doheon Lee. 2008. “Inferring Pathway Activity Toward Precise Disease Classification.” PLoS Computational Biology 4 (11): e1000217.

Liberzon, Arthur, Chet Birger, Helga Thorvaldsdóttir, Mahmoud Ghandi, Jill P Mesirov, and Pablo Tamayo. 2015. “The Molecular Signatures Database Hallmark Gene Set Collection.” Cell Systems 1 (6): 417–25.

Ma, Xiaotu, Yu Liu, Yanling Liu, Ludmil B Alexandrov, Michael N Edmonson, Charles Gawad, Xin Zhou, et al. 2018. “Pan-Cancer Genome and Transcriptome Analyses of 1,699 Paediatric Leukaemias and Solid Tumours.” Nature 555 (7696): 371.

Oleś, Andrzej. 2023. BiocStyle: Standard Styles for Vignettes and Other Bioconductor Documents. https://doi.org/10.18129/B9.bioc.BiocStyle.

Oshlack, Alicia, and Matthew J Wakefield. 2009. “Transcript Length Bias in Rna-Seq Data Confounds Systems Biology.” Biology Direct 4 (1): 14.

Pagès, Hervé, Marc Carlson, Seth Falcon, and Nianhua Li. 2023. AnnotationDbi: Manipulation of Sqlite-Based Annotations in Bioconductor. https://doi.org/10.18129/B9.bioc.AnnotationDbi.

Papaemmanuil, Elli, Moritz Gerstung, Lars Bullinger, Verena I Gaidzik, Peter Paschka, Nicola D Roberts, Nicola E Potter, et al. 2016. “Genomic Classification and Prognosis in Acute Myeloid Leukemia.” New England Journal of Medicine 374 (23): 2209–21.

Parker, Joel S, Michael Mullins, Maggie CU Cheang, Samuel Leung, David Voduc, Tammi Vickery, Sherri Davies, et al. 2009. “Supervised Risk Predictor of Breast Cancer Based on Intrinsic Subtypes.” Journal of Clinical Oncology 27 (8): 1160.

Ross, Mary E, Rami Mahfouz, Mihaela Onciu, Hsi-Che Liu, Xiaodong Zhou, Guangchun Song, Sheila A Shurtleff, et al. 2004. “Gene Expression Profiling of Pediatric Acute Myelogenous Leukemia.” Blood 104 (12): 3679–87.

Thé, Hugues de, Catherine Lavau, Agnès Marchio, Christine Chomienne, Laurent Degos, and Anne Dejean. 1991. “The Pml-Rar\(\alpha\) Fusion mRNA Generated by the T (15; 17) Translocation in Acute Promyelocytic Leukemia Encodes a Functionally Altered Rar.” Cell 66 (4): 675–84.

Tomfohr, John, Jun Lu, and Thomas B Kepler. 2005. “Pathway Level Analysis of Gene Expression Using Singular Value Decomposition.” BMC Bioinformatics 6 (1): 225.

Verhaak, Roel GW, Chantal S Goudswaard, Wim van Putten, Maarten A Bijl, Mathijs A Sanders, Wendy Hugens, André G Uitterlinden, et al. 2005. “Mutations in Nucleophosmin (Npm1) in Acute Myeloid Leukemia (Aml): Association with Other Gene Abnormalities and Previously Established Gene Expression Signatures and Their Favorable Prognostic Significance.” Blood 106 (12): 3747–54.

Wickham, Hadley. 2020. Reshape2: Flexibly Reshape Data: A Reboot of the Reshape Package. https://github.com/hadley/reshape.

———. 2023. Plyr: Tools for Splitting, Applying and Combining Data. http://had.co.nz/plyr.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2023. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://ggplot2.tidyverse.org.

Wu, Po-Yen, John H Phan, and May D Wang. 2013. “Assessing the Impact of Human Genome Annotation Choice on Rna-Seq Expression Estimates.” Bmc Bioinformatics 14 (11): S8.

Xie, Yihui. 2023. Knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.

Zhao, Shanrong, and Baohong Zhang. 2015. “A Comprehensive Evaluation of Ensembl, Refseq, and Ucsc Annotations in the Context of Rna-Seq Read Mapping and Gene Quantification.” BMC Genomics 16 (1): 97.