поиск последовательности в геноме с несоответствием

Question

У меня есть файл fastq с более чем 100 миллионов просмотров в нем и последовательности генома 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 + С: С3 @: C): 1 *):?) ####### #########

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

NATCCCCATCCTCTGCTTGCTTTTCGGGATATGTTGTAGGATTCTCAGC TAGGATTTCAAATGGGTCGAGGTGGTCCGTTAGGTATAGGGGCAACAGG

У меня есть все последовательности в файле, теперь я хочу взять каждую строку последовательности и поиска в последовательности генома с разрешением 3 несовпадений и если он находит напечатать последовательности

пример:

генома файл последовательности:

GGGGAGGAATATGATTTACAGTTTATTTTTCAACTGTGCAAAATAACCTTAACTGCAGACGTTATGACATACATACATTCTATGAATTCCACTATTTTGGAGGACTGGAATTTTGGTCTACAACCTCCCCCAGGAGGCACACTAGAAGATACTTATAGGTTTGTAACCCAGGCAATTGCTTGTCAAAAACATACA

Последовательность поиска файла:

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 символьными длинной genseq ...

Answer 2

Я думаю, вы должны рассмотреть возможность использования инструмента выравнивания предназначенный для этих данных по нескольким причинам:

• Эти инструменты также найдут обратные совпадения (хотя, лет u также может реализовать это).
• Aligners будут правильно обрабатывать сообщения с парным концом и несколько совпадений.
• Большинство выравнивателей написаны на языке C и используют структуры данных и алгоритмы, предназначенные для этого объема данных.

По этим причинам и другим, любой сценарий, который вы придумали, скорее всего, будет не таким быстрым и полным, как инструменты, которые уже существуют. Если вы хотите указать количество несоответствий для сохранения, вместо того, чтобы выравнивать все ваши чтения и затем анализировать выходные данные, вы можете использовать Vmatch, если у вас есть к нему доступ (этот инструмент очень быстр и хорош для многих совпадающих задач).