主讲老师第十五课：TCGA数据的下载与清洗
到此，万事具备，所有的准备工作已经基本完成，就差数据了，当我们想要发一篇生信文章的时候，第一步该想什么呢？巧妇难为无米之炊，要发生信文章需要有自己研究领域的数据集对不对？这就解释了，为什么癌症领域的生信文章特别好发，而一些小众领域的生信文章发的比较少？因为做肿瘤的人多，相应的 GEO 数据库上的数据集也多，而且有了 TCGA数据库及其衍生数据库，有了合适的数据集，
才能进行后续的分析（TCGA.zip）。
1.TCGA 数据库简介
两篇文章题目中明确了这是一项基于 TCGA 数据挖掘的研究，说到 TCGA 数据库，大家绝对是再熟悉不过了，最常用最适合小白的就是 TCGA 数据库，我们总是会用到 TCGA 数据库，然而大家真的了解TCGA 吗？很多 10 分+，乃至 CNS 级别的纯生信文章，TCGA 都是用到的主流数据库，今天，给大家简单总结一下TCGA 存储的数据包括 SNV、转录组分析、生物样本信息、原始测序数据、CNV、DNA甲基化、临床信息等。这些数据可分为三个级别，Level 1: 原始的测序数据（fasta,fastq 等）；Level 2:比对好的 bam 文件；Level 3:为经过处理及标准化的数据；对于这三个层次的数据，我们一般只能下载到 Level 3 的数据，当然，有些大的实验室也可以通过申请得到前两个层次的数据，那么这么多类型的数据该如何选择呢？下面为大家讲解一些 TCGA 数据库必备的知识①Clinical 数据：临床数据 TCGA 数据库是非常全面的，这也是为何 TCGA 数据库用来做生信分析的最主要原因之一，包括病人的基本资料,治疗进程,临床分期(TNM),肿瘤病理,存活情况等
②mRNA 表达矩阵：这是大家十分熟悉的内容
③microRNA 表达矩阵：包含由 microRNA-Sequencing 测序得到的表达量④Copy Number 肿瘤拷贝数变异数据：包含由 SNP microarray 测序得到的肿瘤对比正常组织染色体各片段的比值这一部分对于刚刚生信入门的小伙伴们可能有些陌生，简单介绍一下，这部分数据在很多高分泛癌生信文章中经常被使用，主要用于什么呢？拷贝数变异在肿瘤中发挥着重要作用，肿瘤样品中有着广泛的CNV 变化，这些变化可在一定程度上用于区分肿瘤细胞和正常细胞，肿瘤的形成过程中涉及到了多种类型的基因组变异，比如点突变，拷贝数变异，基因融合等等。肿瘤和遗传病不同，各种基因组变异是后天形成的，所以在肿瘤研究中，关注的是体细胞上的基因组变异，并且在一些单细胞分析研究中，可能通过 CNVs 发现原发和复发肿瘤之间的重叠，肿瘤细胞亚型的系统发育关系，也可以发现哪些突变发生在早期，而哪些具有队列特异性。
⑤Mutation：肿瘤组织测序数据相对于参考基因组序列得到的核苷酸变化，包括 insertion、deletion以及氨基酸的变化，目前所有 TCGA 的 somatic mutation 文件都是可以下载的，里面包含的 somaticmutation 非常多，都是 MAF 格式记录的⑥Protein：这部分数据包含由 reverse-phase protein array 测序得到的约 200 种常见癌症相关蛋白的表达量⑦Methylation 甲基化数据：包含由 methylation microarray 测得的 DNA 甲基化程度以上这些就是比较常用的数据类型，当然，除了这些以外，还可以下载到部分患者对应的病理 HE 图片，影像 CT，MRI 等图，我们还可以选择感兴趣的区域或概要，通过分子异常表达、样本特征、临床变
量等，进行基于两个或多个数据的关联分析
2. TCGA 数据的下载
关于 TCGA 数据库中转录组数据的下载，我们分别通过在线下载和代码下载两种方法来进行讲解。
2.1 在线下载
目前对于TCGA数据下载的途径越来越多，包括TCGA官网提供了专门的数据下载窗口，即 TCGA GDC (https://portal.gdc.cancer.gov/)以及著名的 UCSC Xena 网站，在这里，我们将一一介绍几种常见的下载方式
2.1.1 TCGA GDC 网站
1). TCGA GDC 网站的介绍
1.登入网站，我们可以总体浏览网站所包含的疾病、样本数等基本信息
这里需要提醒一下，在数据下载前，首先确认右上角的 Cart 已经清空，防止与上次下载的数据混淆，点击“repository”按钮，进入下载界面。
2).在 Files 和 Cases 中选取疾病和数据类型https://portal.gdc.cancer.gov/ 如果进不去网，先听讲，这部分快速过一下。
这里我们以胃癌(STAD)的转录组数据为例，TCGA 官网提供了数据储存库，直接进入TCGA GDC 下载。
在 Cases 栏中，Primary Site 中选择 Stomach，Program 中选择 TCGA，Project中选择 TCGA-STAD。
当然后面也可以对疾病的具体分型，性别，年龄，存活状态等进行选择，一般我们不做选择。
在 Files 栏中，Data Category 中选择 transcriptome profiling（转录组分析），Data Type中选择GeneExpression Quantification（基因表达量），Workflow Type中选择HTSeq-Counts。
这里需要注意几点：
1）在Gene Expression Quantification中既包括能编码蛋白的mRNA数据，同时也包含了非编码的lncRNA 数据，在后期的数据清洗中可以进一步的分离；2） 对 于 miRNA 的研究数据，在 miRNA Expression Quantification 中选择进行下载，但是，注意了，不能同时点击两个一起下载；3）在 Workflow Type 中，Counts，FPKM 和 FPKM-UQ 分别代表三个不同的数据呈现形式，对于差异分析，我们一般选择下载 Counts 用于后续的研究分析，因为目前的三个测序相关数据的差异分析 R 包都是基于 counts 数据来进行设计的，而如果对于单基因的研究，我们也可以选择经过标准化处理的FPKM 数据形式或进一步标准化后的 TPM 数据，对于 FPKM 和 TPM 两个的选择，两者都没有错，都可以用，只是目前更倾向于 TPM 数据，对于这三个数据之间的转换和具体差异，我们等等再讲3).在选择好需要的数据集后，点击“Add All Files to Cart”将所有需要的样品加入到 Cart，这操作和淘宝购物实际上是类似的，加到购物车里。点击右上角的Cart，进入下载界面。
下面来解释一下几个文件的用处：
1. Manifest：解释文件，样本信息；
2. Cart：每个样本中的基因表达文件，其中 cart 文件是最大的，是所有基因的表达文件，受网速影
响，如果下载失败的主要也是这个文件；
3. Metadata：提供样本名称对应的 TCGA 的 ID；
4. Clinical：样本对应的临床信息，关于 clinical 可以在这里简单下载，也可以下载原始的临床信息文
件，加入购物车进行下载，区别在于 file 选择问题，选择 xml 文件，然后同样下载这三个文件，一般来讲，这是最简单直接的下载方式。Cart 文件相对比较大，对网速有一定的要求
2.1.2 UCSC XENA 网站
接下来，我们将介绍另外一个在线网站中 TCGA 数据的下载方法。对于新手而言，从UCSC XENA 下载数据是十分友好的，因为与前面方法不同，该方法可以节省出大量整理的时间，但需要注意的是，UCSC数据库下载的数据是处理后的数据，需要看清楚相应的数据说明。
1.首先打开 UCSC Xena 网站
2.点击 Launch Xena，进入下方界面
3.点击 Dataset，进入数据下载页面
4.找到相应的数据集，这里找到文章中的 TCGA LIHC 数据集
5.在这里，建议大家选择 GDC 入口的 TCGA 数据集
可以看到 LIHC 数据条目下提示有 14 个数据集这里包含了 CNV 或者甲基化等数据，找到表达数据的 FPKM 数据。
点击红色方框进入表达数据，当然，如果选择 count 的话，也可以对应的进入其中。
6.点击 HTSeq-FPKM 进入数据页面
在数据介绍页面我们可以看到这份数据是经过了 log2(x+1)处理的结果，所以，如果进行后续部分分析时，需要对其进行还原，尤其是 count 数据做差异分析，必须把它还原回去
7.点击下方红框内容进行数据下载，结果即为 TCGA LIHC 的 FPKM 数据
8.接着是临床数据的下载：返回上一页面，找到临床数据
在 UCSC XENA 中，其把 TCGA LIHC 数据的临床数据整理为两个文档，放在 phenotype条目下，其中Phenotype 里面包含了临床病理参数数据，Survival data 条目是相关样本的生存时间和生存状态数据，
在进行相应分析时候，需要把两部分结果都下载
9.点击 Survival data
网页下方有相应的表格预览，点击红色方框中的链接，即可下载生存数据，同理，Phenotype 的数据也是点击相应位置即可下载
2.2 代码下载
虽然我们可以通过在线下载 TCGA 数据库的数据信息，但是其对网速的要求相对较高，下载过程往往容易中断，因此，在这种情况下，学习一下如何使用代码下载 TCGA 数据则十分重要。
1). GDC.exe 软件的使用首先，通过 dir.create()函数完成新建文件夹的操作，顾名思义，dir 为 directory（目录），dir.create =create directory（新建目录）运行可以发现在工作目录下多了一个名为“01.GDC”的文件夹，用于保存 GDC.exe 软件下载得到的数据结果，随后，我们来看一下如何下载 GDC.exe 软件，打开 TCGA GDC 网站，可以看到下载链接根据自己电脑的类型，下载适合的软件，并将其放到“01.GDC”的目录下，这里，我已经把对应的两个版本的软件进行下载其中，有.exe 后缀的是 windows 电脑，没有后缀的是针对 mac 电脑的，同时把之前下载的关于counts 数据的 manifest 文件复制黏贴到当前文件夹下新建两个文件夹分别用于保存下载下来的表达数据和临床数据信息接下来，我们开始数据下载，用 windows 电脑为例。
1.在下方的搜索框里面找到终端命令
点击图标，从而打开电脑终端
接下来，复制文件夹路径的全称：
通过 cd 工作路径 来设置工作目录注意中间有个空格，然后回车。
可以看到，此时的目录已经变成了 1.GDC 文件夹。
2.输入代码，分别下载 counts 表达数据和临床数据，举个 count 数据的例子
分开来一下看下组成，gdc-client.exe 软件名字，download -m 终端下载命令，STAD.gdc_manifest.counts_20201109_170902.txt 这个是前面下载的 manifest 文件名，下载临床数据时将其进行替换即可，文件名不需要复制，输入几个字母，然后直接 tab 键会补齐。注意，把文件后缀名也要输入进去的，完成以后，回车就可以开始进入下载。
下载完成后，可以看到下载的文件是按样本存放的，这里，关于 count 数据和临床数据已经给大家下载好了，课后可以尝试，课上直接用我给大家的数据来进行后续分析就可以了，同时，通过 length() 函数快速检查两个目录下文件夹的数量，即为患者的数量结果显示，表达数据共包含 408 例样本，临床数据共有 444 例患者，在 TCGA 数据中，一个患者是一个文件夹放置的。
2) TCGA 数据的清洗
1.在正式分析前，我们先读取其中两个数据来提取前两个数据先来看看其中的规律
首先，我们随机读取两个样品的数据，看下规律，这样其他的就可以进行循环完成可以看到，在两份 count 文件中，其行和列的数量是完全一样的，查看两个数据集的第一列。
接着，我们来看看两个是不是相同通过 identical()函数判断得到两者的第一列完全相同，因此，在后续的分析中，我们可以每个文件读取进来后直接使用 cbind()函数将两者进行结合。
2.接着，我们开始批量读取所有的 count 文件，进行整合
首先，定义一下文件名的普遍规律，其中，中间 pattern 部分为正则表达式，关于正则表达式的内容，我们明天晚上专门设置一个专题来讲解，解释一下这个语句最后的含义，匹配名为“1.expdata”文件夹下以“htseq.counts.gz”结尾的所有文件名共得到 407 个变量名，这与文件夹的数量是一致的，接着，我们通过一个循环函数，来批量读取文件里面的内容。
file.path("1.expdata/",x)表示当前文件夹下文件的名字，然后，借助一下 lapply()函数的作用。
将文件名作为输入变量，通过 lapply() 函数将 count_files 中的变量名作为 x，逐个代入自定义函数ex 中，每次循环结束将输出一个与前面单独读取时类似的结果，最终组合成为一个列表（ List ），由于每个输出结果的第一列信息都是一样的，我们直接使用 cbind 将其合并到一起do.call()函数的作用是使用 cbind()函数的作用，对 exp 列表进行整理。
这样，一个新的表格就生成了，点击 exp，打开来看下：
其 中 行 名 为 基 因 的 Ensemble id 号 ， 比 如 KRTAP1-5 的 Ensemble id 号 是ENSG00000221852，但是，但是这样产生出来的表达矩阵没有列名的，列名还需要进一步的进行调整，此时，我们需要一个文件名与样本 ID 一一对应的文件，这就需要用到下载的另外一个文件，即 metadata 文件。
由于 metadata 文件是 json 格式的文件，这里需要借助 jsonlite::fromJSON()函数的作用来进行读取，成功读取后，查看一下该文件的列名信息。
在 metadata 文件中，associated_entities 是个列表。
来看看这个列表里包含了什么？
对列表取第一个子集。
可以看到，在该列表里包含数据框，数据框的第四列内容就是 TCGA 样本的 id 号，对列表的自己，进一步取子集。
取数据框的第四个子集，最终可以得到 TCGA 的 id 号。这里需要注意的是，TCGA 样本的 id 号存在于列表 associated_entities 中，但具体的列数不同肿瘤中有差别，需要具体问题具体分析，找到规律以后，我们进一步进行循环提取，再次使用 sapply() 函数，提取 ids 列表中每个元素的第四列信息，并返回一个新的向量。
随后，将文件名与 id 号进行一一对应。
此时，文件名与 TCGA 样本 ID 的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名，需要找到读取文件的顺序，一一对应修改比较文件名和 count_files 名字，可以发现两者的后半部分是一样的，因此以“/”作为标准进行拆分，取后半部分，并使用 match()函数进行排序。
到此，准备文件基本完成，把 file2id 的 ID 对表达矩阵 exp 的列名进行赋值这样整个表达矩阵就整理好了，接着，我们来对临床信息进行提取。
3.临床信息的提取
临床信息的文件保存在.xml 文件中，当然直接打开也可以看到粗粗查看其中的内容，但是，每个 xml 文件分散在每个单独的文件夹下，我们同样需要将其进行整合，对于临床数据的整理，我们同样从一个中查找规律，然后通过循环函数，最终完成提取。
使用 XML 包的 xmlParse()函数读取一个文件看看效果，经过处理，我们可以通过 rootsize 看到其中文件存在两个节点，使用 print()函数可以分别查看每个节点里面的内容。由于每个里面包含的内容比较大，这里就不做展示了，感兴趣的大家可以自行运行，而我们后期分析所需要的内容主要存在于第 2 个节点中进一步进行提取，最终得到了我们所需要的文件形式，使用 colnames()函数查看文件的列名，一共包含了 55 种不同的临床信息可以看到其中包括生存，随访，以及相关的临床资料等等，到此，单个文件的分析成功后，我们同样将整个过程与表达文件整理一样，改装成为一个能够批量快速读取的函数经过批量读取，自定义函数，以及 cbind()函数的结合后，最终得到了整个的临床信息表格，整个处理的过程相对十分复杂，对初学者而言不是十分的友好，因此，我们另外介绍一个相对简单的方法来进行数据的下载与整理，上面的处理方法对有一定基础的同学可以学习一下，对初学者不做要求，了解一下整个处理过程即可，不要过分纠结里面的内容。
3.TCGAbiolinks 包
接着，我们来介绍一个相对常用且简单的方法，除了以上的方法，还有一个工具，集TCGA 数据下载处理可视化为一身，就是 TCGAbiolinks 包。
3.1 对项目进行探索
首先，我们来看一下其中有哪些项目的数据可以进行下载。
安装完成后读取 TCGAbiolinks 包，然后通过 getGDCprojects()函数看下其可以下载的所有数据内容。
列表中展示了可供下载的项目内容，而每个项目中又包含了多种组学的数据进行下载，而且，除了TCGA 项目外，还有一些其他数据库的项目内容。下面，我们来提取其中的 TCGA数据项目看一下内容可以看到，其中包含了 33 个 TCGA 的项目与肿瘤内容。
3.2 下载单个 TCGA 肿瘤项目数据
接着，我们来下载单个的数据。还是以胃癌为例，下载表达量数据。
在新的分析项目之前，首先新建一个名为“02.TCGAbiolinks”的文件夹，接下来的工作目录设置为“02.TCGAbiolinks”，所有内容都在此进行1).使用 GDCquery()函数来查询下载信息2).使用 GDCdownload()函数来下载数据下载的数据，会默认在当前工作目录创建 GDCdata 文件夹，并存放在里面这一步会需要一段时间来下载 110mb 的文件3).使用 GDCprepare()函数批量读取数据此时如果想把读取的数据保存为 Rdata，下次直接加载，可以把 save 参数设置为 TRUE， 且save.filename 给一个名称。
同时，我们需要从 pre.exp 这个对象中提取表达量的信息。
4).使用 assay()函数提取表达量信息同时，我们查看一下前 3 行前 3 列，最终得到行是基因，列是样本的表达量矩阵（代码都整理好了，包括后面的差异分析，只要改个参数，其他全部一键跑完，绝对的一站服务）。
整个下载整理过程就 4 个函数，唯一要修改的是第一个函数中 cancer_type 的名字，将 其改成自己需要的肿瘤名字，其他直接往下运行，即可得到最终的 count 表达矩阵，相对来说难度小了很多，关于临床信息的下载，后面会有个具体实例来进行演示。
4. RTCGA 包
除了这几种方法外，另外还有一个 R 包的下载。来简单介绍一下 RTCGA 包的使用方法。由于 RTCGA 包是先下载所有已经储存好的数据，因此在使用过程中是直接从中进行提取的加载相应的 RTCGA 包后，查看一下每个项目中包含哪些内容。
可以看到，上面罗列了使用 RTCGA 包可以下载的所有数据，直接输入上方的项目名称，即可调取其中的表达数据和临床数据。
可以快速在环境变量中看到相应的表达矩阵和临床信息，尽管使用过程相对方便，但是 一般不建议使用。原因一是由于使用的是下载后的数据，因此后续的数据已经停止更新第二，其中的表达数据是经过了标准化处理后的数据，而且，其对基因进行了一些不知名的过滤条件，下载到的不论是基因信息，还是临床信息，和目前 TCGA 下载得到的有差别所以，如果只是简单快速看看时可以使用一下，其他时候还是不要用了。目前，主流的 TCGA 数据下载方式是这几种，但是 GDC 软件和 TCGA 网页下载一样，受到 TCGA 网速的影响，TCGA卡它也卡，因此一般来说推荐用 TCGAbiolinks 下载