转载请注明:解螺旋·临床医生科研成长平台
今天是最后一题了,照例回顾一下作业:
用RTCGAToolbox包下载结肠癌的mRNAArray的数据,肿瘤类型代号:COAD;版本号:20140115;基因:IL17A。
图中所标的P值的获取方法在课件中没介绍,暂时不用做。
有兴趣的小伙伴还可以试试该文补充材料的FigS1,即同一套数据的RORC基因。
这题用R语言可能对部分小伙伴来说有点难,但由于在R里下载数据,跟后续很多种分析都能轻松对接,所以我还是希望大家能通过这题作业先熟悉R的下载和简单作图。
RTCGAToolbox这个包对windows不太兼容,开发组也一直致力于改进这些问题。而从生信研究人员们的反映来看,它是TCGA相关的包里相对比较稳定的了。而且……给大家写这段代码的老师据说是用Mac电脑,所以之前没碰到这么多妖蛾子@_@这回为了在windows上找到各种问题的解决办法也花了不少力气。我也用windows,我行你肯定也行~
步骤
加载RTCGAToolbox包,用getFirehoseData()按题目要求下载数据,储存在Data中。
课件中还有个clinical参数,但默认是下载(TRUE)的,没特别需求我就不写了,让它默认吧。当下载遇到提示“path too long”的时候,便去解压mRNAArray文件,找到名字长的那个,改成“20140115-COAD-1.txt”,回到R中再次运行上边那行命令。
这样便下载好了,接下来做箱线图。
课件中给了一大段命令,这段几乎可以当成祖传代码了~
写的时候把里面的PIK3CA都换成IL17A。如果你原来写好了这段命令就更简单了,在RStudio中按Ctrl + F,在上方的搜索栏里搜PIK3CA,点Replace替换成IL17A,或者点旁边的All全部替换也一样的。
替换完成后运行,就得到了一张图:
这里面老师已经把颜色、坐标轴、标题啥的都给咱们设好了。要模仿得更像文章里的图的话,就再改一下最后boxplot那句:
就是肿瘤组和正常组换个位置,ylab是Y坐标轴的标题,names是每个箱子的名字,按顺序填好,删掉颜色。再运行就得到了这样的图,调整一下窗口大小,图也就跟着调整了:
接下来,图上方的Export可以将图片导出成PDF或图像。
若上方点选了Image,则还可在弹窗中选择其他格式,SVG、EPS都是受杂志欢迎的矢量图。这样可以转到Adobe Illustrator中跟其他来源的图片拼成一张Figure。
关于纵坐标轴标签
教程里标记的是“mRNAArray”,文献里标记的是“Normalized IL17A expression”,可以说,后者会更加贴切,原因是我们拿到的TCGA数据,都是经过了标准化处理的。本例中,mRNAArray的数据是经过Lowess标准化(normalization)的。还记得上边关于提示“Path too long”的处理步骤吗?其中提到要解压一个压缩文件,而解压后的文件夹,名字超长的,这个名字里就有数据标准化方法的信息(见下图)。
关于分组
TCGA对样本的分组信息,是隐藏在样本的编号里的。TCGA样本编号解读:
图片来源于https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode
其中第一个粉红色框中的数字便是区分肿瘤和正常的编码,01-09为tumor,10-19为normal,所以我们的教程中就有了这一段命令:
sampleIDs1 9则提取出来命名为normalSample组,
而这两句之前的一段for循环程序,则是为了把粉红框里的数字提取出来,转换为数字类型(原本是字符串),这才能通过比较大小来分组。
稍后我们还会推出视频教程,进一步教授大家使用R语言分析TCGA数据的操作,届时会详细讲解代码的意义,这里就先不展开了。
关于P值
下面讲讲作业中没要求大家做的P值,文章中基因差异分析是用经验贝叶斯方法(empiracal Bayes)。课件中所用的基因表达差异分析函数,getDiffExpressedGenes(),其实也是用的这个方法,但它将分析过程和全部基因的结果隐藏在内部,然后再根据我们设置的条件,把符合条件的结果展示出来。
当我们用getDiffExpressedGenes()函数分析基因差异时,目的是找到具有统计学意义的差异表达基因,因此我们设置了几个阈值(比如,限定了p值
而当我们的目的是计算某个基因的p值时,我们需要修改一下思路。原因是我们需要计算p值的这个基因,它的表达情况可能因为不符合我们设置的阈值而没有被筛选出来,因此,我们可以不设定阈值,即获取所有基因的数据,然后从这些数据里找到我们需要的这个基因的数据。
将getDiffExpressedGenes()函数的参数重新设置一下,不定义筛选阈值:
注意,不定义筛选阈值≠不填参数,因为如果不填,系统就会按默认值筛选(比如logFC,默认值是2)。这里把logFC设为0,意即显示变化倍数大于2的0次方=1,其实就是不管有没有变化,都要呈现出来;P值设为1,即p值在1以下的结果都呈现出来,也就是给出所有结果。
随后,将这些结果赋值给“DG”对象。
在RStudio中点开DG,便可以看到,结果中包含了17814个基因的分析数据:
在DG这个数据框里,每一行的行名就是基因名,所以可以通过提取DG的某一行结果,达到显示我们感兴趣基因的分析结果的目的:
文献中的p值是0.01,显然作者用的是未校正的P值(见图中红框)。
需要注意的是,getDiffExpressedGenes()只能分析RTCGAToolbox下载的数据,用其他来源的数据可能会报错的。
最后可以在boxplot()后面加一行text(),把P值填到图上。
0.7和1.9是自己从图上估摸的坐标位置,后面P = 0.01是将要写上的文字。这句写不写都行,反正你极有可能是要跟其他图片一起放到Adobe Illustrator等软件中排版的,到时候再加上也方便。
领取专属 10元无门槛券
私享最新 技术干货