表达芯片数据分析1

原创

Erics blog

修改于 2023-09-27 12:33:27

4180

修改于 2023-09-27 12:33:27

文章被收录于专栏：R语言数据分析R语言数据分析

肿瘤研究数据库：TCGA、ICGC、CCLE、SEER、GEO、NHANES。

可挖掘的数据类型：基因表达芯片、转录组、单细胞、突变，表观等。

基因筛选策略：

数据下载

差异分析

WGCNA（加权共表达网络）

富集分析（PRA）（GSEA）

PPI网络（蛋白互作网络）

预后分析

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'tinyarray') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T,quietly = T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}
#character.only=T避免将pkg作为报名的歧义；

for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T,quietly = T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。

#library报错，就单独安装。

常规图标介绍

1、热图

2、散点图和箱线图

箱线图用于表达单个基因在多个分组之间的表达量差异。

3、火山图

logFC:处理组平均值/对照组平均值的log2.

芯片差异分析的起点是取过Log2的表达矩阵。

logFC可以调整到0.585，log2(1.5)

4、主成分分析

PCA样本聚类图

样本的相对距离反映了样本之间的相似程度，代表样本之间的差异。

主成分分析旨在利用降维的思想，把多指标转化为少数几个综合指标（即主成分）。

GEO背景知识及基因表达芯片的原理

分析思路：

找到GSE数据

下载数据（表达矩阵、临床信息分组信息、GPL编号）

数据探索（有无差异，PCA，热图）

差异分析及可视化（P值及logFC，火山图，热图）

富集分析（KEGG, GO）

芯片数据的表达矩阵

探针ID需要转化为gene symbol；样本信息需要转化为分组信息

芯片的差异分析需要输入表达矩阵（数据分布0-20，无异常值，如NA，Inf等；无异常样本）、分组信息（一一对应，因子，对照组的levels在前）、探针注释（gpl编号，对应关系）。

下载数据

rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) ##R默认设置，60s下载不完成就会停止
options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
##数据在曾老师公司的服务器上。
#研究一下这个eSet
class(eSet)
length(eSet)

eSet = eSet[[1]] 
class(eSet)
#?ExpressionSet
#[1] "ExpressionSet"
#attr(,"package")#ExpressionSet由R包作者定义的数据存储方式
#[1] "Biobase"
#eSet@experimentData##查看；或者直接点击查看。

#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)
range(exp)#看数据范围决定是否需要log，是否有负值，异常值
exp = log2(exp+1) #需要log才log
boxplot(exp,las = 2) #看是否有异常样本
#las:标签是否平等于或垂直于坐标轴las=0:平行；las=2:垂直
##对待异常样本，可以删除异常样本
#或者用函数处理：exp=limma::normalizeBetweenArrays(exp)
#关于表达矩阵里的负值
#取过log有负值,正常；
#没取过log，有负值，错误数据，光信号值不能为负值；一般弃用数据
#有一半负值，做了标准化；一般弃用数据

#(2)提取临床信息
pd <- pData(eSet)

#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

分组及探针注释

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 第一种方法，有现成的可以用来分组的列
  
}else if(F){
  # 第二种方法，眼睛数，自己生成
  Group = rep(c("Disease","Normal"),each = 10)
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
  k = str_detect(pd$title,"Normal");table(k)
  Group = ifelse(k,"Normal","Disease")
}
###rep(c("T","N"),times=c(3,5))
###rep(c("T","N"),each=3)
###group=c("t",rep("N",3),"t")
#str_split(pd$title," ", simplify=T)[,1]
# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("Normal","Disease"))
Group
##检查分组是否正确
data.frame(pd$title,Group)


#2.探针注释的获取-----------------
#捷径
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL570') #这是从27行运行结果里复制下来的代码，能打出代码就不需要再管其他方法了，不能的话看GPL表格里有没有。

#四种方法，方法1里找不到就从方法2找，以此类推。
#捷径里面包含了全部的R包、一部分表格、一部分自主注释
#方法1 BioconductorR包(最常用，已全部收入find_anno里面，不用看啦)
if(F){
  gpl_number #看看编号是多少
  #http://www.bio-info-trainee.com/1399.html #在这里搜索，找到对应的R包
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db") #列出R包里都有啥
  ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
}
# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格，就不能用这种方法。
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
  colnames(b)
  #下一行代码里的列名是从colnames(b)的输出结果里复制的
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol") #和R包里的列名保持一致
  #下面三句代码是将没有对应到基因的探针和对应多个基因的探针去掉
  k1 = ids2$symbol!="";table(k1) 
  k2 = !str_detect(ids2$symbol,"///");table(k2)
  ids2 = ids2[ k1 & k2,]
  # ids = ids2 
  #使用方法二需要将42行F改为T，55行取消注释
}
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,file = "step2output.Rdata")

作业

利用GSE28345表达芯片数据做差异分析火山图。

library(GEOquery)
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) ##R默认设置，60s下载不完成就会停止
options(scipen = 20)#
a="GSE28345"
eSet = getGEO(a, destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据

class(eSet)
length(eSet)

eSet = eSet[[1]] 
class(eSet)


#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)
range(exp)#看数据范围决定是否需要log，是否有负值，异常值
#exp = log2(exp+1) #需要log才log
boxplot(exp,las = 2)
#(2)临床信息
pd <- pData(eSet)

#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 第一种方法，有现成的可以用来分组的列
  
}else if(F){
  # 第二种方法，眼睛数，自己生成
  Group = rep(c("Disease","Normal"),each = 10)
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
  k = str_detect(pd$title,"normotensive");table(k)
  Group = ifelse(k,"normotensive","hypertensive")
}
###rep(c("T","N"),time=c(3,5))

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("normotensive","hypertensive"))
Group
##检查分组是否正确
data.frame(pd$title,Group)


#2.探针注释的获取-----------------
#捷径
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码；find_anno('GPL6244')
ids <- AnnoProbe::idmap('GPL6244') #这是从27行运行结果里复制下来的代码，能打出代码就不需要再管其他方法了，不能的话看GPL表格里有没有。

#四种方法，方法1里找不到就从方法2找，以此类推。
#捷径里面包含了全部的R包、一部分表格、一部分自主注释
#方法1 BioconductorR包(最常用，已全部收入find_anno里面，不用看啦)
if(F){
  gpl_number #看看编号是多少
  #http://www.bio-info-trainee.com/1399.html #在这里搜索，找到对应的R包
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db") #列出R包里都有啥
  ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
}
# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格，就不能用这种方法。
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
  colnames(b)
  #下一行代码里的列名是从colnames(b)的输出结果里复制的
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol") #和R包里的列名保持一致
  #下面三句代码是将没有对应到基因的探针和对应多个基因的探针去掉
  k1 = ids2$symbol!="";table(k1) 
  k2 = !str_detect(ids2$symbol,"///");table(k2)
  ids2 = ids2[ k1 & k2,]
  # ids = ids2 
  #使用方法二需要将42行F改为T，55行取消注释
}
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,file = "step2output.Rdata")
###########################################
rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
nrow(deg)

#3.加change列,标记上下调基因
logFC_t = 1
p_t = 0.05
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#火山图
library(ggplot2)
ggplot(data = deg, 
       aes(x = logFC, 
           y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()
##结束

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

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

r 语言

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

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

r 语言

登录后参与评论

0 条评论

热度