新GEO

原创

用户10758803

发布于 2024-03-10 10:03:52

1720

发布于 2024-03-10 10:03:52

探针的表达量对应基因的表达量

代码和图片来自生信技能树

拿到的exp：行名：探针ID，转化为gene symbol

列名;样本编号需要转化为分组信息

富集分析指定数据：ENTREZID

1.Entrez gene ID：我们一般说的Gnen ID即Entrez gene ID，是用一串数字表示的（在NCBI里面用）

2.Gene Symbol：可以理解为基因的官方名称，如TP53

3.Ensembl ID：Ensembl ID形式：ENSG00000223972

4。log2（FC）=log2（x）-log2（y）=log2（x/y）

log2（FC）常见阈值，1，2，1.5

作者：mayoneday

链接：https://www.jianshu.com/p/3c290fee634f

来源：简书

著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

富集分析结果

description（通路的描述：重要）

pvalue、p.adjust、qvalue：衡量富集分析是否显著的p值（y叔写的R包里p值默认padjust）

geneID：差异基因有哪些是属于这条通路

count：差异基因中属于这条通路的有多少个，即把geneID里的基因个数

GeneRatio：差异基因中有多少个属于该通路 / 差异基因中有多少个被数据库收录（一个通路有很多基因，不可能所有的基因都被数据库收录。我们只是是借用数据库来评估富集）

BgRatio：该通路共有多少个基因 / 数据库中所有通路共有多少个基因

富集分析的意义：衡量每个通路里的基因在差异基因里是否足够多（衡量每条通路中的差异基因？）

即GeneRatio 和 BgRatio 分子的比值

一些图标解释

上调基因、下调基因分别富集的通路。横坐标为p值。

GO term：有颜色的代表我们富集到的通路，没颜色代表我们没富集到，但是通过挖掘GEO数据库里这些通路关系时，能够发掘出的事情。

实战（代码pipeline）

01数据准备，怎么找数据，数据是否要取log

①认准芯片数据：

网页里experiment type：array

②狭义的对象：R包作者以某种方式组织起来的数据

？对象名称

class(eSet)
[1] "ExpressionSet"#对象名称
attr(,"package")
[1] "Biobase"#该对象出自的R包

③不正常表达矩阵：样本里的大量基因应该不存在差别。所以中位线应该在同一条线上。取过log的数据中纵坐标的值在0-20之间

#处理异常表达矩阵

#第一个办法：删掉异常样本

#第二个办法：exp = limma::normalizeBetweenArrays(exp)

整理好数据exp（一行一个基因探针名，一列一个样本）和 pd（临床信息，主要是获取分组）

02 分组与探针注释（对应代码2）

整理好数据，即获取探针名的注释，并将exp的行名替换成基因名。获取分组信息

GPL：平台信息

ids:探针id对应基因注释

1.biocoductor上有针对GPL的注释包（部分）

library(tinyarray)
find_anno(gpl_number)#这里等同于find_anno("GPL6244") ，即gpl_number ="GPL6244"
[1] "`library(hugene10sttranscriptcluster.db);ids <- toTable(hugene10sttranscriptclusterSYMBOL)` and `ids <- AnnoProbe::idmap('GPL6244')` are both avaliable"#现在就可以用两种包中的其中一种给探针注释
ids <- AnnoProbe::idmap('GPL6244')

2.GPL文件解析：下载表格，read.table读取文件，再数据框取子集ID和gene symboll。直接搜芯片数据，例子：直接搜GPL8542，得到以下表格。

特别情况：无注释，无Ensembl ID，只有.Entrez gene ID

这里本来有张图的，图片丢失，大概是指GEO获得的芯片exp里，横坐标不是Ensembl ID，为.Entrez gene ID，就需要我们用另一种函数将其转化为symbol名

这里第一列为探针名，第二行为extrez ID,用将其转化为symbol。详细函数见3.annoGene（只接受ENSEMBL or SYMBOL找注释，所以我们用另一种函数bitr）

gpl = GEOquery::getGEO(filename = "GPL8542_family.soft.gz",destdir = ".")

ids  = gpl@dataTable@table[,1:2]
library(clusterProfiler)
library(org.Hs.eg.db)
e2s = bitr(ids$ORF,fromType = "ENTREZID",toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
ids = merge(ids,e2s,by.x = "ORF",by.y = "ENTREZID")
ids = ids[,2:3]

3.官方

一个探针对应多个基因（symobl）：非特异探针，直接去除

代码3 直接画热图和PCA图

展示的是整理好的表达矩阵直接画热图，或者自己感兴趣的基因画热图、PCA图

这里展示了热图的详细说明

代码4 DEG差异分析

重点：给基因标记上下调！

得到的数据库样子

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
jlibrary(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
#2.加上探针注释
ids = ids[!duplicated(ids$symbol),]
#其他去重方式在zz.去重方式.R。还可以用distinct函数去重
deg <- inner_oin(deg,ids,by="probe_id")
nrow(deg)

#3.加change列,标记上下调基因
logFC_t=1
P.Value_t = 0.05
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类

代码5 差异分析作图

火山图、热图、感兴趣基因的相关性的热图

代码6 富集分析

只要差异分析

gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)

#(2)富集
#以下步骤耗时很长，设置了存在即跳过
f = paste0(gse_number,"_GO.Rdata")
if(!file.exists(f)){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)
  #readable输入数据是entrizid，输出的结果用symbol展示
  ego_BP <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "BP",
                  readable = TRUE)
  #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,ego_BP,file = f)
}
load(f)

#(3)可视化
#条带图
barplot(ego)
barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

#气泡图
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")

代码7 string

PPI网络

标准流程后续：还可以画出PPI网络，网页工具string。要求输入数据为gene symbol

hub基因：处于中枢位置的基因

子网络：把基因拆分成数个亚组

需要把R语言得出得数据导出

cytoscape输入数据有两个：

1.string_interactions_short.tsv

2.一个为属性表格（一列为symbol，一列为logFC）

多分组数据：也只能两两差异分析，因为logFC就是两个组的比值，火山图也只能分开。代码在GEO_learnmore里2.里的多分组数据

多数据联合分析，思路.：1各自两两差异分析，再将差异基因取交集

2.先合并，后差异分析

原则上选择同一芯片平台的GSE？（合并表达矩阵时会丢掉一些基因）

批次效应：

用limma::removeBatchEffect()、sva::ComBat()去除

不要选择一个组全是对照，另一组全是实验组

文献

差异基因和转录因子取交集：差异的转录因子？和铁死亡基因，得到铁死亡的差异基因

一根柱子代表一个样本，展示22种免疫细胞在样本中的分布

ssgsea还可以计算通路富集分数，每个通路对应的基因，作为输入数据，和exp。

行名为通路的名字？

穿插一些实用函数?

?AnnoProbe包里的函数： idmap,geoChina,annogene

1.AnnoProbe里的idmap

ids=idmap('GPL570',destdir=tempdir()) #行名没意义，此函数可以直接找出探针名和基因ID

2.geoChina

geoChina('GSE1009')#直接从网站上下载GSE文件到工作目录下。,,geoChina('GSE1009',destdir=tempdir())#加上后面的参数则只是下载到电脑上

AnnoProbe:geoChina（a list of ExpressionSet, which contains the expression matrix and phenotype data）

可直接当做GEOquery::getGEO() 后面的步骤和代码01_start_GEO里一样

geo_download {tinyarray}直接提取(输出结果：a list with exp,pd and gpl）最简单！

library(AnnoProbe)
a = geoChina('GSE1009')
class(a)
length(a)
a = a[[1]]
exp <- exprs(a)#提取表达矩阵
dim(exp)
exp[1:4,1:4]
library(GEOquery)
pd = pData(a)]#提取临床信息

两个代码块用途一样

gse_number = "GSE56649"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
class(eSet)
length(eSet)
eSet = eSet[[1]]

第三个函数的代码

library(tinyarray)
geo = geo_download("GSE16011")
library(stringr)
#只要tumor样本
k = str_detect(geo$pd$title,"glioma");table(k)
#展示了如果只要exp里的一部分样本，如何提取出来
geo$exp = geo$exp[,k]
geo$pd = geo$pd[k,]

3.annoGene（只接受ENSEMBL or SYMBOL找注释）/clusterProfiler（接受ENTREZID转化为symbol）

①直接拿到探针注释（列名为探针名和symbol），拿到了后要将exp行名（此时为探针名），转化为symbol

GEO里探针名注释用idmap，TCGA基因注释用annogene和trans_array完成exp行名的替换（从ENSEMBLid转化为symbol）

基因注释的作用，告诉你是miRna，ENSEMBLid，在哪条染色体上

library(tinyarray)
?annoGene
IDs <- c("DDX11L1", "MIR6859-1", "OR4G4P", "OR4F5")
ID_type = "SYMBOL"
annoGene(IDs, ID_type)

下面为输出结果

> annoGene(IDs, ID_type)
     SYMBOL                           biotypes
1   DDX11L1 transcribed_unprocessed_pseudogene
3 MIR6859-1                              miRNA
7    OR4G4P             unprocessed_pseudogene
9     OR4F5                     protein_coding
          ENSEMBL  chr start   end
1 ENSG00000223972 chr1 11869 14409
3 ENSG00000278267 chr1 17369 17436
7 ENSG00000268020 chr1 52473 53312
9 ENSG00000186092 chr1 65419 71585

得到

②只拿到探针，且用函数、曾老师的网站都找不到注释。如果这个时候有ENTREZID，则用另一个函数bitr转化为symbol

这里不是！！ENSEMBLid，注意！！

gpl = GEOquery::getGEO(filename = "GPL8542_family.soft.gz",destdir = ".")

ids  = gpl@dataTable@table[,1:2]
library(clusterProfiler)
library(org.Hs.eg.db)
e2s = bitr(ids$ORF,fromType = "ENTREZID",toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
ids = merge(ids,e2s,by.x = "ORF",by.y = "ENTREZID")
ids = ids[,2:3]
head(ids)
colnames(ids) = c("probe_id","symbol")
exp4 = trans_array(geo$exp,ids)

4.trans_array

既可以用于GEO数据行名转化（探针名转化为symbol，但是要自备探针和symbol对应的表格），也可以用于TCGA数据库行名转化

先用annoGene注释，再用trans_array转化（只接受ENSEMBL or SYMBOL找注释）

annoGene(rownames（exp）,ID_type = "ENSEMBL")
library(tinyarray)
exp = trans_array(exp,ids = re,from = "ENSEMBL",to = "SYMBOL")

#GEO数据行名转化
library(clusterProfiler)
library(org.Hs.eg.db)
e2s = bitr(ids$ORF,fromType = "ENTREZID",toType = "SYMBOL",OrgDb = "org.Hs.eg.db")
ids = merge(ids,e2s,by.x = "ORF",by.y = "ENTREZID")
ids = ids[,2:3]
colnames(ids) = c("probe_id","symbol")
exp4 = trans_array(geo$exp,ids)#ids数据见下列

#ids数据
> head(ids)
  probe_id   symbol
1     1_at     A1BG
2     2_at      A2M
3     9_at     NAT1
4    10_at     NAT2

5.tinyarray::draw_boxplot

tinyarray::draw_boxplot 输入数据的行名为基因，列名为样本的表达矩阵exp，还需要分组信息group,，可以挑自己感兴趣的基因g/直接用差异基因出箱线图.带分组信息的箱线图。

灵活运用

这里分组信息也可以是高低风险组，这里的exp横坐标也可以是任何东西，不一定要是基因，也可以是细胞

library(tinyarray)
draw_boxplot(exp[g,],Group)

6.?get_deg_all

?get_deg_all，（集成函数，直接汇总pipeline里的01和02代码，GEO代码3课)

library(tinyarray)
library(stringr)
gse = "GSE56649"
#01 geo_download代码汇总
geo = geo_download(gse)
pd = geo$pd
geo$exp = log2(geo$exp+1)
#,destdir=tempdir()表示不使用工作目录下的路径,
#by_annopbrobe = FALSE，不写这个参数默认为T，从公司服务器上下载数据

#数据有多少个探针，样本，数据分部范围(看是否要取log),geo是一个列表，里面含exp、pd、gpl。但是分组信息仍然需要自己去写
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
             "control",
             "RA")
Group = factor(Group,levels = c("control","RA"))
Group
find_anno(geo$gpl)
ids <- AnnoProbe::idmap(geo$gpl,destdir = tempdir())
dcp = get_deg_all(geo$exp,Group,ids)
head(dcp$deg)
dcp$plots

可以直接出图（差异基因的热图、PCA、火山图，缺点是不能定制，需要去原数据调整）

01-03

一些笔记

仿制PCA图（把自己的数据弄成仿制数据的样子）

热图

改了矩阵里的一个数字，有个数据非常突出，为离群值，这样会使其他的数据都黯然失色，变成一个颜色。所以要设置色带。（色带范围为大部分数据所在的范围，离群值则变为最深的的那个颜色）。设置色带的意义：避免离群值对整张图的影响

用基因画热图，组内各自聚成一簇，说明画热图的基因存在不同的表达模式，所以聚类才能和分组匹配

复杂热图：借助complexheatmap

一行一个基因，一列一个样本，展示基因在不同样本里的表达量

cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",
         breaks = seq(-3,3,length.out = 100)
#设置色带分布范围，矩阵里有离群值时候，超过色带分别范围的值以最深的颜色展示
         )

04差异分析

limma

需要的数据：exp、ids、group、gse_number

最后的表格要加上symbol、entrizid、change。

library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#1.加probe_id列，把行名变成一列
jlibrary(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
#2.加上探针注释
ids = ids[!duplicated(ids$symbol),]
#其他去重方式在zz.去重方式.R。还可以用distinct函数去重
deg <- inner_oin(deg,ids,by="probe_id")
nrow(deg)
#3.加change列,标记上下调基因
logFC_t=1
P.Value_t = 0.05
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类

最后得到的结果

即我们拿到了所有的差异基因。此时可以用这些差异基因做热图、火山图、感兴趣基因的相关性的热图、富集分析。

感兴趣基因的相关性的热图（直接运行代码5）

输入数据为exp（行名为symbol，列名为样本）

library(corrplot)
g = deg$symbol[1:10] # 换成自己感兴趣的基因
g

M = cor(t(exp[g,]))
#计算相关性函数，可直接接受矩阵作为输入数据，
#计算列名的相关性，eg列名是基因，计算基因的相关性
pheatmap(M)

library(paletteer)#用于配色的R包
my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
#有11个配色，将配色逆转一下，
#使红色表示正相关，蓝色表示负相关
my_color = colorRampPalette(my_color)(10)
#把11个颜色重新分配成10种颜色，相关系数取值-1到1
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

r 语言

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

r 语言

#新GEO

作者已关闭评论

0 条评论

热度

新GEO

新GEO

富集分析结果

一些图标解释

实战（代码pipeline）

01数据准备，怎么找数据，数据是否要取log

02 分组与探针注释（对应代码2）

代码3 直接画热图和PCA图

代码4 DEG差异分析

代码5 差异分析作图

代码6 富集分析

代码7 string

PPI网络

文献

穿插一些实用函数?

1.AnnoProbe里的idmap

2.geoChina

3.annoGene（只接受ENSEMBL or SYMBOL找注释）/clusterProfiler（接受ENTREZID转化为symbol）

4.trans_array

5.tinyarray::draw_boxplot

6.?get_deg_all

01-03

一些笔记

热图

04差异分析

感兴趣基因的相关性的热图（直接运行代码5）

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐