我有一个名为"test.fa"的文件，上面写着：
>Sequence 1
TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
>Sequence 2
CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
>Sequence 3
TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG
>Sequence 4
AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC
这个测试文件只有8行，但可以有更多行。我现在也不知道所有序列的长度。因此，首先我计算行数，并根据行数创建一个矩阵。
import numpy as np
filename = "test.fa"
with open(filename, 'r') as file:
n1 = 0
n2 = 0
for line in file.readlines():
if line[0] != '>':
m1 = 1
n1 = n1 + m1
else:
m2 = 1
n2 = n2 + m2
seq = np.chararray((n1,2),itemsize = 99)
接下来，我将值添加到矩阵中。
with open(filename, 'r') as file:
n3 = -1
n4 = -1
for line in file.readlines():
if line[0] != '>':
m3 = 1
n3 = n3 + m3
seq[n3,1] = line
else:
m4 = 1
n4 = n4 + m4
seq[n4,0] = line
如果我调用"seq"，我得到：
chararray([[b'>Sequence 1', b'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT'],
[b'>Sequence 2', b'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG'],
[b'>Sequence 3', b'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'],
[b'>Sequence 4', b'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC']], dtype='|S99')
seq[2,1]: b'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'
seq[2,1][0] : 84
seq[2,1][0:8] : b'TCGACCCT'
seq[2,8][8] : 67
seq[2,1][8:9] : C
这并不是管理这些数据的最佳方式(为什么[2,1][0]返回84？(。我不确定项目大小为99的np.chararray是最好的方法。
我的问题是：组织/管理这些数据的更好方法是什么？我最终需要计算每个核苷酸的出现次数(即有多少As、Ts、Cs、Gs(，并从每个序列中提取子串。就上下文而言，这与期望最大化和吉布斯采样有关。
它返回84，而类型是S99，因为找到的最大长度是99，所以类型应该是S99。你从文件中读取的方式不太正确，管理你的数据的更好方法是使用dict，我不知道它是否适合你的下一步需求，但它是：
seq={}
with open("file") as f:
for line in f:
seq[line]=next(f)
这里它将一次读取两行，因此第一个line具有序列1，下一个(f(具有seqTGC。。。。。，如果发生文件没有一致排序的情况，请使用if语句
通过使用np.chararray，您可以获得一个具有S数据类型的数组，字节字符串：
In [145]: x = np.chararray((1,),10)
In [146]: x[0]='ACT'
In [147]: x
Out[147]: chararray([b'ACT'], dtype='|S10')
In [148]: x[0]
Out[148]: b'ACT'
In [149]: x[0][0]
Out[149]: 65
In [150]: x[0][1]
Out[150]: 67
In [151]: x[0][2]
Out[151]: 84
这是每个字符一个字节。A的ASCII码为65；CCD_ 6是19个字符。
你看到chararray注释了吗：
chararray类的存在是为了向后兼容Numarray，不建议用于新的开发。从numpy开始1.4，如果需要字符串数组，建议使用dtypeobject_、string_或unicode_，并使用免费函数在CCD_ 13模块中用于快速矢量化字符串操作。
In [159]: x = np.array([b'ACT'])
In [160]: x
Out[160]: array([b'ACT'], dtype='|S3')
In [161]: x[0]
Out[161]: b'ACT'
In [162]: x[0][0]
Out[162]: 65
如果我们使用Py3默认字符串类型(unicode(：
In [154]: x = np.array(['ACT'])
In [155]: x
Out[155]: array(['ACT'], dtype='<U3')
In [156]: x[0]
Out[156]: 'ACT'
In [157]: x[0][0]
Out[157]: 'A'
In [158]: x[0][1]
Out[158]: 'C'
===
我认为在列表或字典中收集值会更简单：
In [165]: txt = """>Sequence 1
...: TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
...: >Sequence 2
...: CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
...: >Sequence 3
...: TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG
...: >Sequence 4
...: AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC"""
In [166]: alist = []
In [167]: for row in txt.splitlines():
...: if row.startswith('>'):
...: x = [row]
...: else:
...: x.append(row)
...: alist.append(x)
...:
In [168]: alist
Out[168]:
[['>Sequence 1', 'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT'],
['>Sequence 2', 'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG'],
['>Sequence 3', 'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG'],
['>Sequence 4', 'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC']]
In [169]: dd = {key:val for key,val in alist}
In [170]: dd
Out[170]:
{'>Sequence 1': 'TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT',
'>Sequence 2': 'CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG',
'>Sequence 3': 'TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG',
'>Sequence 4': 'AAAACACTTGAGGGAGCAGATAACTGGGCCAACCATGACTC'}
所以这些字符串的基本处理：
In [171]: [len(val) for key,val in alist]
Out[171]: [41, 41, 44, 41]
In [172]: [val.count('T') for key,val in alist]
Out[172]: [16, 7, 6, 6]
In [173]: [val.count('C') for key,val in alist]
Out[173]: [12, 21, 15, 10]
我刚刚找到了这篇论文，也许它会帮助你继续寻找。https://pastel.archives-ouvertes.fr/tel-01762479/document
TL；DRhttps://en.wikipedia.org/wiki/LCP_arrayhttps://en.wikipedia.org/wiki/Suffix_array
pythonhttps://louisabraham.github.io/notebooks/suffix_arrays.html