错配基因组中的搜索序列

Question

我有一个fastq文件，里面有超过1亿的读取量和10000的基因组序列。

我想从fastq文件中取出序列，在基因组序列中搜索，允许3次错配。

我用这种方式使用awk获得了fastq文件中的序列：

1.财务报表(几行)

@DH1 269 1:269:C1UKCACXX:1:1101:1207:2171 1:N:0:TTAGGC + 1=ADBDDHD;F>GF@FFEFGGGIAEEI?D9DDHHIGAAF:BG39?BB @DH1 269 1:269:C1UKCACXX:1:1101:1095:2217 1:N:0:TTAGGC + ??AABDD4C:DDDI+C:C3@:C):1?*):?)?################ $ awk 'NR%4==2‘1.fq

我把所有的序列都归档了，现在我想取每一行序列，在基因组序列中搜索，允许3次错配，如果它找到了，就打印这些序列。

示例：

基因组序列文件：

搜索序列文件：

GGGGAGGAATATGAT GGGGAGGAATATGAA GGGGAGGAATATGCC TCAAAAACATAGG TCAAAAACATGGG

输出文件：

GGGGAGGAATATGAT 0#0失配精确序列 GGGGAGGAATATGAA 1#1错配 GGGGAGGAATATGCC 2#2错配 TCAAAAACATAGG 2#2错配 TCAAAAACATGGG 3#3错配

Answer 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字符的长条形.

Answer 2

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

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

由于这些原因和其他原因，您想出的任何脚本都不可能像已经存在的工具那样快速完整。如果希望指定要保留的不匹配数，而不是对所有读取进行对齐，然后解析输出，则可以使用Vmatch (此工具非常快速，适合于许多匹配任务)。