条形码 DNA 的聚类分析答案

【问题标题】：Cluster analysis of barcoded DNA条形码 DNA 的聚类分析
【发布时间】：2016-02-11 02:24:23
【问题描述】：

我在 PCR 之前使用条形码标记线粒体 DNA 链。条码序列是未知的，但它们有 18 个核苷酸长，并直接进行已知序列（CATCAT 或 TACTAC）。每个 DNA 分子都会获得一个唯一的条形码标识符。分子进行 PCR 后，我需要根据其 18 个核苷酸的条形码对序列进行聚类，然后根据条形码对齐序列。
举一个过于简单的例子，假设我有 2 个分子正在进行 PCR 反应：

     CATCATBARCODE1SEQUENCE1
     TACTACBARCODE2SEQUENCE2

放大后我有：

     CATCATBARCODE1SEQUENCE1
     CATCATBARCODE1SEQUENCE1 
     TACTACBARCODE2SEQUENCE2
     TACTACBARCODE2SEQUENCE2

然后我想搜索位置 6-13 的序列部分，并根据该序列窗口对它们进行聚类，而不更改序列的其余部分，这实际上就像我上面的内容。然后我可以对相邻的序列进行比对。
关于如何在不考虑序列的其余部分的情况下完成序列窗口的这种聚类的任何想法？谢谢。

【问题讨论】：

使用dplyr 你可以做类似df %>% mutate(sub_seq = substr(dna_seq, start, end)) %>% group_by(sub_seq) %>% ...的事情
嗨史蒂夫，我刚开始研究 dplyr。这似乎是一个不错的选择。我会尝试你的建议。如果需要，我什至可以将数据格式化为制表符分隔的列，分隔 CATCAT|TACTAC、条形码和序列，然后使用 dplyr 中的一些函数来组织它们。我会继续调查这个...谢谢你的提示。
如果你能给出更多关于输入和预期输出的细节，我可以给你一个更具体的例子。
谢谢史蒂夫。我有一个包含 fast5 文件的文件夹，可以将其转换为 fasta，然后修剪序列，使其两侧有已知的 CATCAT|TACTAC。然后我将所有序列发送到一个 txt 文件。由新行分隔。我不确定输出格式应该是什么，因为我还没有找到聚类/对齐方法，但我想我可以处理一个大的序列文本文件......
您是否从dplyr 的建议中获得了您所需要的东西，还是需要更具体的东西？

标签： r perl cluster-analysis bioinformatics

【解决方案1】：

过度简化的 R 代码，但似乎可以满足您的要求：

seqs <- c('CATCATBARCODE1SEQUENCE1',
         'CATCATBARCODE1SEQUENCE1',
         'TACTACBARCODE2SEQUENCE2',
         'TACTACBARCODE2SEQUENCE2' )

clusters <- list()

for (seq in seqs) {
  barcode <- substr(seq, 7, 14)
  if (!is.null(clusters[[barcode]])) {
    clusters[[barcode]] <- append(clusters[[barcode]], seq)
  } else {
    clusters[[barcode]] <- c(seq)
  }
}

print(clusters)

打印：

$BARCODE1
[1] "CATCATBARCODE1SEQUENCE1" "CATCATBARCODE1SEQUENCE1"

$BARCODE2
[1] "TACTACBARCODE2SEQUENCE2" "TACTACBARCODE2SEQUENCE2"

【讨论】：

您好文斯，感谢您的建议。我尝试了这个并得到了很多错误，例如“除非在第 10 行的 void 上下文中使用数字 lt (
你用什么运行代码... R 或 perl？是 R 代码。

【解决方案2】：

假设您已经可以获取以 [CATBARCODEX] 开头的序列，我要做的就是在 python 中处理它。如果您的序列开始不同，那么您可能需要搜索 CATCAT 并丢弃那些看起来位于错误位置的序列。如果条形码的数量非常大，可能会出现一些问题，但我认为大约 100,000 个简单的方法应该可以工作。

无论如何，一旦你找到了 CATCAT，我要做的就是建立一个条形码字典并开始过滤。然后，您可以撕下序列的第一部分并使用任何方法进行比对（我有一个条形码项目，并且使用带蝴蝶结的自定义基因组很方便）。

假设你需要找到这个序列而不是仅仅从它开始，在 python 中的解决方案就像

my_dict= {}
for seq in seqs:
     idx = seq.find("CATCAT")
     idx2 = seq.find("TACTAC")
     if idx==-1 and idx2==-1:continue
     # here you need to consider the location of idx and idx2, both may be present, sequence needs to be long enough etc
     barcode = seq[idx+6, idx+6+18]
     # you may want to shorten the barcode or encode it to a string
     if barcode in my_dict: 
          my_dict[barcode]=1 
     else : 
          my_dict[barcode]+=1;
     seq=seq[idx+24:]

除了计数之外，您还可以 1) 将序列附加到每个条形码的 fasta 文件中，或 2) 将条形码作为注释分配给大型 fasta 文件。

无论您可能希望剥离序列以简化下游分析。

【讨论】：