R语言实现扩增子数据抽平算法

vlambda
2020-02-17

R语言实现扩增子数据抽平算法

写在前面

现在是20年2月中旬了，今年注定是不平凡的一年，尽然现在都还在家中。但是任务量还是十分繁重。谁让做生物信息的也不用限制什么地点。

为什么我要实现抽平算法呢？16年phyloseq的出现很对我的胃口，因为扩增子数据本身较小，自R语言中的操作也是基本没障碍，本以为phyloseq中实现了抽平的方法，我没有查看源码，直接采用了，但是最近再做一些项目中发现这个抽平方法是有问题的：

library(tidyverse)
ps = readRDS("../ps_liu.rds")

total = min(sample_sums(ps));total
standf = function(x,t = total)(t*(x/sum(x)))

ps11 = transform_sample_counts(ps,standf);ps11

sample_sums(ps11)

简单看一下似乎结果没设么问题？但是当我想想这种抽平写法，却实在想不通这是怎么抽平的。

KO1   KO2   KO3   KO4   KO5   KO6   OE1   OE2   OE3   OE4   OE5   OE6   WT1   WT2   WT3   WT4   WT5   WT6
32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049 32049

于是我打开了源代码；从这里我们就发现，这并不是抽平，只是标准化而已。这也很好的解释了standf = function(x,t = total)(t*(x/sum(x)))这行代码的意义。这只是一中标准化方式，意味着如果序列数量较少，按照所谓phyloseq抽平的方式得到的结果不是整数。

function (physeq, fun, ...)
{
    if (!identical(length(fun(1:10)), 10L)) {
        stop("`fun` not valid function.")
    }
    if (taxa_are_rows(physeq)) {
        newphyseq = apply(as(otu_table(physeq), "matrix"),
            2, fun, ...)
        if (identical(ntaxa(physeq), 1L)) {
            newphyseq <- matrix(newphyseq, 1L, nsamples(physeq),
                TRUE, list(taxa_names(physeq), sample_names(physeq)))
        }
    }
    else {
        newphyseq = apply(t(as(otu_table(physeq), "matrix")),
            2, fun, ...)
        if (identical(ntaxa(physeq), 1L)) {
            newphyseq <- matrix(newphyseq, 1L, nsamples(physeq),
                TRUE, list(taxa_names(physeq), sample_names(physeq)))
        }
        newphyseq = t(newphyseq)
    }
    if (!identical(dim(newphyseq), dim(otu_table(physeq)))) {
        stop("Dimensions of OTU table change after apply-ing function. \n",
            " Please check both function and table")
    }
    otu_table(physeq) <- otu_table(newphyseq, taxa_are_rows = taxa_are_rows(physeq))
    return(physeq)
}

这里我们来验证一下：使用100条序列抽平，

total = 100
standf = function(x,t = total)(t*(x/sum(x)))

ps11 = transform_sample_counts(ps,standf);ps11

sample_sums(ps11)

#KO1 KO2 KO3 KO4 KO5 KO6 OE1 OE2 OE3 OE4 OE5 OE6 WT1 WT2 WT3 WT4 WT5 WT6
#100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

sample_sums函数统计出来序列数量是100，但是我们具体查看OTU表格，就会发现巨大的问题，what？这明明是小数啊？



vegan_otu <- function(physeq){
 OTU <- otu_table(physeq)
 if(taxa_are_rows(OTU)){
 OTU <- t(OTU)
 }
 return(as(OTU,"matrix"))
}
otu_table = as.data.frame(t(vegan_otu(ps11)))

head(otu_table)

这样的结果必然会是的一些alpha多样性的算法无法正常实现。

> head(otu_table)
                 KO1       KO2         KO3         KO4         KO5        KO6         OE1         OE2         OE3         OE4         OE5
ASV_1    3.349181512 5.4462433 2.094939796 3.600482864 2.726572529 2.94213176 3.915281931 4.932505725 3.950103950 2.920284136 3.743058682
ASV_1591 0.036109774 0.0276740 0.002567328 0.036739621 0.007702182 0.03247984 0.046243487 0.048210196 0.011550012 0.048570214 0.093051178
ASV_657  0.015045739 0.0193718 0.015403969 0.023618328 0.020539153 0.01623992 0.049326387 0.078341569 0.025987526 0.018213830 0.027014858
ASV_28   0.761314396 0.8385222 0.138635722 1.152049546 0.338896021 0.49261084 0.474766470 0.792455104 0.441787942 0.780159067 0.726399520
ASV_1717 0.009027443 0.0166044 0.000000000 0.002624259 0.000000000 0.01082661 0.006165798 0.018078824 0.005775006 0.000000000 0.003001651
ASV_2407 0.003009148 0.0055348 0.017971297 0.000000000 0.002567394 0.01894657 0.006165798 0.006026275 0.002887503 0.006071277 0.009004953
                 OE6         WT1         WT2         WT3         WT4         WT5         WT6
ASV_1    3.173265937 6.270221129 6.611613307 4.708391436 5.428246384 3.838971294 4.188622433
ASV_1591 0.056163999 0.045455761 0.151092818 0.046482378 0.092095997 0.264112688 0.029573072
ASV_657  0.034322444 0.005347737 0.002605049 0.008202773 0.021669646 0.048020489 0.037638456
ASV_28   0.555399544 1.045482500 1.172271863 0.642550515 0.869494556 1.197844414 0.650607592
ASV_1717 0.006240444 0.010695473 0.000000000 0.002734258 0.002708706 0.000000000 0.021507689
ASV_2407 0.024961777 0.010695473 0.002605049 0.005468515 0.000000000 0.002667805 0.008065383

所以我需要编写一份抽平代码：同样我们设置N为100，代表按照100条序列进行抽平

# 提取OTU表格
otb = as.data.frame(t(vegan_otu(ps)))
dim(otb)
#设置抽平数量
N = 100
#-------------开始抽平------------
i = 1
ii =1
a = c()
aa = c()
for (i in 1:dim(otb)[1]) {
  
  b= paste(row.names(otb[i,]),1:otb[i,1],sep = "&")
  a = c(a,b)
}

ss = sample(a,N,replace = FALSE)
str_split(ss, pattern = "&", n = Inf, simplify = FALSE)

sample.names <- sapply(str_split(ss, pattern = "&", n = Inf, simplify = FALSE), `[`, 1)


aa = as.data.frame(table(sample.names))
sum(aa$Freq)
colnames(aa)= c("ID",colnames(otb)[ii])
head(aa)
row.names(aa) = aa$ID
aa$ID = NULL


for (ii in 2:dim(otb)[2]) {
  
  for (i in 1:dim(otb)[1]) {
    b= paste(row.names(otb[i,]),1:otb[i,ii],sep = "&")
    a = c(a,b)
  }
  
  ss = sample(a,N,replace = FALSE)
  # str_split(ss, pattern = "&", n = Inf, simplify = FALSE)
  
  sample.names <- sapply(str_split(ss, pattern = "&", n = Inf, simplify = FALSE), `[`, 1)
  bb = as.data.frame(table(sample.names))
  # bb[2]
  colnames(bb)= c("ID",colnames(otb)[ii])
  head(bb)
  row.names(bb) = bb$ID
  bb$ID = NULL
  head(bb)
  aa = merge(aa,bb,by = "row.names",all = TRUE)
  row.names(aa) = aa$Row.names
  aa$Row.names = NULL

  }


dim(aa)
aa[is.na(aa)] = 0

head(aa)

这次查看，是否就是正常额抽平结果呢。欢迎大家使用。

 head(aa)
         KO1 KO2 KO3 KO4 KO5 KO6 OE1 OE2 OE3 OE4 OE5 OE6 WT1 WT2 WT3 WT4 WT5 WT6
ASV_1      7   5   6   2   3   2   2   7   3   3   1   3   2   3   3   1   5   3
ASV_10     0   0   5   2   0   3   2   1   1   4   2   1   2   2   0   1   1   0
ASV_100    0   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
ASV_101    0   1   0   0   0   1   1   0   1   0   0   0   1   0   0   0   0   0
ASV_1017   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0   0
ASV_102    0   0   1   0   1   0   0   0   1   0   0   0   0   0   0   0   0   0

写在后面

这个方法实现了抽平，但是运行效率很低，抽平一次也还可以接受，如果我需要做稀释曲线，需要迭代大量的抽平次数，运行时间将会变得很长，在这里我想博采众长，希望各位同行可以给我一些提高速度的建议。

vlambda博客
学习文章列表