def deseq2ParseResults(fdr, vsd=False):
    '''
    Standardises the output format from deseq2.

    Deseq2 has the following output columns:
    baseMean log2FoldChange lfcSE stat pvalue padj
    described in
    https://bioconductor.org/packages/release/bioc/
                 vignettes/DESeq2/inst/doc/DESeq2.pdf

    Standardised columns are generated from this output as follows:

    test_id - the gene or region tested, row names from raw output

    treatment_name - the first group in this differential expression
        comparison (from the design file)

    treatment_mean - the mean expression value for this treatment from the
    normalised count table generated by deseq2

    treatment_std - the standard deviation of experssion for this treatment
    from the normalised count table generated by deseq2

    control_name - the second group in this differential expression
        comparison (from the design file)

    control_mean - the mean expression value for this treatment from the
    normalised count table generated by deseq2

    control_std - the standard deviation of experssion for this treatment
    from the normalised count table generated by deseq2

    pvalue - the pvalue generated by deseq2 (from the pvalue column)

    qvalue - the adjusted pvalue generated by deseq2 (from the padj column)

    l2fold - the log2fold change between normalised counts generated by
    deseq2 (log2FoldChange column)

    fold - control mean / treatment mean

    transformed_l2fold - not applicable, set to 0 (see deseq2 manual,
    "The variance stabilizing and rlog transformations are provided
    for applications other than differential testing,
    for example clustering of samples or other machine learning applications.
    For differential testing we recommend the DESeqfunction applied to raw
    counts"

    signif = True if the qvalue is less than the FDR set in `term`PARAMS.

    status = OK if a pvalue has been generated, else FAIL
    '''
    R(''' fdr=%s ''' % fdr)
    #  assign standard column names
    R('''cols = c("test_id",
                  "treatment_name",
                  "treatment_mean",
                  "treatment_std",
                  "control_name",
                  "control_mean",
                  "control_std",
                  "pvalue",
                  "qvalue",
                  "l2fold",
                  "fold",
                  "transformed_l2fold",
                  "signif",
                  "status") ''')

    #  extract normalised counts
    R('''normalcounts = counts(dds, normalized=T)''')
    #  build empty dataframe
    R('''res2 = data.frame(matrix(nrow=nrow(res), ncol=length(cols)))''')
    R('''colnames(res2) = cols''')

    #  fill columns with values described above
    R('''res2['test_id'] = rownames(res)''')
    R('''g = unique(groups)''')
    R('''g1 = which(groups == g[1])''')
    R('''g2 = which(groups == g[2])''')
    R('''res2['treatment_name'] = rep(unique(groups)[1], nrow(res2))''')
    R('''res2['treatment_mean'] = rowMeans(normalcounts[, g1])''')
    R('''res2['treatment_std'] = apply(normalcounts[, g1], 1, sd)''')
    R('''res2['control_name'] = rep(unique(groups)[2], nrow(res2))''')
    R('''res2['control_mean'] = rowMeans(normalcounts[,g2])''')
    R('''res2['control_std'] = apply(normalcounts[, g2], 1, sd)''')
    R('''res2['pvalue'] = res$pvalue''')
    R('''res2['qvalue'] = res$padj''')
    R('''res2['l2fold'] = res$log2FoldChange''')
    R('''res2['fold'] = res2$control_mean / res2$treatment_mean''')
    R('''res2['signif'] = as.integer(res2$qvalue <= fdr)''')
    R('''res2['status'] = ifelse(is.na(res2$pvalue), "FAIL", "OK")''')

    #  replace l2fold change and fold for expression levels of 0 in treatment
    #  and control with 0
    R('''z1  = which(res2$treatment_mean == 0)''')
    R('''z2 = which(res2$control_mean == 0)''')
    R('''zeros = intersect(z1, z2)''')
    R('''res2$l2fold[zeros] = 0''')
    R('''res2$fold[zeros] = 0''')

    #  occupy transformed l2fold with 0s
    R('''res2$transformed_l2fold = 0''')

    D = R('res2')
    D.index = D['test_id']
    D = D.drop('test_id', 1)
    return D
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/deseq2ParseResults/original