精华内容
下载资源
问答
  • 差异表达分析之FDR

    千次阅读 2020-03-20 18:54:49
    差异表达分析之FDR 随着测序成本的不断降低,转录组测序分析已逐渐成为一种很常用的分析手段。但对于转录组分析当中的一些概念,很多人还不是很清楚。今天,小编就来谈谈在转录组分析中,经常会遇到的一个概念FDR,...

    差异表达分析之FDR
    随着测序成本的不断降低,转录组测序分析已逐渐成为一种很常用的分析手段。但对于转录组分析当中的一些概念,很多人还不是很清楚。今天,小编就来谈谈在转录组分析中,经常会遇到的一个概念FDR,那什么是FDR?为什么要用FDR呢?一起来学习吧!

    什么是FDR

    FDR (false discovery rate),中文一般译作错误发现率。在转录组分析中,主要用在差异表达基因的分析中,控制最终分析结果中,假阳性结果的比例。

    为什么要用FDR

    在转录组分析中,如何确定某个转录本在不同的样品中表达量是否有差异是分析的核心内容之一。一般来说,我们认为,不同样品中,表达量差异在两倍以上的转录本,是具有表达差异的转录本。为了判断两个样品之间的表达量差异究竟是由于各种误差导致的还是本质差异,我们需要根据所有基因在这两个样本中的表达量数据进行假设检验。常用的假设检验方法有t-检验、卡方检验等。很多刚接触转录组分析的人可能会有这样一个疑问,一个转录本是不是差异表达,做完假设检验看P-value不就可以了么?为什么会有FDR这样一个新的概念出现?这是因为转录组分析并不是针对一个或几个转录本进行分析,转录组分析的是一个样品中所转录表达的所有转录本。所以,一个样品当中有多少转录本,就需要对多少转录本进行假设检验。这会导致一个很严重的问题,在单次假设检验中较低的假阳性比例会累积到一个非常惊人的程度。举个不太严谨的例子。

    假设现在有这样一个项目:

    ● 包含两个样品,共得到10000条转录本的表达量数据,

    ● 其中有100条转录本的表达量在两个样品中是有差异的。

    ● 针对单个基因的差异表达分析有1%的假阳性。

    由于存在1%假阳性的结果,在我们分析完这10000个基因后,我们会得到100个假阳性导致的错误结果,加上100条真实存在的结果,共计200个结果。在这个例子中,一次分析得到的200个差异表达基因中,有50%都是假阳性导致的错误结果,这显然是不可接受的。为了解决这个问题,FDR这个概念被引入,以控制最终得到的分析结果中假阳性的比例。

    如何计算FDR

    FDR的计算是根据假设检验的P-value进行校正而得到的。一般来说,FDR的计算采用Benjamini-Hochberg方法(简称BH法),计算方法如下:

    1. 将所有P-value升序排列.P-value记为P,P-value的序号记为i,P-value的总数记为m

    2. FDR(i)=P(i)*m/i

    3. 根据i的取值从大到小,依次执行FDR(i)=min{FDR(i),FDR(i+1)}

    注:实际上,BH法的原始算法是找到一个最大的i,满足P≤i/m*FDR阈值,此时,所有小于i的数据就都可以认为是显著的。在实践中,为了能够在比较方便的用不同的FDR阈值对数据进行分析,采用了步骤3里的方法。这个方法可以保证,不论FDR阈值选择多少,都可以直接根据FDR的数值来直接找到所有显著的数据。

    下面我们以一个包含10个数据的例子来看一下FDR计算的过程差异表达分析之FDR

    在这个例子中,第一列是原始的P-value,第二列是排序后的序号,第三列是根据P-value校正得到的初始FDR,第四列是最终用于筛选数据的FDR数值。如果我们设定FDR<0.05,那么绿色高亮的两个数据就是最终分析认为显著的数据。

    FDR的阈值选择在转录组分析中是非常重要的一个环节,常用的阈值包括0.01、0.05、0.1等。实践中也可以根据实际的需要来灵活选择。例如,在做真菌或者原核生物的转录组分析时,由于这些物种转录本数量较少,假阳性累积的程度较低,所以可以适当将FDR阈值设置的较高一些,这样可以获得较多的差异表达结果,有利于后续的分析。

    浅谈多重检验校正FDR
    Posted: 四月 12, 2017 Under: Basic By Kai no Comments

    例如,在我们对鉴定到的差异蛋白做GO功能注释后,通常会计算一个p值。当某个蛋白的p值小于0.05(5%)时,我们通常认为这个蛋白在两个样本中的表达是有差异的。但是仍旧有5%的概率,这个蛋白并不是差异蛋白。那么我们就错误地否认了原假设(在两个样本中没有差异表达),导致了假阳性的产生(犯错的概率为5%)。

    如果检验一次,犯错的概率是5%;检测10000次,犯错的次数就是500次,即额外多出了500次差异的结论(即使实际没有差异)。 为了控制假阳性的次数,于是我们需要对p值进行多重检验校正,提高阈值。

    第一种方法Bonferroni,最简单严厉的方法。
    例如,如果检验1000次,我们就讲阈值设定为5% / 1000 = 0.00005;即使检验1000次,犯错误的概率还是保持在N×1000 = 5%。最终使得预期犯错误的次数不到1次,抹杀了一切假阳性的概率。但是该方法虽然简单,但是检验过于严格,导致最后找不到显著表达的蛋白(假阴性)。

    第二种方法FDR(False Discovery Rate)
    相对Bonferroni来说,FDR用比较温和的方法对p值进行了校正。其试图在假阳性和假阴性间达到平衡,将假/真阳性比例控制到一定范围之内。例如,如果检验1000次,我们设定的阈值为0.05(5%),那么无论我们得到多少个差异蛋白,这些差异蛋白出现假阳性的概率保持在5%之内,这就叫FDR<5%。

    那么我们怎么从p value 来估算FDR呢,人们设计了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg方法,简称BH法。虽然这个估算公式并不够完美,但是也能解决大部分的问题,主要还是简单好用!

    FDR的计算方法
    除了可以使用excel的BH计算方法外,对于较大的数据,我们推荐使用R命令p.adjust。 p.adjust(p, method = p.adjust.methods, n = length§)

    p.adjust.methods

    c(“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,

    “fdr”, “none”)

    我们还可以从R命令p.adjust的源代码,了解其运行的机制是什么。

    p.adjust
    function (p, method = p.adjust.methods, n = length§){
    method <- match.arg(method)
    if (method == “fdr”)
    method <- “BH”
    nm <- names§
    p <- as.numeric§
    ……
    BH = {
    i <- lp:1L
    o <- order(p, decreasing = TRUE)
    ro <- order(o)
    pmin(1, cummin(n/i * p[o]))[ro]
    }
    ……
    p0
    }
    其实该函数表达的意思是这样的:

    我们将一系列p值、校正方法(BH)以及所有p值的个数(length§)输入到p.adjust函数中。
    将一系列的p值按照从大到小排序,然后利用下述公式计算每个p值所对应的FDR值。 公式:p * (n/i), p是这一次检验的p value,n是检验的次数,i是排序后的位置ID(如最大的P值的i值肯定为n,第二大则是n-1,依次至最小为1)。
    将计算出来的FDR值赋予给排序后的p值,如果某一个p值所对应的FDR值大于前一位p值(排序的前一位)所对应的FDR值,则放弃公式计算出来的FDR值,选用与它前一位相同的值。因此会产生连续相同FDR值的现象;反之则保留计算的FDR值。
    将FDR值按照最初始的p值的顺序进行重新排序,返回结果。

    p值还是 FDR ?
    差异分析
    如何筛选显著性差异基因,p value, FDR 如何选
    经常有同学询问如何筛选差异的基因(蛋白)。已经计算了表达量和p value值,差异的基因(蛋白)太多了,如何筛选。其中最为关键的是需要对p value进行校正。

    基本概念:

    零假设:在随机条件下的分布。

    p值:在零假设下,观测到某一特定实验结果的概率称为p值。

    假阳性:得到了阳性结果,但这个阳性结果是假的。

    假阴性:得到了阴性结果,但这个阴性结果是假的。

    单次检验:

    针对单个基因(蛋白),采用统计检验,假设采用的p值为小于0.05,我们通常认为这个基因在两个(组)样本中的表达是有显著差异的,但是仍旧有5%的概率,这个基因并不是差异基因。

    单多次检验:

    当两个(组)样本中有10000个基因采用同样的检验方式进行统计检验时,这个时候就有一个问题,单次犯错的概率为0.05, 进行10000次检验的话,那么就有0.05*10000=500 个基因的差异被错误估计了。

    多重检验矫正:

    为了解决多次检验带来的问题,我们需要对多次检验进行校正。那如何校正呢?在此介绍两种方法:

    Bonferroni 校正法
    Bonferroni校正法:如果进行N次检验,那么p值的筛选的阈值设定为p/N。 比如,进行10000次检验的话,如果p值选择为0.05, 那么校正的p值筛选为0.000005。 p值低于此的基因才是显著性差异基因。
    该方法虽然简单,但是过于严格,导致最后找的差异基因很少,甚至找不到差异的基因。

    FDR(False Discovery Rate) 校正法
    FDR错误控制法是Benjamini于1995年提出的一种方法,基本原理是通过控制FDR值来决定p值的值域。相对Bonferroni来说,FDR用比较温和的方法对p值进行了校正。其试图在假阳性和假阴性间达到平衡,将假/真阳性比例控制到一定范围之内。
    那么怎么从p值来估算FDR呢,人们设计了几种不同的估算模型。其中使用最多的是Benjamini and Hochberg方法,简称BH法。该方法分两步完成,具体如下:
    2.1 假设总共有m个候选基因,每个基因对应的p值从小到大排列分别是p(1),p(2),…,p(m)
    2.2 若想控制FDR不能超过q,则只需找到最大的正整数i,使得 p(i)<= (i*q)/m . 然后,挑选对应p(1),p(2),…,p(i)的基因做为差异表达基因,这样就能从统计学上保证FDR不超过q。

    如何实现多重检验:

    如果你了解R语言的话,那么采用p.adjust方法就可以了。

    简单使用DESeq做差异分析
    Posted: 五月 06, 2017 Under: Transcriptomics By Kai no Comments

    DESeq这个R包主要针对count data,其数据来源可以是RNA-Seq或者其他高通量测序数据。类似地,对于CHIP-Seq数据或者质谱肽段数据也是使用的。

    由于DESeq是一个R包,因此使用它需要一点点R基础语法。

    首先需要读入一个数据框,列代表每个sample,行代表每个gene

    database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
    database <- database_all[,1:6]
    这里主要对于两两比较的数据,因此我取了数据的前6列,分别是两组样品,每组3个生物学重复

    设定分组信息,也就是样本分组的名称

    type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
    我这里是样品1是LC_1,样品2是LC_2

    由于DESeq包要求接下来的count data必须要整数型,因此我们需要对数据进行取整,然后将数据database和分组信息type读入到cds对象中

    database <- round(as.matrix(database))
    cds <- newCountDataSet(database,type)
    接下里对于不同类型的数据要进行不同的处理,可以粗略分为有生物学重复数据、有部分生物学重复数据以及无生物学重复数据

    4.1 有生物学重复

    cds <- estimateSizeFactors(cds)
    cds <- estimateDispersions(cds)
    res <- nbinomTest(cds,“LC_1”,“LC_2”)
    其通过estimate the dispersion并对count data进行标准化,然后得到每个gene做T test检验

    4.2 对于部分样品有生物学重复

    cds <- estimateSizeFactors(cds)
    cds <- estimateDispersions(cds)
    res <- nbinomTest(cds,“LC_1”,“LC_2”)
    其步骤跟有上述的一样的,DESeq会根据有生物学重复的样品来estimate the dispersion,当然要保证unreplicated condition does not have larger variation than the replicated one

    4.3 对于没有生物学重复

    cds <- estimateDispersions(cds, method=“blind”, sharingMode=“fit-only” )
    res <- nbinomTest(cds,“LC_1”,“LC_2”)
    注意参数method=”blind” 和 sharingMode=”fit-only”即可

    最后就是查看符合阈值的差异基因有多少个即可,然后将结果输出到csv文件中方便查看

    table(res p a d j < 0.05 ) r e s < − r e s [ o r d e r ( r e s padj <0.05) res <- res[order(res padj<0.05)res<res[order(respadj),]
    sum(res$padj<=0.01,na.rm = T)
    write.csv(resdata,file = “LC_1_vs_LC_2_DESeq.csv”)
    Summary
    DESeq在前几年的文章中经常被使用,但是现在有了其升级版DESeq2,后者相比前者对于犯第一类错误卡的并不是那么严格了,所以在同样的padj的阈值下,筛选到的差异基因的数目也会多一点。

    并且上述中,我都只使用nbinomTest来获得由显著差异的gene,其实DESeq包还提供了其他检验方法,比如nbinomGLMTest函数,具体可以查看DESeq文档,还有其他对应对因素分析的使用方法。

    DESeq2和EdgeR都可用于做基因差异表达分析,主要也是用于RNA-Seq数据,同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。

    这两个都属于R包,其相同点在于都是对count data数据进行处理,都是基于负二项分布模型。因此会发现,用两者处理同一组数据,最后在相同阈值下筛选出的大部分基因都是一样的,但是有一部分不同应该是由于其估计离散度的不同方法所导致的。

    DESeq2的使用方法:
    输入矩阵数据,行名为sample,列名为gene;DESeq2不支持无生物学重复的数据,因此我选择了2个样本,3个生物学重复的数据;并对count data取整(经大神指点,这里需要说明下,我的测试数据readcount是RSEM定量的结果,并不是常见的htseq-count的结果,所以count值会有小数点,而DESeq2包不支持count数有小数点,所以这里需要round取整)。

    database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
    database <- database_all[,1:6]
    #type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
    database <- round(as.matrix(database))
    设置分组信息以及构建dds对象

    condition <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
    coldata <- data.frame(row.names = colnames(database), condition)
    dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
    使用DESeq函数进行估计离散度,然后进行标准的差异表达分析,得到res对象结果

    dds <- DESeq(dds)
    res <- results(dds)
    最后设定阈值,筛选差异基因,导出数据

    table(res p a d j < 0.05 ) r e s < − r e s [ o r d e r ( r e s padj <0.05) res <- res[order(res padj<0.05)res<res[order(respadj),]
    resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
    write.csv(resdata,file = “LC_1_vs_LC_2.csv”)
    EdgeR的使用方法:
    跟DESeq2一样,EdgeR输入矩阵数据,行名为sample,列名为gene;DESeq2不支持无生物学重复的数据,因此我选择了2个样本,3个生物学重复的数据。

    exprSet_all <- read.table(file = “readcount”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
    exprSet <- exprSet_all[,1:6]
    group_list <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
    设置分组信息,去除低表达量的gene以及做TMM标准化

    exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
    exprSet <- DGEList(counts = exprSet, group = group_list)
    exprSet <- calcNormFactors(exprSet)
    使用qCML(quantile-adjusted conditional maximum likelihood)估计离散度(只针对单因素实验设计)

    exprSet <- estimateCommonDisp(exprSet)
    exprSet <- estimateTagwiseDisp(exprSet)
    寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计),然后按照阈值进行筛选即可

    et <- exactTest(exprSet)
    tTag <- topTags(et, n=nrow(exprSet))
    tTag <- as.data.frame(tTag)
    write.csv(tTag,file = “LC_1_vs_LC_2_edgeR.csv”)
    Summary
    以上我主要针对单因素两两比较组进行差异分析,其实DESeq2和EdgeR两个R包都可以对多因素进行差异分析。

    DESeq2修改以上代码的分组信息design参数以及在差异分析results函数中添加所选定的分组因素,其他代码基本一样,具体参照DESeq2手册

    EdgeR则需要用Cox-Reid profile-adjusted likelihood (CR)方法来估算离散度,y <- estimateDisp(y, design)或者分别使用三个函数(y <- estimateGLMCommonDisp(y, design),y <- estimateGLMTrendedDisp(y, design), )y <- estimateGLMTagwiseDisp(y, design);然后差异表达分析也跟单因素分析不同,主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test,具体代码可以参照EdgeR手册。

    简单使用limma做差异分析
    Posted: 五月 12, 2017 Under: Transcriptomics By Kai no Comments

    首先需要说明的是,limma是一个非常全面的用于分析芯片以及RNA-Seq的差异分析,按照其文章所说:

    limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments.

    在这我只是对其中的一种情况进行简单的总结,比如这个包可以处理RNA-Seq数据,我简单的以两个比较组进行分组为例,至于其他分组情况,请看limma说明文档,有非常详细的说明,非常亲民。

    首先我们还是输入count矩阵,这里也跟其他差异分析R包一样,不要输入已经标准化的数据。顺便也加载下edgeR这个R包

    library(limma)
    library(edgeR)
    counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
    接着按照文档的说明以及limma包的习惯,我们需要对count进行标准化以及转化为log2的值,这里标准化的方法为TMM,使用edgeR里面的calcNormFactors函数即可

    dge <- DGEList(counts = counts)
    dge <- calcNormFactors(dge)
    logCPM <- cpm(dge, log=TRUE, prior.count=3)
    这里prior.count值我粗略理解为是为了防止log2()遇到过于小的值

    然后跟其他包一样,设置分组信息

    group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
    design <- model.matrix(~group_list)
    colnames(design) <- levels(group_list)
    rownames(design) <- colnames(counts)
    接下来就是常规的差异分析

    fit <- lmFit(logCPM, design)
    fit <- eBayes(fit, trend=TRUE)
    output <- topTable(fit, coef=2,n=Inf)
    sum(output$adj.P.Val<0.05)
    到这里为止,我们主要是用了limma包里对RNA-Seq差异分析的limma-trend方法,该方法主要适用于样本间测序深度基本保持一致的情况,也就是说两个样本的文库(reads数目)大小相差的不悬殊(说明文档中是默认3倍以内?)

    当文库大小在样本间变化幅度相当大的话,可以使用limma的voom方法,可按照下面的代码流程:

    count数据的输入以及数据标准化还是跟之前的一样

    counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
    dge <- DGEList(counts = counts)
    dge <- calcNormFactors(dge)
    还是一样需要分组信息

    group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
    design <- model.matrix(~group_list)
    colnames(design) <- levels(group_list)
    rownames(design) <- colnames(counts)
    接下来进行voom转化

    v <- voom(dge, design, plot=TRUE)
    其实可以不进行TMM标准化,直接用count数据进行voom转化,如:

    v <- voom(counts, design, plot=TRUE)
    最后就是普通的差异分析过程

    fit <- lmFit(v, design)
    fit <- eBayes(fit)
    output <- topTable(fit, coef=2,n=Inf)
    sum(output$adj.P.Val<0.05)
    Summary
    Limma长久以来就是一个非常流行的差异分析R包,其内容涉及的非常广泛,用于RNA-Seq只是其内容的一小部分,并且使其处理RNA-Seq数据也使用芯片类似线性模型下,并且按照其说法,limma包比其他基于负二项式分布模型的差异分析R包更加的优秀。

    其实差异分析不外乎数据的标准化以及统计模型分析差异两个方面的作用,每个差异分析R包都有其自身的优点,个人理解,取舍在于自己的理解以及想法即可。

    其实,自己对于limma包的理解还是比较粗浅的

    转录组差异表达分析小实战(一)
    Posted: 七月 28, 2017 Under: Transcriptomics By Kai no Comments

    读文献获取数据
    文献名称:AKAP95 regulates splicing through scaffolding
    RNAs and RNA processing factors

    查找数据:Data availability
    The RIP-seq an RNA-seq data have been deposited in the Gene
    Expression Omnibus database, with accession code GSE81916. All other data is
    available from the author upon reasonable request.

    获得GSE号:GSE81916

    下载测序数据
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916获取数据信息,并点击网址下方的ftp,下载测序数据

    从https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA323422可知我们需要的mRNA测序编号为SRR3589956到SRR3589962

    通过Apera下载SRR数据,这里以SRR3589956为例:

    ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.nih.gov:sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
    转化fastq测序数据
    通过sratoolkit工具将SRR文件转化为fastq格式的测序数据(写了个shell循环)

    for i in ( s e q 5662 ) ; d o n o h u p f a s t q − d u m p − − s p l i t − 3 S R R 35899 (seq 56 62);do nohup fastq-dump --split-3 SRR35899 (seq5662);donohupfastqdumpsplit3SRR35899{i} &;done
    通过fastqc对每个fastq文件进行质检,用multiqc查看整体质检报告(对当前目录下的fastq测序结果进行质检,生成每个fq文件的质检报告总multiqc整合后统计查看)

    fastqc *.fastq
    multiqc ./
    点击这个url可以查看我这个multiqc报告:http://www.bioinfo-scrounger.com/data/multiqc_report.html

    如果有接头或者质量值不达标的需要进行过滤,这次的数据质量都不错,因此直接进行比对即可
    序列比对
    安装hisat2软件,下载人类的hiast2索引文件

    hisat2下载并安装:

    ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
    unzip hisat2-2.1.0-Linux_x86_64.zip
    下载hisat2的human索引

    ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
    tar zxvf hg19.tar.gz
    用hisat2进行比对,测序数据放在data目录下,索引文件放在reference/index/hisat2/hg19目录下,SRR3589956-SRR3589958为人的测序数据

    for i in KaTeX parse error: Undefined control sequence: \ at position 28: …do hisat2 -p 4 \̲ ̲-x ~/reference/…{i}_1.fastq -2 ./data/SRR35899KaTeX parse error: Undefined control sequence: \ at position 13: {i}_2.fastq \̲ ̲-S SRR35899i.sam >SRR35899${i}.log;done
    用samtools将sam文件转化为bam文件,并使用默认排序

    for i in ( s e q 5658 ) ; d o s a m t o o l s s o r t − @ 5 − o S R R 35899 (seq 56 58);do samtools sort -@ 5 -o SRR35899 (seq5658);dosamtoolssort@5oSRR35899{i}.bam SRR35899${i}.sam;done
    reads计数
    用htseq对比对产生的bam进行count计数

    htseq安装,使用miniconda,省事!唯一的问题是htseq版本不是最新的,是0.7.2。想要最新版还是要正常安装,可参考http://www.biotrainee.com/thread-1847-1-2.html

    conda install -c bioconda htseq
    用htseq将对比后的结果进行计数

    for i in KaTeX parse error: Undefined control sequence: \ at position 48: …m -r pos -s no \̲ ̲SRR35899{i}.bam ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf
    1>SRR35899 i . c o u n t 2 > S R R 35899 {i}.count 2>SRR35899 i.count2>SRR35899{i}_htseq.log;done
    将3个count文件(SRR3589956.count,SRR3589957.count,SRR3589958.count)合并成一个count矩阵,这是就需要脚本来解决这个问题,不然其他方法会稍微麻烦点

    #!/usr/bin/perl -w
    use strict;

    my $path = shift @ARGV;
    opendir DIR, $path or die;
    my @dir = readdir DIR;

    my $header;
    my @sample;
    my %hash;
    foreach my KaTeX parse error: Expected '}', got 'EOF' at end of input: …dir) { if (file =~ /^\w+.*.count/) {
    push @sample, $file;
    KaTeX parse error: Undefined control sequence: \t at position 12: header .= "\̲t̲file";
    open my $fh, f i l e o r d i e ; w h i l e ( < file or die; while (< fileordie;while(<fh>) {
    chomp;
    next if ($_ =~ /^\W+/);
    my @array = split /\t/, $_;
    KaTeX parse error: Expected '}', got 'EOF' at end of input: hash{array[0]} -> {$file} = $array[1];
    }
    close KaTeX parse error: Expected 'EOF', got '}' at position 9: fh; }̲ } print "header\n";
    map{
    my $gene = ; p r i n t " _; print " ;print"gene";
    foreach my KaTeX parse error: Undefined control sequence: \t at position 33: … print "\̲t̲".hash{KaTeX parse error: Expected 'EOF', got '}' at position 5: gene}̲ -> {file};
    }
    print “\n”;
    }keys %hash;
    按照接下来的剧本,应该讲count_matrix文件导入DESeq进行差异表达分析。但是从这篇文章的Bioinformatic analyses部分可以发现,作者的control组的2组数据是来自2个不同的批次(一个是SRR3589956,另外一个来源GSM1095127 in GSE44976),treat组倒是同一个批次(SRR3589957和SRR3589958)。但是对于Mouse cells来说,倒是满足2个control和2个treat都正常来自同个批次,因此打算重新用SRR3589959-SRR3589962重新做个一个count_matrix进行后续差异分析

    转录组差异表达分析小实战(二)
    Posted: 八月 14, 2017 Under: Transcriptomics By Kai no Comments

    差异基因表达分析
    我按照前面的流程转录组差异表达分析小实战(一),将小鼠的4个样本又重新跑了一遍,从而获得一个新的count文件:mouse_all_count.txt,有需要的话,可以下载下来进行后续的差异分析。

    一般来说,由于普遍认为高通量的read count符合泊松分布,所以一些差异分析的R包都是基于负二项式分布模型的,比如DESeq、DESeq2和edgeR等,所以这3个R包从整体上来说是类似的(但各自标准化算法是不一样的)。

    当然还有一个常用的R包则是Limma包,其中的limma-trend和limma-voom都能用来处理RNA-Seq数据(但对应适用的情况不一样)

    下面准备适用DESeq2和edgeR两个R包分别对小鼠的count数据进行差异表达分析,然后取两者的结果的交集的基因作为差异表达基因。

    DEseq2

    library(DESeq2)
    ##数据预处理
    database <- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = T, row.names = 1)
    database <- round(as.matrix(database))

    ##设置分组信息并构建dds对象
    condition <- factor(c(rep(“control”,2),rep(“Akap95”,2)), levels = c(“control”, “Akap95”))
    coldata <- data.frame(row.names = colnames(database), condition)
    dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
    dds <- dds[ rowSums(counts(dds)) > 1, ]

    ##使用DESeq函数估计离散度,然后差异分析获得res对象
    dds <- DESeq(dds)
    res <- results(dds)

    #最后设定阈值,筛选差异基因,导出数据(全部数据。包括标准化后的count数)
    res <- res[order(res$padj),]
    diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
    diff_gene_DESeq2 <- row.names(diff_gene)
    resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
    write.csv(resdata,file = “control_vs_Akap95.csv”,row.names = F)
    最终获得572个差异基因(筛选标准为padj < 0.05, |log2FoldChange| > 1)

    edgeR

    library(edgeR)
    ##跟DESeq2一样,导入数据,预处理(用了cpm函数)
    exprSet<- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
    group_list <- factor(c(rep(“control”,2),rep(“Akap95”,2)))
    exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]

    ##设置分组信息,并做TMM标准化
    exprSet <- DGEList(counts = exprSet, group = group_list)
    exprSet <- calcNormFactors(exprSet)

    ##使用qCML(quantile-adjusted conditional maximum likelihood)估计离散度(只针对单因素实验设计)
    exprSet <- estimateCommonDisp(exprSet)
    exprSet <- estimateTagwiseDisp(exprSet)

    ##寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计),然后按照阈值进行筛选即可
    et <- exactTest(exprSet)
    tTag <- topTags(et, n=nrow(exprSet))
    diff_gene_edgeR <- subset(tTagKaTeX parse error: Expected 'EOF', got '&' at position 19: …le, FDR < 0.05 &̲ (logFC > 1 | l…table,file = “control_vs_Akap95_edgeR.csv”)
    最终获得688个差异基因(筛选标准为FDR < 0.05, |log2FC| > 1)

    取DESeq2和edgeR两者结果的交集

    diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]
    最终的差异基因数目为545个

    head(diff_gene)
    [1] “ENSMUSG00000003309.14” “ENSMUSG00000046323.8” “ENSMUSG00000001123.15”
    [4] “ENSMUSG00000023906.2” “ENSMUSG00000044279.15” “ENSMUSG00000018569.12”
    其他两个R包(DESeq和limma)就不在这尝试了,我之前做过对于这4个R包的简单使用笔记,可以参考下:
    简单使用DESeq做差异分析
    简单使用DESeq2/EdgeR做差异分析
    简单使用limma做差异分析

    GO&&KEGG富集分析
    以前一直没有机会用Y叔写的clusterProfiler包,这次正好看说明用一下。

    GO富集,加载clusterProfiler包和org.Mm.eg.db包(小鼠嘛),然后将ENSEMBL ID后面的版本号去掉,不然后面不识别这个ID,然后按照clusterProfiler包的教程说明使用函数即可。

    library(clusterProfiler)
    library(org.Mm.eg.db)

    ##去除ID的版本号
    diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, “\.”)[[1]][1]}))
    ##GOid mapping + GO富集
    ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db,
    keytype = “ENSEMBL”, ont = “BP”, pAdjustMethod = “BH”,
    pvalueCutoff = 0.01, qvalueCutoff = 0.05)
    ##查看富集结果数据
    enrich_go <- as.data.frame(ego)
    ##作图
    barplot(ego, showCategory=10)
    dotplot(ego)
    enrichMap(ego)
    plotGOgraph(ego)
    KEGG富集,首先需要将ENSEMBL ID转化为ENTREZ ID,然后使用ENTREZ ID作为kegg id,从而通过enrichKEGG函数从online KEGG上抓取信息,并做富集

    library(clusterProfiler)
    library(org.Mm.eg.db)

    ##ID转化
    ids <- bitr(diff_gene_ENSEMBL, fromType = “ENSEMBL”, toType = “ENTREZID”, OrgDb = “org.Mm.eg.db”)
    kk <- enrichKEGG(gene = ids[,2], organism = “mmu”, keyType = “kegg”,
    pvalueCutoff = 0.05, pAdjustMethod = “BH”, qvalueCutoff = 0.1)
    ##查看富集结果数据
    enrich_kegg <- as.data.frame(kk)
    ##作图
    dotplot(kk)
    到这里为止,转录组的差异表达分析算是做完了,简单的来说,这个过程就是将reads mapping 到reference上,然后计数获得count数,然后做差异分析,最后来个GO KEGG,over了。。。

    对于mapping和计数这两部还有其实还有好多软件,具体可见文献:Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis,有时间可以都尝试下。

    至于GO && KEGG这两步,对于人、小鼠等模式物种来说,不考虑方便因素来说,完全可以自己写脚本来完成,数据可以从gene ontology官网下载,然后就是GO id与gene id相互转化了。KEGG 也是一样,也可以用脚本去KEGG网站上自行抓取,先知道URL规则,然后爬数据即可。

    展开全文
  • [ERP]ERP中标准成本的差异分析控制

    千次阅读 2015-10-12 09:07:37
    本文将着重讨论差异分析控制。 一、 标准成本的差异类型标准成本的差异的来源主要存在三个阶段,也即:采购、存储和生产。采购过程的差异主要围绕订单、发票及收货环节的不配比造成。理论上三者的各要素应当一致,...

    当前大多成熟的ERP系统均集成了标准成本管理系统。由于ERP系统的高度集成的特性,各种因素对标准成本影响的结果需要耗费相关人员大量的时间和精力去分析,毕竟再智能的系统依然需要人的参与。本文将着重讨论差异的分析控制。

    一、 标准成本的差异类型
    
    标准成本的差异的来源主要存在三个阶段,也即:采购、存储和生产。
    
    采购过程的差异主要围绕订单、发票及收货环节的不配比造成。理论上三者的各要素应当一致,但事实上除了采购价格和标准成本的人为设置的原因之外,由于其他各种原因造成的差异因素千差万别。为便于控制主要因素,价格、发票、标准等差异几乎成为标准的设置。
    
    理论上存储是静止的,但由于各种原因,实物与账面的差异始终存在。物流过程中,任何存放物料的地点都存在账实不符的风险。
    
    生产环节可能由于各种原因造成的实际用料和物料单的不符,人工和制造费用也会由于标准费率与实际的偏离造成差异。
    
    二、 差异的分析控制
    
    只有充分把握差异的来源和成因,才能从根本上分析成本的差异,并采取有效的手段控制成本减少异常的差异出现,以保持存货和成本的真实性与合理性。
    
    采购环节的差异始终围绕订单,发票,接收三个环节产生。对于采购环节的差异,采购部门应逐月搜集相关信息,找到问题的根源加以解决。针对采购流程的全部作业活动的关键点实施跟踪控制,加强信息流的准确性和及时性。例如针对采购订单的下达和传递过程,应搜集采购员下达的订单信息和供应商收到的信息进行分析。通常由于控制的薄弱造成下单错误(数量,单价,单位,货币等),或者由于供应商的疏忽造成处理错误等。同时由于供应商用错误的订单信息会造成相应的发票,出货的错误。而事实上供应商财务人员的疏忽,以及发货的错误,同时接收人员的大意,都造成差异的产生。若果追溯采购差异,上述所有环节均可能成为问题的根源。所以相应的控制流程必须清晰而且有效。采购部门应在财务部门的配合下,考核差异到订单,加强对于这些控制点的控制,加重在考核指标中的权重。
    
    存储环节通常会由于公司管理的漏洞、库管人员的疏漏、物料价值的重估、BOM问题等造成差异产生。ERP系统中的记录的物料信息和实际的盘存差异,应由责任人员书面提请报告并批准后相应调整库存。由于物料价值重估产生的差异应用后进先出法在存货成本和销售成本之间进行分配。BOM错误会造成存货数量和价值的错误,制造企业应对BOM结构和版本进行审核和控制,避免由此造成的差异。
    
    制造环节的差异通常由于BOM的错误、库存、订单错误、领料等造成各种异常差异和成本波动。公司应对生产订单的差异原因进行分析,并由生产管理人员牵头组织相关人员对差异原因进行剖析,找出原因并采取措施。
    
    三、 结论及建议
    
    由于ERP系统的高度集成性,任何上游的异常业务均会造成差异的产生。公司应标准化各操作流程,加强对系统录入信息的质量控制,并采取差异负责制由责任人对差异负责,逐步从源头减少差异产生的几率,使得标准成本真正成为企业管理的利器。
    
    展开全文
  • 依据评价结果,分析了我国东中西部地区生态文明建设的区域差异,指出子系统评价值在60分以下省份进一步促进生态文明建设的对策建议。主要结论和政策建议包括:①我国生态文明建设按照西部地区、中部地区、东部地区...
  • 研究了大断面黄土隧道中系统锚杆的受力状况,得出了拱部和边墙锚杆的受力特点,从径向位移的衰减规律和围岩压力分布形式两个角度分析了系统锚杆受力差异较大的原因,相关研究结论对黄土隧道的设计和施工具有参考价值...
  • 本文利用WTO各成员国向WTO通报的TBT和SPS数量,对中国与发达经济体的WTO技术性贸易壁垒的绩效差异进行实证分析,得出以下结论:发达经济体的技术性贸易壁垒阻碍了我国的出口贸易,但对其产品出口到中国有利;而中国的技术...
  • 目的:分析重庆市欠发达地区卫生资源配置的现状和差异,为卫生规划和... 结论:建议增加欠发达地区卫生政策的差异性和针对性,实现卫生资源的动态管理和控制,促进卫生人力资源的可持续发展,并确保提高卫生筹资水平。
  • 以四川省1998~2007年各市(州)人口状况和GDP为基础,运用基尼系数法探讨...结论是:1998~2007年间,基尼系数一直保持在0.3左右。这表明,自西部大开发实施以来,四川省各地区经济发展水平差距不显著,收入分配总体比较均衡。
  • (听别人说的面试中问过的一个问题,记不清了,如有偏差希望大家补充纠正) 问题:完成相同的功能,for循环使用++和–哪个效率更高? 通常我们在编写需要for... 结论:通常来讲第2种用–的方法效率更高 原因如下:

    (听别人说的面试中问过的一个问题,记不清了,如有偏差希望大家补充纠正)
    问题:完成相同的功能,for循环使用++和–哪个效率更高?

    通常我们在编写需要for循环实现的程序时有两种实现方法:

    1. for(int i = 0;i<arr.length;i++)
    2. for(int i = arr.length-1;i>=0;i--)

    分析:

    结论:通常来讲第2种用–的方法效率更高

    原因如下:

    • 在循环中,循环条件会被反复计算,如果不使用复杂表达式,而使循环条件值不变的话,程序将会运行的更快。

    • 第2种方法和0比较作为判读,很多cpu和0比较可以减少1个指令.

    参考:
    i–操作本身会影响CPSR(当前程序状态寄存器),CPSR常见的标志有N(结果为负), Z(结果为0),C(有进位),O(有溢出)。i > 0,可以直接通过Z标志判断出来。

    i++操作也会影响CPSR(当前程序状态寄存器),但只影响O(有溢出)标志,这对于i < n的判断没有任何帮助。所以还需要一条额外的比较指令,也就是说每个循环要多执行一条指令。

    • 考虑使用情况的话:arr.length可能会变化, 所以有可能比较时每次都要从arr里面取length的值.是一个问题。

    所以有:

    1,for(int i = 0;i<arr.length;i++) 最耗时,

    2,for(int i=0,len=list.size();i<len;i++) 次之,

    3, for(int i = arr.length-1;i>=0;i--) 最快

    1,最耗时不需多言,
    2与3 基本上差距在10%之内

    总结:(针对上边三种方法)

    • 如果arr.length需要变化则使用++,否则使用–。
    • 至于把arr.length保存到静态变量,前一种情况下不可能,后一种情况下则没必要

    补充加油:

    如何提高for循环效率

    在我们所写的程序中,几乎没有不用到for循环的,但是,对于for循环,很多人确实效率很低的,包括我看得很多代码,for循环的执行效率非常低,下面我就举个例子来说明:

    for(i=0;i<strlen(string);i++)

    这个上边的程序程序我想大家都明白,那我问问读者,你知道这个程序的效率是多少吗?
    你肯定不屑的说,不就是n吗?其实,你错了,你说的n只是在算法层面上的优化,其实对于底层的优化还没做好,这段代码的效率是n^2(n的平方),为什么?是这样的,我们在每次循环的时候,都会调用strlen函数,这个函数的效率也是n,所以,我们要再加一个变量,比如下边所看到的,

    int k=strlen(string);
    for(i=0;i<k;i++)

    这个程序的效率是2n,当n很大时,我就不说了。
    再看这个例子:

     int i,m,k1,k2,k3,k4,k5,k6,k7,k8;  
     m=10000000;  
     for(i=0;i<m;i++)
     {  
     k1++;  k5++;  
     k2++;  k6++;  
     k3++;  k7++;  
     k4++;  k8++;  
     } 
    

    我们再看这个for循环,我们是让这些数每次循环都加1,这个效率也不是达到了最优,是这样的,在每次for循环的时候,编译器会给变量i,变量m每个独占一个寄存器,因为是循环嘛,编译器给i和m也是为了效率考虑,不用再向寄存器加载每次循环的判断变量了,但是,一般基于Intel的处理器都只有8个通用寄存器,这样,我们就剩下了6个寄存器给for循环里的内容中的变量使用,但是,我们的变量有8个,我们这时会在把其他的变量入栈,这样,我们的效率变低了,每次出战或者入栈都会消耗两个CPU时钟周期,这样,我们就总共满了8个周期,但是,如果更多呢?呵呵!
    我们可以这样该我们的程序;

     int i,m,k1,k2,k3,k4,k5,k6,k7,k8;  
     m=10000000;  
    
    for(i=0;i<m;i++){  
      k1++;  k5++;  
      k2++;  k6++;  
    }  
    for(i=0;i<m;i++){  
      k3++;  k7++;  
      k4++;  k8++;  
    }  
    

    这样的效率就得到了很大改善!


    更多优化技巧用到再补充,也欢迎看到这篇文章结尾同仁留言补充

    展开全文
  • 关键词:基因芯片、R、筛选、预处理、差异分析 NCBI-淋巴癌突变与未突变基因的差异分析 PS:好久没分享生信了,这是一年前做的一次生信task(准确来说是2018年11月了),这里分享一下给大家,有助于一些小伙伴们想...

    关键词:基因芯片、R、筛选、预处理、差异分析

    NCBI-淋巴癌突变与未突变基因的差异分析
    PS:好久没分享生信了,这是一年前做的一次生信task(准确来说是2018年11月了),这里分享一下给大家,有助于一些小伙伴们想通过常规的,使用NCBI科研数据库+R编程语言方式,进行对某种癌的差异分析。最近用心做了一些更棒的生信task,相信不久会分享出来~
    PS2:如果这篇笔记有什么不足,或者疑惑不解的地方,可以提出来,看到会回答(留个TODO给大家

    简单介绍

    通过使用计算机对NCBI下载的基因数据,进行筛选,其中包括质量控制、质量分析,聚类分析,对数据进行预处理,并使用差异分析,获取突变与未突变基因的差异表达基因,进行注释,最后生成研究结果。
    本次研究主要以R作为数据处理的语言,多次对38组细胞进行筛选、其中编写了5个版本的实验脚本,研究结果主要有:基因的差异分析,需要对海量的数据进行十分严谨的筛选与预处理,突变与未突变基因在组间差异较为明显。

    一、研究对象

    淋巴癌突变与未突变基因的差异分析

    二、研究工具

    1. NCBI

    NCBI, The National Center for Biotechnology Information biomedical genomic,美国国家生物技术旨在通过提供生物医学和基因组信息提供给热爱科学与健康事业的人们。NCBI功能强大,提供了海量丰富的基因数据库、计算机分析工具、文献等,研究者可以根据自己的研究方向,选择自己所需要的数据与处理方式,推动我们对基因遗传的理解,从而推进健康和疾病的理解。

    2. R

    R语言是用于统计分析,图形表示和报告的编程语言和软件环境。 R语言由Ross Ihaka和Robert Gentleman在新西兰奥克兰大学创建,目前由R语言开发核心团队开发。 R语言在GNU通用公共许可证下免费提供,并为各种操作系统(如Linux,Windows和Mac)提供预编译的二进制版本。 这种编程语言被命名为R语言,基于两个R语言作者的名字的第一个字母(Robert Gentleman和Ross Ihaka),并且部分是贝尔实验室语言S的名称。

    三、研究流程

    1、 数据源载入

    1) 针对参考的文献,通过基因芯片研究突变/未突变的套细胞淋巴瘤患者的关系
    2) 对应的数据下载源:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36000
    3) 数据下载源:
    数据源下载
    数据源下载到本地

    4) 参考数据集:
    参考数据集

    2、 使用R读取下载的数据

    在使用R进行前期准备的过程中,需要先将应用到的一些常用的lib包加载到本机(如果看不懂,那可以先跳过,继续往下看),不同的库会依赖到其他库,那么也需要提前下载好。
    并为了增加可读性,我们需要将基因芯片重命名,并根据已有原本数据的分类进行简单的分类与命名(例如,本次研究的是未突变与突变淋巴细胞)

    PS:下载package的时候,更新了Bioconductor version 3.8,更新到了R的Version 3.5.1(2018-07-02),所以需要使用到BiocManager::install的方式,bioCLite已经提示废弃。

    具体流程

    1. 设置工作目录
    2. 先下载并加载必要的package,已安装并加载的,可以不用重复这些步骤
    3. 可以通过GEOquery直接读取基因芯片的编号(例如本次差异分析中使用到的是“GSE36000”);也可以直接在官网下载raw包,并通过ReadAffy的方式,这里会有AffyBatch特定的格式(本差异分析采取的是这种方式)。
    4. 查看数据类型、基本信息,并重命名芯片名称
      (温馨提示:R代码部分后面会提供)

    代码参考:

    # 本文字符格式基于UTF-8编码
    # 【废弃的前数据-不合适】:1.下载数据:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE47811
    # 1.下载数据:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36000
    # 2.读取下载的数据
    # 2.1 先下载并加载必要的package,已安装并加载的,可以不用重复这些步骤
    # Bioconductor用于生物芯片数据分析和基因组数据分析的软件包
    source("https://bioconductor.org/biocLite.R")
    # 2.1.1 Bioconductor version 3.7版本以前,使用bioCLite下载
    biocLite("affy")
    biocLite("simpleaffy")
    biocLite("pheatmap")
    biocLite("RColorBrewer")
    biocLite("affyPLM")
    biocLite("GOstats")
    
    # 2.1.2 Bioconductor version 3.8 其中biocLite已在3.8废弃,改用BiocManager::install
    BiocManager::install("affy")
    BiocManager::install("simpleaffy")
    BiocManager::install("pheatmap")
    BiocManager::install("RColorBrewer")
    BiocManager::install("affyPLM")
    BiocManager::install("GOstats")
    BiocManager::install("survival")
    
    # 2.2 设置工作目录——放置源文件
    setwd("C:/libtobrain/R/R_Package/GSE36000_RAW")
    
    # 2.3 将GSE36000_RAW.tar解压到GSE36000_RAW目录下
    # 2.4 读取基因数据
    dir.files <- list.files(path = "C:/libtobrain/R/R_Package/GSE36000_RAW/", pattern = "*.CEL.gz")
    dir.files
    library(parallel)
    library(BiocGenerics)
    library(Biobase)
    library(affy)
    affy.data <- ReadAffy(filenames = dir.files)
    
    # 2.4.1 查看数据类型
    data.class(affy.data)
    class(affy.data)
    # 2.4.2 查看芯片的基本信息
    show(affy.data)
    # 2.4.3 查看芯片的样本名
    sampleNames(affy.data)
    
    # 2.4.4 重命名样品名
    # GSM879118	MCL_IGHV_UNMUT(R1394)       mcl.ighv.unmut.1
    # GSM879119	MCL_IGHV_UNMUT(R1328)       mcl.ighv.unmut.2
    # GSM879120	MCL_IGHV_UNMUT(R1329)       mcl.ighv.unmut.3
    # GSM879121	MCL_IGHV_UNMUT(R1306)       mcl.ighv.unmut.4
    # GSM879122	MCL_IGHV_UNMUT(R1332)       mcl.ighv.unmut.5
    # GSM879123	MCL_IGHV_MUT(R1399_M)       mcl.ighv.mut.1
    # GSM879124	MCL_IGHV_MUT(R1400)         mcl.ighv.mut.2
    # GSM879125	MCL_IGHV_UNMUT(R1587)       mcl.ighv.unmut.6
    # GSM879126	MCL_IGHV_UNMUT(R1680)       mcl.ighv.unmut.7
    # GSM879127	MCL_IGHV_UNMUT(268-01-5TR)  mcl.ighv.unmut.8
    # GSM879128	MCL_IGHV_MUT(269-01-5TR)    mcl.ighv.mut.3
    # GSM879129	MCL_IGHV_UNMUT(043-01-4TR)  mcl.ighv.unmut.9
    # GSM879130	MCL_IGHV_MUT(R1302)         mcl.ighv.mut.4
    # GSM879131	MCL_IGHV_MUT(R1304)         mcl.ighv.mut.5
    # GSM879132	MCL_IGHV_MUT(R1305)         mcl.ighv.mut.6
    # GSM879133	MCL_IGHV_MUT(R1628)         mcl.ighv.mut.7
    # GSM879134	MCL_IGHV_MUT(R1629)         mcl.ighv.mut.8
    # GSM879135	MCL_IGHV_MUT(R1341)         mcl.ighv.mut.9
    # GSM879136	MCL_IGHV_MUT(R1333)         mcl.ighv.mut.10
    # GSM879137	MCL_IGHV_MUT(R1585)         mcl.ighv.mut.11
    # GSM879138	MCL_IGHV_MUT(R1589)         mcl.ighv.mut.12
    # GSM879139	MCL_IGHV_MUT(R1591)         mcl.ighv.mut.13
    # GSM879140	MCL_IGHV_UNMUT(R1595)       mcl.ighv.unmut.10
    # GSM879141	MCL_IGHV_MUT(R1640)         mcl.ighv.mut.14
    # GSM879142	MCL_IGHV_MUT(R1626)         mcl.ighv.mut.15
    # GSM879143	MCL_IGHV_MUT(R1610)         mcl.ighv.mut.16
    # GSM879144	MCL_IGHV_MUT(02_201)        mcl.ighv.mut.17
    # GSM879145	MCL_IGHV_MUT(R1339)         mcl.ighv.mut.18
    # GSM879146	MCL_IGHV_MUT(R1301)         mcl.ighv.mut.19
    # GSM879147	MCL_IGHV_UNMUT(R1338)       mcl.ighv.unmut.11
    # GSM879148	MCL_IGHV_UNMUT(R1762)       mcl.ighv.unmut.12
    # GSM879149	MCL_IGHV_MUT(R1938)         mcl.ighv.mut.20
    # GSM879150	MCL_IGHV_MUT(R1940)         mcl.ighv.mut.21
    # GSM879151	MCL_IGHV_MUT(R1942)         mcl.ighv.mut.22
    # GSM879152	MCL_IGHV_MUT(R1941)         mcl.ighv.mut.23
    # GSM879153	MCL_IGHV_MUT(R1970)         mcl.ighv.mut.24
    # GSM879154	MCL_IGHV_UNMUT(R1897)       mcl.ighv.unmut.13
    # GSM879155	MCL_IGHV_UNMUT(R1388)       mcl.ighv.unmut.14
    
    c.all <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.3","mcl.ighv.unmut.4","mcl.ighv.unmut.5",
               "mcl.ighv.mut.1","mcl.ighv.mut.2",
               "mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.8",
               "mcl.ighv.mut.3",
               "mcl.ighv.unmut.9",
               "mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10","mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13",
               "mcl.ighv.unmut.10",
               "mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
               "mcl.ighv.unmut.11","mcl.ighv.unmut.12",
               "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24",
               "mcl.ighv.unmut.13","mcl.ighv.unmut.14")
    
    c.all.status <- c("unmut","unmut","unmut","unmut","unmut",
                      "mut","mut",
                      "unmut","unmut","unmut",
                      "mut",
                      "unmut",
                      "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                      "unmut",
                      "mut","mut", "mut","mut","mut", "mut",
                      "unmut","unmut",
                      "mut","mut", "mut","mut","mut",
                      "unmut","unmut")
    sampleNames(affy.data) <- c.all
    sampleNames(affy.data)
    

    3、 质量分析

    3.1 芯片灰度图

    Affymetrix芯片在印刷时,会在在左上角印制芯片的名称,在四个角印制实心的花纹,我们可以通过这个特征来了解芯片数据是否可靠。比如说,图像特别黑,说明信号强度低;图像特别亮,说明信号可能过饱和,图像亮度适中,说明信号比较可靠。

    从GSE36000芯片组中,38个芯片,以其中的一个样本GSM879118为例构造芯片灰度图,可以看出,四角实心的花纹名下,并且也有芯片名称,图像亮度适中,我们可以认为:该样本信号比较可靠。

    GSE36000芯片组-GSM879118芯片灰度图
    GSE36000芯片组-GSM879118芯片灰度图2

    代码参考

    # 2.5读取芯片灰度图
    # 2.5.1 设置工作目录——放置结果
    setwd("C:/libtobrain/R/R_Package/workSpaces/Version5/result2/")
    # 2.5.2 将基因芯片进行解析读取,存放为result目录下的pdf文件,即芯片灰度图
    pdf(file = "芯片灰度图.pdf")
    image(affy.data[,1])
    dev.off()
    
    3.2 质量控制

    采取了以上两种方式,来评估一个芯片是否质量可靠,还是不够,那么我们就需要采用更加综合的方式来采取分析。

    首先我们采取了对数据进行线性拟合,并将权重图、残差图、残差符号图保存到pdf上(也可以直接在RStudio上观察,通过笔记本来跑较大数据的分析,例如本次下载的数据就高到165MB,还是相对比较吃力的,个人建议,为了减少资源的占用,直接将图片输出会高效一些)

    以本次实验GSE36000为例,可以看出权重图、残差图、残差符号图,色彩平均,那么我们可以认定该芯片质量均匀,不会因为位置的不同而出现色彩过轻过重。如果图片中出现明显的色差变化,那么就认定这组数据质量不达标。

    1、 权重图
    GSE36000-权重图
    2、 残差图
    GSE36000-残差图
    3、 残差符号图
    GSE36000-残差符号图

    参考代码:

    # 2.6 质量分析与控制
    # 2.6.1 质量控制,化权重图、残差图等
    library(preprocessCore)
    library(gcrma)
    library(affyPLM)
    
    # 对探针数据集做线性拟合
    Pset <- fitPLM(affy.data)
    # 根据计算结果,画权重图
    pdf(file = "权重图.pdf")
    image(Pset, type="weights", which=1, main="Weights")
    dev.off()
    # 根据计算结果,画残差图
    pdf(file = "残差图.pdf")
    image(Pset,  type="resids", which=1, main="Residuals")
    dev.off()
    # 根据计算结果,画残差符号图
    pdf(file = "残差符号图.pdf")
    image(Pset, type="sign.resids", which=1, main="Residuals.sign")
    dev.off()
    
    3.3 质量分析报告

    在获取到了基因芯片的数据,那么为了保证数据质量达标,我们还需要使用到Bioconductor的simpleaffy包,使用qc的方法,对Affy芯片的数据,进行质量评估,如果质量评估太差,那么我们也需要放弃这一组数据。
    例如,本次实验使用到的数据(GSE36000),我们进行质量评估,可以发现,第一列,是我们在"读取下载部分的数据"章节,进行重命名的样本名称,第二列的数据是检出率与平均背景噪音,第三列(QC Stats),由actin(空心三角形)、(空心圆形)gapdh、(实心圆形)尺度因子组成。
    从质量评估图中,我们可以看出,这一芯片质量整体良好。指标出现蓝色表示正常,指标红色表示存在质量问题,比如unmut.8芯片相对来说actin远离尺度因子指标,那么也可以考虑弃用这类数据,如果标识有红色bioB表示该样品未检测到,那么这类数据建议不采用。
    另外,在采用GSE36000之前,本次实验原本是想通过基因芯片,GSE47811(已废弃),来研究健康小鼠与患癌小鼠的唾液间的关系。但是因为采用质量评估,发现大量数据不可用,所以弃用了GSE47811这组芯片。
    实验的芯片:这里我们可以剔除数据:mut3/unmut3,unmut8
    质量分析报告

    这里有个反例,是前实验弃用的芯片(已废弃——研究健康小鼠与患癌小鼠的唾液间,GSE47811)
    PS:这也是一个筛选数据的过程,如果不合适那么最好就抛弃,重新找更加合适的数据。
    反例-质量分析报告

    参考代码:

     2.6.2 将基因芯片进行解析读取,获取质量分析报告
    library(genefilter)
    library(gcrma)
    library(simpleaffy)
    affy.qc <- qc(affy.data)
    pdf(file = "质量分析报告.pdf", width = 10, height =20)
    # 图形化显示分析报告
    plot(affy.qc)
    # 剔除这些质量不合格的数据:mut3/unmut3,unmut8
    dev.off()
    
    3.4 RLE和NUSE箱线图

    RLE和NUSE,是质量检测的两种手段。
    RLE表达的是一个探针组在某个样品在探针组所有样品表达值的中位数取对数,而所有芯片如果都趋近于0值,那么表示这个芯片质量比较可靠;
    NUSE表达的是一个探针组在某个样品表达值的PM值标准差除以这个探针组所有样品表达值的PM值标准差的中位数,而所有芯片如果都趋近于1值,那么表示这个芯片质量比较可靠。
    以本次实验实验对象(GSE36000)为例,RLE趋向于0,NUSE趋向于1,芯片质量达标。但是如果芯片质量有问题,那么会严重偏离0值(RLE)、1值(NUSE),这个时候就需要进行舍弃了。
    RLE箱线图

    NUSE箱线图

    参考代码:

    # 2.6.3质量控制,绘制RLE和NUSE箱线图
    # 载入一组颜色
    library(RColorBrewer)
    library(graph)
    colors <- brewer.pal(12,"Set3")
    # 绘制RLE箱线图,如果数据有部分数据远离0值,那么就需要剔除这些质量不合格的数据,这里无须剔除
    Mbox(Pset,ylim=c(-1,1),col=colors,main="RLE",las=2)
    # 绘制NUSE箱线图,如果数据有部分数据远离1值,那么就需要剔除这些质量不合格的数据,这里无须剔除
    boxplot(Pset,ylim=c(0.95,1.2),col=colors,main="NUSE",las=2)
    
    3.5 质量分析(RNA降解)

    芯片质量分析,还有一种方式,那就是根据RNA降解,如果从左到右,保持一定的斜率(不与X值水平),那么表示RNA降解不严重,芯片质量可靠。以本次实验(GSE36000)为例,RNA降解曲线斜率高,没有呈现下降趋势,那么就表示芯片质量可靠。不需剔除数据。
    芯片质量分析-RNA降解

    参考代码:

    # 2.6.4 RNA降解,质量分析,如果图形中5'趋向3'斜率较低,那么就需要把这些不够明显区别的剔除,这里看出都明显有斜率,这里无须剔除
    affy.deg <- AffyRNAdeg(affy.data)
    plotAffyRNAdeg(affy.deg, col = colors)
    legend("topleft", rownames(pData(affy.data)), col = colors, lwd = 1, inset = 0.05, cex = 0.5)
    
    3.6 初步聚类分析

    从聚类分析的角度来看,未突变与突变组基本上能够比较好的分开,但这也不能判定这个实验是成立的。只能说明我们的研究对象,从未突变到突变有一些因素产生了作用,但是其中具体问题还是需要我们去具体分析。这里我们可以简单得出一个小结论,未突变到突变的基因,组间差异大于组内差异,基因差异较为明显。但这里我们不考虑去除不归类的样品。
    初步聚类分析

    4、 数据预处理

    数据的预处理,将质量分析筛选下来合格的数据,进行背景校正、标准化、汇总三个部分的处理,获取用来做实验分析的矩阵
    可以通过bgcorrect.methods()、normalize.methods(affy.data)、express.summary.stat.methods()的方式,获取能用于进行预处理的方法,如下图:
    数据预处理

    具体流程

    affy包提供了通过expresso方法来实现背景校正、标准化、汇总的预处理,这里需要设定不同的参数进行具体方法实现。
    另外,如果想更快的实现同样的功能,我们也可以直接通过采取mas5,gcrma,rma方法来进行一体化的预处理。这里mas5需要进行对数的转换,获取数字化的表达谱矩阵。

    经过上面那么多步骤,这里需要停下来,对数据进行一次处理(筛选时刻),参考代码:

    # 2.6.5 将以上的筛选重新组成实验组
    # GSM879118	MCL_IGHV_UNMUT(R1394)       mcl.ighv.unmut.1 
    # GSM879119	MCL_IGHV_UNMUT(R1328)       mcl.ighv.unmut.2
    # GSM879120	MCL_IGHV_UNMUT(R1329)       mcl.ighv.unmut.3    x
    # GSM879121	MCL_IGHV_UNMUT(R1306)       mcl.ighv.unmut.4
    # GSM879122	MCL_IGHV_UNMUT(R1332)       mcl.ighv.unmut.5
    # GSM879123	MCL_IGHV_MUT(R1399_M)       mcl.ighv.mut.1
    # GSM879124	MCL_IGHV_MUT(R1400)         mcl.ighv.mut.2
    # GSM879125	MCL_IGHV_UNMUT(R1587)       mcl.ighv.unmut.6
    # GSM879126	MCL_IGHV_UNMUT(R1680)       mcl.ighv.unmut.7 
    # GSM879127	MCL_IGHV_UNMUT(268-01-5TR)  mcl.ighv.unmut.8    x
    # GSM879128	MCL_IGHV_MUT(269-01-5TR)    mcl.ighv.mut.3      x
    # GSM879129	MCL_IGHV_UNMUT(043-01-4TR)  mcl.ighv.unmut.9
    # GSM879130	MCL_IGHV_MUT(R1302)         mcl.ighv.mut.4 
    # GSM879131	MCL_IGHV_MUT(R1304)         mcl.ighv.mut.5
    # GSM879132	MCL_IGHV_MUT(R1305)         mcl.ighv.mut.6
    # GSM879133	MCL_IGHV_MUT(R1628)         mcl.ighv.mut.7
    # GSM879134	MCL_IGHV_MUT(R1629)         mcl.ighv.mut.8
    # GSM879135	MCL_IGHV_MUT(R1341)         mcl.ighv.mut.9
    # GSM879136	MCL_IGHV_MUT(R1333)         mcl.ighv.mut.10
    # GSM879137	MCL_IGHV_MUT(R1585)         mcl.ighv.mut.11
    # GSM879138	MCL_IGHV_MUT(R1589)         mcl.ighv.mut.12
    # GSM879139	MCL_IGHV_MUT(R1591)         mcl.ighv.mut.13
    # GSM879140	MCL_IGHV_UNMUT(R1595)       mcl.ighv.unmut.10
    # GSM879141	MCL_IGHV_MUT(R1640)         mcl.ighv.mut.14
    # GSM879142	MCL_IGHV_MUT(R1626)         mcl.ighv.mut.15
    # GSM879143	MCL_IGHV_MUT(R1610)         mcl.ighv.mut.16
    # GSM879144	MCL_IGHV_MUT(02_201)        mcl.ighv.mut.17
    # GSM879145	MCL_IGHV_MUT(R1339)         mcl.ighv.mut.18
    # GSM879146	MCL_IGHV_MUT(R1301)         mcl.ighv.mut.19
    # GSM879147	MCL_IGHV_UNMUT(R1338)       mcl.ighv.unmut.11
    # GSM879148	MCL_IGHV_UNMUT(R1762)       mcl.ighv.unmut.12
    # GSM879149	MCL_IGHV_MUT(R1938)         mcl.ighv.mut.20
    # GSM879150	MCL_IGHV_MUT(R1940)         mcl.ighv.mut.21
    # GSM879151	MCL_IGHV_MUT(R1942)         mcl.ighv.mut.22
    # GSM879152	MCL_IGHV_MUT(R1941)         mcl.ighv.mut.23
    # GSM879153	MCL_IGHV_MUT(R1970)         mcl.ighv.mut.24
    # GSM879154	MCL_IGHV_UNMUT(R1897)       mcl.ighv.unmut.13
    # GSM879155	MCL_IGHV_UNMUT(R1388)       mcl.ighv.unmut.14
    
    
    c.experiment  <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.4","mcl.ighv.unmut.5","mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.9","mcl.ighv.unmut.10",
                       "mcl.ighv.unmut.11","mcl.ighv.unmut.12", "mcl.ighv.unmut.13","mcl.ighv.unmut.14",
                       "mcl.ighv.mut.1","mcl.ighv.mut.2","mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10",
                       "mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13","mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
                       "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24"
                       )
    
    c.experiment.status <- c("unmut","unmut","unmut","unmut","unmut","unmut",
                      "unmut","unmut","unmut","unmut","unmut","unmut",
                      "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                      "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                      "mut","mut", "mut"
                      )
    length(c.experiment)
    length(c.experiment.status)
    affy.experiment <- affy.data[,c.experiment]
    colnames(affy.experiment)
    
    # 2.7 将原始数据组成数据框,对数据进行一个分类的整理(为差异分析做准备)
    affy.experiment.frame <- data.frame(SampleID = c.experiment, Disease = c.experiment.status)
    
    sampleNames(affy.experiment)
    
    c.unmut <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.4","mcl.ighv.unmut.5","mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.9","mcl.ighv.unmut.10",
                 "mcl.ighv.unmut.11","mcl.ighv.unmut.12", "mcl.ighv.unmut.13","mcl.ighv.unmut.14")
    c.mut <- c("mcl.ighv.mut.1","mcl.ighv.mut.2","mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10",
               "mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13","mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
               "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24")
    length(c.unmut)
    length(c.mut)
    
    
    affy.experiment.unmut <- affy.experiment[,c.unmut]
    affy.experiment.mut <- affy.experiment[,c.mut]
    
    # 2.8 预处理,(预处理)背景校正、标准化、汇总
    # 预处理的方式有ms、expresso、gcrma
    bgcorrect.methods()
    normalize.methods(affy.experiment)
    express.summary.stat.methods()
    
    # 2.8.1 通过聚类分析,查看数据质量,预处理
    # 使用gcrma算法来预处理数据,也可以通过mas5、rma算法
    affy.experiment.gcrma <- gcrma(affy.experiment)
    affy.experiment.gcrma.exprs <- exprs(affy.experiment.gcrma)
    
    # 计算样品两两之间的Pearson相关系数
    pearson_cor <- cor(affy.experiment.gcrma.exprs)
    # 得到Pearson距离的下三角矩阵
    dist.lower <- as.dist(1 - pearson_cor)
    # 聚类分析、画图
    hc <- hclust(dist.lower,"ave")
    
    pdf(file = "聚类分析.pdf", width = 10, height =10)
    plot(hc)
    dev.off()
    
    # PCA
    # samplenames <- sub(pattern = "\\.CEL", replacement = "", colnames(affy.data.gcrma.eset))
    groups <- factor(affy.experiment.frame[,2])
    # BiocManager::install("affycoretools")
    # BiocManager::install("affycoretools")
    library("affycoretools")
    affycoretools::plotPCA(affy.experiment.gcrma.exprs, addtext=c.experiment, groups=groups, groupnames=levels(c.experiment.status))
    
    # 2.8.2 使用expresso,进行背景校正、标准化处理和汇总
    affy.experiment.expresso.eset <- expresso(affy.experiment, bgcorrect.method = "mas", normalize.method = "constant", pmcorrect.method = "mas",
                                              summary.method = "mas")
    affy.experiment.expresso.exprs = exprs(affy.experiment.expresso.eset) 
    
    # 2.8.3 直接采用MAS5算法进行数据预处理
    affy.experiment.mas5.eset = mas5(affy.experiment)  
    # 用exprs()从eset中获取数字化的表达谱矩阵
    affy.experiment.mas5.exprs = exprs(affy.experiment.mas5.eset) 
    

    5、 基因芯片差异分析

    5.1 limma处理差异分析

    limma是一个比较全面的,通过用于芯片、RNA-Seq、PCR进行线性模型和差异表达函数的工具包,它通过贝叶斯方法提供稳定的分析结果,limma支持多重复杂的设计实验,我们可以进行多个差异条件进行分析,差异条件越多,就越复杂。除了limma进行差异分析,也可以通过DESeq来集成处理。

    例如,本次实验(GES36000)主要是对未突变和突变两种差异条件进行差异分析,是比较简单的情况。

    使用limma进行差异分析,主要需要准备好:分组矩阵,差异比较矩阵,表达矩阵,进行线性模拟拟合与差值计算、贝叶斯检验、生成检验结果。

    具体流程
    1)获取实验设计矩阵,即未突变、突变两种情况
    实验设计矩阵-突变、未突变
    实验设计矩阵-突变、未突变2

    2)构建对比模型,比较两个实验条件下表达数据,本次实验只有两种差异条件,如果需要读个差异条件,通过makeContrasts进行多对多特征交互比较,得到差异
    构建对比模型
    3)获取表达矩阵,矩阵要求进行过预处理,并为对数转换的的表达谱矩阵
    4)对表达矩阵与实验设计矩阵进行线性模拟拟合
    5)根据对比模型进行差值计算
    6)对差值计算的结果进行贝叶斯检验
    7)生成检验结果的报告
    8)对P.value进行筛选,得到比较理想的差异表达基因
    未完全处理的差异表达基因

    PS:当然也可以直接使用GEO2R做一些处理,后面再做分享吧。

    参考代码:

    # 基因芯片差异分析
    # 3.1 差异分析方法:limma
    # 3.1.1 选取差异表达基因,使用Bioconductor中的limma包
    # 将(未突变、突变)状态数据框数据读取为因子
    library(limma)
    status <- factor(affy.experiment.frame[, "Disease"])
    # 构建实验设计矩阵(即未突变、突变两种差异对比)
    design <- model.matrix(~-1+status)
    head(design)
    # 构建对比模型,比较两个实验条件下表达数据,其中contrasts可以参考design的列名
    contrast.matrix <- makeContrasts(contrasts = "statusmut - statusunmut", levels=design)
    
    # 3.1.2 线性模型拟合
    fit <- lmFit(affy.experiment.gcrma.exprs , design)
    colnames(affy.experiment.gcrma.exprs)
    # 根据对比模型进行差值计算 
    fit1 <- contrasts.fit(fit, contrast.matrix)
    # 3.1.3 贝叶斯检验
    fit2 <- eBayes(fit1)
    View(fit2)
    
    # 3.1.4 生成所有基因的检验结果报告
    result <- topTable(fit2, coef="statusmut - statusunmut", n=nrow(fit2), lfc=log2(2))
    nrow(result)  # 1446
    # 用P.Value进行筛选,得到全部差异表达基因
    result <- result[result[, "P.Value"]<0.01,]
    result2 <- result[result[, "P.Value"]<0.001,]
    # 显示一部分报告结果
    head(result)
    head(result2)
    # 分析不同p.value值下,解析出来匹配的数据行数
    nrow(result)  # 909
    nrow(result2) # 701
    write.table(result, file="statusmut_statusunmut+0.01.txt", quote=F, sep="\t")
    write.table(result2, file="statusmut_statusunmut+0.001.txt", quote=F, sep="\t")
    
    5.2 寻找在不同条件下的差异表达

    另外我们可以采用倍数法、t-test方法来进行基因的差异表达。
    affy芯片中对探针杂交质量评定有三个值P/M/A(有/临界值/没有),分别是Present/Marginal/Absent。如果探针集被标记位P或者A,表示一个芯片中,基因有表达,PM数值与MM相比,没有显著差异。

    具体流程
    1)对表达谱矩阵进行对数化处理(针对mas5进行预处理的情况)
    2)计算样品组与实验组的平均表达值,计算差异的相对表达量
    3)使用倍数法,筛选出2倍以上差异的基因
    4)使用t-test检验差异基因(注意:根据 ?test方式,读取文档中提及,两个实验组需要保证,相同对比组,即样品组与实验组需要同数量)
    5)计算p值
    6)检验每张芯片中的三个值P/M/A(Present/Marginal/Absent),筛选出至少一个有表达的部分
    7)设置FDR(False Discovery Rate)方法,来校正P值
    8)筛选出2倍以上差异,P值小于0.05的,至少一个有表达的基因
    9)构造差异表达基因的热图
    差异表达基因的热图

    参考代码

    # 3.2 寻找在不同条件下存在差异表达的基因
    # 3.2.1 对表达谱矩阵进行对数化处理
    affy.experiment.mas5.exprs.log <- log(affy.experiment.mas5.exprs, 2)
    probeset.id <- row.names(affy.experiment.mas5.exprs.log)
    # 保存数据
    write.table(affy.experiment.mas5.exprs.log, file="affy.experiment.mas5.exprs.log.txt", quote=F, sep="\t")
    
    # 3.2.2 计算每个基因在样本中的平均表达值
    unmut.mean <- apply(affy.experiment.mas5.exprs.log[, c.unmut], 1, mean)
    mut.mean <- apply(affy.experiment.mas5.exprs.log[, c.mut], 1, mean)
    # 计算差异基因的相对表达量 (log(fold change))
    unmut.to.mut.mean <- unmut.mean - mut.mean
    
    # 使用cbind,将数据表拼接在一起
    unmut.to.mut.mean.data <- cbind(affy.experiment.mas5.exprs.log, unmut.mean, mut.mean, unmut.to.mut.mean)
    head(unmut.to.mut.mean.data)
    # 保存数据
    write.table(unmut.to.mut.mean.data, file="unmut.to.mut.mean.data.txt", quote=F, sep="\t")
    
    # 3.2.3 倍数法,寻找在不同条件下表达量存在2倍以上差异的基因(探针组) (log(fold change)>1 or <-1)
    unmut.to.mut.mean.fc.probesets <- names(unmut.to.mut.mean[unmut.to.mut.mean >1 | unmut.to.mut.mean<(-1)])
    
    # 3.2.4 用t test检验某个基因在突变与未突变基因中的表达是否存在差异
    dataset.unmut <- affy.experiment.mas5.exprs.log[1, c.unmut]
    dataset.mut <- affy.experiment.mas5.exprs.log[1, c.mut]
    test.gene <- t.test(dataset.unmut, dataset.mut, "two.sided")
    test.gene
    test.gene$p.value
    length(dataset.mut)
    length(dataset.unmut)
    ?t.test
    
    # 3.2.5 使用apply函数,对t test检验每个基因在未突变与突变基因中的表达是否存在差异,计算p值
    dim(affy.experiment.mas5.exprs.log)
    p.value.experiment.genes <- apply(affy.experiment.mas5.exprs.log, 1, function(x) { t.test(x[1:12], x[13:35]) $p.value } )
    
    # 3.2.6 用mas5calls函数来检测每张芯片中每个基因是否表达,P/M/A(有/临界值/没有)
    affy.experiment.mas5calls <- mas5calls(affy.experiment)
    affy.experiment.mas5calls.exprs = exprs(affy.experiment.mas5calls)
    # 保存数据
    write.table(affy.experiment.mas5calls.exprs, file="affy.experiment.mas5calls.exprs.txt", quote=F, sep="\t")
    affy.experiment.PMA <- apply(affy.experiment.mas5calls.exprs, 1, paste, collapse="")
    head(affy.experiment.PMA)
    
    # 3.2.7 把至少在一个芯片中有表达的基因选出来
    genes.present <- names(affy.experiment.PMA[affy.experiment.PMA != "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"])
    
    # 54675
    length(affy.experiment.PMA)
    # 40673
    length(genes.present)
    affy.experiment.mas5.exprs.log.present <- affy.experiment.mas5.exprs.log[genes.present,]
    
    # 3.2.8 设置FDR(False Discovery Rate)方法,来校正P值
    p.adjust.methods
    raw.pvals.present <- p.value.experiment.genes[genes.present]
    fdr.pvals.present <- p.adjust(raw.pvals.present, method="fdr")
    
    # 3.2.9 对FDR p值按从小到大排序
    fdr.pvals.present.sorted <- fdr.pvals.present[order(fdr.pvals.present)]
    length(fdr.pvals.present.sorted)
    # 输出FDR p值最小的10个基因
    fdr.pvals.present.sorted[1:10]
    fdr.pvals.present.sorted[40663:40673]
    # 将结果整理成表格
    expression.plus.pvals <- cbind(affy.experiment.mas5.exprs.log.present, raw.pvals.present, fdr.pvals.present)
    write.table(expression.plus.pvals, "mas5_expression.plus.pvals.txt", sep="\t", quote=F)
    
    # 3.2.10 找出那些至少在一个芯片中有表达,且在不同组织中表达有差异(FDR p值<0.23)的基因(探针组)
    DE.fdr.probesets <- names(fdr.pvals.present[fdr.pvals.present < 0.05])
    DE.fc.probesets <- names(unmut.to.mut.mean[unmut.to.mut.mean >1 | unmut.to.mut.mean<(-1)])
    DE.probesets <- intersect(DE.fdr.probesets,DE.fc.probesets)
    
    # 获取那些至少在一个芯片中有表达且在不同组织中表达有差异的基因的表达量
    unmut.to.mut.mean.mas5_DE <- unmut.to.mut.mean.data[DE.probesets, c.experiment]
    
    # 3.2.11 构造热图
    library(pheatmap)
    pdf(file = "差异表达基因热图.pdf",width = 14,height = 35)
    # pheatmap(unmut.to.mut.mean.mas5_DE,color=colorRampPalette(c("green","black","red"))(100),clustering_distance_rows = "correlation",clustering_distance_cols = "euclidean",clustering_method="complete")
    pheatmap(unmut.to.mut.mean.mas5_DE,color=colorRampPalette(c("green","black","red"))(100),
             clustering_distance_rows = "correlation",clustering_distance_cols = "euclidean",clustering_method="complete")
    dev.off()
    
    rnames.mas5_DE <- as.matrix(rownames(unmut.to.mut.mean.mas5_DE))
    colnames(rnames.mas5_DE) <- "probeset.id"
    
    affy.experiment.mas5.DE.mat <- cbind(rnames.mas5_DE,unmut.to.mut.mean.mas5_DE)
    #保存为txt文件
    write.table(affy.experiment.mas5.DE.mat, "unmut.to.mut.mean.mas5_DE.txt", sep="\t", quote=F,row.names=FALSE)
    
    5.3 小结

    本次实验中,对P.value的一些感想:
    P.value表示假设成立的条件下,重复n次试验,获得现有统计量以及极端情况下的概率。在样本量较低,使用P.value值可靠性会比较低,因为差异分为组内差异与组间差异(样品组与实验组),如果样本量少,那么组内差异会比较大,那么就会影响到我们需要关注的组间差异。
    P.value代表着阈值,我们可以提高阈值的方式(即将P.value<0.05转为P.value<0.01),虽然这降低了假阳性的概率,不过通过这种方式设定标准,会让原本可能真的存在表达差异的数据达不到我们设定的0.01的阈值标准,这就会提高了假阴性率。
    芯片数据分析,热图,是直观展示出基因表达变化的二维土层,其中数值以红绿黑为基本色调,能够让繁杂的数据一目了然,毕竟图片比文字更能够直观的表达。
    在本次差异分析中,GES36000样本,右侧文本为基因编号,热图的每一行都是代表这一个基因,下方文本为样本编号,热图的每一列都是代表着一个样本,左侧为树状结构,代表着基因之间的相似度聚落关系,顶部也为树状结构,表示样本之间的相似度聚落关系。其中,红色表示上调,绿色表示下调
    这里我们可以看出,GES36000样本,突变与未突变的样本在一些基因上的差异表达还是比较明显的。

    6、 基因功能分析

    6.1 注释工具包

    通过bioconductor收录的芯片注释包,将对应ID值进行映射,得到Gene symbol和Entrez ID两种探针注释信息,有了这些注释,能够直观的了解每个ID值对应的是什么术语?确定每个探针对应检测哪个基因的表达。

    在添加注释工具包时,可以先通过show的方式来自动下载对应的注释包,例如,本次实验中是引用了hgu133plus2.db注释包

    参考代码:

    # 4.基因功能分析
    # 4.1 注释工具包
    # 4.1.1 加载注释工具包
    # BiocManager::install(affydb)
    library(XML)
    library(annotate)
    library(org.Hs.eg.db)
    
    # 4.1.2 获得基因芯片注释包名称
    affy.experiment@annotation
    # 4.1.3 加载注释包"hgu133plus2.db"
    affy.experiment.affydb <- annPkgName(affy.experiment@annotation, type="db")
    library(affy.experiment.affydb, character.only=TRUE)
    # 根据每个探针组的ID,获取那些至少在一个芯片中有表达且在不同组织中表达有差异的基因名、Entrez ID
    unmut.to.mut.mean.symbols <- getSYMBOL(rownames(unmut.to.mut.mean.mas5_DE),affy.experiment.affydb)
    unmut.to.mut.mean.EntrezID <- getEG(rownames(unmut.to.mut.mean.mas5_DE),affy.experiment.affydb)
    unmut.to.mut.mean.mas5.DE.anno=cbind(unmut.to.mut.mean.EntrezID,unmut.to.mut.mean.symbols,unmut.to.mut.mean.mas5_DE)
    # 保存为txt文件
    write.table(unmut.to.mut.mean.mas5.DE.anno, "unmut.to.mut.mean.mas5.DE.anno.txt", sep="\t", quote=F,row.names = FALSE)
    

    根据挑选出来的差异基因,GO富集分析会对实验结果有补充提示的作用,我们可以通过GO富集分析差异基因,找到富集差异基因的GO分类目录,其中数据结果可以保存为本文,同时也可以保存为html的方式,这提供了较大的便利,我们可以通过这些方式,更加直观的寻找不同样品的差异基因可能和哪一些基因功能的改变有关系
    GO富集分析

    参考代码:

    # 5 GO富集分析
    # 加载所需R包
    library(Matrix)
    library(Category)
    library(graph)
    library(GOstats)
    library("GO.db")
    # 5.1 获取基因芯片所有探针组与差异表达基因的EntrezID
    # 提取芯片中affy.data所有探针组对应的EntrezID,注意保证uniq
    entrezUniverse <- unique(getEG(rownames(affy.experiment),affy.experiment.affydb));
    # 提取所有差异表达基因及其对应的EntrezID,去除na值,注意保证uniq
    entrezSelected <- unique(unmut.to.mut.mean.EntrezID[!is.na(unmut.to.mut.mean.EntrezID)]);
    # 设置GO富集分析的所有参数
    params <- new("GOHyperGParams", geneIds = entrezSelected, universeGeneIds = entrezUniverse, 
                  annotation = affy.experiment.affydb, ontology = "BP", pvalueCutoff = 0.001, conditional = FALSE, testDirection = "over");
    # 5.2 对所有的GO term根据params参数做超几何检验,并进行汇总
    hgOver <- hyperGTest(params);
    bp <- summary(hgOver) ;
    
    # 5.3 同时生成所有GO term的检验结果文件,每个GOterm都有指向官方网站的链接,可以获得其详细信息
    htmlReport(hgOver,  file='unmut.to.mut.mean.mas5_DE_go.html') ;
    head(bp)
    # 保存数据
    write.table(bp, "unmut.to.mut.mean.mas5_DE_go.txt", sep="\t", quote=F,row.names = FALSE)
    
    # 根据每个探针组的ID获取对应基因Gene Symbol,并作为新的一列
    result$symbols <- getSYMBOL(rownames(result), affy.experiment.affydb)
    # 根据探针ID获取对应基因Entrez ID
    result$EntrezID <- getEG(rownames(result), affy.experiment.affydb)
    # 显示前几行
    head(result)
    nrow(result)	
    

    四、结论

    通过在NCBI下载一套人类基因(未突变和突变淋巴),其中包含了38组样品基因,其中未突变14组,突变24组,前前后后对这38组数据编辑了5套实验脚本,进行数据的载入、读取、质量控制、质量分析、聚类分析、预处理、差异分析、注解、生成结果分析报告等。

    整一套流程走了几遍,主要收获如下:基因芯片的差异分析,主要分为数据的处理与分析;大多数人会侧重在分析上,也为了得到研究想要获得的结果;但是数据的筛选与质量分析也是一个很重要的部分。

    在本次研究之前(已废弃),是以研究健康小鼠与患癌小鼠的唾液间的关系,但因为基因质量不达标,选择放弃这组实验,另外在基因筛选上发现,如果样品组与实验组(突变与未突变)统计量少(38组数据,只人工采用20组数据),亦或是在p值设置阈值偏低,那么组内的差异很可能因为我们实验设计的不合理,而升高,导致实验结果十分不理想(前后5组实验,对比分析结果)。

    所以,我们在处理数据,进行基因差异分析,需要保证我们的结果的可靠性,那么就需要我们注重我们在分析时的严谨性,并反复测试不同实验变量,对比结果,得出合理的结论。

    五、参考资料

    [1] NCBI NCBI介绍[EB/OL]https://www.ncbi.nlm.nih.gov/home/about/
    [2] W3Cschool R介绍[EB/OL]https://www.w3cschool.cn/r/
    [3] 朱猛进.差异表达基因的分析[CP/OL]http://blog.sciencenet.cn/blog-295006-403640.html
    [4] 显著检验与多重检验校正[S/OL] http://www.omicshare.com/forum/thread-260-1-12.html
    [5] tommyhechina芯片数据分析(探针注释)[S/OL] https://blog.csdn.net/tommyhechina/article/details/80409983

    六、本文源码

    github的repo地址:https://github.com/PinkSmallFan/Pbioinformation

    对应本文的R源码:https://github.com/PinkSmallFan/Pbioinformation/tree/master/source_code/differential_analysis/Human_unmut_mut

    PS:本文产生的源数据(例如芯片灰度图、箱线图、差异分析、聚集分析、PCA等等的数据与对应的原图这里就不共享啦~ 需要的话可以留言,酌情考虑)
    PS2:如果这篇笔记有什么不足,或者疑惑不解的地方(可能对小白不太友好),可以提出来,看到会回答(留个TODO给大家~)

    展开全文
  • 迁移学习和传统机器学习在手写数字识别问题上的性能差异分析报告 工作概述: 运用matlab读取手写数字识别数据集,对其分别进行传统机器学习建模和迁移学习建模。对比分析两者性能差异。 工作计划: 准备数据集,进行...
  • 答:如果要用泊松分布做差异分析模型的话,必须要用 reads count 的。只有 RPKM值的话,可以用 RPKM 的公式反推 reads count 数,再做检验。 问 2:Deseq 是怎么控制 reads 多重比对的? 答:Deseq 只是一个差异...
  • 资深财务总监先后在三家生产型公司分别实施过用友U8、SAP/R3和Odoo后,总结三个ERP软件的设计理念和主要差异点。先看结论: 集团公司/多公司支持、云办公 用友ERP客户端登录界面 这个在U8上是完全不支持...
  • 差异分析流程(五)三大R包比较

    千次阅读 2019-12-19 20:39:12
    2. limma包进行差异分析 3. edgeR包进行差异分析 4. DESeq2包进行差异分析 5. 三大R包比较 1.输入数据比较 edgeR: 原始的count矩阵,支持单个样品和重复样品 DESeq2: 原始的count矩阵(htseq-count...
  • 用主成分分析法、标准差、变异系数、对数偏差均值指数等对甘肃省近20年来的区域经济差异进行了分析研究,得出以下结论:①从区域经济发展水平差异的空间格局上讲,定西、平凉、庆阳3市经济处于上升趋势,武威、临夏...
  • 通过国与国之间的贸易往来,经济溢出效应,跨国公司上市的逐步增加以及股票公司的全球化,从四个方面探讨了中美股票市场的相关方式,并试图分析中美股票市场的形成机理。中美股市之间的联系有四个渠道:国际贸易,...
  • 样本量大的时候做差异分析容易得到有显著性差异结论,原因如下图, 求p值的过程中,n越大,Z0也越大,相对应的p就小了。当然这里默认方差变化不大的情况下,因为一般来说很多数据经过平均后方差不会变化很大,...
  • 方法 采用野外调查、室内实验和检测分析岩石的地质特征、自然环境特征等进行对比分析。结果 红石峡以风蚀打磨字迹消失、模糊或者字迹多处粗糙为主要病害,千佛崖以彩绘造像表面酥粉和粉化脱落为主;千佛崖岩石吸水率...
  • 分析产品差异化的问题上引入合谋成本,修正了以往认为古诺竞争比伯川德竞争的合谋更容易维持的认识,得出在合谋成本的临界值之上,伯川德竞争比古诺竞争的合谋更容易维持的结论,这为经济参与者和决策制定者分析合谋...
  • 虽然主流的数据库系统中都提供了限制结果集行数的方法,但是无论是语法还是使用方式都存在着很大的差异。即使是同一个数据库系统的不同版本(比如MSSQLServer2000和MSSQLServer2005)也存在着一定的差异
  • 1. 主成分分析(principal components analysis,PCA)是一种常用的数据间差异分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征向量,常用于高维数据的降维。原理推荐...
  • 组间差异的非参数检验 若数据无法满足t检验或ANOVA的参数假设,可使用非参数方法 两组的比较 Mann-Whitney U检验 两组数据独立时使用。 用来判断一个总体中获得更高得分的概率是否比另一个总体要大。 > library...
  • rtklibexplorer的一篇文章该问题做了分析,文章Exploring differences between real-time and post-processed solutions.,这是一个有趣的事情,实验步骤稍后描述,我出现的问题也稍后描述,很有意思
  • 为研究克氏原螯虾在不同庇护所生境下的生长差异及其采取的不同生存策略,本实验在水池中设置洞穴、角落2种庇护所生境,并以无任何庇护所的水池为对照,饲养亚成体克氏原螯虾一个月。计算各组克氏原螯虾存活、生长的...
  • 主成分分析(PCA)原理详解

    万次阅读 多人点赞 2018-06-09 15:08:25
    “微信公众号”本文同步更新在我的微信公众号里,地址:https://mp.weixin.qq.com/s/Xt1vLQfB20rTmtLjiLsmww本文同步更新在我的知乎专栏里面:主成分分析(PCA)原理详解 - Microstrong的文章 - 知乎...
  • [结果和结论]结果表明:不同蛋白含量大豆各器官中GS基因家族不同成员的表达量具有明显的差异。GSβ1在根、茎、叶和根瘤中都能高效表达;叶中GS2表达量显著升高,超过其他器官和根瘤。GSγ1在根、茎、
  • 平时大家使用 epoll 时都知道其事件触发模式有默认的 level-trigger 模式和通过 EPOLLET 启用的 edge-trigger 模式两种。从 epoll 发展历史来看,它刚诞生时只有 edge-trigger 模式,...二者的差异在于 level-trigger
  • 新浪微博、微信朋友圈、...具体对比分析如下表所示: 纯文字输入:微信朋友圈的功能是最少的,qq空间发布时附加功能是最有趣的(如趣味动图,语音输入,个性字体),新浪微博是最偏向于媒体化的(如更多的种类标签...
  • 杰哥带着一种好奇心的想法,结合自身的工作经验与业界全国关于招聘运维工程师的岗位做一个初步型的分析,我的一位好朋友 —— 黄伟呢,帮我爬取了 13966 条关于运维的招聘信息,看看有哪些数据存在相关差异化。...
  • 测试结论

    千次阅读 2019-03-23 11:36:34
    测试结论 1.版本功能基本实现且运行稳定,问题修改及时,在预定日期 内完成开发和测试进度 质量评价:通过,可以发布及系统上线 测试结论: 通过,可以发布及系统上线 不通过,需要进行重大修改更新版本重新...
  • meta分析一般步骤

    万次阅读 多人点赞 2018-07-26 16:42:57
    对一些大样本,多中心临床合作已经得到明确结论的的,没必要做meta分析。 二、文献检索 在制定文献检索策略时,总体的要求就是查全和查准。 需要考虑如下几个方面: 1. 圈定搜索数据库(外文有:MEDLINE、the ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 55,658
精华内容 22,263
关键字:

差异分析结论