TCGA数据下载与整理-TCGAbiolinks包

写在前面，如果不是R语言高手，建议利用网站下载
个人爱好——1、官网：https://portal.gdc.cancer.gov
方法见从TCGA数据库下载并整合清洗高通量肿瘤表达谱-临床性状数据。
2、UCSC—xena:UCSC Xena,强推，压缩版，为tsv文件。

1、下载临床信息

以下代码参考单基因生信分析流程（1）一文解决TCGA数据下载整理问题 - 简书

1.1、加载包，rm(list = ls())清空环境，日常操作，以免出什么幺蛾子，好习惯得养成
rm(list = ls())
setwd("D:\\生\\结")
library(TCGAbiolinks)
library(SummarizedExperiment)
library(dplyr)
library(DT)
library(rtracklayer)
library(readr)

我遇到的问题：记得先安装，如果install.packages()安装出现安装不成功的情况，建议用官网的安装代码，或者BiocManager::install()，可以解决。或者更改镜像路径！

下载整理TCGA数据，

以结肠癌（COAD）为例，下载其他类型，请参考TCGA / 癌症简称 / 缩写 / TCGA癌症中英文对照 - 简书

query=GDCquery(project = "TCGA-COAD",
              data.category = "Clinical",#数据类型
              file.type = "xml")
GDCdownload(query)
clinical=GDCprepare_clinic(query,clinical.info = "patient")
1.2、筛选主要信息

由于我本人用R下载数据的时候经常有问题，所以我在xena下载，为tsv文件，以下均以该文件为例。
首先筛选自己所需要的数据

dat=read_tsv("TCGA-COAD.GDC_phenotype.tsv",col_names = TRUE,col_types = NULL)
clinicaldata=dat%>%
  dplyr::select(submitter_id.samples,age_at_initial_pathologic_diagnosis,gender.demographic,submitter_id,
                days_to_death.demographic,days_to_last_follow_up.diagnoses,race.demographic,
                age_at_diagnosis.diagnoses,person_neoplasm_cancer_status,tumor_stage.diagnoses,
                pathologic_T,pathologic_N,pathologic_M,vital_status.demographic,disease_type,
                sample_type.samples)%>%
  distinct(submitter_id.samples,.keep_all = TRUE)#选择所有行

利用dplyr::select选择所需要列，建议先view或者$命令查看下载的列名根据需要选择
distinct选择所需要的行， %>%为管道符号，表示然后意思，能够加快对代码的阅读理解，不加载dplyr包无法使用，偶尔会报错。

1.3、整理死亡患者信息
dead_patient=clinicaldata%>%
  dplyr::filter(vital_status.demographic=='Dead')%>%
  dplyr::select(-days_to_last_follow_up.diagnoses)%>%
  reshape::rename(c(submitter_id='Barcode',
                    submitter_id.samples='sampleID',
                    gender.demographic='Gender',
                    days_to_death.demographic='OS_time',
                    vital_status.demographic='OS',
                    race.demographic='Race',
                    age_at_initial_pathologic_diagnosis='Age',
                    person_neoplasm_cancer_status='cancer_status',
                    tumor_stage.diagnoses='TNM_stage'))%>%
  mutate(OS=ifelse(OS=="Dead",1,0))

备注
filter过滤函数；reshape::rename重新命名函数
此处应注意选择的列，选择死亡患者，则不要days_to_last_followup这一列，然后days_to_death作为生存时间，挑选生存患者的时候则相反，切记切记！
mutate(),转换数值函数

1.4、整理生存患者信息，合并
alive_patient=clinicaldata%>%
  dplyr::filter(vital_status.demographic=='Alive')%>%
  dplyr::select(-days_to_death.demographic)%>%
  reshape::rename(c(submitter_id='Barcode',
                    submitter_id.samples='sampleID',
                    gender.demographic='Gender',
                    days_to_last_follow_up.diagnoses='OS_time',
                    vital_status.demographic='OS',
                    race.demographic='Race',
                    age_at_initial_pathologic_diagnosis='Age',
                    person_neoplasm_cancer_status='cancer_status',
                    tumor_stage.diagnoses='TNM_stage'))%>%
  mutate(OS=ifelse(OS=="Dead",1,0))
surdata=rbind(dead_patient,alive_patient)
write.csv(surdata,file = "COAD_phenotype.csv")

二、下载表达谱与整理

COAD为例

exp=GDCquery(project = "TCGA-COAD",           
             data.category = "Transcriptome Profiling",
             data.type = "Gene Expression Quantification", 
             workflow.type = "HTSeq - FPKM")
GDCdownload(exp)
expdat=GDCprepare(query = exp)
count_matrix=as.data.frame(assay(expdat))

运行时家里网络报：Error in value[3L] :**** GDC server down, try to use this package later，实验室网络也偶尔遇到。主要是官网api问题，所以建议在线下载安全，个人愚见。以下操作以xena下载数据展示，tsv文件。

library(TCGAbiolinks)
library(SummarizedExperiment)
library(dplyr)
library(DT)
library(readr)
library(rtracklayer)
exp=read_tsv("TCGA-COAD.htseq_fpkm.tsv",col_names = T,col_types = NULL,na=c("","NA"))
exp=exp%>%
  reshape::rename(c(Ensembl_ID="gene_id"))
#按装rtracklayer包，不仅可用于gtf文件的导入导出，还可以操作其他多种数据格式(BED，WIG，bigWig......)
#BiocManager::install("rtracklayer")
#加载注释文件，仅仅所需要得id,name,type列，其他信息可不需要，感兴趣得可以view函数看一下到底还有啥
gtf=rtracklayer::import('gencode.v35.long_noncoding_RNAs.gff3')
gff_df=as.data.frame(gtf)
genegtf=gff_df%>%
  dplyr::select(gene_id,gene_type,gene_name)
#提取表达谱内的lncRNA,记得先把表达谱内的Ensembl_ID改成"gene_id")
lncRNA_exp=gff_df%>%
  dplyr::filter(type=="gene",gene_type=="lncRNA")%>%
  dplyr::select(gene_id,gene_type,gene_name)%>%
  dplyr::inner_join(exp,by="gene_id")%>%
  tidyr::unite(gene_id,gene_type,gene_id,gene_name,sep="|")
#做到这一步时，记得留意一下目前第一列得名称得顺序，接下来是把他们分开，顺序需要跟后面删除得对应起来。我把name放在第三列，所有后面删除第一第二列。
lncRNA_exp=lncRNA_exp%>%
  tidyr::separate(gene_id,c("gene_type","gene_id","gene_name"),
                  sep="\\|")
lncRNA_exp=lncRNA_exp[,-(1:2)]#删除第一第二列，保留gene_name列
index=duplicated(lncRNA_exp$gene_name)
lncRNA_exp=lncRNA_exp[!index,]
row.names(lncRNA_exp)=lncRNA_exp$gene_name
lncRNA_exp$gene_name=NULL

提取完毕，可以保存，接下来删除癌旁样本和二次测序的样本,根据TCGA样本命名规则，第14、15个数字表示样本，其中编号0109表示肿瘤，1019表示正常对照，(https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes)。

mete=data.frame(colnames(lncRNA_exp))#取第一行，即样本ID
for (i in 1:length(mete[,1])) {
  num=as.numeric(as.character(substring(mete[i,1],14,15)))#提取14、15列
  if(num==1) {mete[i,2]="T"}
  if(num!=1){mete[i,2]="N"}
}
names(mete)=c("id","group")
mete$group=as.factor(mete$group)
mete=subset(mete,mete$group=="T")
lncRNA_exp1=lncRNA_exp[,which(colnames(lncRNA_exp)%in% mete$id)]
#讲FPKM转为TPM，看自己需求。
fpkmToTpm=function(fpkm)
{
  exp(log(fpkm)-log(sum(fpkm))+log(1e6))
}
expr=as.data.frame(apply(lncRNA_exp1,2,fpkmToTpm))
#数据导出，TXT文档
write.table(expr, file = "COAD_lncRNA_exp_TPM.txt", append = FALSE, 
quote = TRUE, sep = "\t",eol = "\n", na = "NA", dec = ".", 
row.names = TRUE,col.names = TRUE, qmethod = c("escape", "double"),fileEncoding = "")

至此，得到肿瘤样本（01）的表达谱（TPM或者FPKM），其实01-09均是肿瘤样本，本次仅仅删选了01的肿瘤样本。由于xena下载的表达谱数据不包含二代测序数据，所以以下代码仅供参考。

mete1=data.frame(colnames(lncRNA_exp1))
for (i in 1:length(mete1[,1])) {
  chr=as.character(substring(mete1[i,1],22,25))#第22，25列
  if(chr=="A277") {mete1[i,2]="Rep"}
  if(num!="A277"){mete1[,2]="T"}
}
names(mete1)=c("id","group")
mete1$group=as.factor(mete1$group)
mete1=subset(mete1,mete1$group=="T")
lncRNA_exp2=lncRNA_exp1[,which(colnames(lncRNA_exp1) %in% mete1$id)]
write.csv(lncRNA_exp2,file = "COAD_lncRNA_exp.csv")#以csv文件保存
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,542评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,596评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 158,021评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,682评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,792评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,985评论 1赞 291
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,107评论 3赞 410
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,845评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,299评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,612评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,747评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,441评论 4赞 333
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,072评论 3赞 317
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,828评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,069评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,545评论 2赞 362
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,658评论 2赞 350