精华内容
下载资源
问答
  • 该书详细地介绍了基因芯片数据的处理及分析流程
  • 利用R的bioconductor包进行分析。由于安装的是R3.5+版本所以实际用的是用biomanager指令,其他基本一样,系统走一遍流程,记录网络有用资料。

    利用R的bioconductor包进行分析。由于安装的是R3.5以上版本所以实际用的是用biomanager指令,其他基本一样。
    不同的包有各类坑,具体可以查阅bioconductor官网寻找解决办法。R包安装踩过的各种bug大全见这篇草稿,没有什么是重装不能解决的,有就换环境(暴言)

    一、初始配置

    1. bioconductor包安装
    if (!requireNamespace("BiocManager", quietly = TRUE))
         install.packages("BiocManager")  #3.5之后的新版本安装
    
    1. 所需包的安装和加载
    BiocManager::install()  
    #BiocManager::install(c("annotate","affycoretools","graph","gcrma","CLL","simpleaffy","affyPLM","RColorBrewer","affy"))
    library()
    

    二、数据预处理

    以CLL基因表达数据包(详见官方说明书)为例。

    1.原始数据导入

    library(CLL)              #载入包
    data("CLLbatch")     #读入CLLbatch数据集
    CLLrma <- rma(CLLbatch)   #这里调用RMA算法来对数据标准化
    e <- exprs(CLLrma)  #读取预处理后所有样品基因的表达值矩阵
    

    表达矩阵(局部)如下图所示,行对应Probeset(而不是基因),列代表不同sample。

    效果图

    补充说明:

    • rma函数:
      ①CLL包自动调用affy包,affy包含有多种处理函数如rma函数;
      ②rma函数是一体化的标准化函数,芯片其他标准化方法见3
    • 查看其他CLL数据集信息
      data.class(CLLbatch) #查看数据类型,发现是AffyBatch类(和ExpressionSet类,SnpSet类都属于eSet类)
      data(disease) 、disease #读入、查看所有样本的状态信息(稳定期、发展期)
      help(AffyBatch)、phenoData(CLLbatch)、featureData(CLLbatch)、protocolData(CLLbatch)、annotation(CLLbatch) #查看"AffyBatch"类的详细介绍
    • 表达量矩阵调取:assayData槽与exprs方法
      assayData槽是AffyBatch类必不可少的,其方法返回结果的第一个元素是矩阵类型,用于保存基因表达矩阵。当使用exprs方法时,调取的就是这个基因表达矩阵。→ExpressionSet对象介绍

    2.质量评估与控制
    三个层次的质控:直观观察(image函数)、平均值方法(simpleaffy包)、数据拟合方法(affyPLM包)

    2.1 image函数查看芯片的灰度图像
    直接观察法比较粗糙(非量化),但可以对芯片质量产生总体认识。

    image(CLLbatch[,1])     #对第一列(一个样本)画图
    

    在这里插入图片描述

    • 几个灰度图的判断指导:
      ①如果图像特别黑,说明型号强度低;如果图像特别亮,说明信号强度很有可能过饱和;
      ②如果芯片图像有斑块现象就很可能是坏片;
      ③如果四个点看不到印花,说明数据有问题。

    2.2 用simpleaffy包对芯片数据进行质量评估
    平均方法假设一组实验中的每个芯片数据对于某个平均值指标都相差不大(质量均匀)。

    library(simpleaffy)
    data(CLLbatch)
    Data.qc <- qc(CLLbatch) # 读取芯片数据并获取质量分析报告
    plot(Data.qc)  # 图型化显示报告
    

    在这里插入图片描述

    • 结果图解释
      ①第一列:样品名称。
      ②第二列:检出率、平均背景噪声(往往较高的平均背景噪声都伴随着较低的检出率)
      ③第三列:蓝色正常,红色异常
      尺度因子蓝色实心圆:取值(-3,3)。
      GAPDH 3’/5’ 比值“o”:不大于1.25,否则数据质量不好
      actin 3’/5’比值“△”:不大于3,否则数据质量不好。
      bioB:芯片检测没有达标。

    2.3 用affyPLM包对芯片数据进行质量评估
    affyPLM包可对芯片原始数据进行拟合回归,最后得到芯片权重(Weights)图、残差(Residuals)图、相对对数表达(RLE,Relative log expression)箱线图、相对标准差(NUSE,Normalized unscaled standard errors)箱线图等等,

    library(affyPLM)
    data(CLLbatch)     #读取数据
    Pset=fitPLM(CLLbatch)  #对数据集做回归分析,结果为一个PLMset类型的对象Pset
    
    image(Pset,type="weights",which=1, main="Weights") #根据计算结果,画权重图
    image(Pset,type="resids",which=1, main="Residuals") #根据计算结果,画残差图
    image(Pset,type="sign.resids",which=1, main="Residuals.sign") #根据计算结果,画残差符号图
    

    在这里插入图片描述

    2.4 总体质量评估:RLE和NUSE箱线图
    接上得到PLMset结果后,一个样品的所有探针组的分布可用箱型图表示。

    RLE:反映基因表达量的一致性趋势,定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。
    NUSE:定义为一个探针组在某个样品的PM值的标准差除以该探针值在各样品中的PM标准差的中位数。

    理论上,随机分布应该画出接近0均匀分布的权重图和接近1均匀分布的残差图。

    library(RColorBrewer)  #RColorBrewer包用于提供配色方案。
    colors<-brewer.pal(12,"Set3")  #载入一组颜色
    
    #利用上一步的affyPLM包的回归计算结果Pset画箱线图
    Mbox(Pset, ylim=c(-1,1),col=colors,main="RLE",las=3)# 绘制RLE(相对对数表达)箱线图
    boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3) #绘制NUSE(相对标准差)箱线图
    

    在这里插入图片描述
    在这里插入图片描述
    CLL1,CLL10 的质量明显异于其他样品,考虑去除

    2.5 查看RNA降解曲线

    “因为RNA是从5’端开始降解的,理论上5‘端的荧光强度应该低于3’端的荧光强度。如果RNA降解曲线的斜率太小,甚至接近于0,很可能是RNA降解太严重,需要作为坏数据去除。”

    library(affy)
    # data("CLLbatch")
    data.deg <- AffyRNAdeg(CLLbatch) # 获取芯片数据和降解数据
    #colors<-brewer.pal(12,"Set3") 
    plotAffyRNAdeg(data.deg, col=colors) #利用上面载入过的一组颜色画图
    legend("topleft", rownames(pData(CLLbatch)), col=colors, lwd=1, inset=0.05, cex=0.5) #在左上部位加注图
    

    在这里插入图片描述
    RNA降解图找出了CLL13需要剔除。

    综合上述所有数据,把质量差的3个样本的数据从数据集剔除。

    CLLbatch <- CLLbatch[, -match(c("CLL10.CEL", "CLL1.CEL", "CLL13.CEL"), sampleNames(CLLbatch))]
    

    2.5 从聚类分析的角度看数据质量
    从芯片之间的相互关系来检验芯片的质量,可以画聚类分析树,或做PCA分析。

    library(CLL)
    library(gcrma)
    library(graph)
    
    data(CLLbatch)
    data("disease") #数据加载
    
    CLLgcrma<-gcrma(CLLbatch) #使用gcrma算法预处理数据
    eset<-exprs(CLLgcrma) #提取基因表达矩阵
    pearson_cor<-cor(eset)# 计算样品两两之间的Pearson相关系数
    dist.lower<-as.dist(1-pearson_cor)# 得到Pearson距离的下三角矩阵
    
    hc<-hclust(dist.lower,"ave")# 聚类分析并画图
    plot(hc)
    

    在这里插入图片描述
    从聚类分析的结果来看,稳定组(红框)和恶化组分不开。再做主成分分析(PCA):

    library(affycoretools)
    samplenames<-sub(pattern="\\.CEL", replacement ="",colnames(eset))  #提取样本名
    groups<-factor(disease[,2])  #提取样本状态
    plotPCA(eset,addtext=samplenames,groups=groups,groupnames=levels(groups))  #被PCA的是表达值矩阵
    

    在这里插入图片描述
    结果类似。两者分不开。

    注意:
    “①只有当聚类图中有明显的类别差异时,才适合考虑去除个别不归类的样品;如果类似CLL数据集来自不同个体差异的原因整体分类被打乱,则不能简单判定所有样品都出了问题。
    ②使用主成分分析时,还必须考虑前2个主成分是否具有代表性,这要看前2个主成分的累计贡献率,如果低于60%,可以考虑采用多维尺度分析来构建分类图。”

    3.背景校正、标准化和汇总*
    除了质量控制剔除不合格的样品,还需要关注下面的三步预处理层面:

    背景校正background correction:去除背景噪声
    标准化normalization:使各次/组测量或各种实验条件下的测量可以相互比较
    汇总summarization:使用一定的统计方法将前面得到的荧光强度值从探针水平汇总到探针组水平

    • 各类预处理方法

    • 方法:
      非一体化方法
      有affy包的expresso函数、threestep函数。自定义参数的空间比较大,但是使用比较复杂,详见help()命令。
      ②一体化预处理算法(推荐)
      包括gcrma、rma、和mas5。有时需进行算法比较,进而确定哪种算法最合适

    dmas5 <- mas5(data.raw)                        
    drma <- rma(data.raw)                                 
    dgcrma <- gcrma(data.raw)    
    

    三、基因芯片数据分析

    (一)选取差异表达基因

    六个关键步骤详解见参考资料。

    1.构建基因表达矩阵

    library(gcrma)
    library(CLL)
    data("CLLbatch")
    data(disease)      #载入所需包和数据
    
    CLLbatch <- CLLbatch[, -match(c("CLL10.CEL", "CLL1.CEL", "CLL13.CEL"),sampleNames(CLLbatch))]  #移除质控环节决定剔除的数据
    CLLgcrma <- gcrma(CLLbatch) #用gcRMA算法进行预处理
    #################
    
    sampleNames(CLLgcrma) <- gsub(".CEL$", "",sampleNames(CLLgcrma)) #remove .CEL in sample names
    disease <- disease[match(sampleNames(CLLgcrma), disease[, "SampleID"]), ] # remove record in data disease about 异常数据
    
    eset <- exprs(CLLgcrma) #获得余下21个样品的基因表达矩阵
    

    2.构建实验设计矩阵

    disease <- factor(disease[, "Disease"])#提取实验条件信息
    design <- model.matrix(~-1+disease) #构建实验设计矩阵
    

    3.构建对比模型(对比矩阵)

    contrast.matrix <- makeContrasts(contrasts = "diseaseprogres. - diseasestable",levels=design)#构建对比模型比较两个实验条件下表达数据
    

    4.线性模型拟合

    fit <- lmFit(eset, design)
    

    5.贝叶斯检验
    使用limma包实现贝叶斯检验。limma包对输入数据要求是必须经过对数转换的表达值,所以前面预处理需要经过log2变换的gcRMA算法

     library(limma)
     fit1 <- contrasts.fit(fit, contrast.matrix) #根据对比模型进行差值计算 
     fit2 <- eBayes(fit1)
    

    fit2我遇到一个Warning message:Z ero sample variances detected, have been offset away from zero
    据说是因为至少有一个gene它的表达值在所有的样品中没有任何变化,这个gene会在计算中被忽略。不影响我们跑。

    6.生成结果报表。

    dif <- topTable(fit2, coef="diseaseprogres. - diseasestable", n=nrow(fit2), lfc=log2(1.5))#生成所有基因的检验结果报告
    dif <- dif[dif[, "P.Value"]<0.01,] #用P.Value=0.01作为阈值进行筛选,得到全部差异表达基因.阈值可以自行划定。
    head(dif) #显示一部分报告结果
    

    在这里插入图片描述
    →差异表达基因分析补充资料

    (二)注释

    “注释实质上是一个ID映射过程,这里是将芯片探侦组的ID映射到基因国际标准名称(Gene symbol)和Entrez ID两种ID上。这里映射GI的目的是为了将GI映射到基因本体论(GO),然后做富集分析。”

    library(annotate) #加载注释工具包
    affydb <- annPkgName(CLLbatch@annotation, type="db") #获得基因芯片注释包名称
    #BiocManager::install("affydb") 
    library(affydb, character.only = TRUE) # 加载该包
    dif$symbols <- getSYMBOL(rownames(dif), affydb) #根据每个探针组的ID获取对应基因Gene Symbol,并作为新的一列
    dif$EntrezID <- getEG(rownames(dif), affydb) # 根据探针ID获取对应基因Entrez ID
    head(dif)# 显示前几行
    

    在这里插入图片描述

    (三)统计分析及可视化

    差异表达分析可以用富集分析方法做。这里介绍了GO的富集分析和KEGG通路的富集分析,分别由GOstats包和GeneAnswers包实现。
    关于富集分析(enrichment analysis)

    1.GOterm检验
    (1)结果表生成

    library(GOstats)
    entrezUniverse <- unique(unlist(mget(rownames(eset), hgu95av2ENTREZID))) # 提取HG_U95Av2芯片中所有探针组对应的EntrezID,unique函数去重
    entrezSelected <- unique(dif[!is.na(dif$EntrezID), "EntrezID"])
    # 提取差异表达基因及其对应的EntrezID
    
    #设置GO富集分析的所有参数
    params <- new("GOHyperGParams", geneIds=entrezSelected, universeGeneIds=entrezUniverse,
                  annotation=affydb, ontology="BP", pvalueCutoff=0.001, conditional=FALSE,
                  testDirection="over")
                 
    hgOver <- hyperGTest(params) #对所有的GOterm根据params参数做超几何检验
    
    bp <- summary(hgOver) #生成所有GOterm的检验结果表
    
    htmlReport(hgOver, file="ALL_go.html") # 同时生成所有GOterm的检验结果文件,每个GOterm都有指向官方网站的链接,以获得详细信息
    head(bp) #显示局部结果
    

    在这里插入图片描述

    (2)基因表达热图

    #BiocManager::install("pheatmap")
    library(pheatmap)
    selected <- eset[rownames(dif), ] #从基因表达矩阵中,选取差异表达基因对应的数据
    rownames(selected) <- dif$symbols # 将selected矩阵每行的名称由探针组ID转换为对应的基因symbol
    
    pheatmap(selected[1:20,], color = colorRampPalette(c("green", "black", "red"))(100),fontsize_row = 4, scale = "row", border_color = NA) #这里只画前20个基因的热图
    

    在这里插入图片描述

    (3)DAG关系图

    library(Rgraphviz)
    
    ghandle <- goDag(hgOver) #显著富集的GO term的DAG关系图
    subGHandle <- subGraph(snodes=as.character(summary(hgOver)[,1]), graph = ghandle)
    plot(subGHandle) #选一部分数据构建局部图
    

    在这里插入图片描述

    2.KEGG通路分析
    (1)生成结果

    #BiocManager::install("GeneAnswers")
    library(GeneAnswers)  #加载所需包
    
    # 选取dif中的三列信息构成新的矩阵,新一列必须是EntrezID
    humanGeneInput <- dif[, c("EntrezID", "logFC", "P.Value")]
    
    # 获得humanGeneInput中基因的表达值
    humanExpr <- eset[match(rownames(dif), rownames(eset)), ]
    humanExpr <- cbind(humanGeneInput[,"EntrezID"], humanExpr)
    
    # 去除NA数据
    humanGeneInput <- humanGeneInput[!is.na(humanGeneInput[, 1]),]
    humanExpr <- humanExpr[!is.na(humanExpr[,1]),]
    
    # KEGG通路的超几何检验
    y <- geneAnswersBuilder(humanGeneInput, "org.Hs.eg.db", categoryType = "KEGG",
    testType = "hyperG", pvalueT=0.1, geneExpressionProfile=humanExpr,verbose = FALSE)
    
    getEnrichmentInfo(y)[1:6] #查看部分结果
    

    在这里插入图片描述

    (2)KEGG关系图

    BiocManager::install("KEGG.db") #下载数据
    yy<- geneAnswersReadable(y)
    geneAnswersConceptNet(yy, colorValueColumn = "logFC", centroidSize = "pvalue", output="interactive")
    

    最后画图一行出现了tcl/tk library not available的问题,没解决。
    https://stackoverflow.com/questions/41133847/r-tcltk-not-loading-on-windows

    在这里插入图片描述
    (3)KEGG热图

    yyy<- geneAnswersSort(yy,sortBy = "pvalue")
    geneAnswersHeatmap(yyy)
    

    在这里插入图片描述

    引用资料:
    https://www.jianshu.com/p/50bcf4ba9d8a
    https://www.jianshu.com/p/fb4217512ec0
    https://blog.csdn.net/tommyhechina/article/details/80356110
    实验课课件

    展开全文
  • GSEA在全基因组表达谱芯片数据分析中的应用,比较系统的介绍了GSEA的用法,有利于初学者的学习
  • 基因芯片,分析流程,相关代码,结果文件
  • 芯片数据分析步骤6 探针注释

    万次阅读 多人点赞 2018-05-22 18:46:36
    注释探针 注释探针的原因 为了防止非特异性结合造成的干扰,芯片厂商往往会使用多个探针检测同一...而后续分析如GSEA,只能对基因进行分析,因此也要求对探针进行注释。 注释探针的方法 1 使用芯片厂商的...

    注释探针

    注释探针的原因

    为了防止非特异性结合造成的干扰,芯片厂商往往会使用多个探针检测同一个基因的表达。因此,芯片厂商不会使用基因名作为探针的名称,而是使用自己定义的探针名称。要合并重复探针,我们必须先对探针进行注释,确定每个探针对应检测哪个基因的表达,然后再合并重复探针。而后续分析如GSEA,只能对基因进行分析,因此也要求对探针进行注释。

    注释探针的方法

    1 使用芯片厂商的注释信息注释

    这个方法是金标准,但也是最不常用的方法。为什么呢?你去芯片厂商的网站上搜索一下就知道了。操作界面非常的user-unfriendly,我找了半天都没找到我想要的注释信息。就更别提下载下来对手头的芯片数据进行注释了。

    2 使用 GPL 文件注释

    这个方法比较常见,操作起来很简单,既可以手动下载GPL信息,也可以用GEOquery包下载GPL信息。唯一的问题就是GPL文件一般比较大,下载下来不是很方便,还要求我们有一定的R语言基础才能进行注释。

    以下用一个例子来演示如何使用GPL文件注释探针。

    首先找到你想要分析的芯片数据的信息。这里使用的是GEO数据库GSE49382的芯片数据。

    点开GSE49382的页面,可以看到GPL信息。用红框框出。是031058-Agilent ATH NAT array的芯片。我尝试去bioconductor里寻找相关的注释包,没找到。也就是说这个芯片不能用bioconductor进行注释,只能用GPL进行注释

    点开GPL看看,可以看到GPL的注释信息。

    注意到,第四列的”ORF”是我们需要的gene name信息。

    下载数据。

    library(GEOquery)
    gse382<-getGEO('GSE49382',destdir =".")##根据GSE号来下载数据,下载_series_matrix.txt.gz
    gpl515<-getGEO('GPL17515',destdir =".")  ##根据GPL号下载的是芯片设计的信息, soft文件
    

    查看列名

    > colnames(Table(gpl515))
    [1] "ID"   "PROBE_ID" "SEQUENCE" "ORF"  "COL" 
    [6] "ROW"     
    

    可以看到第1列和第4列是我们需要保留的信息,其他列都需要舍弃。由于GPL文件中还包含有其他的一些关于实验纪录的信息,因此不能直接对GPL文件进行操作,需要使用Table函数提取GPL文件中的探针注释信息。

     #提取探针注释信息,保留第1列和第4列
    genename <- Table(gpl515)[,c(1,4)]
    

    注释探针只能对表达矩阵进行,因此我们需要先拿到表达矩阵。

    注意到,下载的gse382是一个list格式的文件,因此不能直接使用exprs函数提取表达矩阵。应该使用代码exprs(gpl515[[1]])提取表达矩阵。

    注意到,表达矩阵是Matrix对象,而我们接下来要用到的merge函数不能对Matrix对象使用,因此要先将表达矩阵转换为data.frame对象。否则会报错。Error in fix.by(by.x, x) : 'by'必需指定唯一有效的列

    注意到,表达矩阵并未包含探针的ID,探针的ID作为表达矩阵的rownames存在。下图红框标出。因此我们要先赋予表达矩阵探针ID信息,才能使用merge函数合并表达矩阵和genename

    我们可以使用如下代码正确合并表达矩阵和genename,并删除探针ID。

    exprSet <- as.data.frame(exprs(gse382[[1]]))# 得到表达矩阵,行名为ID,需要转换
    #转换ID为gene name
    exprSet$ID <- rownames(exprSet)
    express <- merge(x = genename, y = exprSet, by = "ID")
    express$ID =NULL
    

    express就是我们得到的注释了探针的表达矩阵。合并重复探针之后就可以进行差异基因筛选或是GSEA了。

    3 使用 bioconductor 注释

    这是最常用的方法,也是要重点介绍的方法。基本原理就是使用bioconductor中的对应的芯片注释包,提取其中注释的信息去注释我们自己的芯片数据。这有个前提,那就是bioconductor中已经收录了芯片的注释包。如果没有收录就不能使用这个方法,只能用GPL信息去注释探针。上文就是一个经典例子。

    使用biconductor注释探针有两种方法,一种是自己从注释包中提取注释信息,然后用merge函数合并。另一种是使用annotate包的函数注释。两种方法的本质都是一样的,从芯片注释包中提取探针注释信息。下面分别讲解两种方法如何操作。

    自行提取注释信息注释探针

    这个方法需要事先下载芯片注释包。除了在网页上能够看到芯片种类,还可以使用show函数查看芯片种类。推荐使用show函数查看芯片信息。因为使用这个函数查看芯片信息的时候,如果没有下载对应的注释包的话会自动下载。免掉了自己上bioconductor查找下载注释包的麻烦。

    下面的例子使用了ALL包的数据。

    首先加载ALL数据。

    > library(ALL)
    > data(ALL)
    > show(ALL)
    ExpressionSet (storageMode: lockedEnvironment)
    assayData: 12625 features, 128 samples 
      element names: exprs 
    protocolData: none
    phenoData
      sampleNames: 01005 01010 ... LAL4 (128 total)
      varLabels: cod diagnosis ... date last seen (21 total)
      varMetadata: labelDescription
    featureData: none
    experimentData: use 'experimentData(object)'
      pubMedIds: 14684422 16243790 
    Annotation: hgu95av2 
    

    Annotation显示芯片种类为hgu95av2,对应的注释包为hgu95av2.db。运行以下代码得到探针注释信息probe2symbol和表达矩阵exprSet

    library(hgu95av2.db) #加载芯片注释包
    columns(hgu95av2.db) #查看芯片注释包提供的注释信息
    probe2symbol <- toTable(hgu95av2SYMBOL) #提取注释信息,这里选择GENESYMBOL进行注释。也可以选择其他ID进行注释,详情请见注释包的技术支持文档。
    exprSet <- as.data.frame(exprs(ALL)) #提取表达矩阵,转换为data.frame对象
    head(exprSet) #查看表达矩阵
    

    可以看到,表达矩阵exprSetrownames就是探针名,而缺少probe_id列。我们需要将包含探针信息的probe_id列加入表达矩阵exprSet中,然后用merge函数合并。代码如下。

    exprSet$probe_id <- rownames(exprSet) #将行名转化为probe_id
    exprSet_symbol <- merge(probe2symbol, exprSet, by = "probe_id") #合并探针注释和表达矩阵
    dim(exprSet) #查看合并前的探针数目
    dim(exprSet_symbol) #查看合并后的探针数目。可以见到除去了部分没有注释的探针。
    exprSet_symbol$probe_id <- NULL #删掉不要的probe_id
    

    使用 annotate 包注释

    使用annotate包对探针进行注释比上述方法都要简单不少,只需要几步就可以完成。

    下面的例子承接上文的ALL数据和注释包,依然是对探针进行GENESYMBOL注释。

    library(annotate)
    exprSet <- as.data.frame(exprs(ALL)) #提取表达矩阵,转换为data.frame对象
    probeset <- rownames(exprSet) #提取行名为探针名
    Symbol <- as.character(as.list(hgu95av2SYMBOL[probeset])) #提取探针对应的GENESYMBOL
    exprSet_symbol <- cbind(Symbol,exprSet)
    

    这样就得到了注释好的表达矩阵。

    还可以使用annotate包提供的getSYMBOL函数提取GENESYMBOL信息或是使用lookUp函数提取注释信息。

    getSYMBOL函数只能提取GENESYMBOL信息并对探针进行注释。而lookUp函数还可以提取ENTREZID等对探针进行注释。不过注意,getSYMBOL函数的返回值是character对象,因此可以直接使用exprSet_symbol <- cbind(Symbol,exprSet)命令进行注释。但lookUp函数返回值是list,需要先转换为character对象才能合并。

    getSYMBOL函数的用法如下。

    Symbol <- getSYMBOL(probeset, "hgu95av2.db")
    

    其中hgu95av2是ALL数据使用的芯片,如果分析自己的芯片的话需要换成自己芯片的注释包。然后就可以直接合并了。

    lookUp函数的用法如下。

    先用columns命令查看hgu95av2包所能提供的注释信息。

    > columns(hgu95av2.db)
     [1] "ACCNUM"   "ALIAS""ENSEMBL"  "ENSEMBLPROT" 
     [5] "ENSEMBLTRANS" "ENTREZID" "ENZYME"   "EVIDENCE"
     [9] "EVIDENCEALL"  "GENENAME" "GO"   "GOALL"   
    [13] "IPI"  "MAP"  "OMIM" "ONTOLOGY"
    [17] "ONTOLOGYALL"  "PATH" "PFAM" "PMID"
    [21] "PROBEID"  "PROSITE"  "REFSEQ"   "SYMBOL"  
    [25] "UCSCKG"   "UNIGENE"  "UNIPROT" 
    

    可以看到hgu95av2.db包提供了27种注释方式。我们选用”SYMBOL”注释探针。如果想用其他注释就把”SYMBOL”换成其他就好了,如”ENTREZID”。

    Symbol <- as.character(lookUp(probeset, "hgu95av2", "SYMBOL"))
    

    lookUp函数的返回值转换成character对象,然后就可以直接合并了。

    总结

    本文一共给出了三种方法用于注释芯片(其实只有两种),但实际上用的最多的是第三种方法。只有当bioconductor没有提供对应的注释包的时候才会用第二种方法。

    注释探针之后就能进行合并探针的操作了,合并之后才能进行后续的分析。过滤探针的操作可以在注释前也可以在注释后。findLargest函数可以在注释探针前就进行探针合并,但需要提供芯片注释包。

    参考:

    1、芯片探针注释基因ID或者symbol,并对每个基因挑选最大表达量探针

    2、基因表达芯片数据分析——Agilent

    3、Bioconductor系列之hgu95av2.db

    展开全文
  • 大家手头的芯片数据一般有两个来源,一个是自己做的芯片的数据,一个是从数据库下载的芯片数据。 如果是自己做的芯片的数据,是一定要进行芯片质量控制的。虽然厂家会提供芯片质量分析的结果,但如果有可能的话,...

    affy芯片质量控制

    前言

    大家手头的芯片数据一般有两个来源,一个是自己做的芯片的数据,一个是从数据库下载的芯片数据。

    如果是自己做的芯片的数据,是一定要进行芯片质量控制的。虽然厂家会提供芯片质量分析的结果,但如果有可能的话,最好还是自己也进行质量分析。根据分析的结果,决定排除哪些芯片的数据,甚至重做也是有可能的。一定只能用质量好的芯片数据,否则可能影响实验结果。自己做的芯片数据在质量控制的阶段一定要严格把关,分析过程越详细越好。

    如果手头的芯片数据是从数据库下载的话,那一般是没有质量问题的。因为上传者在上传数据之前就已经进行过质量控制了,使用者一般不需要进行严格的质量控制,分析流程也不需要太详细以加快处理速度。

    方法

    1 使用 arrayanalysis 网站

    1、登录arrayanalysis网站

    2、选择Get started

    3、上传包含.CEL文件的Zip压缩文件

    4、点击Run affyQC

    5、查看并下载输出的芯片质控图。

    arrayanalysis网站能够对AffymetrixIllumina芯片的原始数据进行质量控制。操作非常简单,输出的质量控制结果非常的详细,甚至详细到没有必要的程度。

    如果在网站上运行质量控制可能速度太慢,我自己上传的文件稍微大一点就要运行半天。如果不想在网页上运行,网站也提供了R脚本。把R脚本下载下来,照着运行就可以了。

    2 使用affy包、oligo包、simpleaffy包和arrayQuanlityMetrics包

    1 灰度图

    芯片灰度图能够检测芯片表面是否均一,是否存在spatial artifact。affy包与oligo包的代码相同

    不存在spatial artifact

    存在spatial artifact

    使用affy包生成灰度图的代码如下(data为读取CEL文件得到的AffyBatch对象):

    library(affy)
    for (i in 1:length(sampleNames(data))){
      name = paste("image",sampleNames(data)[i],".jpg",sep="")
      jpeg(name)
      image(data[,i])
      dev.off()
    }
    

    自动生成所有芯片的灰度图,并储存在工作目录下。

    效果如下。

    2 MA 图

    MA图最初是为双色荧光芯片设计的,但现在也能用于单色荧光芯片。用于单色荧光芯片的时候,要先假想一个芯片。假想芯片的每个探针的强度是所有芯片该探针强度的中位数。然后拿每一个芯片跟假想芯片进行比较,生成MA图。

    使用affy包或oligo包生成MA图的方法是一样的。使用以下代码即可自动在工作目录下生成每一张芯片的MA图。

    for (i in 1:length(sampleNames(data))){
    name = paste("MAplot",sampleNames(data)[i],".jpg",sep="")
    jpeg(name)
    MAplot(data,which=i)
    dev.off()
    }
    

    affy包MAplot默认参数plot.method = "normal",即生成散点图。而oligo包MAplot默认参数plot.method = "smoothScatter",即生成smoothScatter图。如下图。

    3 Chip pseudo-images

    Chip pseudo-images也是用来检测spatial artifact的。方法是先拟合探针水平模型(probe-level model, PLM),然后做pseudo images。探针水平模型假设probe set中的probes在所有的样品中表现一致,与目标序列结合良好的探针均结合良好,反之亦然。

    根据探针水平模型的权重或是残差,就可以生成对应的pseudo-images。

    用affyPLM包和oligo包生成pseudo-images的代码有些许不同。

    用affyPLM包生成Chip pseudo-images的代码如下。

    1、拟合探针水平模型

    library(affyPLM)
    pset <- fitPLM(data, output.param = list(varcov="none"))
    

    参数output.param = list(varcov="none")省略一些不必要的计算,节省内存,加快运算速度。

    2、生成权重图

    for (i in 1:length(sampleNames(data))){
      name = paste("pseudoimageweights", sampleNames(data)[i], ".jpg", sep = "")
      jpeg(name)
      image(pset, type = "weights",which = i)
      dev.off()
    }
    

    效果如下。

    可以看出图中几乎不存在spatial artifact。给个存在spatial artifact的例子给你们看看。

    可以看到存在相当明显的spatial artifact。一般这样的芯片是不能要的。

    3、生成残差图

    for (i in 1:length(sampleNames(data))){
      name = paste("pseudoimageresids", sampleNames(data)[i], ".jpg", sep = "")
      jpeg(name)
      image(pset, type = "resids",which = i)
      dev.off()
    }
    

    效果如下。

    oligo包中的代码与affyPLM不同。因为oligo包读取的表达谱芯片格式为ExpressionFeatureSet,并非AffyBatch,所以不能用affyPLM包进行作图,只能用oligo包自带的作图工具。

    代码如下。

    1、拟合探针水平模型

    pset <- fitProbeLevelModel(data)
    

    2、生成pseudo-images

    生成权重图的代码与affy包的代码完全相同。

    生成残差图的代码有些许不同。

    for (i in 1:length(sampleNames(data))){
      name = paste("pseudoimageresids", sampleNames(data)[i], ".jpg", sep = "")
      jpeg(name)
      image(pset, type = "residuals",which = i)
      dev.off()
    }
    

    4、生成RLE,NUSE图

    Relative Log Expression (RLE)图是最常用的QC图之一。RLE反映了基因表达量的一致性趋势,它定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。每个样品的中心应该非常接近纵坐标0的位置。如果个别样品的表现与其他大多数明显不同,那说明可能这个样品可能有问题。

    Normalized Unscaled Standard Errors (NUSE)也是非常常见的QC图。NUSE是一种比RLE更为敏感的质量检测手段。我们可以结合NUSE图来确定是否某个芯片质量有问题。NUSE定义为一个探针组在某个样品的PM值的标准差除以该探针值在各样品中的PM标准差的中位数。如果所有芯片的质量都是非常可靠的话,那么它们的标准差会非常接近,因此它们的NUSE值都会在1附近。

    使用affy包或oligo包生成RLE和NUSE的代码是一样的。

    生成RLE图的代码。

    RLE(pset)
    

    生成NUSE图的代码。

    NUSE(pset)
    

    3 使用 simpleaffy 包进行质量控制

    simpleaffy整合了许多affy包常用的功能,剔除了不少基本用不着的功能,新增了质量控制的功能。这里只讲作质控图的流程。代码如下。

    library(simpleaffy)
    data.qc <- qc(data)
    plot(data.qc)
    

    结果如下。

    左边上方的数字为该芯片中探针calls为present所占的比例,下方数字为探针背景中探针calls为present所占的比例。两者应当差别不大,最好不超过10%。

    蓝色区域为scale factor。带点的线应该在蓝色区域内。三角形和圆圈为内参探针的比值,actin应当不超过1.25,gapdh应当不超过3.质量合格的标为蓝色,不合格的标为红色。

    4 使用 arrayQualityMetrics 包进行质量控制

    一般上传到GEO数据库中的芯片已经经过了一轮质量控制,基本不会有质量太差的芯片。所以如果偷懒的话,前面的3步都不需要做,只要做第4步就行了。

    arrayQualityMetrics包整合了常用的芯片质量控制工具,能够输出 高清大图 以及提供 html 用于手动操作。操作傻瓜,能够对affy芯片、Illumina芯片和双通道芯片进行质控。输入的格式为ExpressionSet, AffyBatch, ExpressionSetIlluminaNChannelSet

    代码如下

    library(arrayQualityMetrics)
    arrayQualityMetrics(expressionset = data,
    outdir = "fig",
    force = TRUE,
    do.logtransform = T)
    

    参数outdir为结果输出目录,force = TRUE 表示如果文件已存在就覆盖, do.logtransform = T 是因为原始数据没有进过标准化处理,因此需要先进行log变换,才有可比性。

    结果如下。

    每种图片都有图片格式和PDF格式。点击index可以进入web界面。下面给出部分结果图。

    由于用于测试的数据并未经过标准化处理,所以即使存在差别显著的芯片也不应该立刻舍去。应该先进行RMA,MAS5之类的标准化处理之后,再用arrayQualityMetrics包进行质量控制。如果芯片质量还是不好,那就应该舍去。

    展开全文
  • 从GEO数据库下载数据的方法 1、在GEO DATASETS中输入关键词,选择符合的GSE,在...包含有芯片的探针信息,如cDNAs,寡核苷酸,ORFs,抗体。 以GPLxxx编号。 一个platform可以包含不同人上传的不同sample。 不同...

    从GEO数据库下载数据的方法

    1、在GEO DATASETS中输入关键词,选择符合的GSE,在ftp中进行手动下载

    2、找到符合的GSE,在R中使用GEOquery包进行下载

    GEO数据库的数据种类

    1、Platforms 平台

    包含有芯片的探针信息,如cDNAs,寡核苷酸,ORFs,抗体。

    以GPLxxx编号。

    一个platform可以包含不同人上传的不同sample。

    不同platform的数据需要分开处理。

    2、Samples 样品

    一个以独立方式处理的样品。

    以GSMxxx编号。

    一个sample只能包含于一个platform,一个sample可以包含于多个series。

    3、Series 系列

    一个Series就是一个study。

    以GSExxx编号。

    一个系列一定包含多个sample,可能包含多个platform。

    不同platform的数据需要分开处理。

    4、Datasets 数据集

    数据集包含有被summiter处理过的数据,可以使用GEO数据库自带的tools进行分析,如differentiated gene expression, cluster, heatmap。

    以GDSxxx编号。

    一个dataset的sample来自同一个platform,因此彼此间具有可比性。

    范例

    • gds858 <- getGEO(‘GDS858’, destdir=“.”) ##根据GDS号来下载数据,下载soft文件

    • gpl96 <- getGEO(‘GPL96’, destdir=“.”) ##根据GPL号下载的是芯片设计的信息!

    • gse1009 <- getGEO(‘GSE1009’, destdir=“.”)##根据GSE号下载数据,下载_series_matrix.txt.gz

    下载GDS返回的对象

    gds858返回的对象很复杂

    用Table(gds858)可以得到表达矩阵!

    用Meta(gds858)可以得到描述信息

    names(Meta(gds858))
    Table(gds858)[1:5,1:5]
    

    可以用 GDS2eSet 函数把它转变为 expressionset

    下载GSE返回的对象

    GPLList函数查看GPL信息

    处理函数有:geneNames/sampleNames/pData/exprs

    用命令

    gsmplatforms <- lapply(GSMList(gse), function(x) {Meta(x)$platform_id})
    head(gsmplatforms)
    

    查看GSM对应的GPL信息

    用命令

    gsmlist = Filter(function(gsm) { Meta(gsm)$platform_id=='GPLXX'},GSMList(gse))
    

    提取GPLXX对应的样本(有些实验涉及到不同平台的样品)。

    下载GPL返回的对象

    根据GPL号下载返回的对象跟GDS一样,也是用Table/Meta处理!

    还可以下载cel原始文件!

    tmp=getGEOSuppFiles(GSE1009)
    if (is.null(tmp)) {
      warning("Supplementary data files not provided!\nyou should check this GEO ID in NCBI\n")
    }
    

    参考:

    1、用GEOquery从GEO数据库下载数据

    2、Using the GEOquery Package

    3、GEOquery Reference Manual

    展开全文
  • 7基因芯片数据分析.pptx
  • 流行学习算法应用于基因芯片数据分析.pdf
  • 基因芯片数据分析

    万次阅读 多人点赞 2017-07-02 16:17:34
    前言...2 实验概述...2 实验流程...2 数据下载...2 ...聚类分析及可视化...7 差异表达基因GO注释...10 K邻近分类器...10 SVM分类...12 留一法检验-K邻近分类器...14 留一法检验-svm..15 SVM-RFE.1...
  • 基于概率模型-gMOS的基因芯片数据分析.pdf
  • 用于基因芯片数据分析的模块性图聚类方法.pdf
  • 基于并行计算的大规模外显子芯片数据分析.pdf
  • 芯片数据分析步骤2 读取数据-affy

    千次阅读 2018-05-12 15:51:33
    读取affy表达谱芯片数据的方法 Affymetrix表达谱芯片数据读取的方法分3种: 1、使用affy包读取。(HGU95/HGU133芯片) 2、使用oligo包读取。(Whole Transcriptome 芯片/ NimbleGen 芯片/ SNP芯片等) 3、使用...
  • 糖尿病性勃起功能障碍大鼠的基因表达谱芯片数据分析.pdf
  • 将分布式计算系统实体方法用作自动化基因芯片数据分析处理.pdf
  • 小鼠烧伤后时序芯片数据分析及免疫标志物的初步筛选.pdf
  • 利用50K芯片数据分析中国11个地方绵羊群体的遗传结构.pdf
  • 人膀胱癌细胞NRP-1基因表达降低后基因表达谱芯片数据分析.pdf
  • 大样本芯片数据分析亚洲人乳腺癌差异表达基因及可能信号通路.pdf
  • 利用GEO结肠癌芯片数据分析FOXC1及其相关基因的表达与意义.pdf
  • 基于芯片数据分析同病异证H22荷瘤小鼠甲状腺生长因子表达特征.pdf
  • 芯片数据分析步骤7 合并重复探针

    千次阅读 2018-05-27 12:16:28
    为了避免非特异性结合等干扰因素影响实验结果,芯片厂商往往采取多个探针检测同一基因表达的策略,从而导致注释探针后发现许多探针被注释为同一个基因。但在后续的分析中,程序往往不能接受表达矩阵中存在多个探针...
  • minfi:甲基化芯片数据分析

    千次阅读 2020-01-19 22:16:18
    minfi包适用于分析450K和850K的甲基化矩阵。 每个样本在一个矩阵上以两个不同的颜色通道(red、green)进行测量;每个矩阵测量4.5E5个CpG,对每个CpG的甲基化和未甲基化进行度量。 样品之间的DNA甲基化差异可以在...
  • 芯片制造商和芯片型号。然后在百度里输入芯片制造商和芯片型号
  • 芯片数据分析步骤4 标准化-affy

    万次阅读 2018-05-17 20:05:37
    标准化 标准化的原因 芯片实验中存在大量干扰因素,标准化可以削弱这些干扰因素,使得实验条件下的测量可以相互比较。 常见干扰因素:芯片杂交的...二,如果要在后续分析中使用limma包,请不要进行基于方差(vari...
  • 芯片数据分析步骤5 过滤探针

    千次阅读 2018-05-22 18:42:11
    表达谱芯片上的探针往往能够覆盖到所有人类基因,也就是说,能够同时检测所有人类基因的表达。但先前的实验表明,一个细胞中不可能所有基因都同时表达,能够同时表达的基因反而是少数。同时表达的基因约占总基因的40...
  • r语言中使用Bioconductor 分析芯片数据

    千次阅读 2018-07-26 15:00:00
    芯片数据分析流程有些复杂,但使用 R 和 Bioconductor 包进行分析就简单多了。本教程将一步一步的展示如何安装 R 和 Bioconductor,通过 GEO 数据库下载芯片数据, 对数据进行标准化,然后对数据进行质控检查,最后...
  • Bioconductor分析基因芯片数据

    千次阅读 2018-02-20 08:18:06
    使读者初步了解使用Bionconductor完成基因芯片预处理的流程接着详细讲解戏弄i按预处理和数据分析等内容最后深入了解实际工作中会遇到的芯片处理问题以及如何用学到的只是解决问题目的:掌握芯片分析的整体框架,自行...
  • Bioconductor分析基因芯片数据第五章

    千次阅读 2018-02-21 20:29:29
    使读者初步了解使用Bionconductor完成基因芯片预处理的流程接着详细讲解戏弄i按预处理和数据分析等内容最后深入了解实际工作中会遇到的芯片处理问题以及如何用学到的只是解决问题目的:掌握芯片分析的整体框架,自行...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 104,214
精华内容 41,685
关键字:

芯片数据分析