Grep, который допускает несоответствия подмножеству.fastq
Я работаю с 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
Однако эта стратегия не терпит никаких несоответствий. Например, чтение, содержащее последовательность GAAATG ATA, будет проигнорировано. Мне нужно это извлечение, чтобы допустить одно несоответствие в любой позиции в запрашиваемой последовательности. Итак, мой вопрос: как мне этого добиться? Существует ли пакет выравнивания последовательностей с аналогичной функциональностью для grep? Есть ли какие-либо пакеты поднабора fastq, которые выполняют этот тип операции? Одно замечание, что скорость очень важна. Спасибо за ваше руководство.
3 ответа
Вот решение с использованием 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
Это должно работать, но IDK, если MATCH.fastq
в вашем вопросе ожидаемый результат или нет, или даже если ваш пример ввода содержит какие-либо случаи, которые нуждаются в рабочем решении, чтобы найти, если он на самом деле работает или нет:
$ 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
Вы можете попробовать файл шаблонов -
$: 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"
Это создаст этот конкретный файл. Надеюсь, это поможет, но помните, что это будет намного медленнее, так как вы теперь используете регулярные выражения вместо простого сопоставления простых строк, и вы вводите целый ряд альтернативных шаблонов для сопоставления...