黄致远,江文豪,杨福贵,朱 磊,董 林,汪令伟,冯明昊,魏 亮,王光学*,李钦传,*
1.同济大学附属东方医院心胸外科(上海 200120)
2.同济大学附属东方医院转化医学研究中心(上海 200120)
非小细胞肺癌(Non-small cell lung cancer, NSCLC)是肺癌的主要类型,约占所有病例的80%[1]。NSCLC包括几个组织学亚型,如腺癌、鳞状细胞癌和大细胞癌。microRNA(miRNAs)是19-25个核苷酸长的单链非编码RNA,通常通过与3’非翻译区域(3’-UTR)结合在转录后水平负调控基因表达[2-3]。人们发现在肝细胞癌和前列腺癌中miR-501-3p可能起抑癌作用[4-5]。此外,miRNA微阵列研究表明,miR-501-3p与肺癌骨转移相关[6]。生存分析是探究microRNA表达情况对肿瘤患者生存时间影响大小的常用方法。目前科研工作者在需要使用生存分析相关操作的时候多使用SPSS等专业的统计软件[7-8],其缺点在于专业性高、操作性差。Excel软件由美国的Microsoft公司开发,在全世界普遍使用,可以实现一般的t检验,F检验,线性回归和方差分析等常见统计方法[9-12]。本文旨在利用Excel软件中自带的Visual Basic(VBA)编程功能,分析肺鳞癌患者组织中的miR-501表达水平与患者生存时间之间的关系。
登录https://kmplot.com/analysis/网站,点击start miRPower for pan-cancer按钮,在新页面的Gene symbol框中输入hsa-miR-501。在Split patients by选项中选择median,使用肿瘤组织中has-miR-501表达量的中位数区分高表达和低表达,在plot options中选择Export plot data as text并在癌症类型中勾选Lung squamous cell carcinoma。其余选项保持默认即可。点击Draw Kaplan-meier plot按钮。在新页面中点击Export plot data as text即可获得原始数据,见图1。
图1 数据格式Figure 1.Data Format
Kaplan-Meier法又称为乘积极限法,用于未分组资料的生存分析,是一个充分考虑了在不同的分组中可能有患者中途退出实验时候的统计学处理方法。它可以充分考虑每个时间点上患者的生存率,因此比较适合用于肿瘤等慢性病的统计学分析[13-14]。生存曲线,即Kaplan-Meier曲线。它是以生存时间tk为横轴,生存率为纵轴,绘制而成的连续型的阶梯形曲线。中位生存时间,是指生存率为0.5时对应的生存时间,本文提供的案例采用内插法估计。
对数秩检验法用于对比两个或者多个生存曲线的位置是否有差异。它的基本思想是假定无效假设成立,即两总体的生存曲线位置相同,那么据此应该得出两总体中不同生存时间的初期观察人数和通过理论死亡概率计算获得的理论死亡数与实际死亡数应当相差不大,否则无效假设不成立,认为两生存曲线的位置不相同。该方法不要求生存时间服从某特定的分布,且是对整个生存曲线的比较,而不是某个特定时间生存率的比较,适合于本案例中按照has-mir-501表达量高低分组的肿瘤患者的生存时间分析。本文提供的数据表明两个分组的生存曲线没有相交,粗略预测可以通过ph检验,可使用对数秩检验法。
1.4.1 用户界面设计
在Excel自带的VBE开发环境中点击菜单栏>插入,插入3个用户窗体并使用VBE中的工具箱插入控件,如图2所示。
1.4.2 模块设计与程序逻辑
向该工程插入三个模块。“KM法统计和作图”模块中主要包含5个VB过程,单击窗体上的“统计描述并绘制生存曲线”按钮将会调用他们,其过程名字和参数及返回值如下所示:
图2 Excel插件设计界面Figure 2.Excel plug-in design int
“录入数据”模块包含2个VB过程,用于实现“删除数据并初始化表格”和“新增加案例”两个按钮的功能,其过程名字和参数及返回值如下所示:
“对数秩检验”模块包含 5个VB过程,用于实现“对数秩检验”窗体中提供的功能。其中LogRankTest过程用于调用其余4个过程。代码如下:
对数秩检验的算法步骤较多,操作步骤即图3所示。
以上诸多步骤和最后的画图功能由代码调用下面的其余四个vb过程实现:
以上各个模块各个函数的代码见附件1,读者可自《医学新知》官网(http://www.jnewmed.com/)相对应文章中数据与资料栏下载使用,同时包含具体操作视频。
1.4.3 插件的使用
图3 对数秩检验程序设计思路Figure 3.Design ideas of log-rank test program
在生存分析窗体上点击“删除数据并初始化表格”,输入案例数量2并点击确定,把附件2(下载查看方式同附件1)中两个案例的数据复制到Excel中,点击统计描述并绘制生存曲线,Excel会自动进行所有案例的生存分析并绘图,如果用户需要添加误差线可以选中生存曲线,点击图表工具>设计>添加图表元素>误差线>其他误差线选项,在设置误差线格式中设置垂直误差线的误差量为自定义,点击指定值,在正错误值和负错误值中选择误差线绘制数据一表当中的数据,可以根据需要选择方法一或者方法二产生的误差线长度。回到原始数据表格,点击“生存分析窗体中的对数秩检验按钮”。在弹出的对数秩检验窗体中选择两个案例然后点击“开始按钮”即可获得对数秩检验结果和生存曲线对比图。
本例生存时间以月为单位,并将t月当作一个时点看待。对分析结果表格中各栏含义解释如下:第(1)列为序号k。第(2)列是将生存时间t由小到大依次排列,如某时间点既有完全数据又有截尾数据,将截尾数据排在后面。第(3)列为t月的死亡人数dk。截尾患者即便是已死亡,也非死于研究变量,所以相应的d=0。第(4)列为初期病例数nk,即恰好在t时点以前尚存活的病例数。第(5)列计算各时点死亡概率qk,即在t时点以前尚存活的患者恰好在t时点上(第t个月)死亡的概率。第(6)列计算各时点生存概率pk,即在t时点以前尚存活的患者在t时点上(第t个月)继续存活的概率。第(7)列计算各时点生存率,即在t时点以前尚存活的患者活过t时点的概率。第(8)列为各时间点生存率的标准误。见表1。
可以采用两种方法估计某时点总体生存率的置信区间。方法一采用正态近似原理,使用生存率±1.96×标准误的方法计算每个生存率的置信区间,不适合于曲线尾部或接近尾部总体生存率的置信区间估计。因为此处的正态性较差,所估计的置信区间的上、下限值可能小于0或大于1。此时可以计算经过对数变换后的生存率值以及相应的标准误,据此来估计其置信区间,即方法二。见表2。
表1 极限乘积法估计miR-501不同表达量水平的生存率及标准误Table 1.The limit product method to estimate the survival rate and standard error of miR-501 at different expression levels
Excel把使用内插法计算获得的中位生存时间写入表格尾部。低表达miR-501的患者中位生存时间为37.484830。高表达miR-501的患者中位生存时间为74.371976。
生存曲线中水平横线的长短代表一个t时点到下一个t时点的距离,当最后一个时点的观察对象全部死亡时,曲线与横轴相交。生存曲线图可直观地比较各组观察对象的生存过程,也可对任意时刻的生存率作出粗略估计。绘制该图首先需要生成绘图专用的数据区域,该数据表格给出了每个时间开始和结束时的生存率。然后根据这些数据绘制xy散点图,将每个散点连线。对于截尾数据而言,其生存率没有变化,代表其时间头和时间尾的散点在图上重合为一点,根据这一特点添加删失记号即可完成绘图。见图4、图5。
最终获得结果见表3。
这两组生存曲线的P值小于0.05,表明两组数据的差异具有统计学意义。miR-501在肺鳞癌组织中表达较高的病人生存时间较长。生存曲线图如图6。
本文采用了乘积极限法和对数秩检验等经典生存分析方法对mir-501的表达情况进行分析,发现mir-501表达量较高的分组中肿瘤患者的生存时间较长。目前已有较多的学者研究了microRNA对肿瘤患者生存时间的影响,Shi ming 等人通过荟萃分析和生存分析研究了miR-30d-5p表达量对卵巢癌患者生存时间的影响,发现mir-30d-5p 表达较高的患者生存时间较长[15]。Hui Wang等人通过TCGA数据库获取生存数据对mir-21调控的 LZTFL1基因进行表达量与生存时间相关性的研究[16]。他们发现LZTFL1表达较高的患者生存时间较长。Rui-Sheng Zhou等人通过TCGA数据库的挖掘,探究了多个microRNA, lncRNA, mRNA对舌鳞状细胞癌患者生存时间的影响[17]。可见乘积极限法已经成为一种广泛使用的用于生存时间分析的统计学工具,配合对数秩检验这样的经典检验方法可以良好的区别某个基因表达情况对肿瘤患者预后的影响。
表2 miR-501不同表达量水平各个时间生存率的95%置信区间计算Table 2.Calculation of 95% confidence interval for survival rate of miR-501 at different expression levels at each time
图4 癌组织低表达miR-501的肺鳞癌患者生存曲线Figure 4.Survival curve of lung squamous cell carcinoma patients with low expression of miR-501 in cancer tissue
图5 癌组织高表达miR-501的肺鳞癌患者生存曲线Figure 5.Survival curve of lung squamous cell carcinoma patients with high expression of miR-501 in cancer tissue
表3 低表达miR-501与高表达miR-501的对数秩检验表Table 3.Log-rank test table of low expression miR-501 and high expression miR-501
图6 两组生存曲线对比绘图Figure 6.Comparison of survival curves between the two groups
本文利用了Excel软件中自带的Visual Basic(VBA)编程功能,以miR-501为研究目标分析肺鳞癌患者组织中的miR-501表达水平与患者生存时间之间的关系,拓展了Excel的使用范围,把生存分析这类医学特有的统计方式带入Excel等大众办公软件中,降低了生存分析的软件门槛,有利于生存分析方法的普及,极大地方便了一般科研人员对统计工具的使用。本程序在极端情况下可能会计算出错,例如在office安装不完整的情况下Excel 可能无法提供VBA二次开发功能。如果用户提供的数据中出现大量的删失数据,此时即使K-M法统计已经没有了统计学意义,但是本程序依然会按照算法得出不可靠的结果。经过测试,此程序需要在Excel2013以上版本才能正常执行。在完成Excel的对数秩检验后在两组生存曲线的绘制中可以手动在图片上加上P值,方便读者阅读。目前统计学已经在许多慢性病的科研项目中发挥了巨大的威力[18],许多专业的统计学软件已经大量运用到了医学科研中。
本文为统计软件的使用者提供了更多的选择,旨在起到抛砖引玉的作用。生存分析中还有一些其他的检验方法如Breslow检验等,本程序没有探索,读者可以自行查阅相关检验的原理并编程实现。
综上所述,使用Excel 结合Visual Basic语言编程可以较好地分析miR-501 表达高低对肺鳞癌患者生存时间的影响,其结果可信可靠。