生信分析、meta分析、数据挖掘
TCGA、GEO、SEER、Oncomine

伸出我的小脚,将TCGA轻轻绊倒,然后叉腰哈哈笑

之前我们有一篇帖子讲解TCGA数据的下载 TCGA数据下载方式小结

可以先看一下,前面的部分,一直看到GDC(Genomic Data Commons)把数据下载下来就可以了,我们下面就开始处理数据

假设我们已经把数据下载到了这个文件夹


  1. data_download_from_gdc

打开后发现是1208个文件夹,为什么是1208个,这跟我当时学的数据有关系,可以不用管让 每个文件夹里面还有个压缩文件,我们现在的任务就是,如何把每个文件夹里面的压缩文件放置在统一的文件夹中 假设我们的文件夹是这个,


  1. 00_data_read_in_one_file

并且他就是创建在 data_download_from_gdc中,所以现在文件夹的总体数目是应该是1208加1就是1209个 我们在R语言中输入

  1. length(list.files())

这个命令中list.files()的意思就是列举当前工作目录中文件,length的意思就是有有多少个,一看,果真是1209个 其实

  1. length(dir())

这个命令也能做这个事情。 为什么我的文件夹要起一个这么诡异的名字? 以00开头??

这是因为我需要让他处于文件中第一个位置,我们来验证一下:

  1. dir()[1]

发现第一个文件夹就是他,这样我就每次循环访问后面的1208个文件夹,每次都把看到的东西复制过来就可以了

你说我为什么这么事无巨细,因为我可以很简洁,也可以很啰嗦,取决于我是否想要别人听的懂。

在R语言里面是可以直接创建文件夹的,鼠标右击也可以

  1. dir.create("00_data_read_in_one_file") #创建新的文件夹,确保文件夹排在第一位

遍历和复制,为什么从2开始,因为第一个文件夹你已经知道了

  1. for (dirname in dir()[2:length(dir())]){  

  2.  file <- list.files(dirname,pattern = "*.counts")  #找到对应文件夹中的内容,pattern可以是正则表达式

  3.  file.copy(paste0(dirname,"/",file),"00_data_read_in_one_file")  #复制内容到新的文件夹

  4. }

运行完了之后,我们可以打开文件夹看一下,确实在里面,这时候可以全部选择,

将文件解压到data_unzip文件夹,解压数据1.42GB,

我们发现即使这个样子,文件的名称也是怪怪的

那我们就来转换,转换的信息藏在metadata文件中,这个要去看开始的那个帖子下载

切换一下工作目录

  1. setwd("~/skill_practice/BRCA_999")

注意一开始我就建立了 BRCA_999这个文件夹, data_download_from_gdc是他的子文件夹, 00_data_read_in_one_file又是 data_download_from_gdc的子文件夹 metadata是json格式的

读入json格式的文件,他是一个1208行,15列的数据框

  1. metadata <- jsonlite::fromJSON("metadata.cart.2017-11-15T09_56_59.722935.json")

我们转换的信息就是两列filename和associatedentities,我们把它选出来

  1. require(dplyr)

  2. metadata_id <- metadata %>%

  3.  dplyr::select(c(file_name,associated_entities))

我需要的是 file_name和样本名称,样本名称藏在了 associated_entities 列表中里面包括了 entity_idcase_identity_submitter_identity_type这四个项目,查看第一个了解一下

metadata$associatedentities[1] [[1]] entityid caseid 1 52033f64-1e6f-4657-a4fb-7cfeffc61951 39de7761-e762-4811-b95c-8216b79ae06b entitysubmitterid entitytype 1 TCGA-AN-A0XW-01A-11R-A109-07 aliquot

实际上这边如果能有一个更强力一点json阅读工具,可能结果还要直观一点

现在的想法是我把filename和 associated_entities中的 entity_submitter_id提取出来,做成一个数据框,然后我批量对应转换

  1. naid_df <- data.frame()

  2. for (i in 1:1208){

  3.  naid_df[i,1] <- substr(metadata_id$file_name[i],1,nchar(metadata_id$file_name[i])-3)

  4.  naid_df[i,2] <- metadata_id$associated_entities[i][[1]]$entity_submitter_id

  5. }

现在把1208个小文件读入一个矩阵文件,并且给每一个文件加上 filename和 entity_submitter_id 论坛有道题目就是处理的这个问题

生信编程直播第四题:多个同样的行列式文件合并起来

http://www.biotrainee.com:8080/thread-603-1-1.html

我自己也给出了R语言的解法

在R语言中将多个同样的行列式文件合并起来

http://guoshipeng.com/2017/11/10/14-R-for-binding-colums/

但是当时不知道,TCGA的单个文件是没有列名的,导致无法合并,所以本次要复杂一点

  1. #读入所有解压的文件 1208个

  2. nameList <- list.files("data_unzip/")

  3. location <- which(naid_df==nameList[1],arr.ind = TRUE) ##which函数有一个已知value返回坐标的功能

  4. TCGA_id <- as.character(naid_df[location[1],2]) ##通过坐标,获取TCGA_id

  5. expr_df<- read.table(paste0("data_unzip/",nameList[1]),stringsAsFactors = F, header = F) #读入第一个文件,保存为data.frame

  6. names(expr_df) <- c("gene_id",TCGA_id) #给刚才数据库命名

这边开始批量作业

  1. for (i in 2:length(nameList)){

  2.  location <- which(naid_df==nameList[i],arr.ind = TRUE)

  3.  TCGA_id <- as.character(naid_df[location[1],2])

  4.  dfnew <- read.table(paste0("data_unzip/",nameList[i]),stringsAsFactors = F,header = F)

  5.  names(dfnew) <- c("gene_id",TCGA_id)

  6.  expr_df <- inner_join(expr_df,dfnew,by="gene_id")

  7. }

晚上走的时候没运行完,早上来的时候已经完毕,限速环节应该是read.table,早上再来尝试运行一次总是说内存不够 我尝试了一下fread来解决这个问题:

  1. require(data.table)

  2. nameList <- list.files("data_unzip/")

  3. location <- which(naid_df==nameList[1],arr.ind = TRUE)

  4. TCGA_id <- as.character(naid_df[location[1],2])

  5. expr_df<- fread(paste0("data_unzip/",nameList[1]))

  6. names(expr_df) <- c("gene_id",TCGA_id)

  7. for (i in 2:length(nameList)){

  8.  location <- which(naid_df==nameList[i],arr.ind = TRUE)

  9.  TCGA_id <- as.character(naid_df[location[1],2])

  10.  dfnew <- fread(paste0("data_unzip/",nameList[i]))

  11.  names(dfnew) <- c("gene_id",TCGA_id)

  12.  expr_df <- inner_join(expr_df,dfnew,by="gene_id")

  13. }

结果大概2分钟搞定,速度喜人!!!! 总共60488行,查看最后几行发现有5行不是我们要的

  1. tail(expr_df$gene_id,10)

去掉最后五行

  1. expr_df <- expr_df[1🙁length(expr_df$gene_id)-5),]

保存数据,大概是3个G左右

  1. save(expr_df,file = "expr_df.Rda")


下面开始id转换,信息在GTF文件中表达矩阵里面的gene_id有小数点, 而GTF文件中没有,调整一下,先以“.”分列,在去掉小数点后的列

  1. require(dplyr)

  2. require(tidyr)

  3. expr_df_nopoint <- expr_df %>%

  4.  tidyr::separate(gene_id,into = c("gene_id","drop"),sep="\.") %>%

  5.  dplyr::select(-drop)

  6. save(expr_df_nopoint,file = "expr_df_nopoint.Rda")

  7. load(file = "expr_df_nopoint.Rda")

下载GTF文件来注释 ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens

安装包:

  1. source("https://bioconductor.org/biocLite.R")

  2. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

  3. biocLite("rtracklayer")

  4. biocLite("SummarizedExperiment")

读入GTF数据

  1. gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.90.chr.gtf')

  2. gtf_df <- as.data.frame(gtf1)

保存数据

  1. save(gtf_df,file = "gtf_df.Rda")

读入27个变量,2612129个观测,测试一下显示的不错

  1. test <- gtf_df[1:5,]

  2. View(test)


进展的很快,我们现在可以提取mRNA的表达矩阵啦, 以gtf文件中的 gene_biotype为标准,里面写 protein_coding的就是编码基因

首先要把这些基因提取出来,然后与表达谱融合,我在这个例子还提取了 gene_namegene_id,所以最后的时候,我把三种表达方式合在了一起 这样,以后我无论用什么方式都可以选出我要的基因了

  1. require(dplyr)

  2. require(tidyr)

  3. mRNA_exprSet <- gtf_df %>%

  4.  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #筛选gene,和编码指标

  5.  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%

  6.  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%

  7.  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")

得到19688行,跟我们的认知很吻合

保存数据

  1. save(mRNA_exprSet,file = "mRNA_exprSet.Rda")

下面我们提取非编码基因的是表达矩阵 这里面有个trick,就是编码的gene虽然是用gene_biotype来锁定的 但是,非编码RNA不能这样,应该用转录本来确定,

一个编码基因也有可能转录出非编码基因的,对么?

首先我定义了一个非编码RNA的集合,这个每个人的标准不一样,但是我的原则是,多多益善,这样出来以后会有个问题,就是编码基因转录出非编码基因会无法从基因名称上区分,可以在运行时把geneid换成转录本id,必须要记在心里。

  1. ncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA","processed_transcript","sense_intronic","bidirectional_promoter_lncRNA","non_coding")

这个太长了判断名称有没有写错

  1. ncRNA %in% unique(gtf_df$transcript_biotype)

这时候就开始运行了

  1. LncRNA_exprSet <- gtf_df %>%

  2.  dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA) %>% #注意这里是transcript_biotype

  3.  dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>%

  4.  dplyr::distinct() %>% #删除多余行????

  5.  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%

  6.  tidyr::unite(gene_id,gene_name,gene_id,transcript_biotype,sep = " | ")

最终得到19910行,也还行吧


现在我们有了两个表达矩阵,就可以做差异分析了,但是如果想要做WGCNA或者GSEA分析,需要对数据进行标准化

下面进行标准化

  1. suppressPackageStartupMessages(library(DESeq2))

  2. suppressPackageStartupMessages(library(edgeR))

本次使用的方法是DESeq2中vst法 至于有哪几种方法,还有该选什么方法,明天会有一个帖子介绍

制作metadata,区别肿瘤和正常组,这时候需要对TCGA的id有个了解,从左往右数数,第14,15上的数字很重要 其中01-09是tumor,也就是癌症样本;其中10-29是normal,也就是癌旁

  1. metadata <- data.frame(names(mRNA_exprSet)[-1])

  2. for (i in 1:length(metadata[,1])) {

  3.  num <- as.numeric(substring(metadata[i,1],14,15))

  4.  if (num %in% seq(1,9)) {metadata[i,2] <- "T"}

  5.  if (num %in% seq(10,29)) {metadata[i,2] <- "N"}

  6. }

加个名称

  1. names(metadata) <- c("TCGA_id","sample")

将sample转化成factor

  1. metadata$sample <- as.factor(metadata$sample)

我们可总结一下,有112例normal,1096例肿瘤组织

  1. metadata %>% dplyr::group_by(sample) %>% summarise(n())

保存数据,我想问一下,为什么你老是要保存数据,因为我眼里常含泪花

  1. save(metadata,file = "metadata.Rda")

制作countData

  1. mycounts <- mRNA_exprSet

  2. keepGene=rowSums(edgeR::cpm(mycounts[-1])>0) >=2

  3. table(keepGene);dim(mycounts)

这是数据维度19668 和1209

  1. dim(mycounts[keepGene,])

剩余19546 1209

  1. mycounts <-mycounts[keepGene,]

构建 DESeq对象

  1. dds <-DataSetFromMatrix(countData=mycounts,

  2.                              colData=metadata,

  3.                              design=~sample,

  4.                              tidy=TRUE)

标准化vst法,5分钟

  1. vsd <- vst(dds, blind = FALSE)

  2. mRNA_exprSet_vst <- as.data.frame(assay(vsd))

  3. save(mRNA_exprSet_vst,file = "mRNA_exprSet_vst.Rda")

再把ncRNA的表达谱标准化一下,再自行下载microRNA的数据,就可以构建ceRNA的网络了

一不小心到了这里,我想DEseq的对象都已经构建好了,为什么不进行差异分析呢?

开始计算,这是最耗时间的部分,因为DESeq会自己标准化,不是采用刚才所说的vst法,如果需要知道DESeq标准化的原理 可以看一下这个视频

你一看youtube这几个字,你就知道该怎么做了 直接用刚才构建好的对象dds,你有没有发现,我刚才的代码里面没有保存dds,所以不要嫌我啰嗦,我是保存了的,只是没说,我加载进来先

  1. load(file = "dds.Rda")

下面就是完整的过程

  1.  dds <- DESeq(dds)

  2.  save(dds,file = "dds_DEseq.Rda")

  3.  res <- results(dds, tidy=TRUE) #获取结果

  4.  res <- as_tibble(res)

  5.  require(dplyr)

  6.  res <- res %>%

  7.    separate(row,c("symbol","ensemble","genetype"),sep = " \| ") %>%

  8.    dplyr::select(- c(ensemble,genetype)) %>%

  9.    arrange(desc(abs(log2FoldChange))) %>% #排序。为了去重

  10.    distinct(symbol,.keep_all = TRUE) %>%

  11.    arrange(desc(log2FoldChange))#再次按照log2FoldChange从大到小排序

  12.  save(res,file = "result.Rda")

既然做了差异分析,那对下游基因注释就是顺手 的事情啦,我们没有选择,就用Y叔的包,clusterprofiler 在这之前我们需要获得差异基因列表 这一步,就是基本数据框的基本操作,按照pvalue和fc来排序,选择自己一定数量的基因,数量你来定,最终得到基因列表gene(我没有演示,需要自己做)

加载包,加载之前确定自己装了

  1. library(DOSE)

  2. library(GO.db)

  3. library(org.Hs.eg.db)

  4. library(topGO)

  5. library(GSEABase)

  6. library(clusterProfiler)

基因名称转换,返回的是数据框

  1. gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")

  2. head(gene)

GO分析 细胞组分

  1. ego_CC <- enrichGO(gene = gene$ENTREZID,

  2.                  OrgDb= org.Hs.eg.db,

  3.                  ont = "CC",

  4.                  pAdjustMethod = "BH",

  5.                  minGSSize = 1,

  6.                  pvalueCutoff = 0.01,

  7.                  qvalueCutoff = 0.01,

  8.                  readable = TRUE)

生物过程

  1. ego_BP <- enrichGO(gene = gene$ENTREZID,

  2.                  OrgDb= org.Hs.eg.db,

  3.                  ont = "BP",

  4.                  pAdjustMethod = "BH",

  5.                  minGSSize = 1,

  6.                  pvalueCutoff = 0.01,

  7.                  qvalueCutoff = 0.01,

  8.                  readable = TRUE)

分子功能:

  1. ego_MF <- enrichGO(gene = gene$ENTREZID,

  2.                  OrgDb= org.Hs.eg.db,

  3.                  ont = "MF",

  4.                  pAdjustMethod = "BH",

  5.                  minGSSize = 1,

  6.                  pvalueCutoff = 0.01,

  7.                  qvalueCutoff = 0.01,

  8.                  readable = TRUE)

作图

  1. barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")#条状图,按p从小到大排的

  2. dotplot(ego_BP,title="EnrichmentGO_BP_dot")#点图,按富集的数从大到小的

KEGG分析

  1. kk <- enrichKEGG(gene = gene$ENTREZID,

  2.                organism ="human",

  3.                pvalueCutoff = 0.01,

  4.                qvalueCutoff = 0.01,

  5.                minGSSize = 1,

  6.                #readable = TRUE ,

  7.                use_internal_data =FALSE)

作图

  1. barplot(kk)

  2. dotplot(kk)

昨天有朋友提出ncRNA的生存分析怎么做呢?今天ncRNA的表达矩阵也有了,也就不难了吧,留给大家探索吧。

我一不小心把明天的东西都写进来了,这可怎么办呢 好了,现在TCGA的数据下载,处理,差异表达,分析都已经会搞了,那么下一步就可以进行ceRNA网络构建,WGCNA分析,GSEA分析啦。


赞(0) 打赏
未经允许不得转载:医学SCI科研之家 » 伸出我的小脚,将TCGA轻轻绊倒,然后叉腰哈哈笑
分享到: 更多 (0)

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

meta分析、生信分析

meta、生信交流群综合科研交流群