伸出我的小脚，将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_id, case_id, entity_submitter_id, entity_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_name, gene_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标准化的原理 可以看一下这个视频

https://www.youtube.com/watch?v=UFB993xufUU

你一看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)

相关推荐