2

我的数据文件结构如下：

OTU1    PIA0    1120
OTU2    PIA1    2
OTU2    PIA3    6
OTU2    PIA4    10
OTU2    PIA5    1078
OTU2    PIN1    24
OTU2    PIN2    45
OTU2    PIN3    261
OTU2    PIN4    102
OTU3    PIA0    16
OTU3    PIA1    59
OTU3    PIA2    27
OTU3    PIA3    180
OTU3    PIA4    200
OTU3    PIA5    251
OTU3    PIN0    36
OTU3    PIN1    61
OTU3    PIN2    156
OTU3    PIN3    590
OTU3    PIN4    277
OTU4    PIA0    401
OTU4    PIN0    2

我想创建一个矩阵，显示来自第二列的数据组合，以第一列作为组合计数的参考（显示多少次，通过第一列号测量每个 -OTU1、OTU2、OTU3、OTU4-第二列中的每个数据在同一个 OTU 中彼此一起出现）。它需要看起来像这样：

    PIA0  PIA1  PIA2  PIA3  PIA4  PIA5  PIN0  PIN1  PIN2  PIN3  PIN4
PIA0  1     1     1     1     1     1     2     1     1     1     1 
PIA1  1     0     1     2     2     2     1     2     2     2     2
PIA2  1     1     0     1     1     1     1     1     1     1     1
PIA3  1     2     1     0     2     2     1     2     2     2     2
PIA4  1     2     1     2     0     2     1     2     2     2     2
PIA5  1     2     1     2     2     0     1     2     2     2     2
PIN0  2     1     1     1     1     1     0     1     1     1     1
PIN1  1     2     1     2     2     2     1     0     2     2     2
PIN2  1     2     1     2     2     2     1     2     0     2     2
PIN3  1     2     1     2     2     2     1     2     2     0     2
PIN4  1     2     1     2     2     2     1     2     2     2     0

在同名的行和列之间共享的数据反映了该数据在 OTU 中单独出现的次数。

有任何想法吗？

我在这里阅读了有关 R 库“reshape2”和命令“acast”的信息，但是我只能更改包含所有数据的矩阵的形状，而不是根据需要进行组合计数。我也一直在想一个Biopython脚本，但是我觉得用我对编程的一点了解，写下来会太大而困难。

目标是建立一个像示例中那样的矩阵，以便我可以使用这些数据运行CIRCOS 在线程序。

4

1 回答 1

4

您可以使用dcast创建一个二进制矩阵来指示每个 OTU 中每个 PI 的存在，然后将其与自身相乘以获得计数。

d <- read.fwf( textConnection("
OTU1    PIA0    1120
OTU2    PIA1    2
OTU2    PIA3    6
OTU2    PIA4    10
OTU2    PIA5    1078
OTU2    PIN1    24
OTU2    PIN2    45
OTU2    PIN3    261
OTU2    PIN4    102
OTU3    PIA0    16
OTU3    PIA1    59
OTU3    PIA2    27
OTU3    PIA3    180
OTU3    PIA4    200
OTU3    PIA5    251
OTU3    PIN0    36
OTU3    PIN1    61
OTU3    PIN2    156
OTU3    PIN3    590
OTU3    PIN4    277
OTU4    PIA0    401
OTU4    PIN0    2"), widths=c(8,8,10), header=FALSE, skip=1 )

library(reshape2)
A <- as.matrix( dcast( V1 ~ V2, data=d, length )[,-1]>0 )
#          PIA0     PIA1     PIA2     PIA3     PIA4     PIA5     PIN0     PIN1     PIN2     PIN3     PIN4    
# [1,]     TRUE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE
# [2,]    FALSE     TRUE    FALSE     TRUE     TRUE     TRUE    FALSE     TRUE     TRUE     TRUE     TRUE
# [3,]     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE
# [4,]     TRUE    FALSE    FALSE    FALSE    FALSE    FALSE     TRUE    FALSE    FALSE    FALSE    FALSE
t(A) %*% A
#              PIA0     PIA1     PIA2     PIA3     PIA4     PIA5     PIN0     PIN1     PIN2     PIN3     PIN4    
# PIA0            3        1        1        1        1        1        2        1        1        1        1
# PIA1            1        2        1        2        2        2        1        2        2        2        2
# PIA2            1        1        1        1        1        1        1        1        1        1        1
# PIA3            1        2        1        2        2        2        1        2        2        2        2
# PIA4            1        2        1        2        2        2        1        2        2        2        2
# PIA5            1        2        1        2        2        2        1        2        2        2        2
# PIN0            2        1        1        1        1        1        2        1        1        1        1
# PIN1            1        2        1        2        2        2        1        2        2        2        2
# PIN2            1        2        1        2        2        2        1        2        2        2        2
# PIN3            1        2        1        2        2        2        1        2        2        2        2
# PIN4            1        2        1        2        2        2        1        2        2        2        2
于 2013-08-21T11:31:26.857 回答