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

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

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

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

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

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

根据生信菜鸟团的讲解:

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

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

这里我分成3个部分来说明,因为下载完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) 
#提取第14,15位字符,规则请看我前一篇笔记关于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个,还有一个我就不想去追究了(我好懒。。。),感觉1个重复测序对500个样品里应该没太大影响。

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,547评论 6 477
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,399评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,428评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,599评论 1 274
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,612评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,577评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,941评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,603评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,852评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,605评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,693评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,375评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,955评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,936评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,172评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 43,970评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,414评论 2 342

推荐阅读更多精彩内容