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

1. 整体规划

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

2. 读入同一个文件夹

2.1 写在最前面的思维模式

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

2.2 读入一个文件夹

# multipling the multi-data into one file
dir.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 id
metadata <- 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 data
test <- 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_id
gene_id <- data.table::fread(paste0('data_in_one/',naid_df$filename[1]))$V1
expr_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 files
load("gtf_df.Rdata")

3.1 筛选需要的类型

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

3.1.1 筛选mRNA

# obtaining mRNA
library(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)
# lncRNA
ncRNA <- 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_19
sample <- ifelse(substring(metadata$TCGA_id, 14,15) == "01", "cancer", "normal")
metadata$sample <- sample
metadata$sample <- as.factor(sample)
metadata <- metadata[order(metadata$sample,decreasing = T),]

mycounts <- mRNA_exprSet
dds <- 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 analysis
res <- 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)

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

整理从TCGA下载的数据