转录组差异分析代码（4）R包TCGAbiolinks下载TCGA数据

生信菜鸟团 · 公众号 · 生物 · 2024-09-23 18:28

正文

学习笔记总结于『生信技能树』马拉松课程

完整代码脚本已分享于公众号内

本文继续学习下载TCGA数据，除了转录组差异分析代码（1）数据读取与整理一文介绍的基础方法与代码，还有更为推荐的代码方法：使用R包TCGAbiolinks

一、如果网络没问题

这段代码感兴趣自行研究，但不能因为网络问题下载不了数据而不分析了。所以详细介绍网络有问题该怎么办，请看下文的第二点

# 1.查看TCGA的33个project
rm(list = ls())
library(TCGAbiolinks)
library(stringr)
library(SummarizedExperiment)

#列出TCGA里面的33中癌症，有可能因为网络原因导致这一步报错
projs <- getGDCprojects()$project_id %>% 
  str_subset("TCGA")
projs

# 2.下载并整理表达矩阵
# 如果网络有问题，连上面的projs都得不到，这步几乎也下载不了
# 如果网络可以，下载时间大概3-5min
proj = "TCGA-CHOL" #不要手敲，从上面projs返回的结果中直接复制
f1 = paste0(proj,"expf.Rdata")
if(!file.exists(f1)){
  query = GDCquery(project = proj, 
                 data.category = "Transcriptome Profiling",
                 data.type = "Gene Expression Quantification", 
                 workflow.type = "STAR - Counts",
)
  GDCdownload(query)
  dat = GDCprepare(query)
  exp = assay(dat)
  #tpm = assay(dat,4)#后续生存分析、作图时可以用tpm
  save(exp,file = f1)
}
load(f1)

二、如果网络有问题

1.查看TCGA的33个project

rm(list = ls())
library(TCGAbiolinks)
library(stringr)
library(SummarizedExperiment)
projs <- getGDCprojects()$project_id %>% 
  str_subset("TCGA")
projs

2.下载并整理表达矩阵

人美心善的小洁老师为我们下载好了数据，只要找到我们要的癌症的count和临床信息数据，下载下来放在工作目录下即可

https://share.weiyun.com/ZMQdPBLC 密码：xjlshh

proj = "TCGA-CHOL"
load("chol_exp.Rdata")
exp = chol

假设我们不是研究CHOL，而是KIRC，需要把全篇代码中的CHOL换成KIRC，chol换成kirc

3.下载并整理临床信息

load("chol_clinical.Rdata")
clinical = chol_clinical

4.表达矩阵行名ID转换

library(tinyarray)
exp = trans_exp_new(exp)
exp[1:4,1:4]

5.基因过滤

#需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一
#过滤之前基因数量：
nrow(exp)

#### 常用过滤标准1：
# 仅去除在所有样本里表达量都为零的基因
exp1 = exp[rowSums(exp)>0,]
nrow(exp1)

#### 常用过滤标准2(推荐)：
#仅保留在一半以上样本里表达的基因
exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
nrow(exp)