r - 计算成对汉明距离并保留较小的一个

Question

我有一个矩阵:

A = matrix( c(1, 0, 0, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow = TRUE)
dimnames(A) = list(c("Taxa1", "Taxa2", "Taxa3"), c("A1", "B1", "C1"))
df <- data.frame("ID" = c("A1", "B1", "C1"), "Triplicate" = c("T1", "T1", "T1"))

我想计算一个类群的 A1 和 B1 之间以及 A1 和 C1 之间的汉明距离，并保留两者中的最小值作为 A1 的指定值。然后计算 B1 和 A1(我们已经这样做了)以及 B1 和 C1 之间的汉明距离，并保留 B1 的最小值，依此类推 C1。注意 sample A1、B1、C1 是 T1 所示的同一三份 sample 的一部分。我将有类似的其他 Triplicate 样本 T2、T3 或 T4，并且我想根据 Triplicate 列值对样本进行分组，以计算成对汉明

结果矩阵为:

df$Hamming <- c(2, 2, 2)

由于A1、B1和A1之间的距离，C1为2，因此保留该值为2。

PS:汉明距离的简单 1 分钟描述: https://www.youtube.com/watch?v=P02mJhS9qQ4

添加我正在使用的确切数据: https://www.dropbox.com/s/wrlwmdipeyhbcok/Hamming.RData?dl=0

Answer 1

位翻转是通过xor()函数来识别的，位翻转的总数是通过将xor()函数的结果相加得出的。我没有优化 hamm_dist_min() 中的代码。

xor(0,0)
# [1] FALSE
xor(1,1)
# [1] FALSE
xor(0,1)
# [1] TRUE
xor(1,0)
# [1] TRUE

根据 OP 的要求，计算两个方向的汉明距离。例如:AB、BA、AC、CA、BC 和 CB，形成三联，即 ABC、BCA 和 CAB。如果您只需要一个方向，例如:AB、AC 和 BC，您可以使用 combn() 函数在 hamm_dist_min() 函数内设置列号。

数据:

A = matrix( c(1, 0, 0, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow = TRUE)
dimnames(A) = list(c("Taxa1", "Taxa2", "Taxa3"), c("A1", "B1", "C1"))
df <- data.frame("ID" = c("A1", "B1", "C1"), "Triplicate" = c("T1", "T1", "T1"))

汉明距离

# minimum of hamming distance
hamm_dist_min <- function(data)
{
  # setup combinations of column numbers
  n_col <- ncol(data)
  x <- expand.grid(seq_len(n_col), seq_len(n_col))
  x <- x[ x[, 1] != x[, 2], ]
  x <- x[order(x[, 1]), ]
  x <- split(x, cut(x[, 1], breaks = c(0, seq_len(n_col)), labels = colnames(data) ))

  # minimum of hamming distance
  h_d <- unlist(lapply(x, function(y){
    min( colSums(apply(y, 1, function(z) xor(data[, z[1]], data[, z[2]]))))
  }))
  return(h_d)
}

hamm_dist_min(data = A)
# A1 B1 C1 
# 2  2  2

df$Hamming <- hamm_dist_min(data = A)
df
#   ID Triplicate Hamming
# 1 A1         T1       2
# 2 B1         T1       2
# 3 C1         T1       2

Youtube example :

df1 <- matrix( c(0,1,1,1,1,1,0,0), ncol = 2, byrow = FALSE)
colnames(df1) <- LETTERS[1:2]
hamm_dist_min(data = df1)
# A B 
# 3 3 

<小时/>

编辑:基于问题中添加的新数据集。

注意:如果某个样本类型只有一列，则汉明距离取 0，因为我们至少需要 2 列来计算汉明距离。查看df 的一式三份列中的T71。您可以返回 NA，而不是值 0，这意味着 不可用

load("Hamming.RData")

# setup unique colnames pattern
col_list <- unique(unlist(lapply( colnames(A), function(x){
  substr(x = x, start = 1, stop = nchar(x) - 1)
} )
))

# get hamming distance
my_results <- lapply( col_list, function(x){
  cols_x <- grep(x, colnames(A) )
  if(length(cols_x) == 1 ){  # return 0 for one column
    return( setNames( object = rep( 0, length(cols_x)), nm = colnames(A)[cols_x]))
  } else{ # return minimum of hamming distance
    return(hamm_dist_min(data = A[, cols_x]))
  }
})

# get triplicate id
triplicate <- paste0( "T", rep(seq_along(my_results), 
                               lengths(my_results)))

# final data
my_results <- unlist(my_results)
df <- data.frame( SampleID = names( my_results ),
                  Hamming = my_results,
                  Triplicate = triplicate,
                  stringsAsFactors = FALSE )

head(df)
#                        SampleID Hamming Triplicate
# Affy22_MDA_1       Affy22_MDA_1       2         T1
# Affy22_MDA_2       Affy22_MDA_2       2         T1
# Affy22_MDA_3       Affy22_MDA_3       3         T1
# GutRef001_MDA_1 GutRef001_MDA_1       4         T2
# GutRef001_MDA_2 GutRef001_MDA_2       4         T2
# GutRef001_MDA_3 GutRef001_MDA_3       6         T2