1

我有一个 fastq 文件，其中包含超过 1 亿次读取和 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

GGGGAGGAATGCC

TCAAAAACATAGG

TCAAAAACATGGG

输出文件：

GGGGAGGAATATGAT 0#0错配精确序列

GGGGAGGAATATGAA 1 # 1 不匹配

GGGGAGGAATATGCC 2 #2 不匹配

TCAAAAACATAGG 2 # 2 不匹配

TCAAAAACATGGG 3 # 3 不匹配

4

2 回答 2

2

something like?

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;
}

for you input produces:

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

Honestly, haven't idea how will work with an 10000 chars long genseq...

于 2013-05-02T23:00:35.733 回答
0

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

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

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

于 2013-05-03T14:11:29.610 回答