Grep, который допускает несоответствия подмножеству .fastq - PullRequest
Новая
       67

Grep, который допускает несоответствия подмножеству .fastq

0 голосов
/ 13 декабря 2018

Я работаю с Bash на кластере Linux.Я пытаюсь извлечь чтения из файла .fastq, если они содержат совпадение с запрашиваемой последовательностью.Ниже приведен пример файла .fastq, содержащего три чтения.

$ cat example.fastq

@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.2 2/1
CTATANTATTCTATATTTATTCTAGATAAAAGCATTCTATATTTAGCATATGTCTAGCAAAAAAAA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6

Я хотел бы извлечь чтения, содержащие последовательность GAAATAATA.Я могу выполнить это извлечение, используя grep, как показано в следующей команде. 

$ grep -F -B 1 -A 2 "GAAATAATA" example.fastq > MATCH.fastq

$ cat MATCH.fastq

@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6

Однако эта стратегия не допускает никаких несовпадений.Например, чтение, содержащее последовательность GAAAT G ATA, будет проигнорировано.Мне нужно это извлечение, чтобы допустить одно несоответствие в любой позиции в запрашиваемой последовательности.Итак, мой вопрос: как мне этого добиться?Существует ли пакет выравнивания последовательностей с аналогичной функциональностью для grep?Есть ли какие-либо пакеты поднабора fastq, которые выполняют этот тип операции?Одно замечание, что скорость очень важна.Спасибо за ваше руководство.

Ответы [ 3 ]

0 голосов
/ 13 декабря 2018

Вы можете попробовать файл шаблонов - 

$: cat GAAATAATA
.AAATAATA
G.AATAATA
GA.ATAATA
GAA.TAATA
GAAA.AATA
GAAAT.ATA
GAAATA.TA
GAAATAA.A
GAAATAAT.

затем 

grep -B 1 -A 2 -f GAAATAATA example.fastq > MATCH.fastq

, но это, вероятно, немного замедлит процесс, добавив как полный анализ регулярных выражений, так и альтернативный шаблондля каждого возможного отдельного изменения ...

, отвечая на вопрос в комментариях:

Для заданного значения $word, например word=GAAATAATA, 

awk '{
  for ( i=1; i<=length($0); i++ ) {
     split($0,tmp,""); tmp[i]=".";
     for ( n=1; n<=length($0); n++ ) { printf tmp[n]; }
     printf "\n";
  }
}' <<< "$word" > "$word"

Это создаст этот конкретный файл.Надеюсь, это поможет, но помните, что это будет намного медленнее, поскольку вы теперь используете регулярные выражения вместо простого сопоставления простых строк, И вы вводите целую серию альтернативных шаблонов для сопоставления ...

0 голосов
/ 14 декабря 2018

Это должно работать, но idk, если MATCH.fastq в вашем вопросе является ожидаемым результатом или нет, или даже если ваш пример ввода содержит какие-либо случаи, которые требуют рабочего решения, чтобы найти, так что idk, если он действительно работает или нет: 

$ cat tst.awk
BEGIN {
    for (i=1; i<=length(seq); i++) {
        regexp = regexp sep substr(seq,1,i-1) "." substr(seq,i+1)
        sep = "|"
    }
}
{ rec = rec $0 ORS }
!(NR % 4) {
    if (rec ~ regexp) {
        printf "%s", rec
    }
    rec = ""
}

$ awk -v seq='GAAATAATA' -f tst.awk example.fastq
@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
0 голосов
/ 13 декабря 2018

Вот решение, использующее agrep для получения номеров записей совпадений и awk, который печатает эти записи с некоторым контекстом (из-за отсутствия -A и -B в agrep): 

$ agrep -1 -n  "GAAATGATA" file | 
  awk -F: 'NR==FNR{for(i=($1-1);i<=($1+2);i++)a[i];next}FNR in a' - file

Выход: 

@SRR1111111.1 1/1
CTGGANAAGTGAAATAATATAAATTTTTCCACTATTGAATAAAAGCAACTTAAATTTTCTAAGTCG
+
AAAAA#EEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEA<AAEEEEE<6
@SRR1111111.3 3/1
CTATANTATTGAAATAATAATGTAGATAAAACTATTGAATAACAGCAACTTAAATTTTCAATAAGA
+
AAAAA#EE6EEEEEEEEEEEEAAEEAEEEEEEEEEEEE/EAE/EAE/EA/EAEAAAE//EEAEAA6
...