原先三个样本的HTSeq-count计数的数据可以在我的GitHub中找到，但是前面已经说过Jimmy失误让我们分析的人类就只有3个样本， 另外一个样本需要从另一批数据获取（请注意batch effect)，所以不能保证每一组都有两个重复。

我一直坚信”你并不孤独“这几个字，遇到这种情况的人肯定不止我一个，于是我找到了几种解决方法

以上方法都会在后续进行介绍，但是我们DESeq2必须得要有重复的问题亟待解决，没办法我只能自己瞎编了。虽然是编，我们也要有模有样，不能直接复制一份，要考虑到高通量测序的read是默认符合泊松分布的。我是这样编的。

这仅仅是一种填坑的方法而已，更好模拟数据的方法需要参阅更加专业的文献， 有生之年 我希望能补上这一个部分。

这部分内容最先在 RNA-Seq Data Analysis 的8.5.3节看到，刚开始一点都不理解，但是学完生物统计之后，我认为这是理解所有差异基因表达分析R包的关键。

基本上，统计课都会介绍如何使用 t检验 用来比较两个样本之间的差异，然后在样本比较多的时候使用 方差分析 确定样本间是否有差异。当然前是样本来自于正态分布的群体，或者随机独立大量抽样。

对于基因芯片的差异表达分析而言，由于普遍认为其数据是服从正态分布，因此差异表达分析无非就是用t检验和或者方差分析应用到每一个基因上。高通量一次性找的基因多，于是就需要对多重试验进行矫正，控制假阳性。目前在基因芯片的分析用的最多的就是 limma 。

但是 ，高通量测序(HTS)的read count普遍认为是服从泊松分布（当然有其他不同意见），不可能直接用正态分布的 t检验 和 方差分析 。 当然我们可以简单粗暴的使用对于的 非参数检验 的方法，但是统计力不够，结果的p值矫正之估计一个差异基因都找不到。老板花了一大笔钱，结果却说没有差异基因，是个负结果，于是好几千经费打了水漂，他肯定是不乐意的。因此，还是得要用参数检验的方法，于是就要说到方差分析和线性模型之间的关系了。

线性回归和方差分析是同一时期发展出的两套方法。在我本科阶段的田间统计学课程中就介绍用 方差分析 （ANOVA）分析不同肥料处理后的产量差异，实验设计如下

这是最简单的单因素方差分析，每一个结果都可以看成 yij = ai + u + eij， 其中u是总体均值，ai是每一个处理的差异，eij是随机误差。

注 ：方差分析(Analysis of Variance, ANAOVA)名字听起来好像是检验方差，但其实是为了判断样本之间的差异是否真实存在，为此需要证明不同处理内的方差显著性大于不同处理间的方差。

线性回归 一般是用于量化的预测变量来预测量化的响应变量。比如说体重与身高的关系建模：

当然线性回归也可用处理名义型或有序型因子（也就是离散变量）作为预测变量，如果要画图的话，就是下面这个情况。

如果我们需要通过一个实验找到不同处理后对照组和控制组的基因变化，那么基因表达可以简单写成， y = a + b · treament + e。 和之前的 yij = ai + u + eij 相比，你会发现公式是如此的一致。 这是因为线性模型和方差分析都是 广义线性模型 (generalizing linear models, GLM)在正态分布的预测变量的特殊形式。而GLM本身只要采用合适的 连接函数 是可以处理对任意类型的变量进行建模的。

目前认为read count之间的差异是符合负二项分布，也叫gamma-Possion分布。那么问题来了，如何用GLM或者LM分析两个处理件的差异呢？其实可以简单的用上图的拟合直线的斜率来解释，如果不同处理之间存在差异，那么这个拟合线的斜率必定不为零，也就是与X轴平行。但是这是一种便于理解的方式（虽然你也未必能理解），实际更加复杂，考虑因素更多。

注1 负二向分布有两个参数，均值(mean)和离散值（dispersion). 离散值描述方差偏离均值的程度。泊松分布可以认为是负二向分布的离散值为1，也就是均值等于方差（mean=variance）的情况。

注2 这部分涉及大量的统计学知识，不懂就用维基百科一个个查清楚。

聊完了线性模型和方差分析，下面的设计矩阵（design matrix）就很好理解了， 其实就是用来告诉不同的差异分析函数应该如何对待变量。比如说我们要研究的KD和control之间变化，设计矩阵就是

那么比较矩阵（contrast matrix)就是告诉差异分析函数应该如何对哪个因素进行比较， 这里就是比较不同处理下表达量的变化。

其实read count如何标准化的方法有很多，最常用的是FPKM和RPKM，虽然它们其实是错的-- FPKM/RPKM是错的 。

我推荐阅读 Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data , 了解不同标准化方法之间的差异。

有一些方法是要求原始数据，有一些则要求经过某类标准化后的数据，记得区分。

关于DESeq2分析差异表达基因，其实在 https://www.bioconductor.org/help/workflows/rnaseqGene/ 里面介绍的非常清楚了。

我们已经准备好了count matrix，接下来就是把数据导入DESeq2。DESeq2导入数据的方式有如下4种，基本覆盖了主流read count软件的结果。

注 DESeq2要求的数据是raw count， 没必要进行FPKM/TPM/RPFKM/TMM标准化。

本来我们是可以用DESeq2为htseq-count专门提供的 DESeqDataSetFromHTSeq ，然而很尴尬数据不够要自己凑数，所以只能改用 DESeqDataSetFromMatrix 了 :cold_sweat:

导入数据，构建 DESeq2 所需的 DESeqDataSet 对象

注 ： 这一步到下一步之间可以过滤掉一些low count数据，节省内存，提高运行速度

使用 DESeq 进行差异表达分析： DESeq 包含三步，estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons)， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest），可以分布运行，也可用一步到位，最后返回 results 可用的DESeqDataSet对象。

用results获取结果： results的参数非常的多，这里不好具体展开 :pensive: 但是你们会自己看的吧

我们可用mcols查看每一项结果的具体含义，比如说 log2FoldChange 表示倍数变化取log2结果，还能画个火山图。一般简单粗暴的用2到3倍作为阈值，但是对于低表达的基因，3倍也是噪音，那些高表达的基因，1.1倍都是生物学显著了。更重要的没有考虑到组内变异，没有统计学意义。 padj 就是用BH对多重试验进行矫正。

用summary看描述性的结果，大致是上调的基因占总体的11%，下调的是7.1%（KD vs control）

画个MA图，还能标注p值最小的基因。

下图是没有经过 statistical moderation平缓log2 fold changes的情况

如果经过 lfcShrink 收缩log2 fold change， 结果会好看很多

当然还有火山图，不过留给其他方法作图，我们先把差异表达的基因找出来。

一般p value 小于0.05就是显著了, 显著性不代表结果正确，只用于给后续的富集分析和GSEA提供排序标准和筛选而已。关于P值的吐槽简直无数， 请多注意。

edgeR在函数说明中称其不但可以分析SAGE， CAGE的RNA-Seq，Tag-RNA，或RNA-seq， 也能分析ChIP-Seq和CRISPR得到的read counts数据。嗯，我信了:confused:！

edgeR使用 DGEList 函数读取count matrix数据，也就说你需要提供一个现成的matrix数据，而不是指望它能读取单独的文件，然后进行合并(当然机智的我发现，其实可以用 tximport 或 DESeqDataSetFromHTSeq 读取单独的文件，然后传递给 DGEList )

第一步： 构建DGEList对象

第二步： 过滤 low counts数据。与DESeq2的预过滤不同，DESeq2的预过滤只是为了改善后续运算性能，在运行过程中依旧会自动处理low count数据，edgeR需要在分析前就要排除那些low count数据，而且非常严格。从生物学角度，有生物学意义的基因的表达量必须高于某一个阈值。从统计学角度上， low count的数据不太可能有显著性差异，而且在多重试验矫正阶段还会拖后腿。 综上所诉，放心大胆的过滤吧。

根据经验（又是经验 :dog: ）， 基因至少在某一些文库的count超过10 ~ 15 才被认为是表达。这一步全靠尝试， 剔除太多就缓缓，剔除太少就严格点。 我们可以简单的对每个基因的raw count进行比较，但是建议用CPM（count-per-million) 标准化 后再比较，避免了 文库大小 的影响。

这里的0.5(即阈值）等于 10/(最小的文库的 read count数 /1000000)，keep.lib.size=FALSE表示重新计算文库大小。

第三步： 根据组成偏好(composition bias)标准化。edgeR的 calcNormFactors 函数使用 TMM算法 对DGEList标准化

注 大部分的mRNA-Seq数据分析用TMM标准化就行了，但是也有例外，比如说single-cell RNA-Seq(Lun, Bach, and Marioni 2016), 还有就是global differential expression， 基因组一半以上的基因都是差异表达的，请尽力避免，(D. Wu et al. 2013)， 不然就需要用到内参进行标准化了(Risso et al. 2014).

第四步： 实验设计矩阵(Design matrix)， 类似于DESeq2中的design参数。 edgeR的线性模型和差异表达分析需要定义一个实验设计矩阵。很直白的就能发现是1vs0

第五步： 估计离散值（Dispersion）。前面已经提到负二项分布（negative binomial，NB)需要均值和离散值两个参数。edgeR对每个基因都估测一个经验贝叶斯稳健离散值（mpirical Bayes moderated dispersion），还有一个公共离散值（common dispersion，所有基因的经验贝叶斯稳健离散值的均值）以及一个趋势离散值

还可以进一步通过quasi-likelihood (QL)拟合NB模型，用于解释生物学和技术性导致的基因特异性变异 (Lund et al. 2012Lun, Chen, and Smyth 2016).

注1 估计离散值这个步骤其实有许多 estimate*Disp 函数。当不存在实验设计矩阵(design matrix)的时候， estimateDisp 等价于 estimateCommonDisp 和 estimateTagwiseDisp 。而当给定实验设计矩阵(design matrix)时， estimateDisp 等价于 estimateGLMCommonDisp , estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp 。 其中tag与gene同义。

注2 其实这里的第三， 四， 五步对应的就是DESeq2的 DESeq 包含的2步，标准化和离散值估测。

第六步: 差异表达检验（1）。这一步主要构建比较矩阵，类似于DESeq2中的 results 函数的 contrast 参数。

这里用的是 glmQLFTest 而不是 glmLRT 是因为前面用了glmQLTFit进行拟合，所以需要用QL F-test进行检验。如果前面用的是 glmFit ，那么对应的就是 glmLRT . 作者称QL F-test更加严格。多重试验矫正用的也是BH方法。

后续就是提取显著性差异的基因用作下游分析，做一些图看看

第六步：差异表达检验（2）。上面找到的显著性差异的基因，没有考虑效应值，也就是具体变化了多少倍。我们也可用找表达量变化比较大的基因，对应的函数是 glmTreat 。

经过上面两个方法的洗礼，基本上套路你也就知道了，我先简单小结一下，然后继续介绍limma包的 voom 。

Limma原先用于处理基因表达芯片数据，可是说是这个领域的老大 :sunglasses: 。如果你仔细看edgeR导入界面，你就会发现，edgeR有一部分功能依赖于limma包。Limma采用经验贝叶斯模型（ Empirical Bayesian model）让结果更稳健。

在处理RNA-Seq数据时，raw read count先被转成log2-counts-per-million (logCPM)，然后对mean-variance关系建模。建模有两种方法：

数据预处理 ： Limma使用edgeR的DGEList对象，并且过滤方法都是一致的，对应edgeR的第一步,第二步， 第三步

差异表达分析 : 使用”limma-trend“

差异表达分析 : 使用”limma-voom“

如果分析基因芯片数据，必须好好读懂LIMMA包。

基本上每一个包，我都提取了各种的显著性基因，比较就需要用韦恩图了，但是我偏不 :stuck_out_tongue: 我要用UpSetR.

感觉limma的结果有点奇怪，有生之年在折腾吧。

好吧，这部分我鸽了

[1] Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data

[2] https://www.bioconductor.org/help/workflows/rnaseqGene/

[3] https://www.bioconductor.org/help/workflows/RnaSeqGeneEdgeRQL/

[4] https://www.bioconductor.org/help/workflows/RNAseq123/

“参考网址1”中提到如果只是对整数运算（运算过程和结果都只使用整数），没有必要使用“double”(8 byte)，而应该用更小的“integer”(4 byte)。使用storage.mode(x)查看对象存数的模式，storage.mode(x) <- 进行赋值；使用format(object.size(a), units = 'auto')查看对象占用的内存空间（此处有疑问，即在R中每个integer到底占用了多大的空间？）。

需要解释gc()函数，可以查看内存使用情况。同样，在清除了大的对象之后，使用gc()以释放内存使用空间。

李航在”参考网址2“中提到，对于大矩阵的操作，尽量避免使用cbind和rbind之类，因为这会让内存不停地分配空间。“对于长度增加的矩阵，尽量先定义一个大矩阵，然后逐步增加”和“注意清除中间对象”。

使用bigmemory家族：bigmemory, biganalytics, synchronicity, bigtabulate and bigalgebra， 同时还有

biglm。

bigmemory package的使用：

1. 建立big.memory对象

bigmemory采用C++的数据格式来“模仿”R中的matrix。

编写大数据格式文件时候，可以先建立filebacked.big.matrix

big.matrix(nrow, ncol, type = options()$bigmemory.default.type, init = NULL, dimnames = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared = TRUE)

filebacked.big.matrix(nrow, ncol, type = options()$bigmemory.default.type, init = NULL, dimnames = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL)

as.big.matrix(x, type = NULL, separated = FALSE, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared=TRUE)

使用注意：

big.matrix采用两种方式储存数据：一种是big.matrix默认的方式，如果内存空间比较大，可以尝试使用；另外一种是filebacked.big.matrix，这种储存方法可能会备份文件（file-backings），而且需要descriptor file；

“init”指矩阵的初始化数值，如果设定，会事先将设定的数值填充到矩阵中；如果不设置，将处理为NA

"type"是指在big.matrix中atomic element的储存格式，默认是“double”(8 byte)，可以改为“integer”(4 byte), "short"(2 byte) or "char"(1 byte)。注意：这个包不支持字符串的储存，type = "char"是指ASCII码字母。

在big.matrix非常大的时候，避免使用rownames和colnames(并且bigmemory禁止用名称访问元素)，因为这种做法非常占用内存。如果一定要改变，使用options(bigmemory.allow.dimnames=TRUE)，之后colnames, rownames设置。

直接在命令提示符后输入x（x是一个big matrix），将返回x的描述，不会出现所有x中所有内容。因此，注意x[ , ](打印出矩阵全部内容)；

如果big.matrix有很多列，那么应该将其转置后储存；（不推荐）或者将参数“separated”设置为TRUE，这样就将每一列分开储存。否则，将用R的传统方式（column major的方式）储存数据。

如果建立一个filebacked.big.matrix，那么需要指定backingfile的名称和路径+descriptorfile。可能多个big.matrix对象对应唯一一个descriptorfile，即如果descriptorfile改变，所以对应的big.matrix随之改变；同样，decriptorfile随着big.matrix的改变而改变；如果想维持一种改变，需要重新建立一个filebacked.big.matrix。attach.big.matrix(descriptorfile or describe(big.matrix))函数用于将一个descriptorfile赋值给一个big.matrix。这个函数很好用，因为每次在创建一个filebacked.big.matrix后，保存R并退出后，先前创建的矩阵会消失，需要再attach.big.matrix以下

2. 对big.matrix的列的特定元素进行条件筛选

对内存没有限制；而且比传统的which更加灵活（赞！）

mwhich(x, cols, vals, comps, op = 'AND')

x既可以是big.matrix，也可以是传统的R对象；

cols：行数

vals：cutoff，可以设定两个比如c(1, 2)

comps：'eq'(==), 'neq'(!=), 'le'(<), 'lt'(<=), 'ge'(>) and 'gt'(>=)

op：“AND”或者是“OR”

可以直接比较NA，Inf和-Inf

3.bigmemory中其他函数

nrow, ncol, dim, dimnames, tail, head, typeof继承base包

big.matrix, is.big.matrix, as.big.matrix, attach.big.matrix, describe, read.big.matrix, write.big.matrix, sub.big.matrix, is.sub.big.matrix为特有的big.matrix文件操作；filebacked.big.matrix, is.filebacked（判断big.matrix是否硬盘备份） , flush(将filebacked的文件刷新到硬盘备份上)是filebacked的big.matrix的操作。

mwhich增强base包中的which， morder增强order，mpermute（对matrix中的一列按照特定序列操作，但是会改变原来对象，这是为了避免内存溢出）

big.matrix对象的copy使用deepcopy(x, cols = NULL, rows = NULL, y = NULL, type = NULL, separated = NULL, backingfile = NULL, backingpath = NULL, descriptorfile = NULL, shared=TRUE)

biganalytics package的使用

biganalytics主要是一些base基本函数的扩展，主要有max, min, prod, sum, range, colmin, colmax, colsum, colprod, colmean, colsd, colvar, summary, apply（只能用于行或者列，不能用行列同时用）等

比较有特色的是bigkmeans的聚类

剩下的biglm.big.matrix和bigglm.big.matrix可以参考Lumley's biglm package。

bigtabulate package的使用

并行计算限制的突破：

使用doMC家族：doMC, doSNOW, doMPI, doRedis, doSMP和foreach packages.

foreach package的使用

foreach(..., .combine, .init, .final=NULL, .inorder=TRUE, .multicombine=FALSE, .maxcombine=if (.multicombine) 100 else 2, .errorhandling=c('stop', 'remove', 'pass'), .packages=NULL, .export=NULL, .noexport=NULL, .verbose=FALSE)

foreach的特点是可以进行并行运算，如在NetWorkSpace和snow？

%do%严格按照顺序执行任务（所以，也就非并行计算），%dopar%并行执行任务

...：指定循环的次数；

.combine：运算之后结果的显示方式，default是list，“c”返回vector， cbind和rbind返回矩阵，"+"和"*"可以返回rbind之后的“+”或者“*”

.init：.combine函数的第一个变量

.final：返回最后结果

.inorder：TRUE则返回和原始输入相同顺序的结果（对结果的顺序要求严格的时候），FALSE返回没有顺序的结果（可以提高运算效率）。这个参数适合于设定对结果顺序没有需求的情况。

.muticombine：设定.combine函数的传递参数，default是FALSE表示其参数是2，TRUE可以设定多个参数

.maxcombine：设定.combine的最大参数

.errorhandling：如果循环中出现错误，对错误的处理方法

.packages：指定在%dopar%运算过程中依赖的package（%do%会忽略这个选项）。

getDoParWorkers( ) ：查看注册了多少个核，配合doMC package中的registerDoMC( )使用

getDoParRegistered( ) ：查看doPar是否注册；如果没有注册返回FALSE

getDoParName( ) ：查看已经注册的doPar的名字

getDoParVersion( )：查看已经注册的doPar的version

===================================================

# foreach的循环次数可以指定多个变量，但是只用其中最少？的

>foreach(a = 1:10, b = rep(10, 3)) %do% (a*b)

[[1]]

[1] 10

[[2]]

[1] 20

[[3]]

[1] 30

# foreach中.combine的“+”或者“*”是cbind之后的操作；这也就是说"expression"返回一个向量，会对向量+或者*

>foreach(i = 1:4, .combine = "+") %do% 2

[1] 8

>foreach(i = 1:4, .combine = "rbind") %do% rep(2, 5)

[,1] [,2] [,3] [,4] [,5]

result.122222

result.222222

result.322222

result.422222

>foreach(i = 1:4, .combine = "+") %do% rep(2, 5)

[1] 8 8 8 8 8

>foreach(i = 1:4, .combine = "*") %do% rep(2, 5)

[1] 16 16 16 16 16

=============================================

iterators package的使用

iterators是为了给foreach提供循环变量，每次定义一个iterator，它都内定了“循环次数”和“每次循环返回的值”，因此非常适合结合foreach的使用。

iter(obj, ...)：可以接受iter, vector, matrix, data.frame, function。

nextElem(obj, ...)：接受iter对象，显示对象数值。

以matrix为例，

iter(obj, by=c('column', 'cell', 'row'), chunksize=1L, checkFunc=function(...) TRUE, recycle=FALSE, ...)

by：按照什么顺序循环；matrix和data.frame都默认是“row”，“cell”是按列依次输出（所以对于“cell”，chunksize只能指定为默认值，即1）

chunksize：每次执行函数nextElem后，按照by的设定返回结果的长度。如果返回结构不够，将取剩余的全部。

checkFunc=function(...) TRUE：执行函数checkFun，如果返回TRUE，则返回；否则，跳过。

recycle：设定在nextElem循环到底（“错误: StopIteration”）是否要循环处理，即从头再来一遍。

以function为例

iter(function()rnorm(1))，使用nextElem可以无限重复；但是iter(rnorm(1))，只能来一下。

更有意思的是对象如果是iter，即test1 <- iter(obj)test2 <- iter(test1)，那么这两个对象是连在一起的，同时变化。

==============================================

>a

[,1] [,2] [,3] [,4] [,5]

[1,]159 13 17

[2,]26 10 14 18

[3,]37 11 15 19

[4,]48 12 16 20

>i2 <- iter(a, by = "row", chunksize=3)

>nextElem(i2)

[,1] [,2] [,3] [,4] [,5]

[1,]159 13 17

[2,]26 10 14 18

[3,]37 11 15 19

>nextElem(i2) #第二次iterate之后，只剩下1行，全部返回

[,1] [,2] [,3] [,4] [,5]

[1,]48 12 16 20

>i2 <- iter(a, by = "column", checkFunc=function(x) sum(x) >50)

>nextElem(i2)

[,1]

[1,] 13

[2,] 14

[3,] 15

[4,] 16

>nextElem(i2)

[,1]

[1,] 17

[2,] 18

[3,] 19

[4,] 20

>nextElem(i2)

错误: StopIteration

>colSums(a)

[1] 10 26 42 58 74

>testFun <- function(x){return(x+2)}

>i2 <- iter(function()testFun(1))

>nextElem(i2)

[1] 3

>nextElem(i2)

[1] 3

>nextElem(i2)

[1] 3

>i2 <- iter(testFun(1))

>nextElem(i2)

[1] 3

>nextElem(i2)

错误: StopIteration

>i2 <- iter(testFun(1))

>i3 <- iter(i2)

>nextElem(i3)

[1] 3

>nextElem(i2)

错误: StopIteration

============================================

iterators package中包括

irnorm(..., count)；irunif(..., count)；irbinom(..., count)；irnbinom(..., count)；irpois(..., count)中内部生成iterator的工具，分别表示从normal，uniform，binomial，negativity binomial和Poisson分布中随机选取N个元素，进行count次。其中，negative binomial分布：其概率积累函数(probability mass function)为掷骰子，每次骰子为3点的概率为p，在第r+k次恰好出现r次的概率。

icount(count)可以生成1:conunt的iterator；如果count不指定，将从无休止生成1:Inf

icountn(vn)比较好玩，vn是指一个数值向量（如果是小数，则向后一个数取整，比如2.3 -->3）。循环次数为prod(vn)，每次返回的向量中每个元素都从1开始，不超过设定 vn，变化速率从左向右依次递增。

idiv(n, ..., chunks, chunkSize)返回截取从1:n的片段长度，“chunks”和“chunkSize”不能同时指定，“chunks”为分多少片段（长度从大到小），“chunkSize”为分段的最大长度（长度由大到小）

iapply(X, MARGIN)：与apply很像，MARGIN中1是row，2是column

isplit(x, f, drop=FALSE, ...)：按照指定的f划分矩阵

=============================================

>i2 <- icountn(c(3.4, 1.2))

>nextElem(i2)

[1] 1 1

>nextElem(i2)

[1] 2 1

>nextElem(i2)

[1] 3 1

>nextElem(i2)

[1] 4 1

>nextElem(i2)

[1] 1 2

>nextElem(i2)

[1] 2 2

>nextElem(i2)

[1] 3 2

>nextElem(i2)

[1] 4 2

>nextElem(i2)

错误: StopIteration

用R语言进行关联分析关联是两个或多个变量取值之间存在的一类重要的可被发现的某种规律性。关联分析目的是寻找给定数据记录集中数据项之间隐藏的关联关系，描述数据之间的密切度。几个基本概念1. 项集这是一个集合的概念，在一篮子商品中的一件消费品即为一项（Item），则若干项的集合为项集，如{啤酒，尿布}构成一个二元项集。2. 关联规则一般记为的形式，X为先决条件，Y为相应的关联结果，用于表示数据内隐含的关联性。如：，表示购买了尿布的消费者往往也会购买啤酒。关联性强度如何，由三个概念——支持度、置信度、提升度来控制和评价。例：有10000个消费者购买了商品，其中购买尿布1000个，购买啤酒2000个，购买面包500个，同时购买尿布和面包800个，同时购买尿布和面包100个。3. 支持度（Support）支持度是指在所有项集中{X, Y}出现的可能性，即项集中同时含有X和Y的概率：该指标作为建立强关联规则的第一个门槛，衡量了所考察关联规则在“量”上的多少。通过设定最小阈值（minsup），剔除“出镜率”较低的无意义规则，保留出现较为频繁的项集所隐含的规则。设定最小阈值为5%，由于{尿布，啤酒}的支持度为800/10000=8%，满足基本输了要求，成为频繁项集，保留规则；而{尿布，面包}的支持度为100/10000=1%，被剔除。4. 置信度（Confidence）置信度表示在先决条件X发生的条件下，关联结果Y发生的概率：这是生成强关联规则的第二个门槛，衡量了所考察的关联规则在“质”上的可靠性。相似的，我们需要对置信度设定最小阈值（mincon）来实现进一步筛选。具体的，当设定置信度的最小阈值为70%时，置信度为800/1000=80%，而的置信度为800/2000=40%，被剔除。5. 提升度（lift）提升度表示在含有X的条件下同时含有Y的可能性与没有X这个条件下项集中含有Y的可能性之比：该指标与置信度同样衡量规则的可靠性，可以看作是置信度的一种互补指标。R中Apriori算法算法步骤：1. 选出满足支持度最小阈值的所有项集，即频繁项集；2. 从频繁项集中找出满足最小置信度的所有规则。>library(arules) #加载arules包>click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1)) #读取txt文档（文档编码为ANSI）>rules <- apriori(click_detail, parameter =list(supp=0.01,conf=0.5,target="rules")) #调用apriori算法>rulesset of419 rules>inspect(rules[1:10]) #查看前十条规则解释1)library(arules)：加载程序包arules，当然如果你前面没有下载过这个包，就要先install.packages(arules)2)click_detail =read.transactions("click_detail.txt",format="basket",sep=",",cols=c(1))：读入数据read.transactions(file, format =c("basket", "single"), sep = NULL,cols = NULL, rm.duplicates =FALSE, encoding = "unknown")file:文件名，对应click_detail中的“click_detail.txt”format:文件格式，可以有两种，分别为“basket”,“single”，click_detail.txt中用的是basket。basket: basket就是篮子，一个顾客买的东西都放到同一个篮子，所有顾客的transactions就是一个个篮子的组合结果。如下形式，每条交易都是独立的。文件形式：item1,item2item1item2,item3读入后：items 1 {item1,item2}2 {item1}3 {item2,item3}single: single的意思，顾名思义，就是单独的交易，简单说，交易记录为：顾客1买了产品1， 顾客1买了产品2，顾客2买了产品3……(产品1，产品2，产品3中可以是单个产品，也可以是多个产品)，如下形式：trans1 item1trans2 item1trans2 item2读入后：items transactionID1 {item1}trans12 {item1, item2}trans2sep:文件中数据是怎么被分隔的，默认为空格，click_detail里面用逗号分隔cols:对basket, col=1,表示第一列是数据的transaction ids(交易号)，如果col=NULL，则表示数据里面没有交易号这一列；对single，col=c(1,2)表示第一列是transaction ids，第二列是item idsrm.duplicates:是否移除重复项，默认为FALSEencoding:写到这里研究了encoding是什么意思，发现前面txt可以不是”ANSI”类型，如果TXT是“UTF-8”，写encoding=”UTF-8”，就OK了.3)rules <- apriori(click_detail,parameter = list(supp=0.01,conf=0.5,target="rules"))：apriori函数apriori(data, parameter = NULL, appearance = NULL, control = NULL)data:数据parameter：设置参数，默认情况下parameter=list(supp=0.1,conf=0.8,maxlen=10,minlen=1,target=”rules”)supp:支持度（support）conf:置信度（confidence）maxlen,minlen:每个项集所含项数的最大最小值target:“rules”或“frequent itemsets”（输出关联规则/频繁项集）apperence:对先决条件X（lhs），关联结果Y（rhs）中具体包含哪些项进行限制，如：设置lhs=beer，将仅输出lhs含有beer这一项的关联规则。默认情况下，所有项都将无限制出现。control：控制函数性能，如可以设定对项集进行升序sort=1或降序sort=-1排序，是否向使用者报告进程（verbose=F/T）补充通过支持度控制：rules.sorted_sup = sort(rules, by=”support”)通过置信度控制：rules.sorted_con = sort(rules, by=”confidence”)通过提升度控制：rules.sorted_lift = sort(rules, by=”lift”)Apriori算法两步法：1. 频繁项集的产生：找出所有满足最小支持度阈值的项集，称为频繁项集；2. 规则的产生：对于每一个频繁项集l，找出其中所有的非空子集；然后，对于每一个这样的子集a，如果support(l)与support(a)的比值大于最小可信度，则存在规则a==>(l-a)。频繁项集产生所需要的计算开销远大于规则产生所需的计算开销频繁项集的产生几个概念：1， 一个包含K个项的数据集，可能产生2^k个候选集 2，先验原理：如果一个项集是频繁的，则它的所有子集也是频繁的（理解了频繁项集的意义，这句话很容易理解的）；相反，如果一个项集是非频繁的，则它所有子集也一定是非频繁的。 3基于支持度（SUPPORT）度量的一个关键性质：一个项集的支持度不会超过它的子集的支持度（很好理解，支持度是共同发生的概率，假设项集{A,B,C}，{A,B}是它的一个自己，A,B,C同时发生的概率肯定不会超过A,B同时发生的概率）。上面这条规则就是Apriori中使用到的，如下图，当寻找频繁项集时，从上往下扫描，当遇到一个项集是非频繁项集（该项集支持度小于Minsup），那么它下面的项集肯定就是非频繁项集，这一部分就剪枝掉了。一个例子（百度到的一个PPT上的）：当我在理解频繁项集的意义时，在R上简单的复现了这个例子，这里采用了eclat算法，跟apriori应该差不多：代码：item <- list(c("bread","milk"),c("bread","diaper","beer","eggs"),c("milk","diaper","beer","coke"),c("bread","milk","diaper","beer"),c("bread","milk","diaper","coke"))names(item) <- paste("tr",c(1:5),sep = "")itemtrans <- as(item,"transactions") #将List转为transactions型rules = eclat(trans,parameter = list(supp = 0.6,target ="frequent itemsets"),control = list(sort=1))inspect(rules) #查看频繁项集运行后结果：>inspect(rules)items support1{beer, diaper}0.62{diaper, milk} 0.63{bread,diaper}0.64{bread,milk} 0.65{beer} 0.66{milk} 0.87{bread} 0.88{diaper} 0.8以上就是该例子的所有频繁项集，然后我发现少了{bread,milk,diaper}这个项集，回到例子一看，这个项集实际上只出现了两次，所以是没有这个项集的。规则的产生每个频繁k项集能产生最多2k-2个关联规则将项集Y划分成两个非空的子集X和Y-X，使得X ->Y-X满足置信度阈值定理：如果规则X->Y-X不满足置信度阈值，则X’->Y-X’的规则一定也不满足置信度阈值，其中X’是X的子集Apriori按下图进行逐层计算，当发现一个不满足置信度的项集后，该项集所有子集的规则都可以剪枝掉了。

---