使用字符串作为HARD中的输入计算矩阵上的比较

Question

我有两个DNA序列串，我想逐个字符比较这两个序列，以便得到一个具有comparisson值的矩阵。总的想法是有三个要点：

1. 如果存在互补AT(A在一个序列中，T在另一个序列(，则2/3
2. 如果存在互补CG(C在一个序列中，G在另一个序列(，则1
3. 否则，将返回0

例如，如果我有两个序列ACTG，那么结果将是：

| A   | C  | T   | G |
A| 0   | 0  | 2/3 | 0 |
C| 0   | 0  | 0   | 1 |
T| 2/3 | 0  | 0   | 0 |
G| 0   | 1  | 0   | 0 |

我看到这篇文章有一些帮助在Python中从等长字符串计算相似性/差异矩阵，如果你只使用4个核苷酸长的序列，它真的很有效-我尝试使用更大的序列，并打印出此错误：ValueError: shapes (5,4) and (5,4) not aligned: 4 (dim 1) != 5 (dim 0)

我有R中的代码，它是

##2.1 Separas los strings
seq <- "ACTG"
seq1 <- unlist(as.matrix(strsplit(seq,""),ncol=nchar(seq),
nrow=nchar(seq)))
a <- matrix(ncol=length(seq),nrow=length(seq))
a[,1] <- seq1
a[1,] <- seq1
b <- matrix(ncol=length(a[1,]),nrow=length(a[1,]))
for (i in seq(nchar(seq))){
for (j in seq(nchar(seq))){
if (a[i,1] == "A" & a[1,j] == "T" | a[i,1] == "T" & a[1,j] == "A"){
b[[i,j]] <- 2/3
} else if (a[i,1] == "C" & a[1,j] == "G" | a[i,1] == "G" & a[1,j] == "C"){
b[[i,j]] <- 1
} else
b[[i,j]] <-  0
}

但我无法在python中获取代码。

Answer 1

我认为你让它变得比需要的更难。

import numpy as np
seq1 = 'AACCTTGG'
seq2 = 'ACGTACGT'
matrix = np.zeros((len(seq1),len(seq2)))
for y,c2 in enumerate(seq2):
for x,c1 in enumerate(seq1):
if c1+c2 in ('TA','AT'):
matrix[x,y] = 1.
elif c1+c2 in ('CG','GC'):
matrix[x,y] = 2/3
print(matrix)