TCGA数据读入及使用DESeq分析

代码主要是向微信公众号 果子学生信， 生信技能树，学习获得，写在简书，主要是方便自己查找代码。

将TCGA数据批量读入

# 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")
}

将文件的编号和文件名称对应

# 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")

metadata 中存放的信息主要如下：

 >colnames(metadata)
 [1] "md5sum"                "experimental_strategy" "data_type"             "data_format"          
 [5] "state"                 "file_id"               "file_name"             "data_category"        
 [9] "access"                "associated_entities"   "submitter_id"          "file_size"            
[13] "analysis"              "annotations"  

把TCGA和临床信息进行合并

# 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")

# omiting the dot of ensemble ID
expr_df_nopoint <- expr_df %>% 
  tidyr::separate(gene_id, into = c("gene_id"), sep = "\\.")

需要注意的是，每次都应当对数据的的基本结构进行检查，比如这里面读入的最后5行就不是表达矩阵的信息。

加载gtf文件，为后续的分子筛选提供靶点。

# obtaining mRNA
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 = " | ")

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

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

ncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA","processed_transcript","sense_intronic","bidirectional_promoter_lncRNA","non_coding")
ncRNA %in% unique(gtf_df$transcript_biotype) #这个太长了判断名称有没有写错

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 = " | ")

到这里我们获得了表达矩阵，包含有mRNA的矩阵和包含非编码lncRNA 的矩阵，是没有经过标准化的数据。

紧接着，我们使用DESeq2 进行标准化处理，具体使用DESeq2中vst法。针对counts数据还可以使用edgeR进行处理

suppressPackageStartupMessages(library(DESeq2))

suppressPackageStartupMessages(library(edgeR))

suppressPackageStartupMessages() 只加载package，不显示相关信息。
制作metadata，区别肿瘤和正常组。需要理解TCGA-id的命名规则，第14,15上的数字很重要 其中01-09是tumor，也就是癌症样本；其中10-29是normal，也就是癌旁。

对分析数据的癌症组，癌旁组进行统计。

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),]

或者通过dplyr package进行统计癌症组和癌旁组。

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

保存数据

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

制作counts数据

mycounts <- mRNA_exprSet

制作DESeq 对象

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")

差异分析

# 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)

来自果子学生信的代码

dds <- DESeq(dds)

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

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

  res <- as_tibble(res)

  require(dplyr)

  res <- res %>%

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

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

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

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

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

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

提取差异基因

table(res$padj < 0.01 & abs(res$log2FoldChange) > 1)
res <- res[order(res$padj),]
#resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized = T)), by="row.names",sort=FALSE)
signres <- na.omit(res[res$padj<0.01 & abs(res$log2FoldChange) > 1,])