主讲老师第十七课：TCGA临床基线表整理
在学会了数据下载和清洗转换后，接下来，我们来看下另外一块新的内容，也就是 TCGA患者临床数据的整理。在生信分析的文章中，往往第一张表格是基于 TCGA 数据库患者临床信息的基线资料表，以展示分析的整体水平在生信分析的文章中，往往第一张表格是基于TCGA 数据库患者临床信息的基线资料表，以展示分析的整体水平，下面，我们来看下如何基于 TCGA 数据的临床信息绘制临床基线资料表（TCGA2table.zip）。
1.TCGA 临床数据下载
首先，使用 TCGAbiolinks 包下载 TCGA-LIHC 患者的相关临床信息，首先第一步还是使用 GDCquery()函数来查询下载信息这里注意一下这是下载临床数据的相关参数设置，其中，参数 project 与昨天的一样，写入肿瘤名字，剩下的则是临床数据下载的默认参数，注意区别，接着，使用 GDCdownload() 函数来下载数据在批量读取数据中，函数略有改变使用 GDCprepare_clinic()数据来进行临床数据的批量读取。下载到这一步，大部分需要联网完成的操作基本已经结束，为了避免后续使用时再次下载，可以将其保存为.rda 文件。
可以看到，在当前目录下多了一个.rda 文件，后续在使用时直接使用 load()函数进行读取即可。
在结果中，其中包含 415 个患者的 82 个不同的临床信息。为了后续使用方便，将临床信息的第一列列名修改为“Tumor_Sample_Barcode”。
查看前 3 列前 3 行，并保存为 csv 文件
2.病人分组筛选
接下来，读取保存好的临床信息，用于后续的临床数据整理和分组筛选。
在这里，为了演示方便，我们通过存活和死亡的方法来对所有患者进行分组作为演示简单解释一下，分别根据临床信息数据框中“vital_status”的信息，筛选“Alive”的患者作为存活组别，筛选“Dead”的患者作为死亡组别，同时去除其中缺失的患者后，作为纳入的所有患者，当然，在实际应用过程中，也可以根据自己的数据需求，对患者进行重新分组可以看到，最终纳入分析对患者中，存活的患者为 315 例，死亡的患者为 100 例，总计 415 例
3.选取性别分组
接着，对临床数据进行整理，首先对当前分组下性别数据内容进行整理这里，使用 group_by()函数结合 summarise()函数，分别统计存活 Alive 的患者中男性和女性数量，结果显示，其中女性 96 个，男性 219 个。
同样，对死亡患者进行整理，结果显示，在死亡 dead 的患者中女性 40 个，男性 60 个接着，使用 merge() 函数对两个统计分析结果进行合并接着，将第一列作为行名，并去除第一列随后，使用 chisq.test()函数进行卡方检验，对统计结果进行分析，计算对应的 P 值结果显示，经过卡方检验，P 值为 0.099，说明性别在死亡和存活患者之间无显著性差异。最后，对细节信息进行整理。
统计添加所有患者数量，分别对行名和列名进行设置，同时添加百分比。
其中，通过 sprintf("%1.1f%%", sex[2,]/sum*100)，对百分比保留一位小数点，到此，关于性别的整理就基本完成了。同理，如果有其他的二分类变量，可以进行一样的延伸整理
4.选取年龄分组
随后，选取一个数值型变量，来看看如何对数值型变量进行整理分析，在此，我们来对患者的年龄情况
进行统计分析。
使用 summarise 函数来进行统计，包括 mean，sd 值，中位值，min 和 max，同时使用round()函数来设置保留小数点后一位将计算得到的结果使用 paste0()函数进行整合成标准形式
同样的，对死亡的患者年龄进行统计分析
最后，对所有患者的年龄进行统计分析，这样，在环境变量中多了三个年龄的计算结果。
使用 cbind() 函数将三个的统计结果进行合并，并将性别 sex 中整理好的列名赋值给年龄 age两种常见的数值型变量的展现形式。在实际使用过程中，任选其中一种方法进行展示即可
5.选取分期分组
除了二分类变量外，接着讲一下如果存在多个分类等级时如何处理。在肿瘤分析过程中，肿瘤分期是个重要信息，接下来，对肿瘤分期进行分组整理筛选-分组-统计，是不是和当时那份作业的内容十分像？去除分期 Stage 中缺失的患者，并按pathologic_stage 分组情况进行分组，然后使用 summarise() 函数进行统计汇总。同样的，我们对另外一个分组以及所有的患者进行一样的整理。
接着，对三组临床分期进行合并，这里介绍两个用于合并的函数，还是之前那句话，哪个用着简单，用
着方便，后续分析用那种即可，因为两个的结果是一样的
方法一： merge() 函数
注意一点的是，记得加上参数 all = TRUE，表示即使另一组为缺失值，这个变量分组依然会保存，举个例子，比如，I 期患者存活的有 90 人，死亡的 0 人，那么不加这个参数的话，I 期会因为死亡组为 NA，自动抹去这个分组，具体的可以尝试下去掉这个参数后，比较前后差异。
方法二： full_join() 函数
区别在于其不需要加 all 参数，会自动保留所有分组。
随后，对行名和列名进行设置，将最终结果打印出来看下。
随后，计算两组之间的统计显著性差异。既然是计算存活和死亡之间组别的差异，那么自然需要将第三组 total 去掉，并将其赋值给一个新变量。
去除包含 NA 的分组，随后使用 chisq.test()函数进行卡方检验。
最后，将最后的结果汇总输出首先，将缺失值以 NA 代替，同时与性别一样，计算其各个值所占的百分比，并保留一位小数，随后，对行名和列名进行相应的设置。
6.临床信息合并
最后，使用 rbind() 函数将性别，年龄，和分期信息进行合并，并使用 as.data.frame() 函数将结果转换为数据框格式。
7.输出结果
所有内容都整理好了，最后将结果以 csv 形式输出，并命名为 TABLE1.csv。
注意，使用参数 quote = FALSE 将引号去除，其余默认，可以看到，在当前工作目录下多了一个名为TABLE1.csv 的结果文件（TABLE1.csv）。
这样，整个临床信息的基线资料表就整理完成了，当然，如果你对其他信息感兴趣的话，也可以进行相应的整理。
biomart包简介：bioMart包本质上其实就是一个链接bioMart数据库的R语言接口，能通过这个软件包自由的连接到bioMart数据库作用这里其实大家可以简单理解就是把任意数据库之间的ID号进行转换（个人认为最主要，可以和Y叔的clusterfile包中的biter函数互补）可以通过在线网站BioMart对bioMart进行更多的了解主要包含以下功能#安装R包source("http://bioconductor.org/biocLite.R")biocLite(“biomaRt”)主讲老师解读：然后既然这个R包是一个功能型R包，那么我们就来探索一下这个R包的更多内容#查看其中包含的数据库资源listMarts()主讲老师解读：当然探索过多对我们的使用并没有太多的作用，接下来我们将通过几个实例帮助大家掌握这个R包#GEO下载芯片数据，从基因表达矩阵中分别提取lncRNA、mRNA和miRNA的表达矩阵#并且根据基因名提取对应的表达矩阵library(GEOquery)gset <- getGEO("GSE36895", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)#这里我们需要查看一下我们下载的数据里包不包含与gene symbol或其他任何跟ensembl共有的ID，就需要用序列比对的方式去找这段探针组所对应的基因名，才能提取lncRNAstr(gset)
主讲老师解读：从这里我们可以知道两个细节
第一：我们getGEO函数 AnnotGPL 函数选择TRUE的时候，会在下载的ExpressionSet中添加注释信息第二：我们如果想要进行Ensembl ID的转换必须要有中间桥梁，也就是Gene Symbol等可以与之有对应关系的#提取表达矩阵exprdf<-data.frame(exprs(gset))dim(exprdf)#write.csv(exprdf,"not_easy_input.csv",quote = F)#这个时候我们输出的行名是探针名，变量名是GSM，所以我们要通过注释信息进行转换exprdf$gsym<-gset@featureData@data$`Gene symbol`#同时我们还可以对表达矩阵进行修饰exprdf<-exprdf[exprdf$gsym!="",]#删除没有Gene Symbol的行#删除同个探针对应多个基因的行exprdf<-exprdf[!grepl("///", exprdf$gsym),]#多个探针对应同一个基因，取中间值#如果想取平均值，就把median改为meanexprdf_uniq<-aggregate(.~gsym,exprdf,median)dim(exprdf_uniq)#整理完毕后，把行名替换为Gene Symbol并删去Gene Symbol列rownames(exprdf_uniq)<-exprdf_uniq$gsymexprdf_uniq<-subset(exprdf_uniq,select = -gsym)#保存#write.csv(exprdf_uniq,file="gene_exp.csv",row.names = T,quote = F)#为了后续方便，这里仅仅保存前四个样本write.csv(exprdf_uniq[,1:4],file="easy_input.csv",row.names = T,quote = F)主讲老师解读：这样的话我们就得到了行名为基因名，列名为样本名的表达矩阵因为有些芯片是表达谱芯片，这里的范围其实就很宽泛，包括mRNA、lncRNA、miRNA等等但是也有的芯片介绍是mRNA，那么就不用提取了其实依照经验提取这些本身就是用来分析，如果你不提取，在最后进行富集分析得到感兴趣的结果后再识别也是可以的exprdf_uniq<-read.csv("easy_input.csv",header = T,row.names = 1)rownames(exprdf_uniq)[1:4]#[1] "A1BG" "A1BG-AS1" "A1CF" "A2M"
主讲老师解读：接下来就涉及到 biomart 包的使用了
#加载R包library("biomaRt")#显示可以连接的数据库#因为是在线的，所以运行速度很大部分取决于网速listMarts()#这里我们就可以选择数据库，这里的操作其实和SingleR包很类似，都是需要选择的#选择完后来查看当前数据库所包含的基因组注释plant = useMart('ensembl')#选择大数据库listDatasets(plant)#查看大数据库下的子数据库# dataset#1 acalliptera_gene_ensembl#2 acarolinensis_gene_ensembl#3 acitrinellus_gene_ensembl#4 amelanoleuca_gene_ensembl#5 amexicanus_gene_ensembl#6 anancymaae_gene_ensembl#7 aocellaris_gene_ensembl#···#···#54 hsapiens_gene_ensembl
主讲老师解读：根据样品物种来选择后面的参数
#你需要哪个物种，就复制它在dataset列里的词，放在下面这行的`dataset = `参数里ensembl <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",dataset = "hsapiens_gene_ensembl", #人#dataset = "mmusculus_gene_ensembl", #小鼠#dataset = "rnorvegicus_gene_ensembl", #大鼠#dataset = "dmelanogaster_gene_ensembl", #果蝇host = "www.ensembl.org")主讲老师解读：这里我们还是需要解释一个这个函数 useMart 函数useMart(biomart,dataset,host="www.ensembl.org")#biomart参数表示BioMart所用的数据库名称#dataset参数表示biomart参数选中的大数据库中的小数据库#所以这里我们是先选择Ensembl数据库，然后在选择人的数据库#使用下面的函数可以搜索特定的小数据库类型searchDatasets(mart=ensembl,pattern="mouse")
主讲老师解读：
我们需要用gene symbol或ensembl ID来提取gene_biotype然后用gene biotype区分lncRNA、miRNA和mRNA#查看Filters和Attributes提供了哪些信息listFilters(ensembl)#获得所有输入结果（即我们已经有的结果）# name description#1 chromosome_name Chromosome/scaffold name#2 start Start#3 end End#4 band_start Band Start#5 band_end Band EndlistAttributes(ensembl)#获得所有输出结果（即我们想要的结果）# name description page#1 ensembl_gene_id Gene stable ID feature_page#2 ensembl_gene_id_version Gene stable ID version feature_page#3 ensembl_transcript_id Transcript stable ID feature_page#4 ensembl_transcript_id_version Transcript stable ID version feature_page#5 ensembl_peptide_id Protein stable ID feature_page
主讲老师解读：这里其实就是我们输入参数和输出参数的选择
#可以看到输出结果Gene Symbol在62行listAttributes(ensembl)[62,]# name description page#62 arrayexpress Expression Atlas ID feature_page主讲老师解读：在使用 getBM() 函数获得注释的时候，我们需要先了解一下这个函数
attributes 参数：我们要获取的数据类型
filters 参数：已知的数据类型
values 参数：想要转换的基因名
mart 参数：检索的数据库
feature_info <- getBM(attributes = c("gene_biotype",#"transcript_biotype",#还可以提取transcript_biotype#如果基因名是gene symbol，就运行下面这行"hgnc_symbol"),#如果基因名是ensembl ID，就运行下面这行#"ensembl_gene_id"),#如果基因名是gene symbol，就运行下面这行filters = "hgnc_symbol", #小鼠是mgi_symbol，大鼠是mgi_symbol#如果基因名是ensembl ID，就运行下面这行#filters = "ensembl_gene_id",values = rownames(exprdf_uniq), mart = ensembl)#这个函数中的参数和前面那两者listxx函数是对应的关系主讲老师解读：这里我们的输出数据的类型选择两个（输出多个结果在 attributes 参数内添加多个参数即可）最后我们会得到这样一个表格#有些gene symbol对应多个ensembl id，因此会有多个biotype，例如feature_info[feature_info$hgnc_symbol == "ARHGEF26-AS1",]#TCGA数据不会遇到这个问题，因为ensembl id跟gene_biotype是一一对应的关系#然后我们统计一下独一无二的数据个数length(unique(feature_info$hgnc_symbol))#把基因的biotype保存到文件write.csv(feature_info[,c(2,1)],"gene_biotype.csv",quote = F,row.names = F)主讲老师解读：接下来我们识别lncRNA、mRNA和miRNA依据是http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html#查看gene biotype的类型unique(feature_info$gene_biotype)#[1] "protein_coding" "lncRNA"#[3] "transcribed_unprocessed_pseudogene" "transcribed_processed_pseudogene"#[5] "transcribed_unitary_pseudogene" "processed_pseudogene"#[7] "unprocessed_pseudogene" "TEC"#[9] "polymorphic_pseudogene" "IG_C_gene"#[11] "IG_V_gene" "IG_V_pseudogene"#[13] "IG_J_gene" "miRNA"#[15] "translated_unprocessed_pseudogene" "snRNA"#[17] "ribozyme" "scaRNA"#[19] "snoRNA" "TR_C_gene"#[21] "TR_V_gene" "unitary_pseudogene"#此处定义protein_coding作为mRNAmRNA <-"protein_coding"#根据实际研究目的，调整定义为lncRNA的gene_biotype，此处根据Vega定义如下8种biotype为lncRNAlncRNA <-paste("non_coding","3prime_overlapping_ncRNA","antisense","lincRNA","sense_intronic","sense_overlapping","macro_lncRNA","bidirectional_promoter_lncRNA",sep ="|")#还可以定义miRNAmiRNA <-"miRNA"#下面就是表达矩阵里的mRNA、lncRNA、miRNA及其数量，把每类基因的基因名保存到相应的文件里。
mRNA.list<-feature_info[grepl(mRNA, feature_info$gene_biotype),]write.table(mRNA.list,"mRNA.list.txt",quote = F,row.names = F, col.names = F)nrow(mRNA.list)主讲老师解读：最后我们会获得纯粹的mRNA类型与Gene Symbol的对应关系#其它类型的RNA操作类似lncRNA.list<-feature_info[grepl(lncRNA, feature_info$gene_biotype),]write.table(lncRNA.list,"lncRNA.list.txt",quote = F,row.names = F, col.names =F)nrow(lncRNA.list)miRNA.list<-feature_info[grepl(miRNA, feature_info$gene_biotype),]write.table(miRNA.list,"miRNA.list.txt",quote = F,row.names = F, col.names = F)nrow(miRNA.list)主讲老师解读：至此我们就成功获得了mRNA、miRNA、lncRNA对应的Gene Symbol接下来我们需要针对基因名来提取表达矩阵exprdf_uniq<-read.csv("easy_input.csv",header = T,row.names = 1)head(exprdf_uniq)#以lncRNA基因列表为例lncRNA.list<-read.table("lncRNA.list.txt")head(lncRNA.list)mRNA_expr <- exprdf_uniq[as.character(mRNA.list[,2]),]write.csv(mRNA_expr,"mRNA_expr.csv",quote = F,row.names = T)head(mRNA_expr) 某学习平台丨lncRNA_expr <- exprdf_uniq[as.character(lncRNA.list[,2]),]write.csv(lncRNA_expr,"lncRNA_expr.csv",quote = F,row.names = T)head(lncRNA_expr)miRNA_expr <- exprdf_uniq[as.character(miRNA.list[,2]),]write.csv(miRNA_expr,"miRNA_expr.csv",quote = F,row.names = T)head(miRNA_expr)主讲老师解读：至此我们最后会获得mRNA、miRNA和lncRNA单独的表达矩阵这个操作可以分为三个步骤，下载、转换、生成每一个步骤都可以单独使用转换那里不光这三种RNA类型，支持的转换还有很多，可以通过 listAttributes 函数查看参考教程：实例讲述bioMart包的用法 - 简书 (jianshu.com)ID转换一文通晓 (360doc.com)