R数据可视化12:曼哈顿图

2023-02-28 01:35:02Python014

R数据可视化12:曼哈顿图,第1张

曼哈顿图是一种散点图，通常用于显示具有大量数据点，许多非零振幅和更高振幅值分布的数据。该图通常用于全基因组关联研究（GWAS）以显示重要的SNP（来源wiki）。

在图中每个点代表一个SNP，纵轴为每个SNP计算出来的Pvalue取-log10，横轴为SNP所在的染色体。基因位点的Pvalue越小即-log10(Pvalue)越大，其与表型性状或疾病等关联程度越强。而且通常来说受到连锁不平衡的影响，强关联位点周围的SNP也会显示出相对较高的信号强度，并依次向两边递减，所以会出现上图中红色部分的现象。一般，在GWAS的研究中，Pvalue的阈值在10^-6 或者10^-8以下。（现在可能要求更高了？好久没看过文章）

用于做曼哈顿图最常用的一个R包叫做qqman——an R package for creating Q-Q and manhattan plots。本文我们直接使用该包中的例子进行讲解（毕竟我也没有可以绘图的GWAS数据哈哈哈）。没有安装的可以先输入 install.packages("qqman") 安装该包。当然qqman包由于是为曼哈顿图服务所以其实有很多限制，如果想要完全DIY我们可以使用ggplot。本文将会介绍使用这两个R包进行绘图。

下述内容来源于 Manhattan plot in R: a review ,我只是一个搬运工。

1）需要什么格式的数据

qqman提供的数据例子很直接就叫做"gwasResults"，数据格式如下：

第一列为SNP的名字，第二列CHR为所在染色体，第三列BP为染色体上所在位置。要注意如果你的CHR中存在X，Y这样的，需要给他们转化为数字如赋予23，24等，其中第一列SNP的名字是可选择的，后三列是必须提供的。

2）如何作图

利用 manhattan 函数即可作出相应的曼哈顿图。

如果你想要标记其中一系列你感兴趣的SNP，怎么办呢？给出你感兴趣的snpsOfInterest列表即可。

如果你想知道每条染色体上pvalue最小的SNP，你可以通过下述方式：

如果不喜欢黑色和灰色的搭配，也可以自行改变颜色。

那么使用ggplot要如何作图呢？

don是用于作图的主要数据表，而axisdf是用于处理x轴，因为我们想要他们按照染色体的位置排布。上述数据处理完成后，我们就可以使用ggplot作图：

我们这里展示一下有无 scale_x_continuous( label = axisdf$CHR, breaks= axisdf$center ) 的差异，可以看到x轴的变化

那么如果想要把某些SNP标记出来呢？那么我们在前期处理数据的时候需要将这些数据标记出来，这个过程和之前火山图标记显著的基因很类似:

然后画图的时候geom_point在颜色上进行区分，并使用geom_label_repel标注出来即可：

那么，关于曼哈顿图的分享就到这里啦。

往期R数据可视化分享

R数据可视化11: 相关性图

R数据可视化10: 蜜蜂图 Beeswarm

R数据可视化9: 棒棒糖图 Lollipop Chart

R数据可视化8: 金字塔图和偏差图

R数据可视化7: 气泡图 Bubble Plot

R数据可视化6: 面积图 Area Chart

R数据可视化5: 热图 Heatmap

R数据可视化4: PCA和PCoA图

R数据可视化3: 直方/条形图

R数据可视化2: 箱形图 Boxplot

R数据可视化1: 火山图

我们经常用随机森林等机器学习又或者是其他数据挖掘的方法寻找某些疾病的biomarker或者候选基因。但是来自临床的数据包括了生存事件等信息，数据的内容有所不同，所以需要一些和之前不太一样分析方法，其中常见的就是通过制作生存曲线图获取结论。

生存曲线可以帮助我们回答许多问题：参与者生存5年的概率是多少？两组之间的生存率是否存在差异（例如，在临床试验中分配给新药还是标准药的两组之间）？某些行为或临床特征如何影响参与者的生存机会？

通常，在这类分析中，我们会关注特定事件（如死亡或疾病复发）的事件，并比较两组或更多组患者发生这些特定事件的事件。

可以看到上图显示了经常玩棋类游戏的老年人和很少玩这类游戏的老年人之间的痴呆风险Kaplan-Meier曲线。纵轴为非痴呆老人的比例，横轴为跟踪的年数，从图中可以看到经常玩棋类游戏的老年人患痴呆的风险较低。

在制作生长曲线之前，我们需要首先了解几个相关的术语

参考： R语言-Survival analysis（生存分析）

Event（事件）： 指在随访过程中发生的某个结果，如癌症研究中，可能为复发（Relapse）、恶化（Progression）、死亡（Death）

Survival time（生存时间）： 指某个事件开始到终止的时间，在癌症研究中经常用到的几个指标：

Overall survival（OS）：

指从开始到任意原因死亡的时间，一般常见的5年生存率、10年生存率都是基于OS计算的

Progression-free Survival（PFS，无进展生存期）：

指从开始到肿瘤发生任意进展或者死亡的时间，可用于评估治疗方法的临床效益

Time to Progress（TTP，疾病进展时间）：

从开始到肿瘤发生任意进展或者进展前死亡的时间，与PFS相比仅包括肿瘤的恶化，而不包括死亡。

Disease-free Survival（DFS，无病生存期）：

指从开始到肿瘤复发或任何原因死亡的时间，常用于根治性手术治疗或放疗后的辅助治疗的评估

Event Free Survival（EFS，无事件生存期）：

指从开始到发生包括肿瘤进展、死亡、治疗方案的改变等各种事件的时间

Censoring（删失）： 一般指不是由于死亡造成的数据丢失，可能是由于失访、非正常原因推出、时间终止而事件未发生等，一般在展示时用“+”表示

生存分析的方法一般可以分为三类：

1、参数法：已知生存时间的分布模型，根据数据估计模型参数，最后以分布模型计算生存率

2、半参数法：不需要知道生存时间的分布，但是仍通过模型来评估影响生存率的因素，常见方法如 Cox回归模型

3、非参数法：不需要知道生存时间的分布，根据样本统计量估计生存率，常见方法如 Kaplan-Meier方法、寿命法

具体地，我们通过同样一个例子介绍常用的Kaplan-Meier方法和寿命法的异同。

例子：一项探究死亡时间的前瞻性队列研究，研究涉及20位65岁以上的参与者，招募时间为5年，整个研究进行长达24年的随访直至死亡、研究结束或退出研究（失访）。因此，如果参与者是在研究开始后加入的，他们的最长随访时间应该少于24年。具体数据如下，其中有6位参与者死亡，3位接受了完整的随访（24年），其余11位由于在研究开始后加入或失访而少于24年随访:

寿命法

寿命法经常用于保险行业中估计预期寿命并设置保费。不过，我们只关注生物领域的使用，我们称为随访生命表，该表记录了参与者在队列研究或临床试验中在预定的随访期内的经历，直到目标事件发生或研究结束为止。

要构建生命表，我们要将随访时间分割成间距相等的几组，上述例子中我们随访的最长时间为24年，所以我们考虑5年一个间隔（0-4，5-9，10-14，15-19和20-24年）。然后统计每个时间间隔开始时活着的参与者人数，和该期间死亡人数和每个时间间隔中删失的人数。

然后，我们来定义几个参数：

N t =在时间间隔t内没有发生目标事件的但处于风险中的人数（如本研究中目标事件为死亡，而参与者都处于可能死亡的风险之中）

D t =在时间间隔t内死亡的人数

C t =在时间间隔t内删失的人数

N t * =在时间间隔t内有风险的参与者的平均数（计算公式为：N t * =N t -C t /2）

q t =时间间隔t内死亡比例，q t =D t /N t *

p t =时间间隔t内生存比例，p t =1-q t

S t ，累计生存概率，S 0 =1,S t+1 =p t+1 *S t

因此，对于第一个间隔0-4年和第二个5-9年的间隔，可以计算出如下数据：

所以完整的随访寿命表为：

Kaplan-Meier

Edward Kaplan和Paul Meier于1958年在《American Statistical Association》共同发表了Kaplan-Meier非参数估计方法，让我们能够估计生存函数。

从寿命表的方法可以看出生存概率会根据不同的间隔改变，尤其是对于小样本而言这种改变可能会很剧烈。Kaplan-Meier通过每次时间发生时重新估计生存概率来解决该问题。

Kaplan-Meier是基于这样的假设进行的：删失与事件发生的可能性无关，且在研究早期和后期被招募的参与者生存率是可比的。这些前提很重要，比如在不同组比较时要保证删失的可能性一致。

Kaplan-Meier与寿命法的计算方式类似，主要区别是时间间隔，寿命法中我们选择的时间间隔相等，而在Kaplan-Meier的方法中我们使用观察到的事件时间和删失时间。

上述的内容原版，以及关于进一步的检验和Cox模型的内容可以阅读Boston大学的教材 Boston Univeristy Suvival Analysis 。在这里暂时就不再解释啦。

今天我们要用到以下几个R包：survival，survminer和dplyr

使用KM方法，通过 ggsurvplot 作图，该函数作图需要两部分数据，具体见下：

1）需要什么格式的数据

我们使用的数据集为ovarian，来自survival包。该数据集来源于文章《Different Chemotherapeutic Sensitivities and Host Factors Affecting Prognosis in Advanced Ovarian Carcinoma vs. Minimal Residual Disease》，主要研究化疗敏感性和宿主因素对晚期卵巢癌和微小残留病变的预后影响，具体含有以下几个指标：

futime: survival or censoring time 生存时间

fustat: censoring status 确定参与者生存时间是否发生缺失

age:in years

resid.ds: residual disease present (1=no,2=yes) 评估肿瘤的消退情况

rx: treatment group 接受两种治疗方案中的一种

ecog.ps:ECOG performance status (1 is better, see reference)依据ECOG评估的患者表现

为了更直观的获取信息，我们根据说明修改一下部分指标的标记方式：

然后我们来看一下年龄的分布 hist(ovarian$age)

然后我们进行生存曲线的分析，使用futime和fustat两列，首先根据是否发生删失对数据进行处理。

可以看到发生删失的都带上了加号。

然后拟合Kaplan-Meier曲线：

2)如何作图