整理从TCGA下载的数据

如果从TCGA官网网页下载数据,或者使用gdc-client工具下载的数据,都是以单个的文件夹形式存储,并且文件夹中的为压缩文件,所以,下载数据后,第一步就是如何把这些文件复制在同一个文件夹中,以利于后续的数据读入。采用R将文件复制并读入同一个文件夹,并制作表达矩阵。

1. 整体规划

对于数据整理,需要解决以下几个问题。
1.将所有的压缩文件存放在一个文件夹中,这些文件不需要解压,因为R可以直接读取这种类型的压缩文件。
2.将这些文件读入R,并与TCGA-id对应
3.与临床信息相对应。

2. 读入同一个文件夹

2.1 写在最前面的思维模式

编程的思维就是把人从重复的劳动中解放出来,如果你知道一件事情如何做,那么和这件事情相似的事情,就可以交给程序。R虽然不是完全的编程语言,但是作为一个统计学语言,最基本的重复、循环这些依然是可以完成的。所以,重点就是,我们要把自己做这件事情的每一步分解细化,告诉计算机需要做什么。

2.2 读入一个文件夹

# multipling the multi-data into one filedir.create('data_in_one')for(dirname in dir('rawdata/')){  file <- list.files(paste0(getwd(),"/rawdata/", dirname), pattern = "*.counts.gz")  file.copy(from = paste0(getwd(),"/rawdata/",dirname,"/",file),to = "data_in_one")}
  新建data_in_one文件夹,用来存放压缩文件。从TCGA下载的文件存放于 rawdata。

2.3 匹配文件名和TCGA-ID

# matching the TCGA idmetadata <- jsonlite::fromJSON("metadata.cart.2020-03-18.json")naid_df <- data.frame()for (i in 1:nrow(metadata)) {  naid_df[i,1] <- metadata$file_name[i]  naid_df[i,2] <- metadata$associated_entities[i][[1]]$entity_submitter_id}colnames(naid_df) <- c("filename", "TCGA_id")

读取JSON文件,使用jsonlite::fromJSON相比较而言,是比较好用的一个函数。

2.4 读入表达数据

# loading the datatest <- data.table::fread(paste0('data_in_one/',naid_df$filename[1]))expr_df <- data.frame(matrix(data = NA, nrow = nrow(test),ncol =  nrow(naid_df)))for (i in 1:nrow(naid_df)) {  print(i)  expr_df[,i] = data.table::fread(paste0("data_in_one/",naid_df$filename[i]))[,2]}colnames(expr_df) <- naid_df$TCGA_idgene_id <- data.table::fread(paste0('data_in_one/',naid_df$filename[1]))$V1expr_df <- cbind(gene_id=gene_id, expr_df)tail(expr_df$gene_id,10)expr_df <- expr_df[1:(nrow(expr_df)-5),]save(expr_df,file = "expr_df.Rdata")

2.5 删去ensemble ID 的点号

expr_df_nopoint <- expr_df %>%   tidyr::separate(gene_id, into = c("gene_id"), sep = "\\.")

之所以删去点,是因为附带点号,不能与gene symbol 对应。

3. 载入gene symbol 与ensemble id对应的文件

# obtaining gtf filesload("gtf_df.Rdata")

3.1 筛选需要的类型

根据TCGA 数据库的数据,可以选取相应的类型,基因表达mRNA,非编码RNA(ncRNA),包括lncRNA。

3.1.1 筛选mRNA

# obtaining mRNAlibrary(dplyr)library(tidyr)mRNA_exprSet <- gtf_df %>%   dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #筛选gene,和编码指标  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%   dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%   tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")

3.1.2 筛选lncRNA

library(dplyr)library(tidyr)# lncRNAncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA",           "processed_transcript","sense_intronic",           "bidirectional_promoter_lncRNA","non_coding",           "antisense_RNA")LncRNA_exprSet <- gtf_df %>%   dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA) %>% #注意这里是transcript_biotype  dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>%   dplyr::distinct() %>% #删除多余行  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>%   tidyr::unite(gene_id,gene_name,gene_id,transcript_biotype,sep = " | ")

3.2 标准化数据

# normalization metadata <- data.frame(TCGA_id = colnames(expr_df)[-1])table(substring(metadata$TCGA_id, 14,15))   ####01- 411, 11-19; cancer_411, normal_19sample <- ifelse(substring(metadata$TCGA_id, 14,15) == "01", "cancer", "normal")metadata$sample <- samplemetadata$sample <- as.factor(sample)metadata <- metadata[order(metadata$sample,decreasing = T),]mycounts <- mRNA_exprSetdds <- DESeqDataSetFromMatrix(countData = mycounts,                              colData = metadata,                              design = ~sample,                              tidy = T)dds <- DESeq(dds)save(dds, file = "mRNA_exprSet_dds_sample.Rdata")vsd <- vst(dds, blind = FALSE)plotPCA(vsd, "sample")mRNA_exprSet_vst <- as.data.frame(assay(vsd))save(mRNA_exprSet_vst,file = "mRNA_exprSet_vst.Rda")

3.3 差异分析

# differential analysisres <- results(dds, tidy=TRUE,contrast = c('sample','cancer','normal'))save(res,file = "DEGs_results.Rda")write.csv(res, file = 'DEGs_results.csv',row.names = res$row)

这一步,需要注意比较那两者进行比较。

(0)

相关推荐