Извлечение нескольких последовательностей из файла fasta размером 7 ГБ с использованием идентификаторов из отдельного текстового файла

После нескольких часов поиска на этом веб-сайте и пробования множества разных вещей, которые не работали, я решил опубликовать свой вопрос. В настоящее время у меня есть текстовый файл (id.txt), который содержит около 100 строк следующих IDS в этой форме:

5377-P3-D5-MSITS2a_R1reads1_1125821

5377-P3-D5-MSITS2a_R1reads1_1126992

У меня есть фаста на 7 гб с записями в виде

>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT  

>5377-P3-D5-MSITS2a_R1reads1_1129826 M00532:203:000000000-BKM3D:1:1110:14480:9405 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAAACTCTCGAGGTTACAGCCTTGCTGAATTATTAACCCTTGTCGTTCGCGTACTTCTTGTTTCCTTGGTGTGTTCGCCCACCACAAGTAAAAACATAAACCTTTTGTAA

Все идентификаторы из id.text можно найти в seq.fasta. Ожидаемый результат найдет соответствующий идентификационный номер в файле fasta из файла id.text и выдаст:

>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA  

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

В настоящее время я могу извлекать одну последовательность из файла fasta за раз, используя grep в bash, просто копируя и вставляя один идентификатор из файла.

Ex: grep 5377-P3-D5-MSITS2a_R1reads1_1126992 seq.fasta -A 1

Результат:

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

Однако у меня есть несколько текстовых файлов, каждый из которых содержит от 50 до 300 идентификаторов, которые я хотел бы использовать для извлечения последовательностей из файла FASTA, а извлечение последовательностей по отдельности кажется излишне трудоемким. Я хотел бы найти способ найти и вывести последовательности из файла fasta для нескольких идентификаторов, расположенных в отдельном текстовом файле. Я в основном экспериментировал с командами awk и grep в bash, основываясь в основном на других ответах на этом сайте, и почти каждая команда, которую я пытаюсь, не дает результата и сообщения об ошибке.

Примеры, которые я пробовал:

awk -F '>' 'NR==FNR{ids[$0]; next} NF>1{f=($2 in ids)}f' id.txt seq.fasta

awk 'NR==FNR{ids[$0];next} /^>/{f=($1 in ids)} f' id.txt seq.fasta

grep -Fwf id.txt seq.fasta

grep -Ff id.txt seq.fasta

Мне кажется, что я пробовал много вариантов этих двух команд (на основе других предложений о переполнении стека и биостаре), а в bash ничего не происходит, нет результата или нет сообщения об ошибке. Я тоже относительно новичок в кодировании, поэтому я не могу точно определить, что идет не так. Я также открыт для любого Python или другого кода, который также может быть использован. Любая помощь или совет будут оценены. Спасибо!


person Emily W    schedule 13.05.2020    source источник


Ответы (2)


Мне кажется, что grep - лучшая идея. Я думаю, вам может потребоваться удалить символы * из строк поиска, поскольку они не соответствуют тому, что находится в файле. С этим изменением, похоже, он работает, когда я пробую ваш экстракт:

$ cat fasta 
*>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0   
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA* 

*>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT* 

*>5377-P3-D5-MSITS2a_R1reads1_1129826 M00532:203:000000000-BKM3D:1:1110:14480:9405 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0    
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAAACTCTCGAGGTTACAGCCTTGCTGAATTATTAACCCTTGTCGTTCGCGTACTTCTTGTTTCCTTGGTGTGTTCGCCCACCACAAGTAAAAACATAAACCTTTTGTAA*
$ cat ids.txt 
5377-P3-D5-MSITS2a_R1reads1_1125821
5377-P3-D5-MSITS2a_R1reads1_1126992
$ grep -Ff ids.txt fasta 
*>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0   
*>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0 
$
person RGoodman    schedule 13.05.2020
comment
Извините, этих звездочек не было, когда я пробовал разные коды, и что-то странно отформатировалось при публикации. Я ищу не только вывод заголовка, но и вывод последовательности. Однако даже когда я использую команду grep с обоими файлами, я не получаю никаких результатов. Есть ли вероятность, что мой файл fasta слишком большой? Спасибо. - person Emily W; 13.05.2020
comment
Если вы добавите параметр -A1 в команду grep, она должна напечатать совпавшую строку плюс одну строку после нее. Если вы опубликуете свой фактический ввод и команду, возможно, станет яснее, почему вы не получаете никакого вывода. Как видите, когда я пробую, он работает нормально. - person RGoodman; 13.05.2020
comment
Если я создам новый файл с несколькими образцами, перечисленными в моем вопросе, он также отлично подойдет мне. Тем не менее, когда я использую полный файл fasta, вывода вообще нет, поэтому меня это беспокоит. Я могу попробовать разделить свой файл fasta на части, чтобы посмотреть, поможет ли это. - person Emily W; 13.05.2020

Я обнаружил, что снятие флага -F и использование просто -f для входного файла шаблонов и -A 1 для получения последовательности отлично работает, с решением:
grep -A 1 -f ids.txt seq.fasta

Кроме того, если вам не нужен разделитель -- между извлеченными записями, добавьте | grep -v "\-\-", чтобы удалить эти строки (обратная косая черта необходима для экранирования дефисов).

Полный вывод:

% cat ids.txt
5377-P3-D5-MSITS2a_R1reads1_1125821
5377-P3-D5-MSITS2a_R1reads1_1126992

% cat seq.fasta
>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA

>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

>5377-P3-D5-MSITS2a_R1reads1_1129826 M00532:203:000000000-BKM3D:1:1110:14480:9405 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAAACTCTCGAGGTTACAGCCTTGCTGAATTATTAACCCTTGTCGTTCGCGTACTTCTTGTTTCCTTGGTGTGTTCGCCCACCACAAGTAAAAACATAAACCTTTTGTAA

% grep -A 1 -f ids.txt seq.fasta
>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA
--
>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT

% grep -A 1 -f ids.txt seq.fasta | grep -v "\-\-"
>5377-P3-D5-MSITS2a_R1reads1_1125821 M00532:203:000000000-BKM3D:1:1101:10654:16493 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGAACCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACTTCTTGTTTCCTTGGTGGGTTCGCCCACCACTAGGACAAACATAAACCTTTTGTATTGGCA
>5377-P3-D5-MSITS2a_R1reads1_1126992 M00532:203:000000000-BKM3D:1:1104:27124:5463 1:N:0:213 orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0
AAGTCGTAACAAGGTCTCCGTAGGTGAACCTGCGGAGGGATCATTACACAAATATGAAGGCGGGCTGGACCCTCTCGGGGTTACAGCCTTGCTGAATTATTCACCCTTGTCTTTTGCGTACATCTTGTTTCCTTTGTTGTTTCTCCCACCCCTAGGACAAACATAAACCTTTAGTAATTTCAATCAGCGT
person Sam    schedule 28.05.2020