如何生成WGCNA分析所需的临床表型信息

根据文章生信菜鸟团：一文学会WGCNA分析所讲：

WGCNA输入数据的准备，主要有两个数据需要提前准备好：
（１）是表达矩阵， 如果是芯片数据，常规的归一化矩阵即可。如果是转录组数据，最好是RPKM值或者其它归一化好的表达量。
（２）临床信息或者其它表型，也就是样本的属性。

在前几篇学习笔记里，我得到了TCGA数据库里的count矩阵（整合gdc-client批量下载的文件），并且把它转化成了FPKM（如何把从TCGA下载的HTseq count转换成FPKM）。现在要做的就是得到临床信息了。

这里你需要进入你之前TCGA下载的页面：

之前，我们点击了Download下载了count矩阵，点击Metadata，下载了与count矩阵相对应的一些信息，比如各种id，各种测序的信息。现在需要另一个东西，就是上图里的“clinical”。点击它，选择“JSON”格式：

根据生信菜鸟团的讲解：

临床信息表型的dataframe需要自己制作，这个是学习WGCNA的基础，本次实例代码都是基于这两个数据。至于如何做出上面代码的两个例子，取决于大家自己的项目。

所以可能每个人所需要的临床信息不一样，比如说肺癌，你可能会需要“smoke”的信息；或者口腔癌，你可能需要“alcohol”的信息。或者你想研究原发和复发的差异，那么你应该关注treatment里的信息。我这里只是举个例子，知道怎么提取信息，就可以为所欲为了。。。

这里我分成３个部分来说明，因为下载完clinical的信息你会发现，临床的样品和测序的count样品数是不一样的，我查了一些文章，说是有个别patient样品重复测序的情况。你需要分析这个样品数的区别到底是什么。

从count矩阵对应的metadata里提取信息（cancer or normal）

对于一个count矩阵来说，可能我们需要关注它的点在于，每一个样品是来自tumor组织还是normal对照。

#载入FPKM矩阵
> a = read.csv("TCGA_HNSCC_count2FPKM.csv",header = TRUE, sep=",")
#养成随时查看你生成的东西的好习惯
> View(a)　
#发现第一列是基因名，把第一列设置成行名
> rownames(a) = a[,1]
#删除第一列　
> a = a[,-1]　

根据metadata.json中的信息，对数据进行分组cancer or normal :

#因为FPKM矩阵的最后一列是基因长度，所以只要前546列样品的信息
> b = a[,1:546]
#group_name: etc: TCGA-BB-4224-01A-01R-1436-07
> group_name=colnames(b) 
#提取第１４，１５位字符，规则请看我前一篇笔记关于TCGA编号的规则
> group_name=substr(group_name,14,15) 
#小于10的是cancer，大于10的是正常或癌旁
> group=ifelse(as.numeric(group_name)<10,1,0) 
> group=factor(group,levels = c(0,1),labels = c('normal','cancer'))
#将分组存为一个data.frame
> group = as.data.frame(group)
> rownames(group) = colnames(b)
> head(group)
> colnames(group)

现在，我的group长这样：

把FPKM矩阵的样品名换成case_id,为了下面和临床信息进行比较：

> library(rjson)
#分别读入两个文件，metadata是和FPKM矩阵对应的，clinical是临床文件
> x = fromJSON(file='metadata.cart.2020-04-17.json')
> trait = fromJSON(file='clinical.cart.2020-04-17.json')
> n = ncol(b)
> case_id=rep(0,n)
> TCGA_id = rep(0,n)
> for(i in 1:n){
  case_id[i]=x[[i]]$associated_entities[[1]]$case_id
  TCGA_id[i]=x[[i]]$associated_entities[[1]]$entity_submitter_id
} #这里我从metadata里提取了case_id和TCGA的标准编号
> sample_matrix = data.frame(case_id = case_id, TCGA_id = TCGA_id)
> head(sample_matrix)

看一下我从metadata里提取的case_id和TCGA编号：

接下来把FPKM矩阵的样品名换成case_id，这样后面可以和临床信息对应上：

> colnames(b)=sample_matrix$case_id
> View(b)

把前面的group和case_id，TCGA_id合并起来：

> sample546_info <- cbind(group,sample_matrix)
> head(sample546_info)

上面就是从count矩阵对应的metadata里提取的信息。

从临床JSON文件里提取需要的信息

上面说了，每个人研究的重点不一样，关注点也不一样，所以一定要根据你的实际情况来选择性的提取，不要盲目copy代码。
比如我研究的cancer，它的原发部位会有不同，那么我就想提取临床信息里的case_id，肿瘤分期，以及肿瘤的起始部位,并把它们存成一个dataframe：

> m = length(trait)
> clinic_case_id=rep(0,m)
> stage=rep(0,m)
> tumor_origin_site = rep(0,m)
> for (i in 1:m){
  clinic_case_id[i] = trait[[i]]$case_id
  stage[i] = trait[[i]]$diagnoses[[1]]$tumor_stage
  tumor_origin_site[i] = trait[[i]]$diagnoses[[1]]$tissue_or_organ_of_origin
}
> clinic_info = data.frame(case_id= clinic_case_id, tumor_stage = stage,tumor_origin = tumor_origin_site)
> head(clinic_info)
> rownames(clinic_info)
> colnames(clinic_info)

看一下从临床JSON文件里提取的信息：

NOTE:我下载的clinical文件只有501个样品，而FPKM样品数是546个。有些文章说可能存在同一个病人的样品重复测序，那么就来看一下是不是：

合并两个dataframe

利用merge函数合并sample546_info和临床信息，这里为什么用merge呢？merge可以合并两个行数不同的dataframe，参数信息见下面：

#通过case_id字段作为连接列，将clinic_info中信息匹配到sample546_info中
#by="case_id"指定连接列（两个数据集均有的列）为case_id字段。
#all.x = TRUE表示以sample546_info数据集为参照，将clinic_info中的信息匹配过来。
> all_info <- merge(sample546_info, clinic_info, by="case_id",all.x = TRUE)
> head(all_info)

看一下合并之后的：

查看是否有重复的case

现在来看看是否有重复的patient测序的信息，查看case_id列的重复项：

# 显示重复项
> all_info[duplicated(all_info$case_id),]

 case_id  group                      TCGA_id  tumor_stage                                       tumor_origin
16  0858c8b7-e2eb-4461-b65e-9d476029ad8d cancer TCGA-CV-7250-01A-11R-2016-07    stage iva                                        Larynx, NOS
18  0871333c-1088-4af1-a365-12256643c5bd cancer TCGA-CV-6935-01A-11R-1915-07    stage iva                                        Larynx, NOS
20  08a45833-16a3-4a57-a310-307ec086e558 normal TCGA-CV-7438-11A-01R-2132-07      stage i                                        Tongue, NOS
63  1dd23cd6-3aa8-4553-8813-04701451846e normal TCGA-CV-7103-11A-01R-2016-07    stage iva                                        Tongue, NOS
78  241d9310-9137-42dd-b28d-0dc50c44cb43 cancer TCGA-CV-7101-01A-11R-2016-07     stage ii                                        Larynx, NOS
..........

看一下有多少个重复项：

> dup = all_info[duplicated(all_info$case_id),]
> nrow(dup)
[1] 45 ＃说明有45个重复的case_id

这里显示有重复的case，那么要删除吗？先来看一下，我选取了上面显示的第一个重复项，去all_info里查看一下区别在哪儿（在Rstudio里View页面ctrl+F搜索）：

从查询的结果可以看到：虽然有两个相同的case_id，测序样品来自同一病人，但是两个测序一个是测的正常/癌旁组织，一个测的是肿瘤组织。这样的情况就不能删除，以防后续我们要进行正常和肿瘤的差异比较。

但是如果你的重复样品确实是同样的组织测了两次，那么就要删除重复的，具体如何操作，请看文章：数据挖掘专题 | TCGA样本命名详解 。

那么我这个FPKM矩阵里有多少个正常/癌旁组织呢？

> sum(all_info$group == "normal")
[1] 44

NOTE:现在我知道FPKM矩阵里的样品数是546个，临床样品数是501个，那么501+44= 545个，还有一个我就不想去追究了（我好懒。。。），感觉１个重复测序对500个样品里应该没太大影响。

现在对于我下载的这个TCGA的数据来讲，有了FPKM矩阵，又有了需要的临床信息，基本上不管是WGCNA还是差异基因的分析应该都可以进行了。