Lefse

最近在折腾lefse分析，linux下直接用conda安装LEfSe没成功（conda install -c biobakery lefse或conda install -c bioconda lefse），基本安装最后都报错。因为刚开始接触Linux，不懂这些错误，也没有深究。因为lefse是调用R包和python进行的分析，因此参照CSDN吴苏博文在Ubuntu中手动配置R和python环境，然后调用函数进行分析

refs: 
https://blog.csdn.net/qq_38854576/article/details/86763650
https://cloud.tencent.com/developer/article/1667023
https://zhuanlan.zhihu.com/p/265681612            # 后面两个用conda安装的博主

Method1-Win10 Ubuntu

## win10下子系统请参照刘永鑫老师博客文章进行安装，这里不做赘述
## ref: https://mp.weixin.qq.com/s/0PfA0bqdvrEbo62zPVq4kQ

## 建立一个环境
(base) ly@LAPTOP-DHF4UD6Q:~$ conda create --name lefse python=2.7 r=3.6 
(base) ly@LAPTOP-DHF4UD6Q:~$ conda activate lefse  #激活lefse环境 
(lefse) ly@LAPTOP-DHF4UD6Q:~$ conda env list

# conda environments:
#
base                     /home/ly/miniconda3
lefse                 *  /home/ly/miniconda3/envs/lefse

## 配置所需R包
## 笨方法，直接进入R环境，install.packages()安装需要用到的R包
(lefse) ly@LAPTOP-DHF4UD6Q:~$ R
>>> install.packages("MASS")

## 简书—灵木er博主也给出了另一种方法， 建立一个R脚本，在Ubuntu里运行
## ref: https://www.jianshu.com/p/fa554a029c0e
## 后面不再添加 (lefse) ly@LAPTOP-DHF4UD6Q:~$
cat > installr.R << EOF
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
package_list <- c("splines","stats4","survival", "mvtnorm", "modeltools", "coin", "MASS")
for(p in package_list){
 if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
 install.packages(p)
 suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
 }
}
EOF
Rscript installr.R  #运行R脚本

## 配置python包
## rpy2 (v. 2.1 or higher), numpy, matplotlib (v. 1.0 or higher), argparse
## 因为我不需要matplotlib画图，就没装这个
conda install numpy
conda install argparse

## rpy2安装参照上面灵木er博文
## 默认下载在/home/ly/目录下
wget -c https://files.pythonhosted.org/packages/32/54/d102eec14f9cabd0df60682a38bd45c36169a1ec8fb8a690bf436cb6d758/rpy2-2.8.6.tar.gz
tar xzvf rpy2-2.8.6.tar.gz
cd rpy2-2.8.6
python2 setup.py install

## github获得lefse源代码文件
## 我下的是2016年的1.0.8版本
wget -c https://github.com/biobakery/lefse/archive/1.0.8.tar.gz
tar xzvf 1.0.8.tar.gz
cd lefse-1.0.8

## 检验主函数调用
## 参数 -c -u见解释
(lefse) ly@LAPTOP-DHF4UD6Q:~/lefse-1.0.8$ format_input.py -h
usage: format_input.py [-h] [--output_table OUTPUT_TABLE] [-f {c,r}]
                       [-c [1..n_feats]] [-s [1..n_feats]] [-o float]
                       [-u [1..n_feats]] [-m {f,s}] [-n int]
                       [-biom_c BIOM_CLASS] [-biom_s BIOM_SUBCLASS]
                       INPUT_FILE OUTPUT_FILE

LEfSe formatting modules

positional arguments:
  INPUT_FILE            the input file, feature hierarchical level can be
                        specified with | or . and those symbols must not be
                        present for other reasons in the input file.
  OUTPUT_FILE           the output file containing the data for LEfSe

optional arguments:
  -h, --help            show this help message and exit
  --output_table OUTPUT_TABLE
                        the formatted table in txt format
  -f {c,r}              set whether the features are on rows (default) or on
                        columns
  -c [1..n_feats]       set which feature use as class (default 1)
  -s [1..n_feats]       set which feature use as subclass (default -1 meaning
                        no subclass)
  -o float              set the normalization value (default -1.0 meaning no
                        normalization)
  -u [1..n_feats]       set which feature use as subject (default -1 meaning
                        no subject)
  -m {f,s}              set the policy to adopt with missin values: f removes
                        the features with missing values, s removes samples
                        with missing values (default f)
  -n int                set the minimum cardinality of each subclass
                        (subclasses with low cardinalities will be grouped
                        together, if the cardinality is still low, no pairwise
                        comparison will be performed with them)
  -biom_c BIOM_CLASS    For biom input files: Set which feature use as class
  -biom_s BIOM_SUBCLASS
                        For biom input files: set which feature use as
                        subclass

## 检验主函数调用
## 参数-l 2定义LDA的检验阈值
(lefse) ly@LAPTOP-DHF4UD6Q:~/lefse-1.0.8$ run_lefse.py -h
usage: run_lefse.py [-h] [-o str] [-a float] [-w float] [-l float]
                    [--nlogs int] [--verbose int] [--wilc int] [-r str]
                    [--svm_norm int] [-b int] [-e int] [-c int] [-f float]
                    [-s {0,1,2}] [--min_c int] [-t str] [-y {0,1}]
                    INPUT_FILE OUTPUT_FILE

LEfSe 1.0

positional arguments:
  INPUT_FILE      the input file
  OUTPUT_FILE     the output file containing the data for the visualization
                  module

optional arguments:
  -h, --help      show this help message and exit
  -o str          set the file for exporting the result (only concise textual
                  form)
  -a float        set the alpha value for the Anova test (default 0.05)
  -w float        set the alpha value for the Wilcoxon test (default 0.05)
  -l float        set the threshold on the absolute value of the logarithmic
                  LDA score (default 2.0)
  --nlogs int     max log ingluence of LDA coeff
  --verbose int   verbose execution (default 0)
  --wilc int      wheter to perform the Wicoxon step (default 1)
  -r str          select LDA or SVM for effect size (default LDA)
  --svm_norm int  whether to normalize the data in [0,1] for SVM feature
                  waiting (default 1 strongly suggested)
  -b int          set the number of bootstrap iteration for LDA (default 30)
  -e int          set whether perform the wilcoxon test only among the
                  subclasses with the same name (default 0)
  -c int          set whether perform the wilcoxon test ing the Curtis's
                  approach [BETA VERSION] (default 0)
  -f float        set the subsampling fraction value for each bootstrap
                  iteration (default 0.66666)
  -s {0,1,2}      set the multiple testing correction options. 0 no correction
                  (more strict, default), 1 correction for independent
                  comparisons, 2 correction for independent comparison
  --min_c int     minimum number of samples per subclass for performing
                  wilcoxon test (default 10)
  -t str          set the title of the analysis (default input file without
                  extension)
  -y {0,1}        (for multiclass tasks) set whether the test is performed in
                  a one-against-one ( 1 - more strict!) or in a one-against-
                  all setting ( 0 - less strict) (default 0)

## 分析运行
## 将win10本地盘文件copy到 /home/ly/lefse-1.0.8分析
cd  /mnt/e/Python/lefse/
cp pathway_sel_lefse_Wagner.txt /home/ly/lefse-1.0.8
cd ~
cd lefse-1.0.8
format_input.py pathway_sel_lefse_Wagner.txt A_lefse_test.in -c 1 -u 2 -o 1000000
run_lefse.py A_lefse_test.in lefse.out.txt -l 3

未标题-1.jpg

Method2-Win10 cmd

具体安装需求参照吴苏博文
ref： https://blog.csdn.net/qq_38854576/article/details/86763650

E:\Python\lefse>format_input.py pathway_sel_lefse_Wagner.txt A_lefse_test.in -c 1 -u 2 -o 1000000
usage: format_input.py [-h] [--output_table OUTPUT_TABLE] [-f {c,r}]
                       [-c [1..n_feats]] [-s [1..n_feats]] [-o float]
                       [-u [1..n_feats]] [-m {f,s}] [-n int]
                       [-biom_c BIOM_CLASS] [-biom_s BIOM_SUBCLASS]
                       INPUT_FILE OUTPUT_FILE
format_input.py: error: too few arguments

## 依然出现这个问题，暂时无解

Method3-R实现

此方法主要参照公众号微生信生物 文涛老师的代码
博文题目：把握LEfSe核心算法使用R实现整套算法
ref：https://mp.weixin.qq.com/s?__biz=MzUzMjYyMDE2OQ==&mid=2247487316&idx=2&sn=0d21754b4c5cdee690a05ff865627037&chksm=fab139b2cdc6b0a4055491e803d67d4dce9b69d2c21098648124a72fcd56dfc6013a98a495e1&mpshare=1&scene=23&srcid=0204YyUGkk7q2tnh6HraCf4Z&sharer_sharetime=1612782715958&sharer_shareid=21a8baa175f8629f0dc4f5b4b45a0f28#rd

## 详细代码不在此贴出，自己搜索
## 此处有个小疑问，此处公式计算的LDA score, 即从lda函数计算的分组均值中用较大的-较小的，然后去对数
## 而Tax4Fun计算得出的功能基因、通路丰度均是相对丰度，值在[0,1]之间，那么导入此处公式计算的LDA score均是小于2
## LDA score = 2, 则 abs(ldamean$max-ldamean$min)/2项的值至少为99，而使用相对丰度肯定不成
ldamean$LDAscore <- signif(log10(1+abs(ldamean$max-ldamean$min)/2),digits=3)

## 而采用PICRUSt得出的功能基因、通路丰度多是绝对丰度，那么对于某个基因或通路，组间均值的差值很容易超过99吧？？
## 毕竟这些丰度值都比较大。那么lefse算法如何处理这个问题。

## format_input.py函数的参数-o float定义了标准化值（set the normalization value）
## 不知此处如何处理，python的源代码不太熟悉，待以后在深挖吧。。。。