Lasso与宏基因组

发表于 | 分类于 |
字数统计: 754 字 | 阅读时长 ≈ 3 分钟

因为项目需要，去了解了一下LASSO【The Least Absolute Shrinkage and Selection Operator】，又译最小绝对值收敛和选择算子、套索算法。鼓捣了几天下来，我也没有弄太懂，但是大概知道了怎么使用，特此记录一下。

  1. 首先是知乎的介绍：
    https://zhuanlan.zhihu.com/p/42122611

  2. 另外查的过程中发现，LASSO并不产生P值，并且不建议去计算P值，虽然原理没弄大懂，但是可以记住结论。同样，附上参考链接：
    Are p values generally reported in LASSO regressions? : statistics

  3. 如果要计算P值的话，得通过另外一种回归计算：
    https://stats.stackexchange.com/questions/84185/lasso-to-identify-important-variables-in-ordered-logistic-regression
    而我看到的一篇文章里是用有序回归（ordinal regression）计算的P值。
    Chen, D. Q. et al. Identification of serum metabolites associating with chronic kidney disease progression and anti-fibrotic effect of 5-methoxytryptophan. Nat. Commun. 10, 1–15 (2019).

  4. 怎么是使用LASSO分析微生物的数据：
    参考了这个脚本，该脚本只适用于二元分组的数据，并且接入的是 Biom 格式的文件：https://github.com/alifar76/MicrobeNets

同样因为项目需要，我基于该脚本按照我的理解修改成接入普通文件以及多组的情况，组别文件需要将字符型变量转化成数值型变量，脚本如下，不一定准确，仅供参考：
Biom 格式的也是类似的，只不过需要处理文件的最后一列的物种注释信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
## 用于非biom 格式文件, Genus
intable <- "Genus_abundance.xls" # 丰度文件
mapfile <- "Sample_Group.txt" # 分组文件，需要转化成数值型
outfile <- “Genus_lasso_output.xls" 
measuretype <- "deviance” #这些设置的说明可以参考上边github的脚本
famaliy <- "multinomial" #同上

otutab <- read.table(intable,header = T, sep = "\t", check.names = T, row.names =1)
group <- read.table(mapfile,header = T, sep = "\t", check.names = T, row.names =1, comment.char= "")
otutab <- subset(as.matrix(otutab),select=rownames(group))
otutab <- t(otutab)
 
fitted = glmnet(as.matrix(otutab), group[,1],standardize=FALSE,alpha=1)# 1 mean lasso
#plot(fitted,xvar = "lambda", label = TRUE)

cv.dat = cv.glmnet(data.matrix(otutab),group[,1],grouped=FALSE,nfolds=10,alpha = 1,parallel=TRUE,type.measure="deviance",family="multinomial") #输出的时候会有 warning, 据我测试是和nfold有关；
#plot(cv.dat,main=paste("cv_",compid,"_",methodtype,sep=""))
coefval <- coef(cv.dat)

### 输出
write.table(as.matrix(t(c("Tax","coefficient"))),file=outfile, sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE,append=TRUE)   

for (i in 1:length(row.names(coefval[[1]]))){ 
    rowName = row.names(coefval[[1]])[i]
    all_dat <- c()
    if (rowName != "(Intercept)"){
        if (coefval[[1]][i] != 0){
            all_dat <- c(rowName, coefval[[1]][i])
            write.table(as.matrix(t(all_dat)),file=outfile, sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE,append=TRUE)   
        }
        
    }
}

  1. 效果如何?
    就我的结果来看，不是很理想。我另外使用MaAsLin2做了测试，发现后者的结果更加好一些。两种方法的重叠物种也不多。个人经验和意见，仅供参考。
    MaAsLin2

  2. Python如何实现LASSO？
    使用了该方法的文章以及脚本：
    Wilmanski, T. /et al./ Blood metabolome predicts gut microbiome α-diversity in humans. /Nat. Biotechnol./ doi:10.1038/s41587-019-0233-9
    GitHub - PriceLab/ShannonMets: ** Blood metabolome predicts gut microbiome α-diversity in humans**

(✪ω✪)