python计算序列列表中存在和不存在子字符串的数量答案

【问题标题】：python counting number of presence and absence of substrings in list of sequencespython计算序列列表中存在和不存在子字符串的数量
【发布时间】：2013-11-29 23:44:37
【问题描述】：

您可以在这里获取数据！ 2shared底部下载

我正在用 python 分析生物数据。

我已经编写了一个代码来在长字符串列表中查找匹配的子字符串。

子串在一个列表中，长度为 7 个核苷酸。

所以在列表中，从 AAAAAAA 到 TTTTTTT，存在 16384 个主题（子字符串），排列 A、C、G、T。

这段代码有一个for循环，用于子字符串列表和嵌套在里面的长字符串列表。

它工作正常，但由于列表有12000行，代码处理很慢。

换句话说，提供有关 AAAAAAA 的信息以及下一个 AAAAAAC 需要 2 分钟。

因此需要 16384 个图案在 2 分钟内通过 12000 行，这将需要 (16384*2 == 32768 分钟 -> 546 小时 -> 22 天...)

我正在使用 scipy 和 numpy 来获取 Pvalues。

我想要的是计算序列列表中存在和不存在子字符串的数量

长字符串列表和代码如下：

list_of_lists_long  =  [
[BGN,    -0.054,     AGGCAGCTGCAGCCACCGCGGGGCCTCAGTGGGGGTCTCTGG....]
[ABCB7,  0.109,      GTCACATAAGACATTTTCTTTTTTTGTTGTTTTGGACTACAT....]
[GPR143, -0.137,     AGGGGATGTGCTGGGGGTCCAGACCCCATATTCCTCAGACTC....]
[PLP2,   -0.535,     GCGAACTTCCCTCATTTCTCTCTGCAATCTGCAAATAACTCC....]
[VSIG4,  0.13,       AAATGCCCCATTAGGCCAGGATCTGCTGACATAATTGCCTAG....]
[CCNB3,  -0.071,     CAGCAGCCACAGGGCTAAGCATGCATGTTAACAGGATCGGGA....]
[TCEAL3, 0.189,      TGCCTTTGGCCTTCCATTCTGATTTCTCTGATGAGAATACGA....]
....] #12000 lines

是否有更快的逻辑可以更快地执行代码？？

我需要你的帮助！

提前谢谢你。

======================================

有没有更简单的方法，无需实现任何其他东西？

我认为模式匹配的迭代是问题...

我想要找到的是，长度为 7 的基序在整个序列列表中出现了多少次，而且也没有出现！！！。所以如果一个主题存在于一个字符串中，它是 TRUE 作为 bool，然后增加一个值 AND FALSE，然后增加另一个值。

不是字符串中图案的数量。

【问题讨论】：

不重要，但您的代码中有一些错误，例如for sMotif in motifs，可能您的意思是motifList 和listoflists_long ，请将字符串放在引号之间。
看看后缀树：daimi.au.dk/~mailund/suffix_tree.html
@elyase，谢谢，错误已修复。
一个快速提示是使用if subs in s 而不是find 检查字符串。这应该可以让你的速度提高 2 倍。
每个单独的字符串有多长？

标签： python string numpy scipy bioinformatics

【解决方案1】：

我们编写了一个名为 Sylamer 的工具。它计算相同给定长度的单词的出现次数。默认情况下，它在排序基因的上下文中计算超几何测试，或者它可以只输出计数。它可以对给定长度的所有可能的单词执行此操作，但也可以指定较小的单词列表。它是用 C 语言编写的，并针对快速处理进行了非常优化。通过链接到 GNU 科学库，我们使超几何计算变得非常快。

【讨论】：

OP 谈论 python。您的库是否有 python 绑定？否则这个答案是在错误的地方。
我在 18 个月前写了这个答案；也许我的意思是暗示 OP 可以查看代码（它使用位移位和位掩码方法）以获得想法，但让它有点太含蓄了。我认为将自己先验地限制在特定工具上通常不是一个好主意，除非它涉及家庭作业（这在其他方面是有问题的）。我在超几何方面做了一些扩展，因为这是分析生物序列数据的常用方法（Sylamer 专门为此编写）。

【解决方案2】：

我猜问题出在 Fisher 的精确测试安装上。如果我在 for 循环之外计算 P 值，那么计算会变得更快。

【讨论】：

【解决方案3】：

一个干净且非常快速的方法（大约 15 秒使用 OP 的数据）是使用 scikits-learn 的 CountVectorizer，因为它在后台使用 numpy，例如：

from sklearn.feature_extraction.text import CountVectorizer

def make_chunks(s):
    width = 2
    return [s[i:i+width] for i in range(len(s)-width+1)]

l = ['ATTGCGGCTCACGAA', 'ACCTAGATACGACGG', 'CCCCTGTCCATGGTA']

vectorizer = CountVectorizer(tokenizer=make_chunks)
X = vectorizer.fit_transform(l)

现在X 是一个稀疏矩阵，将所有可能的块作为列，将序列作为行，其中每个值都是给定块在每个序列中出现的次数：

>>> X.toarray()
# aa ac ag at ca cc cg ...
[[1  1  0  1  1  0  2  1 1 2 1 0 0 1 1 1]      # ATTGCGGCTCACGAA
 [0  3  1  1  0  1  2  1 2 0 1 0 2 0 0 0]      # ACCTAGATACGACGG
 [0  0  0  1  1  4  0  1 0 0 1 2 1 1 2 0]]     # CCCCTGTCCATGGTA

>>> (X.toarray()>0).astype(int)  # the same but counting only once per sequence
[[1 1 0 1 1 0 1 1 1 1 1 0 0 1 1 1]
 [0 1 1 1 0 1 1 1 1 0 1 0 1 0 0 0]
 [0 0 0 1 1 1 0 1 0 0 1 1 1 1 1 0]]

>>> vectorizer.get_feature_names()     # the columns(chunks)
[u'aa', u'ac', u'ag', u'at', u'ca', u'cc', u'cg', u'ct', u'ga', u'gc', u'gg', u'gt', u'ta', u'tc', u'tg', u'tt']

现在您可以对列求和、屏蔽非值或您需要做的任何操作，例如：

>>> X.sum(axis=0)
[[1 4 1 3 2 5 4 3 3 2 3 2 3 2 3 1]]

最后要找出给定motif出现的次数，你必须找到相应motif/chunk的索引，然后在之前的总和中求值：

>>> index = vectorizer.vocabulary_.get('ag')    # 'ag' is your motif
2   # this means third column

在您的情况下，您必须将列表分为两部分（正值和负值）以包含向下条件。我使用 DSM 答案中的列表进行了快速测试，在我的计算机上运行大约需要 3-4 秒。如果我使用 12 000 个长度为 4000 个的序列，则只需不到一分钟。

编辑：使用 OP 数据的整个代码可以在 here 找到。

【讨论】：

感谢@elyase 的帮助，但您能根据我的数据调整代码吗？这样我可以更容易更清楚地得到你的吗？
我认为使用小测试样本的解释更清楚，我用我用来测试你的数据的代码添加了一个编辑。
我已经安装了 scikits-learn 并测试了你的代码。它非常快，但没有返回我所期望的。它只返回一个没有主题信息的数组。我期望的是有多少 AAAAAAA 在下调序列列表中，而不是在列表中以及上调序列中。所以最终打印出来的将是'motif'、'motif_down'、'not_motif_down'、'motif_notdown'、'notmotif_notdown'。非常感谢您的帮助，但我是新手，有点难以理解。
顺便说一句，实际上，下调是-0.5的截止值，而不是负数。不包括-0.5。但是 [down] 和 [~down] 做了什么？他们是否实现了区分位于我的第二个数据位置的负值和正值的功能？
在 python 中你可以逐行运行代码并查看每个变量包含的内容。 down 是一个布尔掩码，~down 只是它的否定。 Numpy 掩码是具有布尔值的简单数组，用于根据条件选择数据。老实说，我不完全理解您想要的最终输出（列表？字典？）。但在我的代码中，您会看到 names 是一个图案数组，counts 是一个大小相同的数组，具有相应的计数，因此 names[0] -> counts[0], names[1] - > 计数[1]。看看变量，可能有一个更短的例子，你会看到它。

【解决方案4】：

基本上你的问题是序列比较。按顺序寻找基序是生物信息学中的一个基本问题。我认为您可以搜索一些现有的算法或软件包。我在 google 中搜索了关键字“motif match”，这是我在第一页找到的： http://biowhat.ucsd.edu/homer/motif/ http://meme.nbcr.net/meme/doc/mast.html http://www.biomedcentral.com/1471-2105/8/189 http://www.benoslab.pitt.edu/stamp/

【讨论】：

【解决方案5】：

您的代码有几个奇怪的地方。

您所说的“排列”看起来更像笛卡尔积，可以使用itertools.product 计算。
您的 OddsRatio, PValue 和 Enrichment 计算在循环内，但计数的归零和 print 都不是，这意味着您正在为每个新行累积计算它们但不做任何事情有了这些信息。
在典型情况下，您会多次重新计算 i[2].find(sMotif)。该结果不会被缓存。

假设我了解您要计算的数字——我很可能是错的，因为你正在做的有几件事我不明白——我会颠倒逻辑。不要循环遍历每个主题并尝试在每一行中计算它，而是遍历每一行并查看那里有什么。这将是大约 7* 行数，而不是图案数 * 行数。

例如：

import random
from itertools import product
from collections import defaultdict, Counter

N = 12000
datalength = 400
listoflists = [[str(i), random.uniform(-1, 1), 
                ''.join([random.choice('AGCT') for c in range(datalength)])]
               for i in range(N)]

def chunk(seq, width):
    for i in range(len(seq)-width+1):
        yield seq[i:i+width]

def count_motifs(datatriples, width=7):
    motif_counts_by_down = defaultdict(Counter)
    nonmotif_counts_by_down = defaultdict(Counter)
    all_motifs = set(''.join(p) for p in product('AGCT',repeat=width))
    for symbol, value, sdata in datatriples:
        down = value < -0.5

        # what did we see?
        motifs_seen = set(chunk(sdata, width))
        # what didn't we see?
        motifs_not_seen = all_motifs - motifs_seen

        # accumulate these
        motif_counts_by_down[down].update(motifs_seen)
        nonmotif_counts_by_down[down].update(motifs_not_seen)

    return motif_counts_by_down, nonmotif_counts_by_down

（我降低行长只是为了使输出更快；如果行长 10 倍，则代码需要 10 倍。）

这会在我的慢速笔记本电脑上产生（插入一些换行符后）：

>>> %time mot, nomot = count_motifs(listoflists, 7)
CPU times: user 1min 50s, sys: 60 ms, total: 1min 50s
Wall time: 1min 50s

所以我认为完整的问题大约需要 20 分钟，这对于这么少的代码来说还不错。（我们可以改为通过算术来加速motifs_not_seen 部分，但无论如何这只会让我们获得两倍。）

在更小的情况下，更容易看到输出：

>>> mot, nomot = count_motifs(listoflists, 2)
>>> mot
defaultdict(<class 'collections.Counter'>, 
{False: Counter({'CG': 61, 'TC': 58, 'AT': 55, 'GT': 54, 'CA': 53, 'GA': 53, 'AC': 52, 'CT': 51, 'CC': 50, 'AG': 49, 'TA': 48, 'GC': 47, 'GG': 45, 'TG': 45, 'AA': 43, 'TT': 40}), 
True: Counter({'CT': 27, 'GT': 26, 'TC': 24, 'GC': 23, 'TA': 23, 'AC': 22, 'AG': 21, 'TG': 21, 'CC': 19, 'CG': 19, 'CA': 19, 'GG': 18, 'TT': 17, 'GA': 17, 'AA': 16, 'AT': 16})})
>>> nomot
defaultdict(<class 'collections.Counter'>, 
{False: Counter({'TT': 31, 'AA': 28, 'GG': 26, 'TG': 26, 'GC': 24, 'TA': 23, 'AG': 22, 'CC': 21, 'CT': 20, 'AC': 19, 'GA': 18, 'CA': 18, 'GT': 17, 'AT': 16, 'TC': 13, 'CG': 10}), 
True: Counter({'AA': 13, 'AT': 13, 'GA': 12, 'TT': 12, 'GG': 11, 'CC': 10, 'CA': 10, 'CG': 10, 'AG': 8, 'TG': 8, 'AC': 7, 'GC': 6, 'TA': 6, 'TC': 5, 'GT': 3, 'CT': 2})})

【讨论】：

感谢您的友好回答，我试图找出长度为 7 的基序在整个序列列表中出现的次数。因此，如果一个主题存在于一个字符串中，这是真的，那么增加这个值。我应该在这个问题上提到它。谢谢帝斯曼让我理解你的。
这就是count_motifs(listoflists, 7) 计算的内容。它返回在 down 案例中看到主题的次数，以及没有看到的次数。它基本上返回nMotif_Down、nMotif_NotDown、nNotMotif_Down和nNotMotif_NotDown。
我一直在使用您的代码，但在打印出宽度为 7 的 defaultdict 时遇到了问题。我应该迭代每个键和值，但它不起作用。你能给我什么建议吗？
我认为您的最接近我想要的，但是您能否将最终结果打印出每个主题，motif_counts_by_down 和 nonmotif_counts_by_down？请。
@Karyo，您介意共享您的文件以对结果进行基准测试吗？

【解决方案6】：

很好的问题。这是一个经典的计算机科学问题。是的，确实有更好的算法。你的处理每个长字符串 16384 次。更好的方法是每个长字符串只处理一次。

与其在每个长字符串中搜索每个主题，不如只记录每个长字符串中出现的主题。例如，如果您在以下字符串中搜索长度为 2 的图案：

s = 'ACGTAC'

然后您可以在长度为 2 的子字符串上运行一个循环，并记录哪些子字符串出现在 dict 中：

motifAppearances = {}
for i in range(len(s)-1):
    motif = s[i:i+2]                   # grab a length=2 substring
    if motif not in motifAppearances:
        motifAppearances[motif] = 0    # initialize the count
    motifAppearances[motif] += 1       # increment the count

现在您已经只处理了整个字符串一次，并找到了其中存在的所有图案。在这种情况下，结果字典看起来像：

motifAppearances = {'AC':2, 'CG':1, 'GT':1, 'TA':1}

为您的案例做类似的事情应该可以将您的运行时间缩短 16384 倍。

【讨论】：

那么，要找到长度为 7 的图案，代码应该是 range(len(s)-6): motif=s[i:i+7] 对吧？试图实现这一点，但这有点简单，让我不知道如何做类似的事情。无论如何谢谢@dg99
@Karyo 是的，听起来不错。我建议先尝试仅使用其中一个长字符串来解决问题。一旦你有了motifAppearances的字典，你就可以遍历每个可能的主题m（你已经有了生成这些的代码）并检查字典以查看m是否存在（if m in motifAppearances ...）并计算根据需要打分。一旦你得到这个工作，然后继续处理多个输入字符串的更大问题。