宏基因组数据分析：差异分析（LEfSe安装使用及LDA score计算）

导读

• LEfSe（Linear discriminant analysis Effect Size）通过将用于统计显着性的标准检验与编码生物一致性和效果相关性的检验相结合，确定最有可能解释类别之间差异的特征。
• LEfSe分析可以实现多个分组之间的比较，还进行分组比较的内部进行亚组比较分析，从而找到组间在丰度上有显著差异的物种（即biomaker）。

原理

1. 首先在多组样本中采用的非参数检验Kruskal-Wallis秩和检验检测不同分组间丰度差异显著的特征；
2. 然后在上一步中获得的显著差异特征，用成组的Wilcoxon秩和检验进行组间差异分析（若没有亚组，该步跳过）；
3. 最后用线性判别分析（LDA）对数据进行分类和评估差异显著的物种的影响力（即LDA score）。

LDA score计算

计算步骤包括（详细计算过程看源码）：

1. 拟合lda模型，获取第一特征向量；
2. 对第一特征向量进行标准化；
3. 根据标准化后的第一特征向量，计算样本新坐标；
4. 根据分组信息，计算组间距离，作为效应系数；
5. LDA score = 效应系数 * 标准化后特征向量 ；
6. 对LDA score进行log转换。

def test_lda_r(cls,feats,cl_sl,boots,fract_sample,lda_th,tol_min,nlogs):
    fk = list(feats.keys())
    means = dict([(k,[]) for k in feats.keys()])
    feats['class'] = list(cls['class'])
    clss = list(set(feats['class']))

    for uu,k in enumerate(fk):
        if k == 'class':
            continue

        ff = [(feats['class'][i],v) for i,v in enumerate(feats[k])]

        for c in clss:
            if len(set([float(v[1]) for v in ff if v[0] == c])) > max(float(feats['class'].count(c))*0.5,4):
                continue

            for i,v in enumerate(feats[k]):
                if feats['class'][i] == c:
                    feats[k][i] = math.fabs(feats[k][i] + lrand.normalvariate(0.0,max(feats[k][i]*0.05,0.01)))

    rdict = {
   }

    for a,b in feats.items():
        if a == 'class' or a == 'subclass' or a == 'subject':
            rdict[a] = robjects.StrVector(b)
        else:
            rdict[a] = robjects.FloatVector(b)

    robjects.globalenv["d"] = robjects.DataFrame(rdict)
    lfk = len(feats[fk[0]])
    rfk = int(float(len(feats[fk[0]]))*fract_sample)
    f = "class ~ "+fk[0]

    for k in fk[1:]:
        f += " + " + k.strip()

    ncl = len(set(cls['class']))
    min_cl = int(float(min([cls['class'].count(c) for c in set