使用singscore从transcriptomic signatures预测AML中的突变
Oct 2019
Abstract
针对生物转录组的RNA测序（RNA-seq）技术的进步革命性地提高了我们掌握转录调控机制的能力，这些转录调控机制可以帮助我们了解癌症等疾病的机理。最近我们发表了singscore，一种基于基因表达量排名的基因集打分（gene set scoring）方法，该方法量化了单个样本的转录谱与特定的基因集所代表的特征之间的一致性。在这里，我们展示了singscore如何应用于研究不同的急性髓细胞白血病（acute myeloid leukemia，AML）转录谱，这些转录谱与特定的几种突变和基因损伤相关。使用TCGA的基因组与转录组数据，我们展示了对恰当的标签（signature）打分可以区分出有不同突变的样本，反映出这些突变具有驱动异常的转录调控以导致白血病发生的能力。我们认为singscore方法对于研究特定癌症亚型内的异质性特别有用，并且我们展示了singscore识别那些驱动异常转录机制的突变的能力。
R version: R version 3.6.1 (2019-07-05)
Bioconductor version: 3.10
Package version: 1.2.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, Cieślik and Chinnaiyan (2018))。许多癌症的亚型分类主要依赖于识别大型患者队列中的复发性突变（recurrent mutations），通过全基因组或全外显子组测序来获得具有临床意义的亚型(Cancer Genome Atlas Research Network 2013, Papaemmanuil et al. (2016))。
流行的“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, Ross et al. (2004))。然而，具有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.org 和 firebrowse.org 找到。GDC的数据使用了STAR的“two-pass”模式进行序列比对，使用了HTSeq进行定量。用户在使用 GDC data transfer tool 下载GDC中的数据时，可以通过GDC portal选择自己感兴趣的特定文件进行下载。下载之后需要读取并整合这些文件之后才可以进行下游分析。以上这些步骤（包括下载）可以通过R包TCGAbiolinks (Colaprico et al. 2015, Colaprico et al. (2019))完成。这个包支持使用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] "b18b2385b1e916597856067dc6437f3c20b46bca"
##
## $data_release
## [1] "Data Release 19.0 - September 17, 2019"
##
## $status
## [1] "OK"
##
## $tag
## [1] "1.22.0"
##
## $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 = 'HTSeq - Counts'
)
#extract results of the query
rnaseq_res = getResults(query_rna)
dim(rnaseq_res)## [1] 151 28colnames(rnaseq_res)## [1] "data_release" "data_type"
## [3] "updated_datetime" "file_name"
## [5] "submitter_id" "file_id"
## [7] "file_size" "cases"
## [9] "id" "created_datetime"
## [11] "md5sum" "data_format"
## [13] "access" "state"
## [15] "version" "data_category"
## [17] "type" "experimental_strategy"
## [19] "project" "analysis_id"
## [21] "analysis_updated_datetime" "analysis_created_datetime"
## [23] "analysis_submitter_id" "analysis_state"
## [25] "analysis_workflow_link" "analysis_workflow_type"
## [27] "analysis_workflow_version" "tissue.definition"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: 39MB)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对象中包含了56925个特征（features）和151个样本。原始数据包含60483个特征，其中3881个特征因无法和ENSEMBL GRCh38.p12匹配而被丢弃。用rowData(se)和colData(se)函数可以获得特征和样本的注释信息，用assay(se)函数可以获得counts数据。TCGA数据往往包含一些福尔马林固定、石蜡包埋（formalin-fixed paraffin-embedded，FFPE）的样本，这些样本需要被过滤掉以避免protocol不同而引入的误差。这一步在本数据集中不需要，因为这一步骤是针对实体瘤（solid tumors）而非白血病（leukemias）的。
aml_se## class: RangedSummarizedExperiment
## dim: 56512 151
## metadata(1): data_release
## assays(1): HTSeq - Counts
## rownames(56512): ENSG00000000003 ENSG00000000005 ... ENSG00000281912
## ENSG00000281920
## rowData names(3): ensembl_gene_id external_gene_name
## original_ensembl_gene_id
## colnames(151): TCGA-AB-3001-03A-01T-0736-13
## TCGA-AB-2853-03A-01T-0734-13 ... TCGA-AB-2977-03B-01T-0760-13
## TCGA-AB-2995-03A-01T-0735-13
## colData names(61): sample patient ... name is_ffpe3.4 过滤掉counts数低的基因
edgeR包提供了过滤所需的数据标准化与转换方法。这些方法要求数据存储在DGEList对象中，因此我们需要将SummarizedExperiment对象转换成DGEList对象。
library(SummarizedExperiment)
library(edgeR)
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)
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进行序列比对并用HTSeq进行定量（流程细节见https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/）。HTSeq统计了比对到每个基因外显子区域的read数量，因此有效基因长度应当是基因的所有外显子区域长度之和。GENCODE v22被用于定量过程中的注释，因此在计算基因长度时仍然应该用这个注释文件。
#download v22 of the GENCODE annotation
gencode_file = 'gencode.v22.annotation.gtf.gz'
gencode_link = paste(
'ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_22',
gencode_file,
sep = '/'
)
download.file(gencode_link, gencode_file, method = 'libcurl') #(size: 39MB)rtracklayer这个R包提供了解析GTF文件的函数。
library(rtracklayer)
library(plyr)
gtf = import.gff(gencode_file, 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_HTSeq_gencodev22.rds')
gene_lengths## DataFrame with 60483 rows and 3 columns
## ensembl_gene_id gene_type gene_length
## <character> <character> <integer>
## 1 ENSG00000000003 protein_coding 4535
## 2 ENSG00000000005 protein_coding 1610
## 3 ENSG00000000419 protein_coding 1207
## 4 ENSG00000000457 protein_coding 6883
## 5 ENSG00000000460 protein_coding 5967
## ... ... ... ...
## 60479 ENSGR0000275287 misc_RNA 290
## 60480 ENSGR0000276543 miRNA 68
## 60481 ENSGR0000277120 miRNA 64
## 60482 ENSGR0000280767 lincRNA 515
## 60483 ENSGR0000281849 antisense 484SummarizedExperiment对象可以存储基因的注释信息，因此我们需要把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: 17412 151
## metadata(1): data_release
## assays(2): HTSeq - Counts logFPKM_TMM
## rownames(17412): ENSG00000000419 ENSG00000000457 ... ENSG00000281772
## ENSG00000281896
## rowData names(5): ensembl_gene_id external_gene_name
## original_ensembl_gene_id gene_length gene_biotype
## colnames(151): TCGA-AB-3001-03A-01T-0736-13
## TCGA-AB-2853-03A-01T-0734-13 ... TCGA-AB-2977-03B-01T-0760-13
## TCGA-AB-2995-03A-01T-0735-13
## colData names(61): sample patient ... name is_ffpe3.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-3001-03A-01T-0736-13 FALSE FALSE FALSE TRUE
## TCGA-AB-2853-03A-01T-0734-13 TRUE FALSE FALSE FALSE
## TCGA-AB-2929-03A-01T-0735-13 FALSE FALSE FALSE FALSE
## TCGA-AB-2939-03A-01T-0740-13 FALSE FALSE FALSE FALSE
## TCGA-AB-2982-03B-01T-0748-13 FALSE FALSE FALSE TRUE
## ... ... ... ... ...
## TCGA-AB-2815-03A-01T-0734-13 FALSE FALSE FALSE FALSE
## TCGA-AB-3000-03A-01T-0736-13 FALSE FALSE FALSE FALSE
## TCGA-AB-2919-03A-01T-0740-13 TRUE FALSE FALSE FALSE
## TCGA-AB-2977-03B-01T-0760-13 FALSE FALSE TRUE FALSE
## TCGA-AB-2995-03A-01T-0735-13 FALSE FALSE FALSE FALSE3.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(Carlson 2019)完成，它提供了更加稳定的注释信息，提高了分析的可重现性。R包AnnotationDbi(Pagès et al. 2019)中的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)##
## processed_pseudogene antisense
## 1014 895
## lincRNA TEC
## 441 263
## sense_intronic protein_coding
## 249 142
## unprocessed_pseudogene transcribed_unprocessed_pseudogene
## 127 95
## transcribed_processed_pseudogene processed_transcript
## 71 69#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(Morgan, Falcon, and Gentleman 2019)读入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')
mapply(download.file, verhaak_links, verhaak_files, method = 'libcurl')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_files, 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 (435 total)
## types in collection:
## geneIdType: SymbolIdentifier (1 total)
## collectionType: BroadCollection (1 total)为了使得分析过程中的索引更加方便，我们把SummarisedExperiment对象的行名改成Entrez IDs。
rownames(aml_se) = rowData(aml_se)$entrezgene4.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）都是这种形式的。这个模式下，上调和下调基因列表应分别传递给upSet和downSet。如果基因标签中的基因都是上调或者下调的，用户只要把基因列表传递给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)): 24 genes missing: 10265,
## 108, 10924, 11025, 11026, 1672, 200315, 2212, 2215, 3215, 3216, 3569, 3627,
## 3759, 50486, 6346, 6364, 6967, 8337, 861, 8843, 9518, 9627, 9997## Warning in checkGenes(downSet, rownames(rankData)): 30 genes missing: 10232,
## 10267, 11217, 2122, 221981, 2258, 23532, 24141, 25907, 2697, 28526, 28638,
## 3047, 3386, 3848, 3934, 4070, 445, 4680, 4681, 5457, 5625, 5790, 6091,
## 653067, 653145, 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-3001-03A-01T-0736-13 -0.09195329 5277.315 -0.0763509783
## TCGA-AB-2853-03A-01T-0734-13 0.19484119 4906.294 0.0354452471
## TCGA-AB-2929-03A-01T-0735-13 -0.07368493 5444.478 0.0008835206
## TCGA-AB-2939-03A-01T-0740-13 -0.09941877 5729.508 -0.0594909905
## TCGA-AB-2982-03B-01T-0748-13 -0.11059695 5479.319 -0.0983332062
## TCGA-AB-2813-03A-01T-0736-13 0.15285753 4675.008 0.1387749302
## UpDispersion DownScore DownDispersion
## TCGA-AB-3001-03A-01T-0736-13 4511.552 -0.01560231 6043.078
## TCGA-AB-2853-03A-01T-0734-13 5876.285 0.15939594 3936.303
## TCGA-AB-2929-03A-01T-0735-13 5166.120 -0.07456845 5722.836
## TCGA-AB-2939-03A-01T-0740-13 4983.019 -0.03992778 6475.997
## TCGA-AB-2982-03B-01T-0748-13 5034.168 -0.01226374 5924.470
## TCGA-AB-2813-03A-01T-0736-13 4335.864 0.01408260 5014.1534.3 对Verhaak标签的诊断
singscoreR包提供了一系列可视化工具来探索基因标签（gene signature）。比如，这些工具可以用来观察双向标签中上调基因列表和下调基因列表的重要程度，观察基因标签中单个基因对不同样本类别的区分能力，以及还可以用来探索最终得分和离差（dispersion）之间的关系。注释可以直接被叠加在图上，Singscore支持连续型和离散型注释。连续型注释可以以向量（vector）形式或者用字符串（string）指定注释位于data frame的某一列。
针对上调基因、下调基因或者所有基因，我们开始探索分数（score）和离差（dispersion）之间的关系。plotDispersion函数可以用来生成带有注释的图。注释可以是离散型或者连续型变量，如果注释被合并在最后的得分data frame中则直接指定列名即可。singscore包中所有的画图函数都可以通过指定isInteractive为TRUE而生成可交互的动态图。
#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