看黄神器免app免vip

你的位置:看黄神器免app免vip-解开她衣服和裤子视频手机官网-美女 > 看黄神器免app免vip > HNSCC数据分析-GSE33205- GPL5175

HNSCC数据分析-GSE33205- GPL5175

发布日期:2021-10-12 17:09    点击次数:165
五月份的学徒专注于GEO数据库里面的表达量芯片数据处理,主要的难点是表达量矩阵获取和探针的基因名字转换,合理的分组后就是标准的差异分析,富集分析。主要是参考我八年前的...

五月份的学徒专注于GEO数据库里面的表达量芯片数据处理,主要的难点是表达量矩阵获取和探针的基因名字转换,合理的分组后就是标准的差异分析,富集分析。主要是参考我八年前的笔记:

解读GEO数据存放规律及下载,一文就够

解读SRA数据库规律一文就够

从GEO数据库下载得到表达矩阵 一文就够

GSEA分析一文就够(单机版+R语言版)

根据分组信息做差异分析- 这个一文不够的

差异分析得到的结果注释一文就够

下面是sophie的投稿数据集介绍

GEO链接:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33205

芯片平台:GPL5175,  [HuEx-1_0-st] Affymetrix Human Exon 1.0 ST Array [transcript (gene) version] 样品列表如下(样本太多并未一一列出,详见链接)

GSM822131 RNA-HNSCC-01GSM822132 RNA-HNSCC-02GSM822133 RNA-HNSCC-03GSM822134 RNA-HNSCC-04GSM822135 RNA-HNSCC-05GSM822136 RNA-HNSCC-06...文章链接是:Activation of the NOTCH pathway in head and neck cancer. Cancer Res 2014 Feb 15;74(4):1091-104. PMID: 24351288 https://www.ncbi.nlm.nih.gov/pubmed/24351288 Toward Signaling-Driven Biomarkers Immune to Normal Tissue Contamination. Cancer Inform 2016;15:15-21. PMID: 26884679 https://www.ncbi.nlm.nih.gov/pubmed/26884679核心步骤R包加载rm(list = ls())library(AnnoProbe)library(GEOquery) library(ggplot2) library(ggstatsplot) library(reshape2)library(patchwork)library(stringr)获取并且检查表达量矩阵主要需要判断是否需要log# 获取表达量矩阵gse_number <- 'GSE33205'gset <- geoChina(gse_number)a=gset[[1]] dat=exprs(a) dim(dat)dat = dat[!duplicated(rownames(dat)),]# 检查,判断需不需要取logdat[1:4,1:4] #dat=log2(dat)boxplot(dat[,1:4],las=2)library(limma)dat=normalizeBetweenArrays(dat)# 画图,使用ggplot需宽数据变长数据class(dat)data <- as.data.frame(dat)data <- melt(data)head(data)title <- paste (gse_number, "/", a@annotation, sep ="")p1 <- ggplot(data,aes(x=variable,y=value))+ geom_boxplot()+ theme_ggstatsplot()+ theme(panel.grid = element_blank(), axis.text=element_text(size=10,face = 'bold'), axis.text.x=element_blank(), plot.title = element_text(hjust = 0.5,size =15))+ xlab('')+ ylab('')+ ggtitle(title)p1

可以看到,处理前后我们的表达量矩阵的表达量范围箱线图如下所示:

图片

根据生物学背景及研究目的人为分组pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData#查看pd发现是配对样本,但分组信息存在于两列中,因此,可以手动生成分组'control','HNSCC'library(stringr) group_list=ifelse(grepl('control',pd$characteristics_ch1.2),'control','HNSCC')group_list = factor(group_list,levels = c('control','HNSCC'));table(group_list)

为了演示方便,我们这里仅仅是区分"HNSCC"和“control”。

probe_id 和symbol的转换并对应至表达矩阵获取芯片注释信息#代码如下gpl_number='GPL5175'a = getGEO(gpl_number,destdir = ".")b = a@dataTable@tablecolnames(b)ids = b[,c("ID","Gene Symbol")]#提取需要的列colnames(ids) = c("probe_id","symbol")#改成与前边一致的列名#将ids与前边的对应起来。ids = ids[ids$symbol!="" & !str_detect(ids$symbol,"///"),]head(ids)

可以看到,探针ID和symbol的对应关系如下:probe_id       symbol 1  2315100      DDX11L2 2  2315125       OR4F17 3  2315147 LOC100288692 4  2315160     FLJ45445 5  2315163 LOC100132062 6  2315191 LOC100507658

探针基因ID对应以及去冗余

代码如下:

library(tidyr)library(dplyr)#接下来,使探针与基因symbol一一对应ids=as.data.frame(ids)table(rownames(dat) %in% ids$probe_id)dat=dat[rownames(dat) %in% ids$probe_id,]ids=ids[match(rownames(dat),ids$probe_id),]ids$probe_id=as.character(ids$probe_id)rownames(dat)=ids$probe_idids=ids[ids$probe_id %in% rownames(dat),]dat=dat[ids$probe_id,]#下一步是:基因symbol去冗余-按照表达量最大值筛选ids$median=apply(dat,1,median)#ids新建median这一列,列名为median,同时对dat这个矩阵按行操作,取每一行的中位数,将结果给到median这一列的每一行ids=ids[order(ids$symbol,ids$median,decreasing = T),] #对ids$symbol按照ids$median中位数从大到小排列的顺序排序,将对应的行赋值为一个新的idsids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项,'!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果s#获得去冗余之后的dat/expdat=dat[ids$probe_id,] #新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat#把ids的symbol这一列中的每一行给dat作为dat的行名rownames(dat)=ids$symboldat[1:4,1:4]

最后得到了表达量矩阵如下所示:

> dat[1:4,1:4]  #保留每个基因ID第一次出现的信息 GSM822131 GSM822132 GSM822133 GSM822134ZZZ3 7.185941  6.817970  7.633744  7.574088ZZEF1   7.316541  7.076985  6.504934  7.228776ZYX  8.497323  8.841857  9.060300  8.059482ZYG11B  6.804453  6.486858  7.805016  7.064768

以及最简单的2分组,如下所示:

> table(group_list)group_listcontrol   HNSCC   25   44

保存为R数据文件:step1-output.Rdata

标准步骤之质控

得到标准的3张图,包括主成分分析,高变基因的表达量热图,样品相关性热图

代码如下:

## 下面是画PCA的必须操作,需要看说明书。exp = datexp=t(exp)#画PCA图时要求是行名是样本名,列名时探针名,因此此时需要转换exp=as.data.frame(exp)#将matrix转换为data.frame library("FactoMineR")#画主成分分析图需要加载这两个包library("factoextra") dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的# 画图,主成分分析图p2this_title <- paste0(gse_number,'_PCA')p2 <- fviz_pca_ind(dat.pca, geom.ind = "point", # show points only (nbut not "text") col.ind = group_list, # color by groups palette = "Dark2", addEllipses = TRUE, # Concentration ellipses legend.title = "Groups")+ ggtitle(this_title)+ theme_ggstatsplot()+ theme(plot.title = element_text(size=15,hjust = 0.5))p2# 下面是1000_sd热图library(pheatmap)cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行('1'是按行取,'2'是按列取)取每一行的方差,从小到大排序,取最大的1000个n=t(scale(t(dat[cg,]))) n[n>2]=2 n[n< -2]= -2n[1:4,1:4]ac=data.frame(Group=group_list)rownames(ac)=colnames(n)# 画图,高变基因的表达量热图p3p3 <- pheatmap::pheatmap(n, show_colnames =F, show_rownames = F, main = gse_number, annotation_col=ac, breaks = seq(-3,3,length.out = 100))#因为已经手动设置了表达量最大值,所以,可以不用设置breakp3# 画图,样品相关性热图p4colD=data.frame(Group=group_list)exprSet=t(exp)rownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的expp4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列 操作 annotation_col = colD, show_rownames = F, show_colnames = F, main = gse_number, border_color = NA)p4

出图如下:

图片

标准步骤之limma差异分析#差异分析library(limma)design=model.matrix(~factor( group_list ))exprSet=as.data.frame(t(exp))dim(exprSet)fit=lmFit(exprSet,design)fit=eBayes(fit)options(digits = 4) #设置全局的数字有效位数为4deg = topTable(fit,coef=2,adjust='BH', n=Inf) #结果检验boxplot(dat[rownames(deg)[1],]~group_list)deg[1,]#核对无误

差异分析结果前10行如下所示:

> deg[1:10,]   logFC AveExpr t   P.Value adj.P.Val  BDSG2 2.1623   7.263  11.851 1.801e-18 2.466e-14 31.51MAL -3.2186   8.692 -11.737 2.847e-18 2.466e-14 31.07TMPRSS11B -3.8970   6.119 -11.557 5.870e-18 3.390e-14 30.37ENAH 1.8779   7.569  11.470 8.362e-18 3.622e-14 30.03KRT78  -2.7949   7.675 -11.039 4.804e-17 1.665e-13 28.34LAMB1   1.3238   7.977  10.831 1.127e-16 3.254e-13 27.51SH3PXD2B   1.3014   7.774  10.701 1.926e-16 4.767e-13 27.00SLC2A1  1.7822  10.411  10.272 1.141e-15 2.471e-12 25.27ITGB7  -0.9061   7.271 -10.240 1.300e-15 2.503e-12 25.14BLNK   -1.3173   6.918  -9.991 3.676e-15 6.369e-12 24.14

有了差异分析就可以进行标准的可视化,包括火山图和上下调的差异基因热图

代码如下:

nrDEG=deghead(nrDEG)attach(nrDEG)plot(logFC,-log10(P.Value))#简单画图看一下df=nrDEGdf$v= -log10(P.Value) #df新增加一列'v',作为新的绘图参数,值为-log10(P.Value) #设定上下调基因df$g=ifelse(df$P.Value>0.05,'stable', ifelse( df$logFC >1,'up', ifelse( df$logFC < -1,'down','stable') ))#统计上下调基因数量table(df$g)#给绘制火山图用的数据新增一列symboldf$name=rownames(df)head(df)logFC_t = 1#设置可循环使用的plot标题this_tile <- paste0('Cutoff for logFC is ',round(logFC_t,3), '\nThe number of up gene is ',nrow(df[df$g == 'up',]) , '\nThe number of down gene is ',nrow(df[df$g == 'down',]))#画图,火山图p5p5 <- ggplot(data = df, aes(x = logFC, y = -log10(P.Value))) + geom_point(alpha=0.6, size=1.5, aes(color=g)) + ylab("-log10(Pvalue)")+ scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+ geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) + xlim(-3, 3)+ theme_classic()+ ggtitle(this_tile )+ theme(plot.title = element_text(size=12,hjust = 0.5), legend.title = element_blank(), )p5#热图library(pheatmap)x=deg$logFC names(x)=rownames(deg) cg=c(names(head(sort(x),100)), names(tail(sort(x),100)))#对x进行从小到大排列,取前100及后100,并取其对应的探针名,作为向量赋值给cgn=t(scale(t(exprSet[cg,])))n[n>2]=2n[n< -2]= -2n[1:4,1:4]pheatmap(n,show_colnames =F,show_rownames = F)ac=data.frame(group=group_list)rownames(ac)=colnames(n) # 画图,上下调的差异基因热图p6p6 <- pheatmap(n,show_colnames =F, show_rownames = F, cluster_cols = T, main = gse_number, annotation_col=ac) p6

出图如下:

图片

标准步骤之生物学功能数据库注释

我们这里不根据任何武断的阈值来区分统计学显著的上下调基因,而是直接根据基因的变化情况排序进行gsea分析,而且仅仅是展示kegg这个生物学功能数据库的注释情况!

symbol 和 ENTREZID转换

gsea分析需要基因的ENTREZID,根据物种进行转换

# 加ENTREZID列,用于富集分析(symbol转entrezid,然后inner_join)deg$symbol=rownames(deg)library(org.Hs.eg.db)library(clusterProfiler)s2e <- bitr(deg$symbol, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)#人#bitr()用于SYMBOL转ENTREZID#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDbdim(deg)dim(s2e)setdiff(deg$symbol,s2e$SYMBOL)DEG <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))进行gsea富集library(dplyr)library(ggplot2) geneList=DEG$logFCnames(geneList)=DEG$ENTREZIDgeneList=sort(geneList,decreasing = T)head(geneList)library(clusterProfiler)kk_gse <- gseKEGG(geneList = geneList, organism = 'hsa',#按需替换 #nPerm = 1000, minGSSize = 10, pvalueCutoff = 0.9, verbose = FALSE)tmp=kk_gse@resultdim(tmp)kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#按需替换#DOSE::setReadable():mapping geneID to gene Symboltmp=kk@resultdim(tmp)pro='comp1'write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))save(kk,file = 'gsea_kk.Rdata')富集可视化

上面的kk这个变量就存储了kegg这个生物学功能数据库的gsea分析结果,我们进行简单可视化,代码如下:

# 展现前6个上调通路和6个下调通路down_k <- kk_gse[tail(order(kk_gse$enrichmentScore,decreasing = F)),];down_k$group=-1up_k <- kk_gse[head(order(kk_gse$enrichmentScore,decreasing = F)),];up_k$group=1dat=rbind(up_k,down_k)colnames(dat)dat$pvalue = -log10(dat$pvalue)dat$pvalue=dat$pvalue*dat$group dat=dat[order(dat$pvalue,decreasing = F),]# gsea分析结果p7p7<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + geom_bar(stat="identity") + scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) + scale_x_discrete(name ="Pathway names") + scale_y_continuous(name ="log10P-value") + coord_flip() + theme_ggstatsplot()+ theme(plot.title = element_text(size = 15,hjust = 0.5), axis.text = element_text(size = 12,face = 'bold'), panel.grid = element_blank())+ ggtitle("Pathway Enrichment") p7#具体看上面条形图里面的每个通路的gsea分布情况p8library(enrichplot)p8 <- gseaplot2(kk, geneSetID = rownames(down_k))+ gseaplot2(kk, geneSetID = rownames(up_k))p8

图片

图片

如果你也有类似的数据分析需求,却苦于不会写代码,可以考虑找我们的工程师帮忙哦!转录组产品线明码标价之转录组常规测序服务(仅需799每个样品)明码标价之普通转录组上游分析明码标价之转录组下游分析仅需800元明码标价之转录组测序数据的可变剪切明码标价之RNA-Seq数据的内含子保留分析公共数据库产品线明码标价之公共数据库的生存分析明码标价之公共数据集的WGCNA明码标价之公共数据库探索明码标价之探索新流程(以MSIpred为例)明码标价之任意科研图表绘制(以氨基酸的位点变异图为例)

上一篇:郭德纲 于魁智 耿其昌 唱京剧
下一篇:教师节诗词专辑
TOP