当替换某些字符位置( R)时，不要计算字符中的某些位置。

Question

所以，我知道我的头衔有点混乱，但我希望你能帮我解决这个问题。

我有一个数据框架df，其中一个列是一个RNA序列比对。该列的类是一个字符。

然后我还有其他的列："Allele_1"，"Allele_2“，它代表RNA序列中单个位置的变体(第1列)，该位置由第3列(”位置“)表示。然而，这些位置不考虑"-"，例如，在第2行中，等位基因的位置是U--ACCGU--G----UAUUUGAU--CTAD 而不是 U--ACCGU--G----UAUUUGAU--CTAD.

sequence                         Allele_1   Allele_2     Position
UAAGGCUCA----UAGGCAGAU--AUaa     A          U            3
U--ACCGU--G----UAUUUGAU--CTAD    C          G            5
cctaACCGU-UUAGCC---------T       U          C            2

列1中序列的长度可以是可变的。

我想要做的是在特定位置用“位置”替换字符的特定字母，用"Allele_1“和"Allele_2”替换字符。例如，如果该位置与"Allele_2“匹配，则我希望将其替换为"Allele_2"，反之亦然。

我试过：

substr(df[,"sequence"], 
  start = df[,"Position"], 
  stop = df[,"Position"]) <- df[,"Allele_1"]

但是，由于我的职位栏没有考虑到"-"，所以它在错误的位置替换。例如，从第2行到第2行，它将替换这里的U--ACCGU--G----UAUUUGAU--CTADinstead of here U--ACCGU--G----UAUUUGAU--CTAD.

此外，我还没有弄清楚如何做“位置匹配”"Allele_2"，然后我想用"Allele_2“代替它，反之亦然。

sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

真希望你能帮我解决这个问题！

干杯!

更新:对不起，如果该位置与"Allele_1“匹配，那么我希望将其替换为"Allele_2"，反之亦然”而不是“Allele_2，然后我想用"Allele_2”替换它。

Answer 1

这里有两个选择。两者都是区分大小写的，因此在第三个序列中不替换任何东西。如果您不希望它们是这样的，请将适当的变量包装在ifelse中的toupper中。

strsplit

您可以将每个序列拆分为一个字母向量，然后可以根据该向量直接检查相等性。在mapply (多变量版本的sapply )中实现

df$new_seq <- mapply(function(seq, a1, a2, pos){
    seq <- strsplit(seq, '')[[1]]    # split into letters
    to_replace <- seq[seq != '-'][pos]    # identify allele to replace
    # assign appropriate replacement to subset
    seq[seq != '-'][pos] <- ifelse(a1 == to_replace, 
                                   a2, ifelse(a2 == to_replace, 
                                              a1, to_replace))
    paste(seq, collapse = '')    # reassemble vector to string
}, df$sequence, df$Allele_1, df$Allele_2, df$Position)

df
##                        sequence Allele_1 Allele_2 Position                       new_seq
## 1  UAAGGCUCA----UAGGCAGAU--AUaa        A        U        3  UAUGGCUCA----UAGGCAGAU--AUaa
## 2 U--ACCGU--G----UAUUUGAU--CTAD        C        G        5 U--ACCCU--G----UAUUUGAU--CTAD
## 3    cctaACCGU-UUAGCC---------T        U        C        2    cctaACCGU-UUAGCC---------T

如果愿意，可以将操作分解为多个步骤，将每个步骤的结果分配给一个变量。

sub (regex)

如果您对regex感到满意，您可以组装表达式来提取所述等位基因，然后用适当的替换替换它：

df$to_replace <- mapply(function(seq, pos){
    sub(paste0('(?:-*(?:\\w)-*){', pos - 1, '}(\\w).*'), '\\1', seq)
}, df$sequence, df$Position)

df$new_seq <- mapply(function(seq, pos, a1, a2, to_rpl){
    replacement <- ifelse(to_rpl == a1, a2, ifelse(to_rpl == a2, a1, to_rpl))
    sub(paste0('((?:-*(?:\\w)-*){', pos - 1, '})\\w(.*)'), 
        paste0('\\1', replacement, '\\2'), 
        seq)
}, df$sequence, df$Position, df$Allele_1, df$Allele_2, df$to_replace)

df[-5]
##                        sequence Allele_1 Allele_2 Position                       new_seq
## 1  UAAGGCUCA----UAGGCAGAU--AUaa        A        U        3  UAUGGCUCA----UAGGCAGAU--AUaa
## 2 U--ACCGU--G----UAUUUGAU--CTAD        C        G        5 U--ACCCU--G----UAUUUGAU--CTAD
## 3    cctaACCGU-UUAGCC---------T        U        C        2    cctaACCGU-UUAGCC---------T

Answer 2

数据

我创建了第一个小数据框架和一个更大更真实的数据框架。

df <- data.frame(sequence = c('UAAGGCUCA----UAGGCAGAU--AUaa',
                                  'U--ACCGU--G----UAUUUGAU--CTAD',
                                  'cctaACCGU-UUAGCC---------T'),
                     Allele_1 = c('A', 'C', 'U'),
                     Allele_2 = c('U', 'G', 'C'),
                     Position = c(3, 5, 2))

df_big <- data.frame(sequence = rep(c(paste(rep('UAAGGCUCA----UAGGCAGAU--AUaa', 100), collapse=''),
                                  paste(rep('U--ACCGU--G----UAUUUGAU--CTAD', 100), collapse=''),
                                  paste(rep('cctaACCGU-UUAGCC---------T', 100), collapse='')), 100),
                 Allele_1 = rep(c('A', 'C', 'U'), 100),
                 Allele_2 = rep(c('U', 'G', 'C'), 100),
                 Position = rep(c(1000, 1000, 1000), 100))

函数

我创建了一个函数，用于在多重对齐后返回“新”位置(即不包括-)，它适用于向量，如下所示；也是一个更快但不太安全的函数：

find_pos <- function(string, pos) {
  vectorized <- data.frame(string=string, pos=pos)
  sapply(seq_len(nrow(vectorized)), function(i) {
    guess <- vectorized[i,'pos']; oldguess <- 0
    offset <- -1; trynum <- 0
    while(offset != 0 & trynum < 500) {
      offset <- nchar(gsub('[^-]', '', substr(vectorized[i,'string'], oldguess+1, guess)))
      oldguess <- guess
      guess <- oldguess + offset
      trynum <- trynum+1
    }
    return(guess)
  })
}

find_pos_unsafe <- function(string, pos) {
  sapply(seq_along(string), function(i) {
    guess <- pos[i]; oldguess <- 0
    offset <- -1
    while(offset != 0) {
      offset <- nchar(gsub('[^-]', '', substr(string[i], oldguess+1, guess)))
      oldguess <- guess
      guess <- oldguess + offset
    }
    return(guess)
  })
}

第一个变量可以与不同的长度变量一起使用，如下所示(但这会带来一定的开销)：

> find_pos(string, 1:5)
[1] 1 4 5 6 7

出于基准测试的目的，我已经包装了在函数中获得解决方案所需的其他代码。同样，有两个表单，一个调用更快的匹配函数，另一个调用更安全的版本：

ds440 <- function(df) {
pos <- find_pos(df$sequence, df$Position)
toswap <- ifelse(substr(df$sequence, pos, pos)==df$Allele_1, as.character(df$Allele_2), #if A1, A2
                 ifelse(substr(df$sequence, pos, pos)==df$Allele_2, as.character(df$Allele_1), #If A2, A1
                        substr(df$sequence, pos, pos))) # else keep same

df$replaced <- as.character(df$sequence)
substr(df$replaced, pos, pos) <- as.character(toswap)

df
}

ds440_quick <- function(df) {
  pos <- find_pos_unsafe(df$sequence, df$Position)
  toswap <- ifelse(substr(df$sequence, pos, pos)==df$Allele_1, as.character(df$Allele_2), #if A1, A2
                   ifelse(substr(df$sequence, pos, pos)==df$Allele_2, as.character(df$Allele_1), #If A2, A1
                          substr(df$sequence, pos, pos))) # else keep same

  df$replaced <- as.character(df$sequence)
  substr(df$replaced, pos, pos) <- toswap

  df
}

改编自@alistaire的其他职能

alistaire_split <- function(df) {
  df$new_seq <- mapply(function(seq, a1, a2, pos){
    seq <- strsplit(seq, '')[[1]]    # split into letters
    to_replace <- seq[seq != '-'][pos]    # identify allele to replace
    # assign appropriate replacement to subset
    seq[seq != '-'][pos] <- ifelse(a1 == to_replace, 
                                   a2, ifelse(a2 == to_replace, 
                                              a1, to_replace))
    paste(seq, collapse = '')    # reassemble vector to string
  }, as.character(df$sequence), as.character(df$Allele_1), as.character(df$Allele_2), df$Position)

  df
}

alistaire_sub <- function(df) {
  df$to_replace <- mapply(function(seq, pos){
    sub(paste0('(?:-*(?:\\w)-*){', pos - 1, '}(\\w).*'), '\\1', seq)
  }, df$sequence, df$Position)

  df$new_seq <- mapply(function(seq, pos, a1, a2, to_rpl){
    replacement <- ifelse(to_rpl == a1, a2, ifelse(to_rpl == a2, a1, to_rpl))
    sub(paste0('((?:-*(?:\\w)-*){', pos - 1, '})\\w(.*)'), 
        paste0('\\1', replacement, '\\2'), 
        seq)
  }, df$sequence, df$Position, df$Allele_1, df$Allele_2, df$to_replace)

  df[-5]
}

结果

请注意上述匹配中的大小写敏感性。如果您不关心大小写，在检查是否相等之前，请使用toupper或类似的方法。

ds440(df)
                       sequence Allele_1 Allele_2 Position                      replaced
1  UAAGGCUCA----UAGGCAGAU--AUaa        A        U        3  UAUGGCUCA----UAGGCAGAU--AUaa
2 U--ACCGU--G----UAUUUGAU--CTAD        C        G        5 U--ACCCU--G----UAUUUGAU--CTAD
3    cctaACCGU-UUAGCC---------T        U        C        2    cctaACCGU-UUAGCC---------T

基准测试

上述职能已予界定。

library(microbenchmark)

microbenchmark(ds440(df), ds440_quick(df), alistaire_split(df), alistaire_sub(df))
Unit: microseconds
                expr     min       lq     mean   median       uq      max neval cld
           ds440(df) 727.157 747.0065 801.3529 764.1475 781.6400 3658.410   100   c
     ds440_quick(df) 339.784 354.7140 364.4484 364.1440 370.7955  421.400   100  b 
 alistaire_split(df) 138.929 144.9500 150.6806 148.5890 153.8570  261.136   100 a  
   alistaire_sub(df) 815.853 833.6630 855.2414 844.0770 857.2130 1499.370   100   c

microbenchmark(ds440(df_big), ds440_quick(df_big), alistaire_split(df_big))
Unit: milliseconds
                    expr      min       lq     mean   median       uq      max neval cld
           ds440(df_big) 195.4233 199.2827 204.4030 204.9039 208.0195 216.1879   100   c
     ds440_quick(df_big) 136.5985 139.3442 143.7216 145.1585 146.7837 153.9201   100 a  
 alistaire_split(df_big) 138.9117 146.3977 150.8772 148.2299 151.0723 278.7308   100  b 

在小的例子中，明显的赢家是@alistaire的拆分函数，但是当df变得更大时，alistaire_sub函数会中断(regex不会处理>999)，而我的ds440_quick函数实际上计算得更快一些。