在基因组中搜索具有错配的序列答案

【问题标题】：search sequence in genome with mismatches在基因组中搜索具有错配的序列
【发布时间】：2023-03-20 04:33:01
【问题描述】：

我有一个 fastq 文件，里面有超过 1 亿个 reads 和一个长度为 10000 的基因组序列

我想从 fastq 文件中取出序列并在允许 3 个不匹配的基因组序列中搜索

我尝试使用 awk 以这种方式从 fastq 文件中获取序列：

1.fq（几行）

@DH1DQQN1:269:C1UKCACXX:1:1101:1207:2171 1:N:0:TTAGGC NATCCCCATCCTCTGCTTGCTTTTCGGGATATGTTGTAGGATTCTCAGC

+

1=ADBDDHD;F>GF@FFEFGGGIAEEI?D9DDHHIGAAF:BG39?BB

@DH1DQQN1:269:C1UKCACXX:1:1101:1095:2217 1:N:0:TTAGGC TAGGATTTCAAATGGGTCGAGGTGGTCCGTTAGGTATAGGGGCAACAGG

+

??AABDD4C:DDDI+C:C3@:C):1?*):?)?###############

$ awk 'NR%4==2' 1.fq

NATCCCCATCCTCTGCTTGCTTTTCGGGATATGTTGTAGGATTCTCAGC TAGGATTTCAAATGGGTCGAGGTGGTCCGTTAGGTATAGGGGCAACAGG

我有文件中的所有序列，现在我想获取每一行序列并在允许 3 个不匹配的情况下在基因组序列中搜索，如果找到则打印序列

示例：

基因组序列文件：

GGGGAGGAATATGATTTACAGTTTATTTTTCAACTGTGCAAAATAACCTTAACTGCAGACGTTATGACATACATACATTCTATGAATTCCACTATTTTGGAGGACTGGAATTTTGGTCTACAACCTCCCCCAGGAGGCACACTAGAAGATACTTATAGGTTGTAACCCAGGCAATTGCTTGTCAAAAACATACA

搜索序列文件：

GGGGAGGAATATGAT

GGGGAGGAATATGAA

GGGGAGGAATATGCC

TCAAAAACATAGG

TCAAAAAACATGGG

输出文件：

GGGGAGGAATATGAT 0#0错配精确序列

GGGGAGGAATATGAA 1 #1 不匹配

GGGGAGGAATATGCC 2 #2 不匹配

TCAAAAACATAGG 2 #2 不匹配

TCAAAAACATGGG 3 #3 不匹配

【问题讨论】：

多少个搜索序列（您展示的示例长度是否具有代表性？）
这只是一个示例来显示它映射序列中的任何位置（允许不匹配）。在我的实际文件中有大约序列长度（从 25 到 100）@ysth
我的搜索序列文件包含超过 1 亿行，这些行确实使用 awk 从 fastq 文件中提取（如上所示）@ysth
只是为了确保我理解...您想在长度约为 10000 的单个基因组序列中搜索 100000000 个长度为 25-100 的序列？
是的，这就是我要找的东西.. 将每一行逐一提取并在基因组序列中搜索不匹配的@ysth

标签： python perl awk biopython bioperl

【解决方案1】：

类似的东西？

use 5.012;
use strict;
use warnings;
use String::Approx qw(aslice);
use File::Slurp;
use Data::Dumper;

my $genseq = "gseq.txt"; #the long sequence

$_ = read_file($genseq);

#read small patterns from stdin
while(my $patt = <>) {
    chomp $patt;
    my $len = length($patt);
    my($index, $size, $distance) = aslice($patt, ["3D0S3", "minimal_distance"]);
    say "$patt matched approx. at $index with mismatch $distance" if $distance <= 3;
}

为您输入产生：

GGGGAGGAATATGAT matched approx. at 0 with mismatch 0
GGGGAGGAATATGAA matched approx. at 0 with mismatch 1
GGGGAGGAATATGCC matched approx. at 0 with mismatch 2
TCAAAAACATAGG matched approx. at 179 with mismatch 2
TCAAAAACATGGG matched approx. at 179 with mismatch 3

老实说，不知道如何使用 10000 个字符长的 genseq...

【讨论】：

谢谢@jm666 我现在正在使用 BOWTIE 来映射序列。让您知道它是否有效

【解决方案2】：

我认为您应该考虑使用为这些数据设计的对齐工具，原因如下：

这些工具还可以找到反向互补匹配（不过，您也可以实现这一点）。
对齐器将正确处理双末端读取和多个匹配。
大多数对齐器都是用 C 语言编写的，并使用为这种数据量设计的数据结构和算法。

出于这些原因以及其他原因，您想出的任何脚本都可能不会像现有工具那样快速和完整。如果您想指定要保留的不匹配数，而不是对齐所有读取然后解析输出，您可以使用 Vmatch（如果您可以访问它）（此工具非常快速且适用于许多匹配任务）。

【讨论】：