|
致 歉 信 本文早在2020年3月13日就已经在本公众号发布，截止到2020年5月4日，已经有60人付费，由于本人的操作不当，意外把文章删除了。对已经付费的粉丝们，深感抱歉。其实，我早已将视频上传到B站(未提供代码)： https://www.bilibili.com/video/BV117411Z7Cq 相信付费的您已经拿到代码和相关文件，如果有什么问题可以联系我。 邮箱：bioinfocloud@aliyun.com DoubleHelix 2020.05.04   关于TCGA数据库的教程，前期我们已经推出了一些文章： 【2】R语言TCGA-Assembler包下载TCGA数据前面的这些文章虽然介绍了TCGA数据库的差异表达，但部分用了perl脚本，如果不懂perl的同学，数据库更新，数据格式可能会变换，运行脚本就得不到想要的结果，所以这里我们就只用R语言进行系统性的讲解，包括数据下载、整理、融合、基因ID转换以及表达差异分析，最后通过火山图进行可视化。所以你需要有一定的R基础，能看的懂代码，能改一些参数。 关于差异表达分析我们利用DESeq2和EdgeR包，其实在我们前面基因芯片数据挖掘序列文章中都已介绍： 可以先通过五~八这4篇文章了解这2个包的原理和使用教程。其实，你只要能从TCGA数据库下载的数据整理得到表达矩阵，参照七、八这2篇文章就可以得到差异表达结果。下面我们步入正题........... 一.数据下载 数据下载有3中方式，官网在线下载；官放下载工具下载；R语言包下载，比如：TCGAbiolinks。TCGA-Assembler等。我们推荐使用TCGAbiolinks，个人觉得这是挖掘TCGA数据比较好用的包。 ########################## 这里下载的是Counts数据，不是FPKM数据####################setwd('.')##包的安装if (!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager')BiocManager::install('TCGAbiolinks')# 加载相应的包，可能会需要其他包，提示错误就安装缺少的包。# 因为每个人已经安装的包不一样。library(TCGAbiolinks)# 请求数据。query <- GDCquery(project = 'TCGA-LUAD', data.category = 'Transcriptome Profiling', data.type = 'Gene Expression Quantification', workflow.type = 'HTSeq - Counts')# 从query中获取结果表，它可以选择带有cols参数的列，并使用rows参数返回若干行。# 594个barcode samplesDown <- getResults(query,cols=c('cases'))# 533个barcodedataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = 'TP')# 59个barcodedataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = 'NT')# 59个正常组织和533个肿瘤组织样本作为研究dataSmTP_short <- dataSmTP[1:533]dataSmNT_short <- dataSmNT[1:59]# 根据前面的筛选，再次请求数据queryDown <- GDCquery(project = 'TCGA-LUAD', data.category = 'Transcriptome Profiling', data.type = 'Gene Expression Quantification', workflow.type = 'HTSeq - Counts', barcode = c(dataSmTP_short, dataSmNT_short)) 网页筛选后添加购物车就可以直接下载了。最后一种下载方式是通过官方工具Data Transfer Tool： https://gdc./access-data/gdc-data-transfer-tool。 关于网页下载需要下载几个文件。  如果直接使用TCGAbiolinks包下载，并用该包进行数据分析，就没有这么麻烦。我们这里主要是讲解通过网页下载的数据，直接整理，分析。让大家知道更加详细的流程。只是下载数据的话，上面三种方式都可以，不过现在TCGAbiolinks包下载速度可慢了，而且这个包现在好像因为R版本的更新，部分函数好像会出错。网页版下载数据是把筛选好的样本数据打包下载，文件太大，网速不好可能会中断。但一般也没什么问题。 二. 数据的整理 我们网页上下载的数据是一个压缩包，需要解压。  解压后我们会看到很多文件夹：   解压以后是一个txt的文本文件，打开文本文件就是该文件对于病人的RNA-Seq数据，第一列是Ensembl ID ，第二列就是Counts数。   三 .ID转换 得到表达矩阵后，Ensembl ID我们看不懂，需要转换成我们能看的懂的基因名称，比如TP53。我们就需要通过gtf文件进行基因转换，人的基因注释文件下载地址可通过gencode：https://www./human/ 下载，NCBI和Ensembl 也可以下载。  但是我们需要注意的是多个Ensembl ID可能对于一个基因。所以我们需要去重,去重可以使用下面函数，取平均值。 当然也可以选择其中一个，删掉重复的数据，建议取平均值。 四，表达差异分析 得到这样的表达矩阵，我们就可以进行差异表达分析了。  DESeq2和EdgeR包进行表达差异分析需要原始的Counts表达矩阵文件，我们前面下载的也是HTSeq - Counts数据，如果是HTSeq - FPKM的数据可以用limma包直接分析，这2类数据在下载和处理得到矩阵文件过程是一样的。后面也会附代码。 进行差异分析以后，我们绘制火山图进行可视化。 |
|
|
