TCGA任意基因任意肿瘤,随意分析
当我们拿到一个基因，肯定是想知道它的各种基本信息，今天的教程可以做到
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
