---
title: "使用singscore从transcriptomic signatures预测AML中的突变"
author:
- name: Dharmesh D Bhuva
affiliation:
- Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, VIC 3052, Australia
- School of Mathematics and Statistics, University of Melbourne, Parkville, VIC 3010, Australia
email: bhuva.d@wehi.edu.au
- name: Momeneh Foroutan
affiliation:
- Department of Clinical Pathology, The University of Melbourne Centre for Cancer Research, Victorian Comprehensive Cancer Centre, Melbourne, Victoria 3000, Australia
email: momeneh.foroutan@unimelb.edu.au
- name: Yi Xie
affiliation:
- Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, VIC 3052, Australia
email: xie.y@wehi.edu.au
- name: Ruqian Lyu
affiliation:
- Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, VIC 3052, Australia
email: lyu.r@wehi.edu.au
- name: Joseph Cursons
affiliation:
- Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, VIC 3052, Australia
- Department of Medical Biology, University of Melbourne, Parkville, VIC 3010, Australia
email: cursons.j@wehi.edu.au
- name: Melissa J Davis
affiliation:
- Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, Parkville, VIC 3052, Australia
- Department of Medical Biology, University of Melbourne, Parkville, VIC 3010, Australia
- Department of Biochemistry and Molecular Biology, Faculty of Medicine, Dentistry and Health Sciences, University of Melbourne, Parkville, VIC, 3010, Australia
email: davis.m@wehi.edu.au
date: "`r format(Sys.time(), '%b %Y')`"
output:
BiocStyle::html_document:
toc_float: true
fig_caption: true
number_sections: true
nocite: |
@R-ggplot2, @R-plyr, @R-reshape2, @R-gridExtra, @R-BiocStyle, @R-knitr, @R-BiocWorkflowTools, @R-rmarkdown
bibliography: [bibliography.bib, packages.bib]
vignette: >
%\VignetteIndexEntry{Using singscore to predict mutations in AML from transcriptomic signatures (Chinese version)}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
keywords: single sample, gene set scoring, signature scoring, AML mutations, NPM1c mutation, mutation prediction, TCGA
abstract: >
针对生物转录组的RNA测序（RNA-seq）技术的进步革命性地提高了我们掌握转录调控机制的能力，这些转录调控机制可以帮助我们了解癌症等疾病的机理。最近我们发表了singscore，一种基于基因表达量排名的基因集打分（gene set scoring）方法，该方法量化了单个样本的转录谱与特定的基因集所代表的特征之间的一致性。在这里，我们展示了singscore如何应用于研究不同的急性髓细胞白血病（acute myeloid leukemia，AML）转录谱，这些转录谱与特定的几种突变和基因损伤相关。使用TCGA的基因组与转录组数据，我们展示了对恰当的标签（signature）打分可以区分出有不同突变的样本，反映出这些突变具有驱动异常的转录调控以导致白血病发生的能力。我们认为singscore方法对于研究特定癌症亚型内的异质性特别有用，并且我们展示了singscore识别那些驱动异常转录机制的突变的能力。
---
**R version**: `r R.version.string`
**Bioconductor version**: `r BiocManager::version()`
**Package version**: `r packageVersion("SingscoreAMLMutations")`
```{r setup, include=FALSE}
#set knitr chunk options
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#load packages to avoid startup messages later in the code
suppressPackageStartupMessages({library(SingscoreAMLMutations)})
library(ggplot2)
library(TCGAbiolinks)
library(SummarizedExperiment)
library(edgeR)
library(rtracklayer)
library(plyr)
library(org.Hs.eg.db)
library(GSEABase)
library(singscore)
library(reshape2)
library(gridExtra)
library(dcanr)
#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'
)
#automatically create a bib database for R packages
knitr::write_bib(c(
.packages(), 'knitr', 'rmarkdown', 'BiocStyle'
), 'packages.bib')
```
介绍 {#intro}
============
微阵列（microarray）技术的发展和RNA测序技术的快速应用为检测生物样品的转录谱（或转录组）提供了平台[@cieslik18]。转录组学分析通常侧重于样品不同分组之间基因的“差异表达（differential expression）”，然而随着公开可用RNA数据的快速增长，越来越多的人开始使用相对的方法来量化样品与特定的基因标签（gene signatures）代表的特征之间的一致性[@cieslik18]。对于用突变或融合基因进行肿瘤亚型分类与鉴定，以及寻找驱动癌症进展的基因组损伤而言，基因组的测序十分重要，而转录组学分析则可以提供携带这些突变的细胞的形态或表型信息。
癌症是一种异质性疾病，其具有多种临床和病理亚型。以乳腺癌为例，临床上主要依靠激素受体（雌激素受体：ER；孕酮受体：PR）的表达或Erb-B2受体酪氨酸激酶（*HER2*）的过度表达进行分型，这些特征的癌症分型有可直接靶向治疗的药物。一个常见的例子就是临床上使用PAM50（prediction analysis of microarray 50）标签来区分不同的乳腺癌亚型[@parker09, @cieslik18]。许多癌症的亚型分类主要依赖于识别大型患者队列中的复发性突变（recurrent mutations），通过全基因组或全外显子组测序来获得具有临床意义的亚型[@cancer13, @papaemmanuil16]。
流行的“relative approach”有single-sample gene set enrichment analysis (ssGSEA) [@barbie09]，用户可以通过GenePattern web-tool使用该方法（http://software.broadinstitute.org/cancer/software/genepattern）。另一种常见的方法是gene set variation analysis (GSVA) [@hanzelmann13]，以R/Bioconductor package的形式实现(https://bioconductor.org/packages/release/bioc/html/GSVA.html)。该package中还包括了ssGSEA，PLAGE [@tomfohr05]和z-score [@lee08]方法。ssGSEA和GSVA都使用了Kolmogorov-Smirnov like random-walk statistic来将归一化（normalised）的基因排名转换成样本得分。这一“归一化（normalised）”的步骤使得这些方法不适用于单一样本（single-sample），而且样本组成带来的variations（比如样本由不同的肿瘤亚型组成）会对样本得分造成影响。其次，这些方法产生的样本得分具有不同的数值范围（range），为阐释结果带来一些困难。为了克服这些问题，我们开发了一种单样本的基因集打分方法*singscore* [@foroutan18] ()， 它利用了给定基因集中基因的排名并相对这些基因排名的最大值与最小值进行归一化（normalization）。
The Cancer Genome Atlas（TCGA）提供了上千的病人转录组数据，通常这些数据还有配对的基因组或表观基因组数据（通常为DNA甲基化数据）。这些数据可以帮助我们找到突变所对应的功能上的变化，探索由于表观因素或者转录调控而引起的肿瘤异质性。在这里，我们展示了singscore方法 [@foroutan18] 可以利用NPM1c mutation, *KMT2A* (*MLL*) 基因融合（gene fusions）, 和 *PML-RARA* 基因融合（gene fusions）的转录组gene signatures对TCGA中的AML样本进行分类。在不需要参数拟合或估计的情况下，用singscore对基因集打分可以区分出携带有这些突变的样本。在这里，我们将介绍基因集打分不仅可以用于识别差异，还可以衡量临床结果不同的AML亚型之间的相对相似性。
生物问题简述 {#biol-problem}
=====================================
与大多数癌症一样，急性髓细胞白血病（acute myeloid leukemia，AML）是一种具有许多亚型的异质性疾病。通过对TCGA AML基因组数据的分析，有人根据特定“驱动突变（driver mutations）”的存在与否对AML进行了分类，总结并扩展了之前已定义的临床上的亚型 [@cancer13]。最近一项主要针对基因组数据的研究进一步完善了具有临床意义的AML亚型 [@papaemmanuil16]，其中囊括了许多共同发生以及相互排斥的突变。
值得注意的是，在临床AML样本中最常见的突变之一是*NPM1*基因第12外显子上的移码突变（frameshift mutation）[@papaemmanuil16]。这种突变导致核磷蛋白（nucleophosmin）的异常定位与胞质积累（cytoplasmic accumulation），因此这种突变经常被称作NPM1c突变 [@brunetti18]。如 @verhaak05 所述，NPM1c突变与同源框结构域（Hox）转录因子家族的活性失调相关，而该转录因子家族对于发育模式（developmental patterning）是必需的。最近的研究进一步证实了该突变在疾病进展中的作用，NPM1c的缺失会导致AML细胞的分化[@brunetti18]。
AML中的复发性遗传损伤还包括赖氨酸甲基转移酶2A（*KMT2A*，以前被称为*MLL*）融合基因，*KMT2A*基因（*KMT2A* -PTD）内的部分串联重复，以及早幼粒细胞白血病蛋白（*PML*）和视黄酸受体α（*RARA*）之间的基因融合（*PML-RARA*）。鉴于NPM1c突变对Hox基因家族的失调作用，具有MLL基因融合的AML样本显示出Hox家族基因的表达失调成为了一个值得探究的问题 [@hess04, @ross04]。然而，具有*MLL*-PTD的样本似乎显示出与MLL-融合样本相对不同的表型 [@ross04]。尽管有充分的证据证明NPM1c突变和其它遗传损伤在阻碍AML细胞分化中的作用，但具有*PML-RARA*融合的样本被诊断为一种叫做急性早幼粒细胞白血病（acute promyelocytic leukemia，APL)的AML亚型。这种临床AML亚型与French-American-British (FAB)分类系统中的FAB-M3相关。由于早幼粒细胞阶段的分化阻滞，该亚型的细胞显示出特有的形态 [@dethe91]。
在本分析流程中，我们证明了singscore方法 [@foroutan18]从转录组数据中区分不同肿瘤“驱动突变（driver mutations）”的能力。我们使用了已定义的NPM1c突变的gene signature[@verhaak05]，以及*PML-RARA*基因融合和MLL融合的gene signatures。后两个gene signatures来自于儿童AML样本，尽管儿童AML和成人AML之间的突变谱差别较大[@ma18]，但研究表明该signatures能够很好地区分具有类似基因损伤的成人AML样本[@ross04]。这些gene signatures已存在于MSigDB数据库（molecular signatures database）中[@liberzon15]，通过这些gene signatures，我们证明了singscore的双向打分方法可以根据不同的突变将TCGA AML样本分类，且结果具有良好的精确度（precision）和召回率（recall）。该方法的优点之一在于它是能够根据表型特征（phenotypic signatures）将样本映射到二维或者更高维度的空间中。通过比较NPM1c和*KMT2A*-/*MLL*- 融合signatures的得分，我们展示了这种分类可能是由于Hox家族失调产生的共同的下游效应而导致的。我们还将NPM1c突变signature与*PML-RARA* signature进行比较，这两个亚型的明显分离反映了它们显著不同的表型和相互排斥的突变。
数据下载与处理 {#download-and-prepare}
==================================
我们可以通过Genomic Data Commons（GDC）获得多种不同预处理的TCGA数据。转录组数据有count值和FPKM值，有用上四分位数标准化（upper quantile normalisation）前后的值。经过其它方式预处理后的数据可以在 [www.cbioportal.org](http://www.cbioportal.org/) 和 [firebrowse.org](http://firebrowse.org) 找到。GDC的数据使用了STAR的“two-pass”模式进行序列比对，使用了HTSeq进行定量。用户在使用 [GDC data transfer tool](https://gdc.cancer.gov/access-data/gdc-data-transfer-tool) 下载GDC中的数据时，可以通过GDC portal选择自己感兴趣的特定文件进行下载。下载之后需要读取并整合这些文件之后才可以进行下游分析。以上这些步骤（包括下载）可以通过R包`TCGAbiolinks` [@colaprico15, @R-TCGAbiolinks]完成。这个包支持使用GDC API和GDC client进行数据下载，我们会使用该包完成下载、注释并将数据整合成SummarizedExperiment R对象。
在开始分析之前需要完成一下几步数据处理：
1. 创建查询列表以下载特定的文件；
2. 执行查询以下载数据；
3. 将下载好的数据读入R中；
4. 过滤掉表达量低的基因；
5. 校正样本组成带来的偏差并对基因长度进行标准化 [@foroutan18]
查询GDC数据库 {#gdc-query}
-------------------------
任何数据分析的第一步是先确定好数据的版本和下载数据的方式。`getGDCInfo()`函数可以返回GDC数据库中数据的版本号和发布时间。
```{r gdc_query}
library(SingscoreAMLMutations)
library(TCGAbiolinks)
#get GDC version information
gdc_info = getGDCInfo()
gdc_info
```
接下来我们需要创建一个查询命令，从GDC找到并下载特定的文件。这一步和使用GDC portal创建*MANIFEST*文件很类似。查询命令的第一个参数（project）指定了项目的名称（使用`getGDCprojects()`函数或者进入可以获得所有项目名称），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值也是合理的。
```{r gdc_results}
#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)
colnames(rnaseq_res)
```
下载TCGA AML RNA-seq数据（read counts） {#download-data}
--------------------------------------------
`GDCdownload`函数可以执行查询并使用GDC API下载数据。如果需要下载的文件很大（如RNA-seq的read数据或者甲基化数据），应该把`GDCdownload`函数里的下载方法切换成“client”。这里我们暂时将数据存储在“GDCdata”文件夹中，用户应该指定自己的存储路径以妥善保存好数据。`GDCdownload`函数通过这种参数设定的方法使我们能够将不同类型的数据存储在同一个文件夹内进行管理与使用。这里，下载后我们可以看到count-level数据被存储在*TEMPDIR/GDCdata/TCGA-LAML/harmonized/Transcriptome_Profiling/Gene_Expression_Quantification/*路径下。
```{r gdc_download, results='hide'}
datapath = file.path(tempdir(), 'GDCdata')
GDCdownload(query_rna, directory = datapath) #(size: 39MB)
```
将count-level数据读入R {#read-data}
-------------------------------
`GDCprepare`函数可以将下载好的数据读入R并处理成`RangedSummarizedExperiment`对象（来自`SummarizedExperiment`包），该对象可以同时储存read count、基因注释和临床信息。在调用该函数的时候临床信息自动下载并整合入`RangedSummarizedExperiment`对象中。RangedSummarizedExperiment对象和ExpressionSet很类似，但它还可以使用基因坐标进行索引以及存储具有相同结构的多个数据矩阵。关于数据特征（feature）的注释被存储在一个RDA/RDATA文件中。
```{r gdc_prepare, results='hide'}
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）的。
```{r show_se}
aml_se
```
过滤掉counts数低的基因 {#filter-data}
--------------------------------
`edgeR`包提供了过滤所需的数据标准化与转换方法。这些方法要求数据存储在DGEList对象中，因此我们需要将SummarizedExperiment对象转换成DGEList对象。
```{r remove_dups}
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）的基因。在特定情况下可能其它过滤低表达基因的方法更加适用。 @chen16 和 @law16 根据实验设计来过滤基因,他们先在各个组内确定表达量过低的基因，然后在所有样本中过滤掉这些基因。这种过滤策略适合那些样本较少的数据集，AML数据有足够多的样本因此我们直接对所有样本进行过滤。在图\@ref(fig:plot-hist-filtering)中我们可以看到过滤后logCPMs值的分布更加接近理想的log-normal分布。
```{r plot-hist-filtering, fig.wide=TRUE, fig.cap="过滤前后AML数据logCPM值的直方图. 过滤之后数据中的零值减少了。在大部分样本里CPM值小于1（logCPM < 0）的基因被过滤掉，使得最后得到了近似log-normal分布的logCPM值。"}
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)
par(op)
```
```{r remove_low_counts}
#subset the data
aml_dge = aml_dge[keep, , keep.lib.sizes = FALSE]
aml_se = aml_se[keep, ]
```
转换成FPKM值并归一化（normalisation） {#calc-fpkm}
-----------------------------------------------
Singscore要求同一个样本内的基因表达量是可以互相比较的，因此我们需要对基因长度进行归一化[@oshlack09]。数据可以被转换成transcripts per million (TPM)或者reads/fragments per kilobase per million (RPKM/FPKM)，它们都对基因长度进行了归一化。因此针对这两种转换方法，只要library size足够大，singscore得到的结果应当是类似的。edgeR包中的`calcNormFactors`函数提供了三种对基因长度归一化的方法，在不指定的情况下TMM normalisation是默认的方法。 @chen16 和 @law16 讨论了在下游分析如差异表达分析之前进行数据归一化的意义。通常归一化的目的是使得样本之间的比较是有意义的。Singscore的分析是针对单个样本内的因此无需对整体样本进行归一化。同理，该方法也不受其它转换方式（如log transformation）的影响。这里我们仅出于可视化目的使用TMM normalisation。
将数据转换成TPM或者RPKM/FPKM值需要先计算所有基因的长度。基因长度的计算依赖于序列比对过程和定量参数。TCGA转录组数据使用了STAR进行序列比对并用HTSeq进行定量（流程细节见）。HTSeq统计了比对到每个基因外显子区域的read数量，因此有效基因长度应当是基因的所有外显子区域长度之和。GENCODE v22被用于定量过程中的注释，因此在计算基因长度时仍然应该用这个注释文件。
```{r download_gencode, results='hide'}
#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文件的函数。
```{r compute_gene_lengths}
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。
```{r add_biotype}
#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
```
SummarizedExperiment对象可以存储基因的注释信息，因此我们需要把Ensembl IDs、gene length和gene type加进去。类似的，我们也要把这些注释信息加到DGEList对象中。
```{r add_length_annotation}
#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`中。
```{r compute_fpkm}
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
```
对样本注释突变信息 {#annotate-mutations}
-----------------------------------
读者们需要注意的是，本分析流程中我们用了从原始TCGA AML文献中[@cancer13]提取的突变列表(Supplemental Table 01 at )而不是标准TCGA流程中的突变信息（可以在[National Cancer Institute Genomic Data Commons](https://gdc.cancer.gov/)获得），两者存在一定差异。针对我们所关心的遗传病变（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
```{r preproc_mutations, eval=FALSE, include=FALSE}
#preprocessing - read tsv (Joe's version) and convert save as RDS
mut_info = read.csv('PatientMutations.tsv', sep = '\t', colClasses = c('character', rep('logical', 4)))
rownames(mut_info) = mut_info$Barcode
mut_info = mut_info[, -1]
patients = mut_info[substring(colnames(aml_se), 1, 12), ]
rownames(patients) = colnames(aml_se)
saveRDS(patients, file = 'AMLPatientMutationsTCGA.rds')
```
```{r annotate_mutations}
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]
```
匹配Entrez IDs {#map-gene-ids}
-----------------------------
Ensembl注释（Ensembl IDs）在基因组上有更高的覆盖度，因此适合用于识别突变（variant calling）等类似的分析过程[@zhao15]。而RefSeq注释（Entrez IDs）对RNA-seq分析更适用，通常RNA-seq分析需要稳定的注释信息以便于未来的比较[@wu13]。因此我们选择将Entrez IDs和Ensembl IDs匹配后去掉那些没有匹配到Entrez ID的基因。
匹配可以通过biomaRt的bioconductor R包完成，它可以提供最新的注释信息。匹配还可以通过一年更新两次的R包`org.Hs.eg.db`[@R-org.Hs.eg.db]完成，它提供了更加稳定的注释信息，提高了分析的可重现性。R包`AnnotationDbi`[@R-AnnotationDbi]中的`mapIds`函数可以实现匹配功能。
```{r map_ensembl_entrez}
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一对一匹配的数据。
```{r discard_multimapped_genes}
#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)
#subset the data
aml_se = aml_se[keep, ]
```
转录组标签（signatures）预测突变状态 {#transcriptional-mut-sig}
=====================================================
在这里，来自 @verhaak05 的标签（signature）被用来预测NPM1c突变的状态。我们需要量化标签（signature）内的基因和它们在样本中表达水平的一致性（concordance）。如标签（signature）内的上调基因（up-regulated genes）在样本中高表达，下调基因（down-regulated genes）在样本中低表达就会产生较高的分数。接下来，这个分数就可以被用来预测样本的突变状态。
首先，我们先将标签（signature）从MSigDB中下载下来并使用R包`GSEABase`[@R-GSEABase]读入`GeneSet`对象中。接着我们用R/Bioconductor包`singscore`来给Verhaak标签（signature）打分，`singscore`包内的一些函数可用于后续的可视化与诊断。最后，我们对最终分数使用了逻辑回归模型（logistic regression model）来预测突变状态。
下载标签并载入 {#prepare-signature}
----------------------------------
@verhaak05 的标签包含有上调和下调的基因列表。许多标签都是这种形式以增强其对样本的区分力。MSigDB将这样的标签（signatures）分成两部分并分别用后缀“_UP”和“_DN”标记。这里，我们直接将标签名“VERHAAK_AML_WITH_NPM1_MUTATED”和"_UP"或"_DN"连接以形成下载链接。
```{r signature_names}
#create signature names
verhaak_names = paste('VERHAAK_AML_WITH_NPM1_MUTATED', c('UP', 'DN'), sep = '_')
verhaak_names
```
利用刚才创造的链接下载后，"_UP"和"_DN"分别可以得到一个XML文件，在这里我们指定好输出XML文件的文件名（参数`verhaak_files`）。`mapply`函数可以用来循环下载类似的“链接-输出”参数对。
```{r download_signatures, results='hide'}
#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`对象。
```{r read_signatures}
library(GSEABase)
verhaak_sigs = getBroadSets(verhaak_files, membersId = 'MEMBERS_EZID')
verhaak_sigs
```
为了使得分析过程中的索引更加方便，我们把`SummarisedExperiment`对象的行名改成Entrez IDs。
```{r rows_to_entrez}
rownames(aml_se) = rowData(aml_se)$entrezgene
```
使用标签给样本打分 {#score-samples}
--------------------------------------------------
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。这样可以减少零值对排名带来的影响，不过在数据质量控制时过滤掉低表达基因仍然是必不可少的步骤（见\@ref(filter-data)部分）。
```{r compute_ranks_aml}
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），包括这里我们用的 @verhaak05 的标签（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。
```{r compute_scores_verhaak, warning=TRUE}
#apply the scoring function
verhaak_scores = simpleScore(aml_ranked,
upSet = verhaak_sigs[[1]],
downSet = verhaak_sigs[[2]])
```
应该注意的是，singscores由两部分组成，即代表富集程度的分数（score）和对排名的离散程度的估计。在理想情况下，所有预期上调的基因都具有高表达，因此将基因从低表达到高表达排序后得到的值越大，这样的值应当位于分布的右端。Singscore旨在量化这种排名的分布，因此它计算了标签（signature）中基因排名的平均值和离差（dispersion）。平均值的计算和其它单样本打分方法类似，但是我们认为用平均值和离差（dispersion）两个数据来观察标签（signature）内基因排名的分布更为合理。默认且推荐的离差（dispersion）统计方法是具有非参数统计特性的绝对中位差（median absolute deviation，MAD）。其它方法还有四分差（inter-quartile range，IQR），用户可以将`IQR`函数作为参数传递给`dispersionFun`。
```{r}
head(verhaak_scores)
```
对Verhaak标签的诊断{#diagnostics-verhaak}
------------------------------------
`singscore`R包提供了一系列可视化工具来探索基因标签（gene signature）。比如，这些工具可以用来观察双向标签中上调基因列表和下调基因列表的重要程度，观察基因标签中单个基因对不同样本类别的区分能力，以及还可以用来探索最终得分和离差（dispersion）之间的关系。注释可以直接被叠加在图上，Singscore支持连续型和离散型注释。连续型注释可以以向量（vector）形式或者用字符串（string）指定注释位于data frame的某一列。
针对上调基因、下调基因或者所有基因，我们开始探索分数（score）和离差（dispersion）之间的关系。`plotDispersion`函数可以用来生成带有注释的图。注释可以是离散型或者连续型变量，如果注释被合并在最后的得分data frame中则直接指定列名即可。`singscore`包中所有的画图函数都可以通过指定`isInteractive`为`TRUE`而生成可交互的动态图。
```{r plot-dispersion, fig.wide=TRUE, fig.height=4, fig.cap="使用NPM1c signature上调和/或下调基因列表得到的分数可以将携带有NPM1c mutations或*MLL* fusions/PTDs的样本区分出来. 以Scores和基因集内基因排名的median absolute deviation (MAD)为轴作图，图上可以看出scores值接近0的样本点拥有更高的MADs值。"}
#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
```
图\@ref(fig:plot-dispersion)显示下调的基因更能区分出具有NPM1c突变的样本。而且我们可以看到该标签（signature）还能够将*MLL*（*KMT2A*）fusions和PTDs从其它样本中区分开来。在下调基因标签的打分中，具有NPM1c突变的样本分数最高，其次是*MLL*fusions和PTDs样本。在上调基因标签的打分中，分数也有这样的趋势，而且尽管分数的范围变大了，它区分不同类型样本的能力只发生了略微减弱。事实上，上调基因能够更好地区分出那些没有我们感兴趣的突变的样本（标记为“Other”），这些样本大多数具有负的得分。上调基因标签的打分为负分意味着这些基因在样本中的表达量低于中位数，即可能在样本中是下调的基因。观察离差（MAD）的趋势我们发现，分数越接近0值则离差越大。有三种情况会导致得分接近0值：1）标签中基因的表达量在中位数附近；2）标签中基因的表达量平均分布在表达谱两端；3）标签中基因的表达量均匀分布在整个表达谱中（可能性最大）。后两种情况会导致离差（dispersion）值很大。为了进一步弄清这个问题，我们选择了三个样本并对上调和下调基因标签中基因的排名分布进行作图。三个样本分别为总体得分最高的样本、总体得分最低的样本和总体离差（dispersion）值最高的样本。
```{r plot-rank-density, fig.wide=TRUE, fig.height=4, fig.cap="基因排名的分布显示了应该上调的基因在高分样本中上调，在低分样本中下调；应该下调的基因在高分样本中下调，在低分样本中上调. 从左到右分别是总体得分最高的样本、总体得分最低的样本和总体离差（dispersion）值最高的样本。barcode图显示了标签内每个基因的排名及其密度函数。上调基因的颜色为绿色而下调基因为紫色。在最高分和最低分样本中，排名两端呈现正态分布（Gaussian distributions）。总体离差（dispersion）最高的样本得分接近0，因为下调基因的排名呈现均匀分布而上调基因的排名呈现双峰分布。"}
library(gridExtra)
#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)
```
图\@ref(fig:plot-rank-density)显示上调和下调基因列表分别对总体得分最高的样本作出了贡献，这些基因分别位于排名分布的两端。类似地，在得分最低的样本中，上调基因位于排名分布的左端而下调基因位于排名分布的右端从而产生了很低的得分。正如前一张图中看到的那样，上调基因列表提高了NPM1c突变样本和other样本之间的区分度，可能相比下调基因列表能更好地指示野生型NPM1c样本。最后，离差（dispersion）值最高的样本中，下调基因的排名呈均匀分布而上调基因的排名呈双峰分布，双峰分别位于谱的两端。
综上，这些图显示了上调基因和下调基因列表在区分NPM1c突变、*MLL*fusion/PTD和野生型样本时都很重要。这些图能帮助用户在应用之前先确定标签（signature）的重要性，有时候标签（signature）是在特定情境下产生的，本身带有偏差因而只能在特定情况下使用。这些图还能帮助验证这些标签（signature）的适用情境。
用Verhaak标签预测突变 {#predict-verhaak}
------------------------------------------------
突变状态可以通过logit函数的逻辑回归模型（logistic regression model）进行预测。这个模型相对于将标签（signature）中每一个基因作为变量的方法而言更加简便， @verhaak05 的标签（signature）中的429个基因在逻辑回归中可以产生429个预测变量（predictor）。由于样本数少于预测变量（predictor）数，一些特征选取（feature selection）方法可能导致信息的丢失。此外，基因水平的模型在选择更多基因标签（gene signatures）上存在限制。而singscore的模型在某种程度上会继承它非参数的特性。比如该模型对于任何保留了基因排序的数据转换都是稳健的。该模型主要的不足在于丢失了一定的准确性，因为它将429个基因的信息限制到两个分类中。
不管怎样，我们的目的不是讨论哪个模型应该被用于预测突变状态，而是展示singscore和转录组标签（signature）区分样本的能力。因此我们在这里用训练模型的数据评价该模型的表现。我们先把分数和突变注释合并到一起。
```{r create_plotdf_boxplot}
#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
```
在训练模型之前，我们可以先可视化一下，看 @verhaak05 的标签（signature）得到的分数区分突变型和野生型样本的能力。图\@ref(fig:plot-boxplot-mutscores)展示了分数可以区分开NPM1c突变型样本和野生型样本，*MLL* (*KMT2A*) fusions和*PML-RARA* fusions的情况也是如此。
```{r plot-boxplot-mutscores, fig.height=5, fig.cap="NPM1c标签的得分可以区分突变样本与野生型样本. NPM1c标签得分的Boxplot（按照不同突变类型分类）。"}
#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
```
为了量化以上的结果，我们将每个基因损伤（genomic lesion）作为因变量（response variable），样本的分数作为预测变量（predictor）进行逻辑回归并列出回归系数估计，标准差，z值和p值。
```{r predict_mutation_glm}
#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')
```
NPM1c突变和分数显著相关，*MLL* gene fusion也和分数显著相关，这显示了它们在Hox基因失调中共有的作用。有趣的是，*PML-RARA* fusion的系数是负值，这很可能反映了急性早幼粒细胞白血病(acute promyelocytic leukemia, APL)相对于其它AML亚型而言独特的细胞形态或表型特征。
评估预测模型的表现{#prediction-performance}
-------------------------------------
以上统计数字只提供了模型本身的信息而非模型的表现好坏。精确率（precision）和召回率（recall）可以衡量模型表现的好坏。`dcanr`包[@R-dcanr]提供了函数进行这类计算。
```{r eval_prediction_performance}
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
```
这里我们只能计算对NPM1c突变预测的精确率（precision），具有其它突变类型的样本在此情景下的预测结果是野生型（wild-types），因此无法计算它们的精确率。同理，具有其他突变类型的样本的召回率（recall）是0。NPM1c突变的模型的精确率（precision）、召回率（recall）、F1-score和预期的一样高。就像\@ref(score-samples)部分提到的，singscores是由两部分组成的分数，因此在考虑排名的离差（dispersion）后，模型表现应该更加优秀（图\@ref(fig:plot-dispersion)）。
以下的模型显示，使用了singscores中的score及离差（dispersion）后，模型的预测能力显著提升。
```{r predict_with_mads}
#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)
)
```
多个基因标签（signature）的landscapes {#signature-landscape}
=============================================
通常我们会对两个独立的或者相互依赖的表型之间的关系感兴趣，比如EMT过程中epithelial和mesenchymal表型。大多数标签（signature）的作用是为了便于量化一些难以衡量的分子表型。因此，我们可以用对应的标签（signature）来探索两个表型之间的关系。 @foroutan17 首次介绍了使用molecular signature landscapes来探索标签（signature）之间的关系，这些标签（signature）均和EMT及TGF$\beta$导致的EMT相关。之后， @cursons18 计算了EMT表型标签（signature）的singscores得分并用signature landscape证明了将miR-200c转染到mesenchymal细胞系之后发生了向epithelial表型的转变。 @foroutan18 使用signature landscapes在epithelial-mesenchymal尺度上划分乳腺癌亚型并把它包含在了`singscore`包中。这里，我们展示了signature landscapes的进一步应用，即使用代表不同突变的转录组标签（signature）来划分AML样本。
Ross MLL fusion signature vs. Verhaak signature landscape {#sig-landscape-mll}
---------------------------------------------------------
我们现在使用 @ross04 的MLL-fusion标签（signatures）来给TCGA AML样本打分。不像NPM1c标签，这个标签中的基因可以用来区分具有MLL-fusion基因的样本。我们按照\@ref(prepare-signature)中的流程下载并处理标签。
```{r process_rossmll_sigs}
#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')
download.file(rossmll_link, rossmll_file, method = 'libcurl')
rossmll_sig = getBroadSets(rossmll_file, membersId = 'MEMBERS_EZID')
rossmll_sig
```
该标签包含有上调和下调基因以区分突变型和野生型。我们用\@ref(score-samples)部分提到过的singscore的第三种模式进行打分，该模式不需要提前知道或确定基因是上调还是下调的。`simpleScore`函数里的`knownDirection`参数应当设定成`FALSE`。\@ref(score-samples)部分已经计算过的基因排名可以在这里再次使用，以计算新的标签的分数。
```{r score_rossmll}
rossmll_scores = simpleScore(aml_ranked, rossmll_sig[[1]], knownDirection = FALSE)
```
我们可以用`plotScoreLandscape`函数画hexbin图来可视化分数的分布。使用 @verhaak05 标签及 @ross04 *MLL*-fusions标签得到的分数被传递给了该函数。两个分数应该基于相同的样本计算后得到。这两个分数所对应的名字应当传递给`scorenames`参数。`textSize`参数可以用来指定图中文字相对于图片的大小，这在应对不同poster、publications和presentations的要求时非常有用。
```{r plot-signature-landscape, fig.height=7, fig.width=7, fig.cap="MLL fusion和NPM1c标签的landscape. AML中两个标签存在正相关。"}
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
```
图\@ref(fig:plot-signature-landscape)显示，尽管两个标签之间只有16个共同基因（标签大小分别为78和429个基因），两个标签的分数之间仍具有很强的正相关（Spearman's $\rho$ = `r round(cor(verhaak_scores$TotalScore, rossmll_scores$TotalScore, method = 'spearman'), digits = 3)`）。在这样的分析中，我们可能想像 @cursons18 文章里那样把新的数据点投影在图上。或者我们可能希望将已有的数据点投影在图上以观察样本的分类情况。这里我们将*MLL* fusions、*MLL* PTDs、*PML-RARA* fusions和NPM1c mutations四类样本的数据点进行投影。首先我们需要先建立注释。
```{r new_annotations}
#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`参数用于选择想要投影的数据点子集，这里我们选择了感兴趣的四种样本。
```{r plot-project-mll, fig.height=7.5, fig.width=7, fig.cap="Landscape展示了MLL fusions和PTDs之间有很大不同. 不同的突变在landscape上占据了不同的位置，提示了它们具有不同的分子特征。"}
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')
```
图\@ref(fig:plot-project-mll)显示*MLL* fusions和*MLL* PTDs两组样本在两个维度上被分开。在*MLL* fusion标签的方向上，*MLL* PTDs相比*MLL* fusions的得分更低。这些样本在两个维度的平面上没有形成cluster，这和 @ross04 的结论是一致的。在图上进行样本注释的投射有助于帮助我们解释landscape不同部分的意义。
Ross PML-RARA fusion signature vs. Verhaak signature landscape {#sig-landscape-pmlrara}
--------------------------------------------------------------
之前的分析显示*PML-RARA*样本和其它的样本有明显不同。这里我们用来自 @ross04 的*PML-RARA*标签重复了刚才的分析以验证这种不同。该标签在MSigDB中的名字为“ROSS_AML_WITH_PML_RARA_FUSION”。我们下载该标签并对所有样本进行打分，用该分数和 @verhaak05 标签得到的分数画landscape图。最后我们将样本投影到图上（见\@ref(sig-landscape-mll)部分）。该标签的和*MLL* fusions标签来自同一个研究结果因此我们仍然使用相同的参数进行打分。
```{r plot-project-pmlrara, fig.height=7.5, fig.width=7, echo=FALSE, fig.cap="PML-RARA标签和NPM1c标签之间没有关联. L型的landscape提示了这两个突变背后的分子机制是互斥的。"}
newsig_name = 'ROSS_AML_WITH_PML_RARA_FUSION'
#generate URL
newsig_link = paste0(
'http://software.broadinstitute.org/gsea/msigdb/download_geneset.jsp?geneSetName=',
newsig_name,
'&fileType=xml'
)
#download files
newsig_file = paste0(newsig_name, '.xml')
download.file(newsig_link, newsig_file, method = 'libcurl')
newsig_sig = getBroadSets(newsig_file, membersId = 'MEMBERS_EZID')
#score
newsig_scores = simpleScore(aml_ranked, newsig_sig[[1]], knownDirection = FALSE)
#plot and project
p1 = plotScoreLandscape(
verhaak_scores,
newsig_scores,
scorenames = c('VERHAAK_AML_WITH_NPM1_MUTATED', newsig_name),
textSize = relSize
)
#project NPM1 mutations onto the landscape
mutated_gene = rep('Other', ncol(aml_se))
mutated_gene[aml_se$NPM1c.Mut] = 'NPM1c Mut'
mutated_gene[aml_se$KMT2A.Fusion] = 'MLL Fusion'
mutated_gene[aml_se$KMT2A.PTD] = 'MLL PTD'
mutated_gene[aml_se$PML.RARA] = 'PML-RARA'
select_aml = !mutated_gene %in% 'Other'
p2 = projectScoreLandscape(p1,
verhaak_scores[, 1:2],
newsig_scores,
subSamples = select_aml,
annot = as.factor(mutated_gene[select_aml]))
p2 + theme(legend.box = 'vertical')
```
图\@ref(fig:plot-project-pmlrara)展现了和*MLL* fusion标签完全不同的landscape。*PML-RARA*标签使得*PML-RARA*和NPM1c样本完全分开，*PML-RARA*样本是唯一在*PML-RARA*标签中得到高分的样本，而且两个标签之间没有显著相关。在样本背景介绍部分，*PML-RARA* fusion是AML中一种已知的亚型——急性早幼粒细胞白血病（acute prolmyelocytic leukemia，APL）的诊断标志，这种亚型独特的细胞表型反映了早幼粒细胞分化的阻滞[@dethe91]。这种特性反映在了L型的landscape中，以及这些基因组损伤驱动白血病不同的机制。
总结
=======
singscore包通过R/Bioconductor环境给使用者提供了一个用户友好的界面来进行基因集打分（gene set scoring）。 TCGABiolinks包使用户相对容易地访问大型的临床相关数据集，如TCGA及其相关注释。 singscore中包含的诊断和绘图功能允许用户调查感兴趣的基因集（gene set）以确定它们区分样本之间差异的能力。然后不同的基因集（gene set）可以组合在一起以探索不同的细胞表型在大型研究中的关系。以上分析展示了当恰当的基因集被使用时，singscore计算得到的结果可以用于样本分类，结果可以进一步被用到缺少基因组数据的大型转录组数据中。
使用到的R包 {.unnumbered}
=============
此工作流程用到了许多Bioconductor项目`r BiocManager::version()`版本中的包，适用于`r version$version.string`或更高版本。本工作流程中用到的所有包已列在下表中：
```{r session_info}
sessionInfo()
```
致谢 {.unnumbered}
===============
工作流程中的结果全部或部分基于TCGA Research Network（）产生的数据。感谢来自The Walter and Eliza Hall Institute的Christoffer Flensburg在TCGA AML mutation data上的帮助。
参考文献 {.unnumbered}
==========