R:拆分一个“:”分隔的 VCF 文件,使用 'for-loop'(迭代多个列)创建多个矩阵

R: Splitting up a ":" delimited VCF file, using a 'for-loop' (iterating over several columns) to create multiple matrices

我为什么要问这个?

似乎很多人对拆分 VCF 文件和使用 for 循环遍历列都有问题,但我还没有遇到任何以与使用相关的方式解决这两个问题的问题包含许多示例的 VCF 文件 - 将进行解释。

这里是数据结构的例子:

Loci    Sample1
[1]     0/1:15:55:54:49:5:9.26%:2.8371E-2:37:36:49:0:5:0
[2]     0/1:42:55:53:40:13:24.53%:5.2873E-5:34:37:40:0:13:0
[3]     0/1:15:54:54:49:5:9.26%:2.8371E-2:35:33:49:0:5:0

问题是如何在许多基因座(行)和多个样本(列)上创建一个眼睛友好的 table 并具有大量输出统计信息(每个由“:”分隔)?

这个问题我已经解决了一半:

我开发了一个 R 脚本,它可以从单个样本列中获取信息并输出一个矩阵来分隔每个单独的统计数据。代码如下:

data <- vcf.small

# First, create a list representing each row (locus) and separate the
# statistics; second, breakdown the list's structure but maintain data order.
split1 <-strsplit(as.character(data$Sample1),":")
split2 <- unlist(split1)

# Create a matrix: here, there are 14 values by 3 loci.
mtx1a <- matrix(split2, ncol=14, nrow=3, dimnames=list(NULL,c("GT","GQ","SDP","DP","RD","AD","FREQ","PVAL","RBQ","ABQ","RDF","RDR","ADF","ADR")), byrow=TRUE)

# Create some additional variables (columns) to add to the matrix.
sample <- matrix(rep(1,3), ncol=1, nrow=3, dimnames=list(NULL,c("SAMPLE")))
locus <- matrix(1:3, ncol=1, nrow=3, dimnames=list(NULL,c("LOCUS")))

# Add them to the matrix.
mtx1b <- cbind(mtx1a,sample)
mtx1b <- cbind(mtx1b,locus)

瞧,输出:

     GT    GQ   SDP  DP   RD   AD   FREQ     PVAL        RBQ  ABQ  RDF  RDR ADF  ADR SAMPLE LOCUS
[1,] "0/1" "15" "55" "54" "49" "5"  "9.26%"  "2.8371E-2" "37" "36" "49" "0" "5"  "0" "1"    "1"  
[2,] "0/1" "42" "55" "53" "40" "13" "24.53%" "5.2873E-5" "34" "37" "40" "0" "13" "0" "1"    "2"  
[3,] "0/1" "15" "54" "54" "49" "5"  "9.26%"  "2.8371E-2" "35" "33" "49" "0" "5"  "0" "1"    "3" 

'for-loop'问题:

输出是完美的,但现在我无法弄清楚如何制作一个包含上述代码的 for 循环来为每个样本创建一个单独的矩阵。我推理:

for(i in names(data){
    split[i] <-strsplit(as.character(data$[i]),":")
    split[i] <- unlist(split[i])
    mtx[i]a <- matrix(split2, ncol=14, nrow=3,  
[etc etc..]
}       

问题是我需要创建自定义的单独变量来为每个样本(即列)设置每个矩阵。但是,R 不会将 [i] 作为占位符,其中 i = 样本(/列)名称。

理想情况下,每个样本(/列)特定变量将读作:"splitSample1"、"splitSample2"、"splitSample3" 等。这主要是为了允许 for 循环处理所有列,而不必为每个列名称重新创建特定代码。我想我想做的是从 Linux 重新创建“$i”语法,但显然这在这里不起作用。

解决这个问题将使处理非常大的数据集变得更易于管理,我真的尝试过寻找解决方法。非常感谢任何帮助!

我认为最好将结果存储在 data.framedata.table 中,因为每个拆分列的 class 类型都不同。 matrix 只能存储一个 class。如果只有一个 character 列,则所有 columns.

的 class 将为 character

使用 data.tabledevel 版本,我们可以使用 tstrsplit 拆分成列,并将 class 更改为 type.convert=TRUE。开发版可以从here

安装
library(data.table)#v1.9.5+
nm1 <- c('GT', 'GQ', 'SDP', 'DP', 'RD', 'AD', 'FREQ', 'PVAL', 'RBQ',
   'ABQ', 'RDF', 'RDR', 'ADF', 'ADR')

setDT(data)[, (nm1):=tstrsplit(Sample1, ':', type.convert=TRUE)][,
         Sample1:=NULL][, c('sample', 'locus'):= list(1, 1:3)][]
#    GT GQ SDP DP RD AD   FREQ       PVAL RBQ ABQ RDF RDR ADF ADR sample locus
#1: 0/1 15  55 54 49  5  9.26% 2.8371e-02  37  36  49   0   5   0      1     1
#2: 0/1 42  55 53 40 13 24.53% 5.2873e-05  34  37  40   0  13   0      1     2
#3: 0/1 15  54 54 49  5  9.26% 2.8371e-02  35  33  49   0   5   0      1     3

如果数据集中有多个 'Sample' 列,我们可以使用 lapply 循环列并在列表中创建拆分数据集 ('lst')。

nm2 <- paste0('splitSample', 1:ncol(data2))
lst <- setNames(
       lapply(seq_len(ncol(data2)), function(i)
          setDT(list(data2[,i]))[, (nm1) := tstrsplit(V1, ":", 
             type.convert=TRUE)][, V1:=NULL][,
               c('sample', 'locus'):= list(i, 1:.N)]), 
                 nm2)

在'list'中工作会更容易,但如果我们需要在全局环境中有单独的数据集对象(不推荐),我们可以使用list2env.

list2env(lst, envir=.GlobalEnv)
splitSample1
#    GT GQ SDP DP RD AD   FREQ      PVAL RBQ ABQ RDF RDR ADF ADR sample locus
#1: 0/1 15  55 54 49  5  9.26% 2.8371E-2  37  36  49   0   5   0      1     1
#2: 0/1 42  55 53 40 13 24.53% 5.2873E-5  34  37  40   0  13   0      1     2
#3: 0/1 15  54 54 49  5  9.26% 2.8371E-2  35  33  49   0   5   0      1     3

splitSample2
#    GT GQ SDP DP RD AD   FREQ      PVAL RBQ ABQ RDF RDR ADF ADR sample locus
#1: 0/2 15  55 55 49  5 10.26%  2.971E-2  37  32  49   0   5   0      2     1
#2: 0/2 52  55 53 40 13 22.53% 1.2873E-5  34  37  12   0  13   0      2     2
#3: 0/2 17  54 54 49 18  9.29% 3.8371E-2  42  33  49   0   5   0      2     3

注意:在这里,我将输入数据集用作 data.frame.

数据

data <- structure(list(Sample1 =
   c("0/1:15:55:54:49:5:9.26%:2.8371E-2:37:36:49:0:5:0", 
 "0/1:42:55:53:40:13:24.53%:5.2873E-5:34:37:40:0:13:0",
  "0/1:15:54:54:49:5:9.26%:2.8371E-2:35:33:49:0:5:0"
 )), .Names = "Sample1", class = "data.frame", row.names = c(NA, -3L))


 data2 <- structure(list(Sample1 =
   c("0/1:15:55:54:49:5:9.26%:2.8371E-2:37:36:49:0:5:0", 
  "0/1:42:55:53:40:13:24.53%:5.2873E-5:34:37:40:0:13:0",
  "0/1:15:54:54:49:5:9.26%:2.8371E-2:35:33:49:0:5:0"
 ), Sample2 = c("0/2:15:55:55:49:5:10.26%:2.971E-2:37:32:49:0:5:0", 
 "0/2:52:55:53:40:13:22.53%:1.2873E-5:34:37:12:0:13:0",
 "0/2:17:54:54:49:18:9.29%:3.8371E-2:42:33:49:0:5:0")),
.Names = c("Sample1", "Sample2"), class = "data.frame",
row.names = c(NA, -3L))