TCGA任意基因任意肿瘤,随意分析

治疗白癜风第一的医院 http://pf.39.net/bdfyy/bdfzj/

当我们拿到一个基因,肯定是想知道它的各种基本信息,今天的教程可以做到

1任意一个肿瘤在泛癌中的表达情况

2任意一个基因在肿瘤和正常组织中的表达情况

任意一个基因在肿瘤不同分期,不同性别等临床特性的表达情况

4任意一个基因在任意一个肿瘤,或肿瘤的某种特征中的生存分析

5这个数据能做到的太多,只要充分发挥想象力,所有的数据获取方式见文末

首先我们从UCSCXena数据框下载pancancer的标准化后的表达谱和临床资料

放到我们的工作目录上,解压

读取数据

这一步耗时较长(如果这个过程你的电脑hold不住了,可以直接用后面整理好的数据,开始作图)

rm(list=ls())library(tidyverse)#读取数据ALLdata-data.table::fread("tcga_RSEM_gene_tpm",data.table=F)ALLdata[1:5,1:5]

读取临床信息

#读取临床信息clin-data.table::fread("Survival_SupplementalTable_S1__xena_sp",data.table=F)

筛选有用的临床资料,并将复杂的名字改成简单的名字

clin1-clin%%select(sample,cancertypeabbreviation,gender,ajcc_pathologic_tumor_stage,OS,OS.time)%%#选取这些列rename(Type=cancertypeabbreviation,Stage=ajcc_pathologic_tumor_stage,status=OS,time=OS.time)

读取基因注释文件(详见TCGA数据下载与ID转换)

gtf-rtracklayer::import(Homo_sapiens.GRCh8.99.chr.gtf.gz)#转换为数据框gtf-as.data.frame(gtf)#保存save(gtf,file="gtf.Rdata")

去掉基因名小数点及后面的数字方便下一步转换

colnames(ALLdata)[1]-gene_idALLdata1-separate(ALLdata,gene_id,into=c("gene_id"),sep="\\.")ALLdata1[1:5,1:5]

提取编码蛋白,与表达谱进行合并以转换基因名。(也可以取miRNA或者lncRNA等)

mRNA-dplyr::filter(gtf,type=="gene",gene_biotype=="protein_coding")%%#选择编码蛋白select(gene_name,gene_id,gene_biotype)%%#选择有用的三列inner_join(ALLdata1,by="gene_id")%%#与表达谱合并select(-gene_id,-gene_biotype)%%distinct(gene_name,.keep_all=T)mRNA[1:5,1:5]

将gene_name这一列作为行名

row.names(mRNA)-mRNA[,1]mRNA-mRNA[,-1]mRNA[1:5,1:5]

这个时候继续下去可能电脑hold不住了,先保存下文件

save(mRNA,file=pancancer_mRNA.Rdata)save(clin1,file=pancancer_clin.Rdata)

重新开始

rm(list=ls())load(pancancer_mRNA.Rdata)

将表达谱倒置,方便后续与临床资料的合并

mRNA-as.data.frame(t(mRNA))mRNA[1:5,1:5]

接下来是设计Normal和Cancer的分组,根据TCGA数据ID的特性可区分Normal和CancerTCGA差异分析及ggplot作图验证

group_list=ifelse(as.numeric(substr(rownames(mRNA),14,15))10,tumor,normal)design-model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))rownames(design)=rownames(mRNA)head(design)

这里我用了之前做差异分析时分组的代码,应该可以写成更简洁的形式

给design数据数据变成数据框,添加一列ID

class(design)design-as.data.frame(design)designID-row.names(design)design[1:,1:]

将design分组中tumor这一列的1换成Tumor,0换成Normal

designtumor[designtumor==1]-"Tumor"#1换成Tumordesigntumor[designtumor==0]-"Normal"#0换成Normaldesign-design[,-1]#去掉normal这一列colnames(design)[1]-Type#cancer这一列列名改为Typehead(design)

合并数据

mRNAID-row.names(mRNA)mRNA=inner_join(mRNA,design,by="ID",copy=T)mRNA[1:5,1:5]

调整列的顺序便于观察

mRNA-select(mRNA,ID,Type,everything())mRNA[1:5,1:5]

加载临床资料

load(pancancer_clin.Rdata)clin1[1:5,1:5]

为了合并,将列名sample换成ID,由于Type列名重复了,将Type换成Cancer

clin1-rename(clin1,ID=sample,Cancer=Type)clin1[1:5,1:5]

合并

drawdata-dplyr::inner_join(mRNA,clin1,by="ID",copy=T)drawdata[1:5,1:5]

调整列名便于观察

drawdata-select(drawdata,ID,Cancer,gender,Stage,status,time,everything())drawdata[1:5,1:10]

保存数据

save(drawdata,file="pancancer_drawdata.Rdata")

重新加载数据(以后就可以直接加载这个数据画图了,但是别忘了加载包)

rm(list=ls())load(pancancer_drawdata.Rdata)drawdata[1:5,1:10]

看看每个肿瘤和正常组织的个数

table(filter(drawdata,Type==Tumor)Cancer)table(filter(drawdata,Type==Normal)Cancer)

接下来就可以任意画图啦!

看一看单基因在泛癌中的表达

p-ggplot(filter(drawdata,Type==Tumor),aes(x=Cancer,y=SPP1,color=Cancer))+geom_boxplot()+geom_jitter()+theme_bw()+theme(legend.position=none)p

加上正常组织

library(ggpubr)p-ggboxplot(drawdata,x="Cancer",y="SPP1",fill="Type",legend=F,palette=c("#00AFBB","#E7B"))+theme_bw()p

单个肿瘤差异表达

p-ggboxplot(filter(drawdata,Cancer==LIHC),x="Type",y="SPP1",fill="Type",legend=F,palette=c("#00AFBB","#E7B"),bxp.errorbar=T)+theme(legend.position=none)+ylab(label=SPP1expression)p

加上P值

my_

转载请注明:http://www.ningbohuodai.com/yyjj/9602.html

网站简介| 发布优势| 服务条款| 隐私保护| 广告合作| 网站地图| 版权申明

当前时间: