精华内容
下载资源
问答
  • Bioconductor分析基因芯片数据

    千次阅读 2018-02-20 08:18:06
    使读者初步了解使用Bionconductor完成基因芯片预处理的流程接着详细讲解戏弄...chip芯片)的芯片处理方法5.1快速入门例5-1 从数据包CLL中载入芯片数据,完成预处理,最后获得基因(探针组)表达矩阵。注意,探针组表...

    使读者初步了解使用Bionconductor完成基因芯片预处理的流程
    接着详细讲解戏弄i按预处理和数据分析等内容

    最后深入了解实际工作中会遇到的芯片处理问题以及如何用学到的只是解决问题
    目的:掌握芯片分析的整体框架,自行学习其他厂商或种类(例如SNP芯片或CHIP-chip芯片)的芯片处理方法

    5.1快速入门
    例5-1 从数据包CLL中载入芯片数据,完成预处理,最后获得基因(探针组)表达矩阵。注意,探针组表达矩阵的行对应的行对应的探针组,而不是基因,基因和探针组的关系见5.2.1.这段程序从载入原始数据(CEL文件)开始,通过预处理得到基因表达矩阵,是芯片数据处理的一个必须步骤

    CLL 数据是慢性淋巴白血病(Chronic lymphocytic leukemia,CLL)数据集,采用了Affymetrix公司的HG_U95Av表达噗芯片,测量了24个样本,12625个探针
    例5-1采用的实验设计方式:两组之间是对照试验(control test),每组内都是平行实验(parallel test),对照实验,简单来说就是为了阐明某种单一因素的效应或者影响,在保持其他因素不变的前提下,测试一定数据的实验组样本呢和对照组样本,并对结果进行比较。平行实验,简单来说就是对同样的一组样本取两个以上相同的样品,以完全一致的条件下进行试验,测试结果的稳定性
    5.2基因芯片基础知识
    5.2.1探针组
    一张基因芯片(以affymetrix表达谱芯片为例)可以包含上万个的探针(通常由25个碱基组成),他们整齐有序地印刷在芯片上。一组探针或者探针组(probe set),来自于一个基因,通常由20对或者11对探针组成,每一对探针都由匹配探针(perfect match,PM)和错配探针(Mi是match,MM)组成,成为探针对(probe pair),MM与PM的序列只有正中央的那个碱基不同,其余的都一致。但是,在一些高密度芯片中,例如外显子芯片(Exon array),每个探针组只有4个PM探针,没有MM探针。

    探针序列的来源叫做参考序列,通常来自于公开的核酸数据库(例如NCBIGeneBank或RefSeq)对于不同的芯片类型,探针组在参考序列中的分布不同,3’表达谱芯片的探针组排布在参考序列3‘末端附近的一至两外显子上,外显子芯片中,每个长度大于25个碱基的外显子都有针对他的探针组:铺瓦芯片(Tilling array)中,探针组覆盖了几乎所有的外显子和内含子
    需要强调的是,芯片数据领域提到的基因表达矩阵往往是以探针组而忽视以基因为单位的,即每行都对应一个探针组的表达量。后面将要降到的差异基因分析也是找打显著性差异的表达的探针组,然后通过ID映射才对应到探针组代表的基因,探针组与基于的关系往往是多个探针组对应一个基因。但是在实际应用中,经常不太注意区分,探针组有时也会被叫做基因

    5.2.2主要的芯片文件格式
    主要的是CEL文件
    affymetrixi芯片原始数据最常用格式为CEL文件,也是芯片预处理和分析的出发点。CEL文件的主要内容就是每个“cell"的灰度信息,"cell"是整个芯片图像划分后得到的小网格,每个小网格中的图像被看作来自一个探针,自CellHeader开始,每行数据对应芯片上的一个”cell"位点,包含5列信息,依次为X坐标,Y坐标,灰度的平均值,灰度的标准差以及用了多少个像素来求这个平均值
    CEL文件只提供了每个探针的灰度信息,还需要基因芯片探针排布的信息(即哪个探针来自哪个探针组),才可以得到芯片上每个探针组对于的表达数据,这就需要CDF文件。另一个重要的是probe文件,他提供了探针的序列信息。afftymetrix公司为每种型号的芯片都提供了对应的CDF文件和Probe文件。CDF文件中的对应关系用户可以自行更改,例如为了应多多个探针组 的ID对应到同一基因ID的现象,有些研究季候就把对应到同一个基因的多个探针合并为一个探针组,并提供修改后的CDF和Probe文件
    图5-3B是affymerixHG-U133A芯片的Probe文件的部分内容,他只包括了一个探针组(名称是“200688_at")的所有探针,共11条序列,文件中第2和3列是对应探针所在的X和Y 坐标,第4列是序列的第13个奸计(中心)位置对齐到一致性序列的相对位置,第5列是对应探针的序列,最后是样品与探针杂交的方向


    5.3基因芯片数据预处理
    基因芯片数据预处理的目的是将探针水平的数据(杂交信号)转换成基因表达数据,主要的数据结构有Affybatch类和Expressionset类,前者用于存储探针水平的数据(相当于CEL文件的内容),而后者用于存储表达水平的数据(相当于基因表达矩阵的内容)。预处理通过质量控制,剔除不合格的芯片(数据),只保留合格的进入下一步处理。然后通过标准化,将所有芯片数据中的基因表达只变换到一个可以比较的水平,用于后续分析
    5.3.1数据输入
    例5-1中,芯片数据的输入是从数据包中得到的,但是在实际应用中,更常见的情况是从CEL中获得数据,无论是数据包还是文件输入,读入的数据会存入一个“affybatch"类型的对象中,可以通过执行help(affybatch)获得更详细的介绍
    头文件:用于描述实验样本、平台等相关信息,其中包括phenoData,featureData,protocolData以及annotation等几个类

    assayData:这是affybatch类必不可少的,他的第一个元素是矩阵类型,用于保存基因表达矩阵。该矩阵的行对应不同的探针组(probe sets),用一个无重复的索引值表示,列对应不同的样品。当使用exprs方法时,调取的就是这个基因表达矩阵
    experimentData:一个MIAME类型的数据,设计这个MIAM类的目的就是用于保存MIAME原则建议的注释信息.MIAME原则是一组指导方针,他建议了一组标准来记录与基因芯片实验设计相关的资料


    5.3.2质量控制
    质量控制对于后续的分析至关重要,原始图像(DAT文件)级别的质量控制一般用个芯片公司自带的软件(如affymetrix公司的GCOS)完成。本节中,质量控制主要集中在CEL文件级别的处理,从最简单的直观观察,到平均值方法,再到比较高级的数据拟合方法。这三个层次的质量控制分别功能分别用image函数simpleaffy包和affyPLM包实现
    直观的查看一下芯片上所有位点的灰度图像
    image函数表示选取的CLLbatch中的第一个基因芯片(即“CLL10.CEL")的数据,然后调用image函数根据CEL文件中的灰度信息画图,affymetrix芯片在印刷时会在四个角印刷特俗的花纹,并且在左上角印刷芯片的名称,花纹与芯片名称可以帮助我们借助这个图像分辨率来了解芯片数据是否可靠。如果无法分辨四角花纹或芯片名称,很可能数据有问题
    根据image函数的图像信息,可以对芯片的信号强度产生一个总体认识:如果图像特别黑,说明信号强度低;如果图像特别亮,说明信号强度很可能过饱和
    尺度因子affymetrix公司规定,用于比较的芯片之间的尺度因子的比例必须小于三
    检测值(detection call)和检出率(percent present):一组探针能否被检测到,用检测值有(present,简称R)、无(Absent,简称A)和不确定(Marginal presen,简称M)来表示检测范围的上下边界(a1及a2)选用了默认值0.04和0.06.检出率,是用所有检测值为p的探针数量除以芯片所有探针组数控得出的百分比。如果检出率过低,表示大部分的基因都未被检测到,很难说明是该芯片实验有问题,还是这个样品的大多数基因本身就很难检测到,有原因是表达量基地或是其他。因此,需要看多个样品之间的相对差别,如果有的样品的检出率与其他的有比较大的差别,那很可能该样品出现了问题
    平均背景噪声(average background):对于每一块芯片,根据所有的MM值作出统计,可以得到背景噪声的平均值、最小值和最大值。往往较高的背景噪声都伴随着最低的检出率,因此这两个指标可以结合使用
    标准内参(internal control genes):mRNA是按照5‘端到3’端的顺序来降解的,芯片探针组也是根据这个顺序来设计的,因此探针组的测量结果可以体现这一趋势。因为大部分的细胞都有β-action和GAPDH基因,所有affymetrix在大部分的芯片里都将他们设置为一组观察RNA降解成都的内参基因。根据这两个基于设计的探针组很好的涵盖了他们3‘端和5’端的每一个区段。通过比较他们3‘端相对于中间或者5’端的信号强度,可以很好地指示出实验质量。affymetrix建议这个比值对于β-action不大于3,对于GAPDH不大于1.25,即可以说明这个芯片的质量可以接受。如果这个比值很高,表明不完整的β-action或者GAPDH的存在,可能源于体外转录不好或者降解非常严重。如果使用的是affymetrix的小样本实验流程(small sample protocol)而不是常用的标准流程(standard protocol),建议使用3’端相对于中间的比值。原因是小样本流程有更扩增次数,有可能产生更多较短的转录序列,不可避免的带来3‘端的偏倚。为了验证杂交的质量
    根据上述标准,可以使用Bionconductor的simpleaffy包对affymetrix芯片数据进行质量评估,最后得到质量控制总览图(图5-8)

    qc图的看法,图5-8是CLL数据集中全部24个芯片数据的质量控制总览图。图5-8中从左至右,第一列是所有样品的名称;第2列是两个数字(对应每个样品),上面是以百分比形式出现的检出率,下面的数字表明平均背景噪音;第3列("QQ stats")最下面的横轴是尺度因子等指标对应的坐标,取值范围从-3到3,用浅蓝色虚线作为边界。第3列用到了三项指标:尺度因子、GAPDH3'/5'比值和action3比值(记做graph3/graph5和action3/action5),分别用实心圆、空心圆和三角标志表示出来。另外,如果第三列中出现红色的”bioB"字样,说明该样品中未能检测到BioB
    简单地讲,所有指标出现蓝色表示正常,红色表示可能存在质量问题。但是根据实际情况不同,还要进一步分析。一般来讲,如果有一个芯片各项指标都不太正常,尤其是BioB无法检测到,建议判定为该芯片实验失败。如图5-3中的样品”CLL15.CEL",这个数据的检出率(38.89%)明显低于其他样品,action3/action5远大于3,而且没有检测到BioB,因此可以判定此数据无效。如果多个芯片都出现了相同的问题,原因则可能是多方面的;如左侧第2列24个芯片的检出率和背景噪声都很高,原因是阈值设定过高,如果降低阈值,大部分就会变蓝;再如,全部芯片都不能检测到BioB,有可能是嵌入探针所针对的DNA溶液加入比例不对
    基于平均值家建设的评价指标都有一个,默认的假设,那就是对于每一块新片,质量是均匀的,不会随着位置变化发生较大的变化。但如果关注芯片的每个小格(Grid),就会发现格与格之间的质量也是有差异的,这可能由于芯片印刷的问题,也可能是杂交过程中出现的问题。那么如何才能得到比较可靠的质量评估,这需要设计多种能反映芯片数据全貌的指标综合分析从而得出最终的结论。这些指标要在对原始数据拟合(回归)的基础上计算得到,然后以图的形式显示,包括:权重(weights)&(residuals)图、相对对数表达(relative log expression,RLE)箱线图、相对标准差(normalized unscaled standard errors,NUSE)箱线图、RNA降解曲线、聚类分析(cluster analysis)图、主成分分析(principal component analysis,PCA)图、信号强度分布图及MA图等,以上功能由Bionconductor中的affyPLM包实现
    一般情况下,在权重图中,绿色代表较低的权重(接近0),白色、灰色代表较高的权重(接近1);在残差图中,红色代表正的高残差,蓝色代表负残差;在残差符号中,红色代表正的残差,蓝色代表负的残差。如果权重和残差都是随机分布的,应该看到绿色均匀分布的权重图和红蓝均匀分布的残差图。图5-9中,左上为原始图像,右上为权重图,左下为残差图,右下为残差符号图。另外,还可以看到,图中左上部出现了一些白色的条块,这是正常的现象,因为有些时候,探针会按照GC比率(GC ratio)排布从而导致白斑的,那什么样的权重和残差图是不可接受的呢
    在对比实验中,即使是相互比较的对照组与实验组之间,大部分基因的表达量还是应该保持一致的,平行实验之间一致性更强。相对对数表达(RLE)箱线图可以反映上述趋势,它定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。一个样品的所有探针组的RLE的分布可以用一个统计学中常用的箱型图形表示。如果使用RLE箱线图来控制CLL数据集的实验质量,每个样品的中心应该非常接近纵坐标0的位置(图5-11)。如果个别样品的表现与其他样品的表现与其他大多数明显不同,那说明这样品有问题

    NUSE是一种比RLE更为敏感的质量检测手段。如果根RLE箱线图对某个芯片的质量产生怀疑,那么再结合NUSE图,这种怀疑就可以确定下来。NUSE定义为一个探针组在某个样品的PM值的标准差除以该探针组在各样品中PM值标准差的中位数,如果所有芯片的质量就是非常可靠的话,那么他们的标准差会十分接近,因此他们的NUSE值会都在1附近。然而,如果有某些芯片质量有问题的话,就会严重地偏离1,进而导致 其他芯片的NUSE值偏向相反的方向。当然,还有一中非常极端的情况,那就是大部分芯片有质量问题,但是他们的标准差却比较接近,反而会显得没有质量问题的NUSE值明显偏离1,所以必须结合RLE及NUSE两个图才能作出更可靠的判断。例如结合图5-11和6-12,可以看出CLL1和CLL10的质量明显有其他yan



    展开全文
  • 芯片质量分析芯片数据预处理获取差异表达基因GO和KEGG分析聚类分析 (本文于2013.09.04更新) 基因芯片技术的特点是使用寡聚核苷酸探针检测基因。前一节使用ReadAffy函数读取CEL文件获得的数据是探针水平的...

    (本文于2013.09.04更新)

    基因芯片技术的特点是使用寡聚核苷酸探针检测基因。前一节使用ReadAffy函数读取CEL文件获得的数据是探针水平的(probe level),即杂交信号,而芯片数据预处理的目的是将杂交信号转成表达数据(即表达水平数据,expression level data)。存储探针水平数据的是AffyBatch类对象,而表达水平数据为ExpressionSet类对象。基因芯片探针水平数据处理的R软件包有affy, affyPLM, affycomp, gcrma等,这些软件包都很有用。如果没有安装可以通过运行下面R语句安装:

    library(BiocInstaller)
    biocLite(c("affy","gcrma", "affyPLM", "affycomp"))
    

    Affy芯片数据的预处理一般有三个步骤:

    • 背景处理(background adjustment)
    • 归一化处理(normalization,或称为“标准化处理”)
    • 汇总(summarization)。

    最后一步获取表达水平数据。需要说明的是,每个步骤都有很多不同的处理方法(算法),选择不同的处理方法对最终结果有非常大的影响。选择哪种方法是仁者见仁智者见智,不同档次的杂志或编辑可能有不同的偏好。

    1 需要了解的一点Affy芯片基础知识

    Affy基因芯片的探针长度为25个碱基,每个mRNA用11~20个探针去检测,检测同一个mRNA的一组探针称为probe sets。由于探针长度较短,为保证杂交的特异性,affy公司为每个基因设计了两类探针,一类探针的序列与基因完全匹配,称为perfect match(PM)probes,另一类为不匹配的探针,称为mismatch (MM)probes。PM和MM探针序列除第13个碱基外完全一样,在MM中把PM的第13个碱基换成了互补碱基。PM和MM探针成对出现。

    我们先使用前一节的方法载入数据并修改芯片名称:

    library(affy)
    library(tcltk)
    filters <- matrix(c("CEL file", ".[Cc][Ee][Ll]", "All", ".*"), ncol = 2, byrow = T)
    cel.files <- tk_choose.files(caption = "Select CELs", multi = TRUE,
                                 filters = filters, index = 1)
    # 最好查看一下文件名称
    basename(cel.files)
    
    ## [1] "NRID9780_Zarka_2-1_MT-0HCA(SOIL)_Rep1_ATH1.CEL" 
    ## [2] "NRID9781_Zarka_2-2_MT-0HCB(SOIL)_Rep2_ATH1.CEL" 
    ## [3] "NRID9782_Zarka_2-3_MT-1HCA(SOIL)_Rep1_ATH1.CEL" 
    ## [4] "NRID9783_Zarka_2-4_MT-1HCB(SOIL)_Rep2_ATH1.CEL" 
    ## [5] "NRID9784_Zarka_2-5_MT-24HCA(SOIL)_Rep1_ATH1.CEL"
    ## [6] "NRID9785_Zarka_2-6_MT-24HCB(SOIL)_Rep2_ATH1.CEL"
    ## [7] "NRID9786_Zarka_2-7_MT-7DCA(SOIL)_Rep1_ATH1.CEL" 
    ## [8] "NRID9787_Zarka_2-8_MT-7DCB(SOIL)_Rep2_ATH1.CEL"
    
    data.raw <- ReadAffy(filenames = cel.files)
    sampleNames(data.raw) <- paste("CHIP",1:length(cel.files),sep="")
    

    用pm和mm函数可查看每个探针的检测情况:

    pm.data.raw <- pm(data.raw)
    head(pm.data.raw, 2)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8
    ## 501131 127.0 166.3   112 139.8 111.3  85.5 126.3 102.8
    ## 251604 118.5 105.0    82 101.5  94.0  81.3 103.8 103.0
    
    mm.data.raw <- mm(data.raw)
    head(mm.data.raw, 2)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8
    ## 501843  89.0  88.0  80.5  91.0  77.0    75  79.0  72.0
    ## 252316 134.3  77.3  77.0 107.8  98.5    75  99.5  71.3
    

    上面显示的列名称就是探针的名称。而基因名称用probeset名称表示:

    head(geneNames(data.raw))
    
    ## [1] "244901_at" "244902_at" "244903_at" "244904_at" "244905_at" "244906_at"
    

    probeset不一定和实际基因一一对应,这些后面对探针进行基因名称映射时会看到。

    2 背景处理(background adjustment)

    虽然说是背景处理,但是这一步既处理背景值,又处理噪声信号。注意背景和噪声是两个概念,比如说乡间夜晚的蛙叫声虽嘈杂但很稳定,算是背景,如果突然来一声狗叫,那就是噪声。这两者在统计上可以区分。

    芯片的背景处理理论上很简单,因为Affy公司设计MM的目的就是检测非特异杂交信号,PM -MM 不就结了?但是研究发现居然有多达30%的MM探针获得的信号强度比相应PM探针的还强(嘿嘿),所以啊,研究者的饭碗就出来了,整些看起来还合理的方法吧。

    R软件包affy用于芯片背景噪声消减的函数是bg.correct(),而MAS和RMA方法是最常用的两种方法。

    MAS方法将芯片分为k(默认值为16)个网格区域,用每个区域使用信号强度最低的2%探针去计算背景值和噪声。RMA方法的原理比较复杂,可以参看文献:

    R. A. Irizarry, B. Hobbs, F. Collin, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics, 4:249–64, 2003b. 11, 18, 27, 232, 241, 432, 443。

    data.rma <- bg.correct(data.raw, method="rma")
    data.mas <- bg.correct(data.raw, method="mas")
    class(data.rma)
    
    ## [1] "AffyBatch"
    ## attr(,"package")
    ## [1] "affy"
    
    class(data.mas)
    
    ## [1] "AffyBatch"
    ## attr(,"package")
    ## [1] "affy"
    
    class(data.raw)
    
    ## [1] "AffyBatch"
    ## attr(,"package")
    ## [1] "affy"
    

    可以看到:ReadAffy()读入的CEL芯片数据以AffyBatch类数据形式存储,而背景消减后得到的依然是AffyBatch类数据。

    MAS方法应用后PM和MM的信号强度都被重新计算。RMA方法仅使用PM探针数据,背景调整后MM的信号值不变。

    head(pm(data.raw)-pm(data.mas),2)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8
    ## 501131 79.34 62.93 63.23 72.87 62.48 64.43 62.97 60.86
    ## 251604 77.57 63.06 63.73 80.69 63.07 66.53 63.33 60.34
    
    head(pm(data.raw)-pm(data.rma),2)
    
    ##        CHIP1  CHIP2 CHIP3  CHIP4 CHIP5 CHIP6  CHIP7 CHIP8
    ## 501131 111.2 102.21 93.23 116.36 93.75 76.15 102.82 85.21
    ## 251604 103.9  88.69 72.57  90.75 82.23 72.64  87.75 85.32
    
    head(mm(data.raw)-mm(data.mas),2)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8
    ## 501843 79.34 62.93 63.24 72.90 62.48 64.44 62.97 60.85
    ## 252316 77.56 63.06 63.73 80.66 63.07 66.51 63.33 60.34
    
    #差值为全部为0,说明rma方法没有对mm数值进行处理
    head(mm(data.raw)-mm(data.rma),2)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8
    ## 501843     0     0     0     0     0     0     0     0
    ## 252316     0     0     0     0     0     0     0     0
    
    identical(mm(data.raw), mm(data.rma))
    
    ## [1] TRUE
    

    3 归一化处理(normalization)

    同一个RNA样品用相同类型的几块芯片进行杂交,获得的结果(信号强度等)都不可能完全相同,甚至差别很大。为了使不同芯片获得的结果具有可比性,必需进行归一化处理。这一步的方法也很多。归一化处理的affy函数为normalize(),以Affybatch对象和处理方法为参数。

    3.1 线性缩放方法

    这是Affy公司在其软件(4.0和5.0版本)中使用的方法。这种方法先选择一个芯片作为参考,将其他芯片和参考芯片逐一做线性回归分析,用回归直线(没有截距)对其他芯片的信号值做缩放。Affy公司的软件做回归分析前去除了2%最强和最弱信号。使用很简单,mas方法做背景消减后做归一化分析的R语句是:

    data.mas.ln <- normalize(data.mas, method = "constant")
    head(pm(data.mas.ln)/pm(data.mas), 5)
    
    ##        CHIP1  CHIP2 CHIP3 CHIP4  CHIP5 CHIP6  CHIP7  CHIP8
    ## 501131     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 251604     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 261891     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 230387     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 217334     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    
    head(mm(data.mas.ln)/mm(data.mas), 5)
    
    ##        CHIP1  CHIP2 CHIP3 CHIP4  CHIP5 CHIP6  CHIP7  CHIP8
    ## 501843     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 252316     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 262603     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 231099     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    ## 218046     1 0.8788 1.002 1.141 0.9929 1.029 0.7578 0.8834
    

    可以看出,线性缩放方法以第一块芯片为参考,它的数值没有被处理,而其他芯片都被缩放了。对同一块芯片,不同探针的缩放倍数是一个常数。PM和MM的缩放方法完全一样。

    3.2 非线性缩放方法

    此方法获得的结果比线性方法要好,做非线性拟合时不是取整张芯片而仅取部分(一列)作为基线。

    data.mas.nl <- normalize(data.mas, method = "invariantset")
    head(pm(data.mas.nl)/pm(data.mas), 5)
    
    ##        CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6  CHIP7  CHIP8
    ## 501131     1 1.050 1.417 1.359 1.399 2.023 1.0086 1.3124
    ## 251604     1 1.348 2.279 1.838 1.587 2.431 1.1126 1.3059
    ## 261891     1 1.564 1.397 1.675 1.497 1.556 1.3167 1.5509
    ## 230387     1 1.259 1.683 1.390 1.360 1.504 1.0145 1.2239
    ## 217334     1 1.009 1.127 1.241 1.229 1.263 0.8417 0.9934
    

    可以看到,同一芯片不同探针的信号值的缩放倍数是不一样的。

    3.3 分位数(quantile)方法

    这种方法认为(或假设)每张芯片探针信号的经验分布函数应完全一样,使用任两张芯片的数据做QQ图应该得到一条斜率为1截距为0的直线。

    data.mas.qt <- normalize(data.mas, method = "quantiles")
    head(pm(data.mas.qt)/pm(data.mas), 5)
    
    ##         CHIP1  CHIP2 CHIP3 CHIP4 CHIP5 CHIP6  CHIP7  CHIP8
    ## 501131 0.7176 0.9602 1.140 1.181 1.112 1.187 0.8145 1.0156
    ## 251604 0.6984 0.9814 1.199 1.209 1.112 1.172 0.8287 1.0143
    ## 261891 0.6939 0.9956 1.137 1.205 1.111 1.186 0.8389 1.0579
    ## 230387 0.7653 0.9777 1.171 1.183 1.115 1.185 0.8153 0.9921
    ## 217334 0.9508 0.9547 1.063 1.162 1.130 1.173 0.7945 0.9266
    

    3.4 其他

    如循环局部加权回归法(Cyclic loess)和 Contrasts方法。

    data.rma.lo <- normalize(data.rma, method = "loess")
    
    ## Done with 1 vs 2 in iteration 1 
    ## Done with 1 vs 3 in iteration 1 
    ## Done with 1 vs 4 in iteration 1 
    ## Done with 1 vs 5 in iteration 1 
    ## Done with 1 vs 6 in iteration 1 
    ## Done with 1 vs 7 in iteration 1 
    ## Done with 1 vs 8 in iteration 1 
    ## Done with 2 vs 3 in iteration 1 
    ## Done with 2 vs 4 in iteration 1 
    ## Done with 2 vs 5 in iteration 1 
    ## Done with 2 vs 6 in iteration 1 
    ## Done with 2 vs 7 in iteration 1 
    ## Done with 2 vs 8 in iteration 1 
    ## Done with 3 vs 4 in iteration 1 
    ## Done with 3 vs 5 in iteration 1 
    ## Done with 3 vs 6 in iteration 1 
    ## Done with 3 vs 7 in iteration 1 
    ## Done with 3 vs 8 in iteration 1 
    ## Done with 4 vs 5 in iteration 1 
    ## Done with 4 vs 6 in iteration 1 
    ## Done with 4 vs 7 in iteration 1 
    ## Done with 4 vs 8 in iteration 1 
    ## Done with 5 vs 6 in iteration 1 
    ## Done with 5 vs 7 in iteration 1 
    ## Done with 5 vs 8 in iteration 1 
    ## Done with 6 vs 7 in iteration 1 
    ## Done with 6 vs 8 in iteration 1 
    ## Done with 7 vs 8 in iteration 1 
    ## 1 0.05775
    
    data.rma.ct <- normalize(data.rma, method = "contrasts")
    
    ## Data size 502156 x 8 Desired subset size 5000 +- 50 
    ## Comuting ranks of old subset....Size of new subset:  109873 
    ## Comuting ranks of old subset....Size of new subset:  64185 
    ## Comuting ranks of old subset....Size of new subset:  27795 
    ## Comuting ranks of old subset....Size of new subset:  14886 
    ## Comuting ranks of old subset....Size of new subset:  12234 
    ## Comuting ranks of old subset....Size of new subset:  7835 
    ## Comuting ranks of old subset....Size of new subset:  7001 
    ## Comuting ranks of old subset....Size of new subset:  5289 
    ## Comuting ranks of old subset....Size of new subset:  5230 
    ## Comuting ranks of old subset....Size of new subset:  5164 
    ## Comuting ranks of old subset....Size of new subset:  5124 
    ## Comuting ranks of old subset....Size of new subset:  5065 
    ## Comuting ranks of old subset....Size of new subset:  5041 
    ## Fitting normalizing curve
    ## Normalizing chip  1 
    ## Normalizing chip  2 
    ## Normalizing chip  3 
    ## Normalizing chip  4 
    ## Normalizing chip  5 
    ## Normalizing chip  6 
    ## Normalizing chip  7 
    ## Normalizing chip  8
    

    4 汇总(Summarization):

    最后一步汇总是使用合适的统计方法通过probeset(包含多个探针)的杂交信号计算出计算表达量。affy包的函数为computeExprSet。需要注意的是computeExprSet函数除需要指定统计方法外还需要指定PM校正的方式:computeExprSet(x, pmcorrect.method, summary.method, …)

    两个参数可以设定的值可以通过下面函数获得:

    pmcorrect.methods()
    
    ## [1] "mas"        "methods"    "pmonly"     "subtractmm"
    
    generateExprSet.methods()
    
    ## [1] "avgdiff"      "liwong"       "mas"          "medianpolish"
    ## [5] "playerout"
    

    常用的汇总方法是medianpolish, liwong和mas。liwong方法仅使用PM做背景校正(pmcorrect.method="pmonly")。例如:

    eset.rma.liwong <- computeExprSet(data.rma.lo, pmcorrect.method="pmonly",
                                      summary.method="liwong")
    
    ## 22810 ids to be processed
    ## |                    |
    ## |####################|
    

    计算过程需要一定的时间,耐心等候完成。最后的结果 ExpressionSet 类型数据:

    class(eset.rma.liwong)
    
    ## [1] "ExpressionSet"
    ## attr(,"package")
    ## [1] "Biobase"
    
    class(data.raw)
    
    ## [1] "AffyBatch"
    ## attr(,"package")
    ## [1] "affy"
    

    5 三合一处理和“自动化”函数:

    了解芯片预处理的原理和步骤后你完全可以用一个R函数完成三步处理。affy包提供的三合一处理函数为 expresso(),用法为:

    # NOT RUN. 函数说明,不可运行。
    expresso(
            afbatch,
            bg.correct = TRUE,
            bgcorrect.method = NULL,
            bgcorrect.param = list(),
            normalize = TRUE,
            normalize.method = NULL,
            normalize.param = list(),
            pmcorrect.method = NULL,
            pmcorrect.param = list(),
            summary.method = NULL,
            summary.param = list(),
            summary.subset = NULL,
            verbose = TRUE,
            widget = FALSE)
    

    例如:

    # NOT RUN:
    eset.mas <- expresso(data.raw, bgcorrect.method="mas", normalize.method="constant",
                         pmcorrect.method="mas", summary.method="mas")
    

    使用的各类方法可用bgcorrect.methods,pmcorrect.methods 和 express.summary.stat.methods函数了解。

    expresso速度较慢,一些R软件包提供速度较快的预编译三合一函数,如affyPLM软件包的threestep函数。affyPLM提供的处理方法很丰富,有兴趣可以自学。

    下面介绍介3个“自动化”函数,这些函数已经预定义了预处理各步骤所采用的方法和参数,不再需要设定。

    5.1 mas5函数

    由affy包提供:

    eset.mas5 <- mas5(data.raw)
    
    ## background correction: mas 
    ## PM/MM correction : mas 
    ## expression values: mas 
    ## background correcting...done.
    ## 22810 ids to be processed
    ## |                    |
    ## |####################|
    

    mas5据说是expresso函数和mas方法的封装。但使用下面语句获得的结果似乎与 mas5(data.raw)的结果不一样(自己去验证看看):

    # NOT RUN:
    expresso(data.raw, bgcorrect.method="mas", normalize=FALSE,
             pmcorrect.method="mas", summary.method="mas")
    

    5.2 rma函数

    也是由affy包提供,其背景处理方法为rma法,归一化处理使用分位数法,而汇总方法使用medianpolish:

    eset.rma <- rma(data.raw)
    
    ## Background correcting
    ## Normalizing
    ## Calculating Expression
    

    它等价于:

    # NOT RUN:
    eset.rma <- expresso(data.raw, bgcorrect.method="rma", normalize.method="quantiles",
                         pmcorrect.method="pmonly", summary.method="medianpolish")
    

    但rma函数是预编译好的C语言函数,速度非常快!

    5.3 gcrma函数

    由gcrma包提供。 Affy的软件(比如mas方法)使用MM数据做背景处理,但由于MM出现的问题(上面提到过),这些方法可能高估了背景值。而rma方法在做背景处理时没有使用MM数据,这可能又低估了背景值。MM序列公布后有人对其特异性进行了评估,并使用这些评估结果建立了新方法。gcrma就是这样的方法,也是封装好的三合一方法。

    library(gcrma)
    eset.gcrma <- gcrma(data.raw)
    
    ## Adjusting for optical effect........Done.
    ## Computing affinities.Done.
    ## Adjusting for non-specific binding........Done.
    ## Normalizing
    ## Calculating Expression
    

    6 芯片处理方法的优劣评估

    芯片预处理的方法这么多,哪个好?我选哪个?知道得越多越迷惑。幸好这些已经有人做了,牛人Rafael Irizarry 和 Leslie Cope专门写了一个R软件包affycomp用于方法评估。

    评估方法的优劣必需有数据,而且是包含已知因素的数据。affycomp需要两个系列的数据,一个是RNA稀释系列芯片数据,称为Dilution data,另一个是使用了内参/外标RNA的芯片,称为Spike-in data。Spike-in RNA是在目标物种中不存在、但在芯片上含有相应检测探针的RNA,比如Affy的拟南芥芯片上有几个人或细菌的基因检测探针。由于稀释倍数已知,内参/外标的RNA量和杂交特异性也已知,所以结果可以预测,也就可以用做方法评估。对于严格的芯片实验来说,这些实验都是必须的。但是绝大多数人不做,因为成本很高,尤其是只做几张芯片的时候,一般直接使用别人认可的方法如RMA或MAS。

    基因芯片(Affymetrix)分析2:芯片数据预处理 - xxx - xxx的博客

    7 Session Info

    sessionInfo()
    
    ## R version 3.0.1 (2013-05-16)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## 
    ## locale:
    ##  [1] LC_CTYPE=zh_CN.UTF-8       LC_NUMERIC=C              
    ##  [3] LC_TIME=zh_CN.UTF-8        LC_COLLATE=zh_CN.UTF-8    
    ##  [5] LC_MONETARY=zh_CN.UTF-8    LC_MESSAGES=zh_CN.UTF-8   
    ##  [7] LC_PAPER=C                 LC_NAME=C                 
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
    ## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C       
    ## 
    ## attached base packages:
    ## [1] parallel  tcltk     stats     graphics  grDevices utils     datasets 
    ## [8] methods   base     
    ## 
    ## other attached packages:
    ## [1] ath1121501probe_2.12.0 gcrma_2.33.1           ath1121501cdf_2.12.0  
    ## [4] AnnotationDbi_1.23.18  affy_1.39.2            Biobase_2.21.6        
    ## [7] BiocGenerics_0.7.4     zblog_0.0.1            knitr_1.4.1           
    ## 
    ## loaded via a namespace (and not attached):
    ##  [1] affyio_1.29.0         BiocInstaller_1.11.4  Biostrings_2.29.15   
    ##  [4] DBI_0.2-7             digest_0.6.3          evaluate_0.4.7       
    ##  [7] formatR_0.9           highr_0.2.1           IRanges_1.19.28      
    ## [10] preprocessCore_1.23.0 RSQLite_0.11.4        splines_3.0.1        
    ## [13] stats4_3.0.1          stringr_0.6.2         tools_3.0.1          
    ## [16] XVector_0.1.0         zlibbioc_1.7.0
    展开全文
  • 芯片数据分析步骤5 过滤探针

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

    过滤探针

    过滤探针的原因

    表达谱芯片上的探针往往能够覆盖到所有人类基因,也就是说,能够同时检测所有人类基因的表达。但先前的实验表明,一个细胞中不可能所有基因都同时表达,能够同时表达的基因反而是少数。同时表达的基因约占总基因的40%左右。

    由于探针与目标之间一定存在着非特异性结合,所以所有的探针均会产生信号。如果不加以过滤,认为这些探针对应的基因都表达,即不符合事实,也会对后续的分析产生影响。因此,过滤掉表达量低的探针是十分必要的。

    注意,limma包的说明里面提供了两点建议。一,如果要进行探针过滤(filter),最好在进行标准化之后再过滤。二,如果要在后续分析中使用limma包,请不要进行基于方差(variance)的过滤,否则会影响方差分布,从而导致limma包处理产生糟糕的结果(poor results)。

    过滤探针的方法

    bioconductor提供了两种过滤探针的方法,一种是使用专门进行探针过滤的genefilter包进行过滤,另一种是使用affy包中的mas5calls函数进行探针过滤。下面会详细解释如何进行探针过滤。

    1 使用genefilter包进行探针过滤

    根据genefilter包的技术文档介绍,genefilter包设计的出发点是进行independent filtering。我们知道,检测差异基因表达是表达谱芯片最重要的功能之一。而检测差异基因表达的方法是进行对每一个基因进行统计检验。

    所谓的independent filtering,意思是在进行统计检验之前,筛掉那些不能或是很可能不能通过显著性检验的探针。independent filtering不依赖于统计检验,在保持一类错误率(即假阳性率)不变的条件下,能够增加检验效能。

    一个好的filtering criterion拥有下面三个特点:

    1、在零假设的条件下,独立于统计检验

    2、在备择假设的条件下,与统计检验相关联

    3、在联合检验中不显著改变依赖结构

    genefilter包中有两个常用的函数能够进行探针过滤,一个是genefilter函数,另一个是nsFilter函数。

    使用genefilter函数进行探针过滤

    genefilter函数需要用户指定筛选检验方式。用法如下。

    genefilter(expr, flist)
    

    expr是表达矩阵,需要用expr(ExpressionSet)expr(AffyBatch)对ExpressionSet或AffyBatch进行提取。flist是包含检验方法的list,可以同时包含多个检验方法,同时进行多个检验,取交集。

    举个例子。例子来源于genefilter包的技术支持文档。sample.ExpressionSet来自Biobase包。

    library("Biobase") 
    library("genefilter")
    data(sample.ExpressionSet)
    f1 <- kOverA(5, 200)
    ffun <- filterfun(f1)
    wh1 <- genefilter(exprs(sample.ExpressionSet), ffun)
    sum(wh1)
    [1] 159
    

    f1函数的意思是指筛选出至少在5个样本中表达量超过200的基因。最终筛选到了159个基因。

    再给一个例子。

    f2 <- ttest(sample.ExpressionSet$type, p=0.1)
    wh2 <- genefilter(exprs(sample.ExpressionSet), filterfun(f2))
    sum(wh2)
    [1] 88
    

    f2函数的意思是进行 t.test,显著性水平 p=0.1,筛选出 t.test 中 p < 0.1 的基因,最终筛选到88个基因。

    将两个筛选方法合并。

    ffun_combined <- filterfun(f1, f2)
    wh3 <- genefilter(exprs(sample.ExpressionSet), ffun_combined)
    sum(wh3)
    [1] 35
    

    满足f1和f2的基因只有35个。

    过滤探针

    sample.filter <- sample.ExpressionSet[wh3,]
    
    nrows(sample.filter)
    [1] 35
    

    用户可以按照上述例子的方式,用genefilter包内整合的检验方法,构造自己的检验标准进行探针过滤。详情请见bioconductor网页的genefilter包技术支持文档

    使用nsFilter函数进行探针过滤

    虽然使用genefilter函数能够自由地指定筛选标准,但很多时候我们并不需要对筛选标准进行改变,,而是希望使用同一套标准,或是懒得去调参数。

    nsFilter函数提供了一站式的探针过滤,能够一步对ExpressionSet的探针进行过滤,返回一个listlist中的eset为过滤后的ExpressionSetfilter.log为每一步过滤到多少探针的记录。nsFilter函数也可以用于oligo包读取的对象。

    使用方法如下。

    nsFilter(eset, require.entrez=TRUE,
    remove.dupEntrez=TRUE, var.func=IQR,
    var.cutoff=0.5)
    

    一般只对以下几个参数进行调整。

    参数 注释
    require.entrez 若为TRUE,过滤掉没有ENTREZ ID的探针。默认为TRUE。
    remove.dupEntrez 若为TRUE,当几个探针对应统同一ENTREZ ID的时候,留下 var.func 值最大的探针,其余过滤。默认为TRUE。
    var.func 用于过滤的统计参数。默认为IQR。
    var.cutoff 截断值。默认为0.5,即过滤到50%的基因。

    这里要说明一下。由于探针过滤只能对 ExpressionSet 进行,不能对 AffyBatch 进行,因此要先进行标准化,再过滤。默认 var.cutoff = 0.5, 即过滤掉50%的探针。

    nsFilter 函数还提供了合并探针的功能,不过只针对 EntrezID。因此,如果要使用 nsFilter 函数提供的合并探针功能,就必须先对探针进行EntrezID注释。如果不使用 nsFilter 函数的合并探针功能,那先注释或后注释甚至不注释都没关系。

    我喜欢用 gene symbol 注释探针而不喜欢用 EntrezID,因此也就不使用 nsFilter 函数合并重复探针的功能。而且,nsFilter函数只能对ExpressionSet对象进行,而我们一般不会对ExpressionSet对象进行注释,所以合并重复探针的功能用不着。我一般习惯使用如下代码。

    eset.filter <- nsFilter(eset, require.entrez=F, remove.dupEntrez=F) #不使用函数的合并探针功能
    
    eset.filter$filter.log #查看每一步筛掉多少探针
    
    eset.filter <- eset.filter$eset #只留下ExpressionSet对象
    

    这样就能得到过滤掉低表达探针后的ExpressionSet了。

    如果想要使用nsFilter函数的合并探针功能,需设置require.entrez=TRUE, remove.dupEntrez=TRUE 然后设置var.func 为自己想要用来合并探针的统计参数即可,默认为 IQR。

    2 使用 mas5calls 函数过滤探针

    mas5calls 函数整合在 affy 包中。 mas5calls 函数的算法是 Wilcoxon signed rank-based gene expression presence/absence detection(不翻译反而意思更清楚)。这个算法整合与 mas5 算法中。这个函数的功能是对 AffyBatch 的探针状态进行评估,分类为”P”,”M”和”A”。我们可以只提取其中分类为”P”的探针,认为其他探针不表达;亦或是只删去分类为”A”的探针。由此过滤掉大部分表达量低的探针。

    注意,由于 mas5calls 函数只能针对 AffyBatch 对象,因此只能用于标准化前的数据。而且不能用于oligo包读取的数据。

    以下代码是使用 mas5calls 函数过滤探针的例子,数据来自于 affydata 包的 Dilution 数据。代码的功能是过滤掉在所有样本中均不表达的探针。

    library(affy)
    library(affydata)
    data("Dilution")
    eset <- rma(Dilution)
    calls <- mas5calls(Dilution) # 获得探针的PMA信息
    calls <- exprs(calls) # 提取包含探针的PMA信息的表达矩阵
    head(calls)
    absent <- rowSums(calls == "A") # 统计每个探针在所有样本中分类为"A"的个数
    absent <- which(absent == ncol(calls)) # 提取在所有样本中分类均为"A"的探针,也就是在所有样本中都不表达的探针。
    rmaFiltered <- eset[-absent,] # 删去不表达的探针
    nrow(exprs(eset)) # 过滤前的总探针数目
    nrow(exprs(rmaFiltered)) # 过滤后的总探针数目
    

    稍微更改一下条件。提取在所有样本中均表达的探针。代码如下。

    present <- rowSums(calls == "P") # 统计每个探针在所有样本中分类为"P"的个数
    present <- which(present == ncol(calls)) # 提取在所有样本中分类均为"P"的探针,也就是在所有样本中都表达的探针。
    rmaFiltered <- eset[present,] # 只留下在所有样本中均表达的探针
    

    再更改一下条件,提取至少在一个样本中表达的探针。代码如下。

    present <- rowSums(calls == "P") # 统计每个探针在所有样本中分类为"P"的个数
    present <- which(present != 0) # 提取至少在一个样本中表达的探针。
    rmaFiltered <- eset[present,] # 留下至少在一个样本中表达的探针
    

    大家可以按照自己的需求修改上述代码,达到过滤探针的目的。

    使用 genefilter 函数过滤探针的好处在于能够自定义筛选条件,但多数情况下用不到。而使用 nsFilter 函数或 mas5calls 函数过滤探针的方法虽然固定,但操作方便。大家根据自己的需求使用相对应的函数就可以了。

    3 使用 paCalls 函数过滤探针

    mas5calls函数只能针对AffyBatch对象进行,因此无法用于oligo包读取的。oligo包读取的数据对象如下。可以看到无论哪种类型的对象都不能用于mas5calls函数。

    oligo包提供了paCalls函数用于过滤探针。对于Whole Transcript arrays (Exon/Gene),可选择的method有”DABG” (p-values for
    each probe) 和 “PSDABG” (p-values for each probeset)。对于Expression arrays 来说,只能选”MAS5”。

    使用oligoData包中的affyExpressionFS数据作为例子。

    首先载入数据。

    library(oligo)
    library(oligoData)
    data("affyExpressionFS")
    

    使用oligo包提供的rma函数进行标准化。

    注意,affy包中也有rma函数。但affy包提供的rma函数不能处理oligo包读取的对象。因此若同时加载了affy包,直接使用rma函数可能会报错。可以选择不加载affy包,或是使用oligo::rma()命令使用oligo包中的rma函数。

    > data.eset <- rma(affyExpressionFS)
    Background correcting
    Normalizing
    Calculating Expression
    

    使用paCalls函数。

    calls <- paCalls(affyExpressionFS) #返回calls
    view(calls)
    

    观察calls的结构,发现calls是一个list,包含callsp两个元素。calls是包含探针PMA信息的matrix,而p是包含p值信息的matrix。所谓的p值,是指探针表达量与背景相同的概率,越小说明探针与背景的差异越显著,也就是指探针表达的可能性越大。

    接下来的步骤与方法2相同。详情请看方法2.

    pmas <- calls$calls
    absent <- rowSums(pmas == "A")
    absent <- which(absent == ncol(pmas))
    rmaFiltered <- data.eset[-absent,]
    nrow(data.eset)
    nrow(rmaFiltered)
    

    这里提供另一个有趣的例子。GSE72154 。这个例子中paCalls函数的返回值比较特别。下面将会详细介绍使用paCalls函数过滤该研究。为了省时间和突出重点,芯片质量控制就不做了。

    首先在ftp下载GSE72154_RAW.tar,解压缩到工作目录。然后用oligo包读取数据。

    celfiles <- list.files(path = ".", pattern = ".CEL", full.names = T)
    data.raw <- read.celfiles(celfiles) #读取数据
    sampleNames(data.raw) <- c("ler1","ler2","wrky1","wrky2")
    data.eset <- rma(data.raw) #标准化
    calls <- paCalls(data.raw) #获得calls
    class(calls) #查看calls类型
    head(calls) #查看calls信息
    

    注意到,和上文不同,paCalls函数在这个例子中返回值callsmatrix对象。同时要注意到,探针名称发生了改变。

    data.raw和calls的探针名称一致,而data.eset的探针名称发生了改变。如果要对data.eset进行探针过滤,就势必要构建两种探针名称的对应关系。

    观察calls的结构,发现这个例子中的calls,其实相当于上个例子中的calls$p,也就是包含p值的矩阵。因此,我们可以根据p值来过滤探针。

    P <- apply(calls, 1, function(x) any(x < 0.05)) #选择在所有样本中至少有一个探针p值<0.05的探针名
    pids <- as.numeric(names(P[P])) #将筛选后的转换探针名称转化为numeric对象,因为后续分析提取的转换探针名是numeric格式。
    head(pids)
    

    获得了符合条件”在所有样本中至少有一个探针p值<0.05”的探针名。接下来要将转换探针名称与原始探针名称对应起来。

    pinfo <- getProbeInfo(data.raw) #获取转换探针名称与原始探针名称的对应关系
    head(pinfo)
    fids <- pinfo[pinfo$fid %in% pids, 2] #获取筛选后的原始探针名称。第二列为原始探针名称。
    

    过滤探针,返回ExpressionSet对象。

    data.filter <- data.eset[rownames(data.eset) %in% fids, ] #保留符合筛选条件的探针,返回ExpressionSet对象
    nrow(data.eset) #查看过滤前探针数目
    nrow(data.filter) #查看过滤后探针数目
    

    其实方法3的第一个例子也可以用第2个例子的方法来做,只是比较麻烦而已,所以如果能用第一个方法做就不要用第二个方法。如果都嫌麻烦的话,可以用nsFilter函数来做,一步到位。

    上面的例子过滤前有38408个探针,过滤后有35697个探针。只过滤了10%的探针。可以考虑更加严格的条件,比如apply(calls, 1, function(x) any(x<0.01))(即p < 0.01), 或是 apply(calls, 1, function(x) sum(x < 0.05) > 2)(即至少在3个样本中p < 0.05).

    参考:

    1.用affy包读取affymetix的基因表达芯片数据-CEL格式数据

    2.使用oligo软件包处理芯片数据

    展开全文
  • Bioconductor分析基因芯片数据第五章

    千次阅读 2018-02-21 20:29:29
    使读者初步了解使用Bionconductor完成基因芯片预处理的流程接着详细讲解戏弄...chip芯片)的芯片处理方法5.1快速入门例5-1 从数据包CLL中载入芯片数据,完成预处理,最后获得基因(探针组)表达矩阵。注意,探针组表...

    使读者初步了解使用Bionconductor完成基因芯片预处理的流程
    接着详细讲解戏弄i按预处理和数据分析等内容

    最后深入了解实际工作中会遇到的芯片处理问题以及如何用学到的只是解决问题
    目的:掌握芯片分析的整体框架,自行学习其他厂商或种类(例如SNP芯片或CHIP-chip芯片)的芯片处理方法

    5.1快速入门
    例5-1 从数据包CLL中载入芯片数据,完成预处理,最后获得基因(探针组)表达矩阵。注意,探针组表达矩阵的行对应的行对应的探针组,而不是基因,基因和探针组的关系见5.2.1.这段程序从载入原始数据(CEL文件)开始,通过预处理得到基因表达矩阵,是芯片数据处理的一个必须步骤

    CLL 数据是慢性淋巴白血病(Chronic lymphocytic leukemia,CLL)数据集,采用了Affymetrix公司的HG_U95Av表达噗芯片,测量了24个样本,12625个探针
    例5-1采用的实验设计方式:两组之间是对照试验(control test),每组内都是平行实验(parallel test),对照实验,简单来说就是为了阐明某种单一因素的效应或者影响,在保持其他因素不变的前提下,测试一定数据的实验组样本呢和对照组样本,并对结果进行比较。平行实验,简单来说就是对同样的一组样本取两个以上相同的样品,以完全一致的条件下进行试验,测试结果的稳定性
    5.2基因芯片基础知识
    5.2.1探针组
    一张基因芯片(以affymetrix表达谱芯片为例)可以包含上万个的探针(通常由25个碱基组成),他们整齐有序地印刷在芯片上。一组探针或者探针组(probe set),来自于一个基因,通常由20对或者11对探针组成,每一对探针都由匹配探针(perfect match,PM)和错配探针(Mi是match,MM)组成,成为探针对(probe pair),MM与PM的序列只有正中央的那个碱基不同,其余的都一致。但是,在一些高密度芯片中,例如外显子芯片(Exon array),每个探针组只有4个PM探针,没有MM探针。

    探针序列的来源叫做参考序列,通常来自于公开的核酸数据库(例如NCBIGeneBank或RefSeq)对于不同的芯片类型,探针组在参考序列中的分布不同,3’表达谱芯片的探针组排布在参考序列3‘末端附近的一至两外显子上,外显子芯片中,每个长度大于25个碱基的外显子都有针对他的探针组:铺瓦芯片(Tilling array)中,探针组覆盖了几乎所有的外显子和内含子
    需要强调的是,芯片数据领域提到的基因表达矩阵往往是以探针组而忽视以基因为单位的,即每行都对应一个探针组的表达量。后面将要降到的差异基因分析也是找打显著性差异的表达的探针组,然后通过ID映射才对应到探针组代表的基因,探针组与基于的关系往往是多个探针组对应一个基因。但是在实际应用中,经常不太注意区分,探针组有时也会被叫做基因

    5.2.2主要的芯片文件格式
    主要的是CEL文件
    affymetrixi芯片原始数据最常用格式为CEL文件,也是芯片预处理和分析的出发点。CEL文件的主要内容就是每个“cell"的灰度信息,"cell"是整个芯片图像划分后得到的小网格,每个小网格中的图像被看作来自一个探针,自CellHeader开始,每行数据对应芯片上的一个”cell"位点,包含5列信息,依次为X坐标,Y坐标,灰度的平均值,灰度的标准差以及用了多少个像素来求这个平均值
    CEL文件只提供了每个探针的灰度信息,还需要基因芯片探针排布的信息(即哪个探针来自哪个探针组),才可以得到芯片上每个探针组对于的表达数据,这就需要CDF文件。另一个重要的是probe文件,他提供了探针的序列信息。afftymetrix公司为每种型号的芯片都提供了对应的CDF文件和Probe文件。CDF文件中的对应关系用户可以自行更改,例如为了应多多个探针组 的ID对应到同一基因ID的现象,有些研究季候就把对应到同一个基因的多个探针合并为一个探针组,并提供修改后的CDF和Probe文件
    图5-3B是affymerixHG-U133A芯片的Probe文件的部分内容,他只包括了一个探针组(名称是“200688_at")的所有探针,共11条序列,文件中第2和3列是对应探针所在的X和Y 坐标,第4列是序列的第13个奸计(中心)位置对齐到一致性序列的相对位置,第5列是对应探针的序列,最后是样品与探针杂交的方向


    5.3基因芯片数据预处理
    基因芯片数据预处理的目的是将探针水平的数据(杂交信号)转换成基因表达数据,主要的数据结构有Affybatch类和Expressionset类,前者用于存储探针水平的数据(相当于CEL文件的内容),而后者用于存储表达水平的数据(相当于基因表达矩阵的内容)。预处理通过质量控制,剔除不合格的芯片(数据),只保留合格的进入下一步处理。然后通过标准化,将所有芯片数据中的基因表达只变换到一个可以比较的水平,用于后续分析
    5.3.1数据输入
    例5-1中,芯片数据的输入是从数据包中得到的,但是在实际应用中,更常见的情况是从CEL中获得数据,无论是数据包还是文件输入,读入的数据会存入一个“affybatch"类型的对象中,可以通过执行help(affybatch)获得更详细的介绍
    头文件:用于描述实验样本、平台等相关信息,其中包括phenoData,featureData,protocolData以及annotation等几个类

    assayData:这是affybatch类必不可少的,他的第一个元素是矩阵类型,用于保存基因表达矩阵。该矩阵的行对应不同的探针组(probe sets),用一个无重复的索引值表示,列对应不同的样品。当使用exprs方法时,调取的就是这个基因表达矩阵
    experimentData:一个MIAME类型的数据,设计这个MIAM类的目的就是用于保存MIAME原则建议的注释信息.MIAME原则是一组指导方针,他建议了一组标准来记录与基因芯片实验设计相关的资料


    5.3.2质量控制
    质量控制对于后续的分析至关重要,原始图像(DAT文件)级别的质量控制一般用个芯片公司自带的软件(如affymetrix公司的GCOS)完成。本节中,质量控制主要集中在CEL文件级别的处理,从最简单的直观观察,到平均值方法,再到比较高级的数据拟合方法。这三个层次的质量控制分别功能分别用image函数simpleaffy包和affyPLM包实现
    直观的查看一下芯片上所有位点的灰度图像
    image函数表示选取的CLLbatch中的第一个基因芯片(即“CLL10.CEL")的数据,然后调用image函数根据CEL文件中的灰度信息画图,affymetrix芯片在印刷时会在四个角印刷特俗的花纹,并且在左上角印刷芯片的名称,花纹与芯片名称可以帮助我们借助这个图像分辨率来了解芯片数据是否可靠。如果无法分辨四角花纹或芯片名称,很可能数据有问题
    根据image函数的图像信息,可以对芯片的信号强度产生一个总体认识:如果图像特别黑,说明信号强度低;如果图像特别亮,说明信号强度很可能过饱和
    尺度因子affymetrix公司规定,用于比较的芯片之间的尺度因子的比例必须小于三
    检测值(detection call)和检出率(percent present):一组探针能否被检测到,用检测值有(present,简称R)、无(Absent,简称A)和不确定(Marginal presen,简称M)来表示检测范围的上下边界(a1及a2)选用了默认值0.04和0.06.检出率,是用所有检测值为p的探针数量除以芯片所有探针组数控得出的百分比。如果检出率过低,表示大部分的基因都未被检测到,很难说明是该芯片实验有问题,还是这个样品的大多数基因本身就很难检测到,有原因是表达量基地或是其他。因此,需要看多个样品之间的相对差别,如果有的样品的检出率与其他的有比较大的差别,那很可能该样品出现了问题
    平均背景噪声(average background):对于每一块芯片,根据所有的MM值作出统计,可以得到背景噪声的平均值、最小值和最大值。往往较高的背景噪声都伴随着最低的检出率,因此这两个指标可以结合使用
    标准内参(internal control genes):mRNA是按照5‘端到3’端的顺序来降解的,芯片探针组也是根据这个顺序来设计的,因此探针组的测量结果可以体现这一趋势。因为大部分的细胞都有β-action和GAPDH基因,所有affymetrix在大部分的芯片里都将他们设置为一组观察RNA降解成都的内参基因。根据这两个基于设计的探针组很好的涵盖了他们3‘端和5’端的每一个区段。通过比较他们3‘端相对于中间或者5’端的信号强度,可以很好地指示出实验质量。affymetrix建议这个比值对于β-action不大于3,对于GAPDH不大于1.25,即可以说明这个芯片的质量可以接受。如果这个比值很高,表明不完整的β-action或者GAPDH的存在,可能源于体外转录不好或者降解非常严重。如果使用的是affymetrix的小样本实验流程(small sample protocol)而不是常用的标准流程(standard protocol),建议使用3’端相对于中间的比值。原因是小样本流程有更扩增次数,有可能产生更多较短的转录序列,不可避免的带来3‘端的偏倚。为了验证杂交的质量
    根据上述标准,可以使用Bionconductor的simpleaffy包对affymetrix芯片数据进行质量评估,最后得到质量控制总览图(图5-8)

    qc图的看法,图5-8是CLL数据集中全部24个芯片数据的质量控制总览图。图5-8中从左至右,第一列是所有样品的名称;第2列是两个数字(对应每个样品),上面是以百分比形式出现的检出率,下面的数字表明平均背景噪音;第3列("QQ stats")最下面的横轴是尺度因子等指标对应的坐标,取值范围从-3到3,用浅蓝色虚线作为边界。第3列用到了三项指标:尺度因子、GAPDH3'/5'比值和action3比值(记做graph3/graph5和action3/action5),分别用实心圆、空心圆和三角标志表示出来。另外,如果第三列中出现红色的”bioB"字样,说明该样品中未能检测到BioB
    简单地讲,所有指标出现蓝色表示正常,红色表示可能存在质量问题。但是根据实际情况不同,还要进一步分析。一般来讲,如果有一个芯片各项指标都不太正常,尤其是BioB无法检测到,建议判定为该芯片实验失败。如图5-3中的样品”CLL15.CEL",这个数据的检出率(38.89%)明显低于其他样品,action3/action5远大于3,而且没有检测到BioB,因此可以判定此数据无效。如果多个芯片都出现了相同的问题,原因则可能是多方面的;如左侧第2列24个芯片的检出率和背景噪声都很高,原因是阈值设定过高,如果降低阈值,大部分就会变蓝;再如,全部芯片都不能检测到BioB,有可能是嵌入探针所针对的DNA溶液加入比例不对
    基于平均值家建设的评价指标都有一个,默认的假设,那就是对于每一块新片,质量是均匀的,不会随着位置变化发生较大的变化。但如果关注芯片的每个小格(Grid),就会发现格与格之间的质量也是有差异的,这可能由于芯片印刷的问题,也可能是杂交过程中出现的问题。那么如何才能得到比较可靠的质量评估,这需要设计多种能反映芯片数据全貌的指标综合分析从而得出最终的结论。这些指标要在对原始数据拟合(回归)的基础上计算得到,然后以图的形式显示,包括:权重(weights)&(residuals)图、相对对数表达(relative log expression,RLE)箱线图、相对标准差(normalized unscaled standard errors,NUSE)箱线图、RNA降解曲线、聚类分析(cluster analysis)图、主成分分析(principal component analysis,PCA)图、信号强度分布图及MA图等,以上功能由Bionconductor中的affyPLM包实现
    一般情况下,在权重图中,绿色代表较低的权重(接近0),白色、灰色代表较高的权重(接近1);在残差图中,红色代表正的高残差,蓝色代表负残差;在残差符号中,红色代表正的残差,蓝色代表负的残差。如果权重和残差都是随机分布的,应该看到绿色均匀分布的权重图和红蓝均匀分布的残差图。图5-9中,左上为原始图像,右上为权重图,左下为残差图,右下为残差符号图。另外,还可以看到,图中左上部出现了一些白色的条块,这是正常的现象,因为有些时候,探针会按照GC比率(GC ratio)排布从而导致白斑的,那什么样的权重和残差图是不可接受的呢
    在对比实验中,即使是相互比较的对照组与实验组之间,大部分基因的表达量还是应该保持一致的,平行实验之间一致性更强。相对对数表达(RLE)箱线图可以反映上述趋势,它定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。一个样品的所有探针组的RLE的分布可以用一个统计学中常用的箱型图形表示。如果使用RLE箱线图来控制CLL数据集的实验质量,每个样品的中心应该非常接近纵坐标0的位置(图5-11)。如果个别样品的表现与其他样品的表现与其他大多数明显不同,那说明这样品有问题

    NUSE是一种比RLE更为敏感的质量检测手段。如果根RLE箱线图对某个芯片的质量产生怀疑,那么再结合NUSE图,这种怀疑就可以确定下来。NUSE定义为一个探针组在某个样品的PM值的标准差除以该探针组在各样品中PM值标准差的中位数,如果所有芯片的质量就是非常可靠的话,那么他们的标准差会十分接近,因此他们的NUSE值会都在1附近。然而,如果有某些芯片质量有问题的话,就会严重地偏离1,进而导致 其他芯片的NUSE值偏向相反的方向。当然,还有一中非常极端的情况,那就是大部分芯片有质量问题,但是他们的标准差却比较接近,反而会显得没有质量问题的NUSE值明显偏离1,所以必须结合RLE及NUSE两个图才能作出更可靠的判断。例如结合图5-11和6-12,可以看出CLL1和CLL10的质量明显有其他样品,所以需要舍弃。
    RNA降解是影响芯片数据质量的一个很重要因素,因为RNA是从5‘端开始降解的,所以理论上探针5’端的荧光强度应该低于3‘端的荧光强度。RNA降解曲线的斜率表示了这种变化趋势,斜率越小,说明降解较少;反之,则降解越多。但是,如果斜率太小,甚至接近0,就要特别注意,这不仅不代表基本没降解,而且可能全部被降解。因为,在实际实验中国,基本没降解是不可能的,很可能是因为RNA降解太严重,才导致计算值接近,从图5-13中,可以看出CLL13对应的曲线几乎平行于横轴,因此判断很可能降解严重,需要作为坏数据去除
    最后经过上面的综合分析,需要去除的三个样品数据:CLL1、CLL10和CLL13
    前面讲到的几种质量控制放大都是基于“平均值”思想的。其实,还可以从另外一个角度来对芯片质量进行检验。这就是利用芯片之间的相互关系,例如在对照试验中,理论上组内同种类型的芯片数据应该聚拢在一起,两个组之间应该明显地分离。这个思想是非常合理的,需要做的就是找到一种指标来刻画芯片数据之间的相似度或距离,Pearson线性相关系数就是最常用的这类指标。基于“相互关系”的方法,其核心是相关系数矩阵,它包括了全部关系信息。计算相关系数矩阵,苦役使用预处理子琪娜的芯片数据,也可以使用标准化之后的数据(见例5-8)。例5-8中,通过查看相关系数矩阵“pearson_cor”,可以看到组内(稳定组和恶化组)和组件相似度差异不大。在实际应用中,往往不是直接查看相关关系矩阵,而是根据有相关系数矩阵,而是根据由相关系数矩阵导出的距离矩阵,进行聚类分析或主成分分析以对样品归类并图形化显示(见例5-8)
    从聚类分析的整体结果看(上图),稳定组和恶化组根本不能很好地分开,这样还不能简单判定实验完全失败,所有样品数据都不能用。理论上讲,如果总体上两组数据是分开的,那么说明我们关心的导致癌症从稳定到恶化的因素起主导作用;如果不是,很可能其他因素起主导昨天,因此导致聚类被整体打乱,则不能简单判定所有样品出了问题。芯片分析往往采用两个主成分来构建分类图,从图5-15也可以看出稳定组(矩形)和恶化组(菱形)根本就不能很好分开。使用主成分分析时,还必须考虑前2个主成分是否具有代表性,这要看前2个主成分的累计贡献率,如果低于60%,可以考虑另外一种类似的方法来构建分类图,即多维尺度分析(metric multi-dimensional scaling method)

    5.5.3背景矫正、标准化和汇总
    芯片数据通过质量控制,剔除不合格的样品,留下的样品往往需要通过散步处理(背景矫正、标准化和汇总)才能得到下一步分析所需要的基因表达矩阵
    首先,讲一下背景矫正,前面提到的芯片中MM探针的作用是检测非特异性杂交信号。理论上,MM只有非特异性杂交,而不会有特异性杂交,MM的信号值永远小于其对应的PM信号值,那么可以用简单的书学方法处理一下,做一个PM-MM或者PM/MM即可去除背景噪声的影响。但实际中,经常发现大量的MM信号值比PM信号值还要高。因此,需要应用更为复杂的统计模型来去除背景噪声,这个过程叫做背景矫正
    其次,介绍一下标准化。标准化的目的是使各组/次测量或各种实验条件下的测量可以相互比较,相处测量间的非实验差异,非实验差异可能来源于样品制备、杂交过程或杂交信号处理等,芯片数据标准化,根据其基本假设总体上分为两种:“bulk normalization”和“control-based normalization”。前者假定仅有一小部分基因表达值在不同条件下有差异,而绝大部分基于表达值不变,因此使用所有的基因表达值作为参考进行标准化;而后者使用表达值被认为是恒定不变的参考基因(通常为芯片制造商提供的外源参考基因)作为标准进行标准化。在实际应用中,芯片数据标准化只采用第一种方法
    ’最后,使用一定的统计方法将前面得到的荧光强度从探针(probe)水平汇总到探针组(probe set)水平,这个过程被称为汇总(summarization)
    上述散步处理过程可由一个函数实现,它就是affy软件包中的expresso函数,通过控制这个函数的参数,就可以分别制定三步处理具体应该采用的方法
    expresso参数复杂,可以通过help(expresso)命令获得他的全部参数说明


    芯片内标准化方法针对双通道(见2.3.1),又可分为全句话方法(global normalization)和荧光强度依赖的方法(intensity-depent normalization),前一种方法假设红色染料的信号强度是正比例关系的,即R=kG(R:红色信号强度;G:绿色信号强度,k:假设为常数)。差异表达值(log2(R/G)在标准化之后相当于平移了一个常量c=log2(k),数学上表示为log2(R/G)-c=log2(R/kG)=0。但实际上,c并不是一个常数,而是另外一个变量的A的倍数c(A),这里A=1/2*log(R/G),这一点可以从MA图(图5-16A)中看到M的总趋势不是平行于X轴的
    MA(M代表Minus,A代表Average)图的英文全称是:The distribution of the red/green intensity ratio plotted by the average intensity .MA图中,定义M=log2(R/G)


    5.3.4 预处理的一体化方法
    前面5.3.3讲到了通过设置参数,expresso函数可以自动化实现整个预处理过程(背景矫正、标准化和汇总)。除了expresso函数,affyPLM软件包提供了three step函数可以更快的实现同样的功能
    例5-11
    从信号强度分布图来看(图5-17),MASS算法处理后的数据出现了很多负数,从图5-17中还可以看出原本不重合的多条分布曲线(图5-17A)在经过了RMA算法处理后重合到了一起(图5-17C),有利于下一步的差异表达分你想。但是他却出现了两个峰值,这并不符合高斯正太分布。如果采用gcRMA算法处理,不但所有的曲线很好地重合到了一起,而且他们的分布也更加近似高斯分布(图5-17D)。因此,gcRMA算法对RMA算法的改进在这一组数据上表现得很明显。然而,这并不意味着gcRMA算法总是优于RMA算法,对于不同的数据进行算法比较,才能进一步确定哪种算法最合适
    通过箱线图(图5-18)可以看到三种算法处理后的个样品的中值十分接近。MASS算法总体而言还是不错的,只是有一定的拖尾现象。而gcRMA的拖尾现象比RMA要明显的多。这说明针对地表达量的基因,RMA算法比gcRMA算法表现的更好
    还可以通过MA图(图5-19)来查看标准化处理的效果,从例5-12中(只示例CLL中一部分数据)可以看出:在原始数据中,中值(红色曲线)偏离0,经过gcRMA预处理之后,中值基本保持在零线上。注意,运行例5-12最后一行代码时,MAplot函数不支持ExpressionSet类型的数据CLLgcrma,读者可以将其转换为Affybatch类型后再运行
    例5-12



    5.4基因芯片数据分析
    本书2.2.3提到了基因芯片表达差异的显著性分析在基因表达数据分析中的特殊地位,而且这个地位很大程度上都是基于芯片领域的经验得来的。尽管研究人员不断改进芯片试验和统计学方法,并不断寻求一些新的方法(例如机器学习)来分析芯片数据,当前最主要的应用依然还是基因表达差异的显著性分析。本书从例5-13到5-17的程序涵盖了一个显著性分析的完整流程,读者可以一次运行全部代码,也可分开运行以便于逐个掌握,但是必须连续运行所有程序,因为后面的程序依赖前面程序的输出。通过这几个实例,读者可以清晰地把握Bioconductor处理芯片数据的整个流程

    5.4.1选取差异表达基因
    基因表达差异的显著性分析的第一步就是选取表达具有显著性差异的基因。总体来说,这类分析的基本假设是标准化的芯片数据符合正太分布,因此所用的统计方法基本上就是T/F检验和方差分析。当前,常用的分析方法主要有:T检验、SAM(significance analysis of microarrays)方法、CyberT方法、经验贝叶斯(Empirical Bayes)方法、方差分析(The Analysis of avriance,anova)和RP(Rank produces)方法
    RP方法通过计算基因表达值的集合平均值及其排序的变化来比较两组间的差异。SAM、CyberT和经验贝叶斯都是调整后的T检验,而且后两种方法都采用了贝叶斯方法进行调整。CyberT将标准差及信号强度的关系使用线性模型进一步强化,提高了准确率,有研究指出,它的计算结果要好于SAM算法。经验贝叶斯又在CyberT基础上进行了改进,首先,经验贝叶斯在计算标准差时考虑的全部基因,而不是排序后相近的(人为设定的同一个窗口范围内)基因;其次,经验贝叶斯不在局限于两组数据,可以通过设计实验对比矩阵,计算多种复杂条件下的差异表达。因此,经验贝叶斯是当前最为常用的分析方法,他已经完整地由Bioconductor中的limma包实现。但是,总体来睡哦,现在没有仍和理论或者经验能够证明哪种算法是最好的
    limma包是基于R和Bioconductor平台的分析芯片数据的综合软件包,例5-13是应用limma包计算CLL数据集中差异表达基因的整个流程


    首先,可以从最终结果(即变量“dif")中查看所有的两组数据(即恶化期与稳定期)之间差异表达基因的信息。每行数据对应一个探针组,包括8列信息:第1列是探针组在基因表达矩阵eset中的行号;第2列“ID”是探针组的AffymetrixID;第3列“log FC”是两组表达值间已以2为底对数化的变话倍数(Fold change,FC),注意由于基因表达矩阵eset本身已经取得了对数值,因此这里实际上只是两组基因表达值均值之差;第4列“AveExpr"是该探针组在所有样品中的平均表达值(average exprssion value);第5列”t"贝叶斯得到的调整后的俩组表达值T检验中的t值;第6列“P.Value"是贝叶斯经验得到的P值;第7列“adj.P.Value"是调整后的p值;第8列”B“是经验贝叶斯得到的标准差的对数化值,由于涉及较深的数学基础,为了加深limma的计算过程,可以用简答函数来得到探针组”39400_at"的行号、“AveExpr"和“log FC“
    然后逐次介绍这个分析过程的六个关键步骤:构建基因表达矩阵、构建实验设计矩阵、构建对比模型(也叫对比矩阵)、线性模型拟合、贝叶斯检验和生成结果报表(dif)
    构建基因表达矩阵时,需要注意的是,limma对输入数据的要求是必须是经过对数转换的表达值。例5-13调用了gcRNA算法来对数据进行预处理,得到标准化的基因表达矩阵eset,这个矩阵是经过对数变转换的。但是,如果是从其他算法(例如MASS)得到的数据,还需要自行编程进行对数转换
    实验设计矩阵需要调用model.matrix函数构建,该函数需要用户指定一个公式,构建好的实验设计矩阵design要提供给下一步的拟合函数lmFit。通过查看design变量,可以看到下面内容:实验设计矩阵的每一行对应一个样品的编号,每一列对应样品的一个特征,每个特征实际上形成了一个包含若干表达基因。比如,在例5-13中,一共有21个样品,他只考虑了一个因素,即疾病状态(disease),这个因素有两个水平,即恶化(progressive)和稳定(stable),最后实验矩阵中出现了diseaseprogres和diseasestable两个特征。多因素和多水平的实验设计,会产生更多的特征
    比较模型需要调用makecontrasts函数构建,该函数需要用户制定一个公式,这个公式表明用户对实验矩阵design中的哪一个特征和哪一列特征进行比较,以得到差异。例5-13指定的是在恶化和稳定两个水平之间进行比较,以寻找这两个水平之间的差异表达给予你,因此,公式表示为contrasts='diseaseprogress. - diseasetable",zhuyi "diseaseprogress."中的“.",是CLL数据集中对:progress”简写带来的,不是运算符号
    接下来是根据实验设计矩阵条用函数对基因表达矩阵做线性拟合lmFit(eset,design),根据对比模型进行插值计算,最后是贝叶斯检验(见5.4.1)。由于这些射击较深的统计学背景
    最后,重点讲一下topTbale函数,它的主要功能有三项:1.对贝叶斯检验得到的“P.Value"进行调整得到”adj.P.Value",调整的算法默认是BH(Benjamini-hochberg)算法2.生成全部基因的检验结果报表3.还可以通过某个参数来筛选具有显著性差异表达的基因,通常使用“adj.P.Value,常用的阈值一般是0.05或者0.01,也可以使用”P.Value"(见例5-13)这里有三点需要注意:1.topTable提供了多种方法可以做基因筛选,例5-13就通过对数化的变化倍数“lfc"去掉了一些在两组条件下变化不大的基因,但是这样做的理由并不是很充分,因为变化倍数不大的不一定就是没有显著变化2.topTable还提供了参数可以对基因进行排序,比如使用“adj.P.Value"从小到大排序,可以清楚地看到变化最显著的基因3.显著基因的选取具有一定的主观性,阈值设定是0.01还是0.05并没有严格的规定


    5.4.2注释
    找到了差异表达基因,接下来是使用注释包对差异表达基因进行注释。在4.2.3中的注释一部分讲解中提到过Bionconductor的几种注释方式,对affymetrix芯片产生的差异表达基因的注释就采用第一类注释方式,即下载对应具体平台的注释包,进行本地注释(这部分代码需要在例5-13后执行)。例5-14只用两种基因ID 来对探针组进行注释,有关用基因本体路(GO)和通路(pathway)注释的内容与5.4.3的GO和通路富集分析一起讲解


    例5-14注释实质上就是一个ID映射的过程(见2.3.1),也就是把芯片探针组的ID映射到基因国际标准名称(gene symbol)和Gene symbol是由人类基因命名委员会(the hugo gene nomenclature committee,hgnc)为每个人类基因提供的唯一命名,一般是大写拉丁字母缩写形式,后面可加数字。NCBI对于每一条提交的序列,根据其存入的NCBI数据库时的先后顺序赋给一个整数,这就是GI。这里增加了一列GI的目的,就是为了下一步通过GI映射到基因本体论(GO),然后做GO 的富集分析
    “GO:0022904"


    5.4.3统计分析及可视化
    差异基因注释后的下一步工作就是统计分析和可视化(见2.3)。对于差异表达分析,最主要的两种统计分析就是GO的富集分析(见2.4.4)和KEGG通路的富集分析(见2.4.5)。这两种分析方法分别由Bioconductor的Gostats包(见例5-15)和geneanswer包(见例5-16)实现


    从例5-15最终结果(即变量"bp")可以看到每个显著性富集的GOterm含有六列信息(不包括行号):第1列是GO term的ID,该ID对应的内容在后面列出,如”GO:0022900",对应后面的“respiratory electron transport chain";第2列”p.value"是超几何检验的p值;第3列“oddsratio”是超几何分布中的比值;第5列“count”是差异表达基因中世纪属于这个GOterm的基因数量;低6列“size"是总基因中属于这个GO term的基因数量。以“GO:0022904"为例,此次分析的总基因数量为8804,差异表达基因数量是1138804个基因中有75个基因(即”size")属于”GO:0022904",如果从8804个基因中随机抽取113个基因,那么113个基因中期望属于“GO:0022904"的基因数量应该是2.25(即”ExpCount"),而实际上12个(即‘Count"),根据这个情况,计算出来的P值应该是1.506871e-10(远远小于0.01),因此可以说差异显著基因在“GO:0022904"上是显著富集的。为了加深理解GO富集分析的计算过程,读者可以用简单函数来计算P值。另外5-15还通过函数htmlReport输出了HTML的报告文件,它在前面六列的基础上,多加了一列GO term的描述,并且链接到GO的官方网站
    值得注意的是,对比例5-15和5-14的结果报表,可以看到例5-15的报表bp没有根据p值来筛选统计上显著富集的GO term,因此包括了全部的GO term

    例5-16调用了GeneAnswers包实现了KEGG通路的注释、统计和可视化的给你。而且GeneAnswers功能强大,除了KEGG,还可以支持GO\REACTOME和CABIO等多个数据库,可以通过设定参数categoryType分别指定注释类型。从例5-16最终结果可以看到每个显著性富集的通路;第2列”precent in the observed List"表示在观察到的基因列表中的比例;第3列“percent in他和genome”是在基因组中的比例;第4列“fold of overrepresents"是基因过表达的倍数;第5列”OddRatio"是超几何分布中的比值比;第6列“P.Value"是超几何检验的P值
    可视化可以直观显示统计结果,帮助研究人员进一步理解实验结果并找到下一步工作的思路,因此可视化和统计分析密不可分。Bioncoductor的所有统计分析包几乎提供了相应的函数来显示数据分析结果。这里根据前面的分析结果,调用pheatmap包来绘制差异表达热谱(图5-20),调用Rgraphviz包来绘制显著富集的GO term的关系图;最后绘制显著富集的KEGG 通路的关系图和热图

    实例3例5-24主要完成找到对比1和4之间、对比2和3之间共同表达的差异表达基因,对比5的差异表达基因,并对三组差异基因座注释和GO富集分析。每组数据输出三个HTML格式的报告文件,分别对应GO三个领域的富集分析的结果,该文件内容请看例5-15有详细介绍。由于pathway分析,特别是pathway的显示要占用较大内存资源,必须在Linux服务器上运行



    展开全文
  • 使用oligo软件包处理芯片数据

    万次阅读 2017-03-27 18:15:10
    Affy芯片的处理方法 ,其中所使用的软件包有一定的局限性,无法读取和分析一些新版Affy芯片。本文介绍oligo软件包的处理方法以解决这些问题。oligo软件包并不是新出现的软件包,只因新类型芯片的不断推出,关注它的...
  • 常用的基因表达数据来自基因芯片或高通量测序。虽然矩阵看起来差不多,但是由于服从不同的分布,因此在进行差异表达的时候需要用不同的方法。对于一般的生命科学领域科研人员来说,了解晦涩的算法并没有太大价值。...
  • W5500 是一款全硬件 TCP/IP 嵌入式以太网控制器,为嵌入式系统提供了更加简易的互联网连接方案.W5500 集成了 TCP/IP 协议栈,10/100M 以太网数据链路层(MAC) 及物理层(PHY),使得用户使用单芯片就能够在他们的应用...
  • 如何利用廉价的51单片机来控制网卡芯片进行数据传输,加载TCP/IP协议连接到互联网,实现网络通信成了众多设计者的目标。但由于指令及资源的限制,实施过程会有许多困难。我们在设计方案中舍弃了耗费资源的高级协议,...
  • BCM芯片FP原理及相关SDK数据结构介绍

    万次阅读 2014-02-26 21:20:11
    BCM芯片有几个大的模块: VLAN、L2、L3和FP等几个,其中FP的使用也最为灵活,能解析匹配数据包文的前128字节比特级的内容,动作包括转发、丢弃、结合qos修改相应字段、分配vid、流镜像、流重定向、指定端口转发...
  • TYPEC-CC逻辑芯片-E-MARK数据线-浅析

    千次阅读 2020-07-10 00:18:36
    所以硬件工程师的位置上宁愿有很多的线缆或者USB-HUB,不为别的,就为了电脑一开机只要连一个HUB线,其它的USB数据线、U盘、USB设备等等,都统一接上,不需要再逐个去连接,不然一个个连(而且又不能一下插上的话)...
  • W7100A iMCU是一个单片机以太网嵌入式控制芯片,它的的结构是:内嵌8051单片机 + TCP/IP协议栈 + 10/100 高速以太网络MAC/PHY W7100A 是 W7100 的升级版。它增加了一些新功能,如记忆锁定功能(Memory Lock)、休眠...
  • STM32+MFRC522完成IC卡号读取、密码修改数据读写

    千次阅读 多人点赞 2021-05-20 14:36:48
    使用MFRC522模块完成对IC卡卡号读取、卡类型区分、IC卡扇区密码修改、扇区数据读写等功能;底层采用SPI模拟时序,可以很方便的移植到其他设备,完成项目开发。 现在很多嵌入式方向的毕业设计经常使用到该模块,比如:...
  • 芯片复位

    千次阅读 2014-12-19 10:16:17
     机器在为芯片上电或程序运行出错是,内部是随机的混乱状态,地址指针不正确,使得程序的获取不了预期的数据。 3.冷复位于热复位  区别在于.冷复位在复位时有加电,而热复位在复位时没有加电。加电的后果使得冷...
  • 芯片质量分析芯片数据预处理获取差异表达基因GO和KEGG分析聚类分析 (本文于2013.09.04更新) TAIR,NASCarray 和 EBI 都有一些公开的免费芯片数据可以下载。本专题使用的数据(Exp350)来自NASCarray,也可以...
  • 本文主要讨论在更换了SDRAM芯片后,初始化代码中内存相关参数应该如何修改。这里以ISSI的IS42S16400在“龙芯1c库”中的配置为例(pmon中类似),参考SDRAM芯片手册中,修改SDRAM相关参数。龙芯1c库是把龙芯1c的常用...
  • 续代码走读(七)的步骤,抓取到的数据的下载方式在文末。...代码其实只加了两个赋值,对程序没有影响,修改的目的一是为了方便实时看数据(毕竟在memory里看数据不是那么直观);二是为了在末尾处加断点。
  • 人工智能芯片研究报告

    千次阅读 2018-12-14 23:30:27
    2010年以来,由于大数据产业的发展,数据量呈现爆炸性增长态势,而传统的计算架构又无法支撑深度学习的大规模并行计算需求,于是研究界对 AI 芯片进行了新一轮的技术研发与应用研究。AI 芯片是人工智能时代的技术...
  • 芯片设计流程 芯片的设计原理图

    万次阅读 多人点赞 2019-05-04 10:09:50
    本文探讨的就是芯片在字面以外的意义,以及芯片是怎么被设计成的。 芯片 芯片,又称微电路(microcircuit)、微芯片(microchip)、集成电路(英语:integrated circuit, IC)。是指内含集成电路的硅片,体积很小...
  • DSP芯片介绍

    千次阅读 2006-03-18 14:36:00
    签于最近汉芯造假丑闻非常流行,这里介绍一下dsp芯片借此机会,给大家一个dsp芯片的概念...DSP芯片的内部采用程序和数据分开的哈佛结构,具有专门的硬件乘法器,广泛采用流水线操作,提供特殊的DSP 指令,可以用来快速
  • 将其嵌入到水、电、气、暖智能(卡)表、机顶盒、智能电器或其它专用设备中,可以完成数据的加密解密、双向身份认证、访问权限控制、通信线路保护、临时密钥导出、数据文件存储等多种功能。ESAM模块使用的安全芯片...
  • 蓝牙芯片

    千次阅读 2014-03-14 11:14:46
    基本上,读过datasheet都应该知道,CSR拥有或者曾经拥有如下三种芯片,HCI ROM,HCI RFCOMM ROM,Full embeded solution. 这三种芯片的区别在于支持的蓝牙协议栈的层次不同。  HCI ROM: 这种芯片支持到HCI 接口...
  • 芯片手册学习总结

    2019-09-29 16:02:07
    芯片手册学习,总结
  • 摄像机DSP芯片介绍

    千次阅读 2019-04-25 13:43:38
    DSP 芯片的内部采用程序和数据分开的哈佛结构,具有专门的硬件乘法器,可以用来快速的实现各种数字信号处理算法。 在当今的数字化时代背景下, DSP 己成为通信、计算机、消费类电子产品等领域的基础器件。 什么是...
  • 目前国内做AI芯片的公司可谓不少了,AI芯片已然成为了当下芯片行业最热领域。但是大部分人对AI芯片的架构应该都不是太了解。那么AI 芯片和传统芯片有何区别?AI芯片的架构到底是怎么样的? (更多相关资讯搜索圈T...
  • 写这篇博客的主要原因是因为公司的产品涉及到电池充放电管理,而且充电电压和电池电压可能会有多种组合,针对这种设计需求,发现目前流行的PD快充协议正好是多级电压,所以在TI支持PD快充的芯片中选择支持多节电池...
  • 自动驾驶芯片之——FPGA和ASIC介绍

    千次阅读 2019-01-10 08:08:25
    当前阶段,GPU 配合 CPU 仍然是 AI 芯片的主流,而后随着视觉、语音、深度学习的算法在 FPGA以及 ASIC芯片上的不断优化,此两者也将逐步占有更多的市场份额,从而与GPU达成长期共存的局面。从长远看,人工智能类脑...
  • MTK芯片

    千次阅读 2008-07-19 09:40:00
    目前联发科技已开发出MT6205、MT6217、MT6218、MT6219、MT6226、MT6227、MT6228等系列平台,其中 MT6205、MT6217、MT6218、MT6219、MT6226、MT6227、MT6228均为基带芯片,所有芯片均采用ARM7的核。MT6305、MT6305B为...
  • TI针对DM6467提供的UBOOT和内核默认都是串口0作为调试串口输出的,但现在我需要使用DM6467的UART0的modem功能,所以修改代码,改变调试串口为串口2。 需要修改的主要有几部分内容: 1. UBL 代码(这部分代码在刚上

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 79,618
精华内容 31,847
关键字:

如何修改芯片数据