整理从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)