用于处理基因型数据的 R 代码答案

【问题标题】：R code to work on genotype data用于处理基因型数据的 R 代码
【发布时间】：2015-07-21 03:11:29
【问题描述】：

我有这个名为mydf 的数据。

我需要将REF 和ALT 列中的字母（DNA 字母）与colnames(x) ("A","T","G","C") 匹配，并将对应的数值粘贴在一起作为"REF,ALT"。

但是，在 TYPE 列中有一些行我有 "snp:+[0-9]" 和 "flat$"。

现在对于"flat$" 我想要的行：

sum 来自对应 "start" id 的多达 "snp:+[0-9]" 的 ALT 值，包括扁平线本身，如果 ALT 字母是唯一的（请参阅用大括号括起来的脚本一个扁平线的括号）
粘贴 ALT 值再次作为"REF,ALT" （REF 值将是 "snp:+[0-9]" 和 "flat$" 的起始 ID 相同）
得到输出，如结果所示。

我已经为一条扁平线做了这个，但我需要帮助为flatcase 制作函数，以便它对所有扁平线都做同样的事情。

如何为flatcase 创建一个函数来执行此操作？

代码

normalCase <- function(x, ns) {
      ref.idx <- which(ns == "REF")
      ref.allele <- x[ref.idx]
      ref.count <- x[which(ns == ref.allele)]

      alt.idx <- which(ns == "ALT")
      alt.allele <- x[alt.idx]
      alt.count <- x[which(ns == alt.allele)]

      paste(ref.count, alt.count, sep=",")
    }



    flatcase??{

     g<-x[,"start"]=="chr16:2530921"& grepl("snp:+[0-9]",x[,"TYPE"])
     myt<-x[g,]
     x[g,"ALT"]
     unique(x[g,"ALT"])
     c<-unique(x[g,"ALT"])
     flat<-myt[grepl("flat$",myt[,"TYPE"]),]
     c<-unique(x[g,"ALT"])
    alt.count<- sum(as.numeric(flat[c]))
    }

    calculateAD <- function(x, mat, ns) {
      if (grepl("flat$", x[which(ns == 'TYPE')])) {
        flatCase(x, mat, ns)
      } else {
        normalCase(x, ns)
      }
    }


    bamAD <- function(x) {
       new.x <- cbind(x, apply(x, 1, calculateAD, x, colnames(x)))
      colnames(new.x)[ncol(new.x)] <- "bam.AD"
      new.x      
    }

我为 flatCase 尝试过的功能是：

flatCase <- function(x, mat, ns) {
  id.idx <- which(ns == 'start')
  type.idx <- which(ns == 'TYPE')
  ref.idx <- which(ns == 'REF')
  alt.idx <- which(ns == 'ALT')


  id <- x[id.idx]
  #m <- mat[mat[, id.idx] == id & mat[, type.idx] == "snp", ]
  #m <- mat[mat[, id.idx] == id & mat[, type.idx] == "snp", ]
  m<-mat[grepl(id,mat[, id.idx]) & grepl("snp:+[0-9]",mat[, type.idx]),]
  #flat<-mat[grepl("flat$",mat[, type.idx]),]
  ref.allele <- x[ref.idx]
  ref.count<-x[which(ns == ref.allele)]


  alt.count <- sum(apply(m, 1, function(x) as.numeric(x[which(ns == x[alt.idx])])))
  paste(ref.count, alt.count, sep=",") 
}

mydf

x <- as.matrix(read.csv(text="start,A,T,G,C,REF,ALT,TYPE 
chr20:5363934,95,29,14,59,C,T,snp
chr5:8529759,24,1,28,41,G,C,snp
chr14:9620689,65,49,41,96,T,G,snp
chr18:547375,94,1,51,67,G,C,snp
chr8:5952145,27,80,25,96,T,T,snp
chr14:8694382,68,94,26,30,A,A,snp
chr16:2530921,49,15,79,72,A,T,snp:2530921
chr16:2530921,49,15,79,72,A,G,snp:2530921
chr16:2530921,49,15,79,72,A,T,snp:2530921flat
chr16:2533924,42,13,19,52,G,T,snp:2533924flat
chr16:2543344,4,13,13,42,G,T,snp:2543344flat
chr16:2543344,4,23,13,42,G,A,snp:2543344
chr14:4214117,73,49,18,77,G,A,snp
chr4:7799768,36,28,1,16,C,A,snp
chr3:9141263,27,41,93,90,A,A,snp", stringsAsFactors=FALSE))

结果：

       start           A    T    G    C    REF ALT TYPE              bam.AD  
 [1,] "chr20:5363934" "95" "29" "14" "59" "C" "T" "snp"             "59,29" 
 [2,] "chr5:8529759"  "24" " 1" "28" "41" "G" "C" "snp"             "28,41" 
 [3,] "chr14:9620689" "65" "49" "41" "96" "T" "G" "snp"             "49,41" 
 [4,] "chr18:547375"  "94" " 1" "51" "67" "G" "C" "snp"             "51,67" 
 [5,] "chr8:5952145"  "27" "80" "25" "96" "T" "T" "snp"             "80,80" 
 [6,] "chr14:8694382" "68" "94" "26" "30" "A" "A" "snp"             "68,68" 
 [7,] "chr16:2530921" "49" "15" "79" "72" "A" "T" "snp:2530921"     "49,15" 
 [8,] "chr16:2530921" "49" "15" "79" "72" "A" "G" "snp:2530921"     "49,79" 
 [9,] "chr16:2530921" "49" "15" "79" "72" "A" "T" "snp:2530921flat" "49,94"
[10,] "chr16:2533924" "42" "13" "19" "52" "G" "T" "snp:2533924flat" "19,13"  
[11,] "chr16:2543344" "42" "13" "13" "42" "G" "T" "snp:2543344flat" "13,55" 
[12,] "chr16:2543344" "42" "23" "13" "42" "G" "A" "snp:2543344"     "13,42" 
[13,] "chr14:4214117" "73" "49" "18" "77" "G" "A" "snp"             "18,73" 
[14,] "chr4:7799768"  "36" "28" " 1" "16" "C" "A" "snp"             "16,36" 
[15,] "chr3:9141263"  "27" "41" "93" "90" "A" "A" "snp"             "27,27"

【问题讨论】：

你能解释一下第 10 行“snp:2533924flat”：这里的 ALT 是 T，这是与 start 唯一的一行，所以“REF，ALT”不应该是“19,13” ?此外，如果 ALT 不是唯一的会发生什么，例如假设“snp:2543344flat”行（第 11 行）具有 ALT“A”（即 4），第 12 行（具有 A=42）也是如此 - 您使用哪个 A？
@mathematical.coffee 这些是基因型，在人类样本或小鼠身上进行。 SNP是单核苷酸多态性。长数字是该染色体上看到突变的碱基位置。 REF 是在该位置测序的规范碱基，而 ALT 是在抽样人群中看到的 SNP。输出显示 REF 碱基在总体中出现的次数，ALT 是“N”的替代碱基出现在该位置的次数。扁平线只是他们在该位置查看多个 ALT 碱基的摘要线。
@AMR 很好，但我不是遗传学家，不知道你在说什么，所有这些信息对我没有帮助。我确实看到您已经编辑了您的问题，但现在 mydf 定义与您想要的输出不匹配 - 在定义 x 的 mydf 代码中，您已将第 11 行和第 12 行更改为在“A”列，但在你想要的结果中你没有。请确保您的输入和期望的输出是一致的。
@mathematical.coffee ...我没有发布问题。我的评论是向您解释发布的数据。我了解数据，所以我了解他的问题与数据的相关性。虽然我没有用 r 编码，所以我无法提供 r 答案。
@AMR 抱歉。但是，这些信息对我的具体问题没有帮助（OP 在其中一行上犯了一个错误，我需要澄清它。他们已经修复了它，因为虽然引入了更多......）。

标签： r bioinformatics genetics

【解决方案1】：

这里有一种方法可以做到这一切，矢量化。

首先，请注意，无论类型如何，REF 都是相同的。我们可以通过使用 REF 作为矩阵中的坐标来快速查找它，例如第 1 行有 REF C，所以如果我们查找坐标 (1, "C")，我们会得到该行的 REF 值。

# the REFs are the same regardless of TYPE
rownames(x) <- 1:nrow(x)
ref <- x[cbind(1:nrow(x), x[, 'REF'])]

看看cbind(1:nrow(x), x[, 'REF'])：这只是一个坐标列表(row number, REF)，我们用它来查找REF号。

然后我们对 ALT 做同样的事情：

alt <- x[cbind(1:nrow(x), x[, 'ALT'])]

但是我们必须确保如果类型是“flat”，我们会将所有其他 ALT 添加到“flat”行的 ALT（如您所说，只有唯一的）。

首先，找出哪些行是平的：

which.flat <- grep('flat$', x[, 'TYPE'])

接下来，对于每个扁平行，查找具有相同“开始”（即x[, 'start'] == x[i, 'start'] 位）的其他行的 ALT，并排除具有重复 ALT 的行（即x[, 'ALT'] != x[i, 'ALT'] 位）。这里i是当前平线的索引。将它们全部添加到扁平线的 ALT 中。 sapply 只是将这一切矢量化为每条扁平线。

# add the other alts to the alt of the 'flat' line.
alt[which.flat] <- as.numeric(alt[which.flat]) + sapply(which.flat,
    function (i) {
        sum(as.numeric(alt[ x[, 'start'] == x[i, 'start'] &
             x[, 'ALT'] != x[i, 'ALT'] ]))
    })

现在我们只是粘贴在一起：

x <- cbind(x, bam.AD=paste(ref, alt, sep=','))

结果与您的结果相同，除了第 10 行，我认为您犯了一个错误 - 只有一行带有“chr16:2533924”并且它的 ALT 是“T”（值 13），所以 bam.AD 是"19,13"（你有 "19,42" 好像 ALT 是 "A"，但它不是）。

如果您必须在问题中坚持使用函数形式（非常缓慢且效率低下！），这与我所做的基本相同（因此，为什么您可以在没有 apply 调用的情况下完成它并完全跳过循环）：

flatCase

# get the other rows with the same 'start' and different 'ALT'
xx <- mat[mat[, 'start'] == x['start'] & mat[, 'ALT'] != x['ALT'], ,drop=F]
if (nrow(xx) > 0) {
  # grab all the alts as done before
  rownames(xx) <- 1:nrow(xx)
  alt <- alt + sum(as.numeric(xx[cbind(1:nrow(xx), xx[, 'ALT'])]))
 }

ref <- x[x['REF']]
return(paste(ref, alt, sep=','))
}

但是，如前所述，如果将其矢量化，则上面的整个代码将减少到几行，并且速度更快：

newBamAD <- function (x) {
    # the version above
    rownames(x) <- 1:nrow(x)
    ref <- x[cbind(1:nrow(x), x[, 'REF'])]
    alt <- x[cbind(1:nrow(x), x[, 'ALT'])]
    which.flat <- grep('flat$', x[, 'TYPE'])
    alt[which.flat] <- as.numeric(alt[which.flat]) + sapply(which.flat,
        function (i) {
            sum(as.numeric(alt[ x[, 'start'] == x[i, 'start'] &
                 x[, 'ALT'] != x[i, 'ALT'] ]))
        })
    cbind(x, bam.AD=paste(ref, alt, sep=','))
}

library(rbenchmark)
benchmark(
  bamAD=bamAD(x),
  newBamAD=newBamAD(x)
)
#       test replications elapsed relative user.self sys.self user.child sys.child
# 1    bamAD          100   0.082    3.905     0.072    0.004          0         0
# 2 newBamAD          100   0.021    1.000     0.020    0.000          0         0

矢量化版本几乎快 4 倍。

【讨论】：

谢谢，不过我上面贴的功能需要扩展一下。
它也不会将值“snp:+[0-9]”和扁平 ALT 字母添加到扁平线。
原始问题中您想要的输出不会将值“snp:+[0-9]”添加到 bamAD，那我为什么要这样做？我正在匹配您想要的输出。
对于 snp:+[0-9] 线我不需要添加，它们必须保持原样，但对于扁平线 snp:+[0-9] 和扁平线如果 ALT 在任何 snp:+[0-9] 中与扁平线不同，则必须添加线。如果每个 snp:+[0-9] ALT (T) 只有一个 snp:+[0-9]flat 如果它也有 ALT (T) 即相同的字母，则不需要添加。但是如果有两个 snp:+[0-9] 你只需将它们添加为 flat 的 ALT 值。无论如何，我非常感谢你的帮助！
这不正是我所做的吗？我的输出与您的输出完全匹配，除了您犯了错误的那一行（好吧，由于您编辑了问题，您的输入不再与您的输出匹配，因为您已将其中两行更改为不同的值，但我的代码 确实做你所要求的）。

【解决方案2】：

另一种方法：

# create dataframe
mydf <- as.data.frame(x, stringsAsFactors=FALSE)
# create temporary values based on REF and ALT
mydf$REFval <- diag(as.matrix(mydf[, mydf$REF]))
mydf$ALTval <- diag(as.matrix(mydf[, mydf$ALT]))

在下一步中，您说“如果 ALT 字母是唯一的”对 ALT 求和，但没有指定在 ALT 相同但值不同时使用哪个值。由于值相同，因此在您的示例数据集中并不重要，因此在下面的代码中，我假设要使用最后一个 ALT 值。

# sum up ALT values for all start ID
require(dplyr)
mydfs <- mydf %>% group_by(start, ALT) %>%
  summarize(ALTkeep=last(ALTval)) %>%  # assume keep last one if same ALT
  group_by(start) %>%
  summarize(ALTflat=sum(as.numeric(ALTkeep)))

# merge back into main dataframe
mydf <- left_join(mydf, mydfs)
# select ALT value for bam.AD depending on "flat$" in TYPE
mydf$bam.AD <- with(mydf,
  paste(REFval, ifelse(grepl("flat$", TYPE), ALTflat, ALTval), sep=","))

# optional clean up of temporary values
mydf <- mydf[, !(names(mydf) %in% c("REFval", "ALTval", "ALTflat"))]

你想要的输出

                                   start  A  T  G  C REF ALT            TYPE bam.AD
1                          chr20:5363934 95 29 14 59   C   T             snp  59,29
2                           chr5:8529759 24  1 28 41   G   C             snp  28,41
3                          chr14:9620689 65 49 41 96   T   G             snp  49,41
4                           chr18:547375 94  1 51 67   G   C             snp  51,67
5                           chr8:5952145 27 80 25 96   T   T             snp  80,80
6                          chr14:8694382 68 94 26 30   A   A             snp  68,68
7                          chr16:2530921 49 15 79 72   A   T     snp:2530921  49,15
8                          chr16:2530921 49 15 79 72   A   G     snp:2530921  49,79
9                          chr16:2530921 49 15 79 72   A   T snp:2530921flat  49,94
10                         chr16:2533924 42 13 19 52   G   T snp:2533924flat  19,13
11                         chr16:2543344  4 13 13 42   G   T snp:2543344flat  13,55
12                         chr16:2543344 42 23 13 42   G   A     snp:2543344  13,42
13                         chr14:4214117 73 49 18 77   G   A             snp  18,73
14                          chr4:7799768 36 28  1 16   C   A             snp  16,36
15                          chr3:9141263 27 41 93 90   A   A             snp  27,27

【讨论】：

谢谢，但是为什么会出现这个错误：Error: unknown column 'start'
您使用的是 x（您提供的样本数据，它是一个矩阵）还是 mydf（我使用我的代码的第一行将 x 转换为数据框）？
我正在使用 mydf，因为您已更改。我只是想复制你所做的。
嗯，很奇怪。你能做一个names(mydf) 看看第一列的名字是什么吗？应该是start。如果没有，可以重命名，或者将代码中对start的引用改为第一列的名称。
啊，很高兴听到这个消息。我认为这让我失去了你最初给我的“公认答案”，不是吗；/