今天是生信星球陪你的第614天
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
接着差异分析往下做生存分析和富集分析~
当然，生存分析只是肿瘤数据可以做，需要对应的病人生存信息。不是肿瘤的数据就直接跳过生存分析。目录等下次发咯~
大家有没有发现，蓝色排版的那个豆豆消失好几天了，请大家一起到文末投票，召唤豆豆。
4.将差异甲基化位点拿来做批量生存分析
rm(list=ls())
library(data.table)
library(stringr)
library(survival)
library(survminer)
load('./Rdata/step2_filtered_pd_myNorm.Rdata')
load('./Rdata/step3.df_DMP.Rdata')
cg<-rownames(df_DMP[df_DMP$change!='NOT',])
myNorm_tumor<-myNorm[cg,]
logrank test批量生存分析
这里和TCGA系列里的内容是一样的，详见：两种方法批量做生存分析
suv_dat<-data.table::fread('./raw_data/TCGA-HNSC.survival.tsv.gz')
suv_dat$sample=str_sub(suv_dat$sample,1,15)
suv_dat<-suv_dat[suv_dat$sample%in%colnames(myNorm_tumor),]
suv_dat<-suv_dat[str_sub(suv_dat$sample,14,15)=='01',]
suv_dat=merge(pd,suv_dat,by.x='sampleID',by.y='sample')
myNorm_tumor<-myNorm_tumor[,suv_dat$sample]
identical(colnames(myNorm_tumor),suv_dat$sample)
#>[1]TRUE
library(survival)
logrankP<-apply(myNorm_tumor,1,function(x){
#x<-myNorm_tumor[1,]
suv_dat$group<-ifelse(x>mean(x),'High','Low')
res<-coxph(Surv(OS.time,OS)~group,data=suv_dat)
beta<-coef(res)
se<-sqrt(diag(vcov(res)))
p.val<-1-pchisq((beta/se)^2,1)
})
table(logrankP<0.05)#110个CpG位点
#>
#>FALSETRUE
#>1676110
得到110个对生存影响显着的差异甲基化位点，取前20个画图。
surv_gene<-names(sort(logrankP))[1:20]
choose_matrix<-myNorm[surv_gene,]
annotation_col<-data.frame(Sample=pd$group_list)
rownames(annotation_col)<-colnames(choose_matrix)
ann_colors=list(Sample=c(Normal='#4DAF4A',Tumor='#E41A1C'))
library(pheatmap)
pheatmap(choose_matrix,show_colnames=T,
annotation_col=annotation_col,
border_color=NA,
color=colorRampPalette(colors=c('white','navy'))(50),
annotation_colors=ann_colors)
image.png也可以画画生存分析的图
gs=head(surv_gene,4)
exprSet=myNorm_tumor
meta=suv_dat
splots<-lapply(gs,function(g){
meta$gene=ifelse(exprSet[g,]>median(exprSet[g,]),'high','low')
sfit1=survfit(Surv(OS.time,OS)~gene,data=meta)
ggsurvplot(sfit1,pval=TRUE,data=meta)
})
arrange_ggsurvplots(splots,print=TRUE,
ncol=2,nrow=2)
5.富集分析
利用ChAMP包对过滤后的数据做了差异甲基化位点分析。如果是肿瘤数据的话，可以加一步生存分析。
rm(list=ls())
load(file='Rdata/step3.df_DMP.Rdata')
library(ggplot2)
library(stringr)
library(clusterProfiler)
library(org.Hs.eg.db)
ID转换
length(unique(df_DMP$gene))
#>[1]16338
s2e<-bitr(unique(df_DMP$gene),fromType='SYMBOL',
toType=c('ENTREZID'),
OrgDb=org.Hs.eg.db)
df_DMP=merge(df_DMP,s2e,by.y='SYMBOL',by.x='gene')
table(!duplicated(df_DMP$ENTREZID))
#>
#>FALSETRUE
#>8377814188
gene_up=unique(df_DMP[df_DMP$change=='UP','ENTREZID'])
gene_down=unique(df_DMP[df_DMP$change=='DOWN','ENTREZID'])
gene_diff=c(gene_up,gene_down)
gene_all=unique(df_DMP$ENTREZID)
富集分析及其可视化
之前完整的富集分析会把上下调基因分开、合并都做。GO富集分析会将MF、CC、BP三个部分分开做。代码太多了，在这简化一下。
kkgo_file='./Rdata/kkgo_file.Rdata'
if(!file.exists(kkgo_file)){
kk<-enrichKEGG(gene=gene_diff,
universe=gene_all,
organism='hsa',
pvalueCutoff=0.05)
go<-enrichGO(gene_diff,OrgDb='org.Hs.eg.db',ont='all')
save(kk,go,file=kkgo_file)
}
load(kkgo_file)
dotplot(kk)
image.pngbarplot(go,split='ONTOLOGY',font.size=10)+
facet_grid(ONTOLOGY~.,scale='free')+
scale_x_discrete(labels=function(x)str_wrap(x,width=45))
