Как я могу устранить дублирующиеся последовательности в файле FASTA - PullRequest
1 голос
/ 22 апреля 2020

Я пытаюсь построить базу данных бактерий жанра, используя все опубликованные последовательности, чтобы вычислить охват моих операций чтения с этой базой данных, используя bowtie2 для картирования, для этого я объединяю все последовательности геномов, которые я скачал из ncbi, в один файл fasta_library (я объединение 74 файлов в файл fasta), проблема в том, что в этом файле fasta (созданной мной библиотеке) у меня много дублированных последовательностей, и это сильно повлияло на покрытие, поэтому я спрашиваю, есть ли способ чтобы устранить дублирование, которое есть в моем Library_File, или если есть какой-либо способ объединить последовательности без дублирования, или также если есть какой-либо другой способ рассчитать охват моих операций чтения по ссылочным последовательностям

Я надеюсь, что я ' достаточно ясно, пожалуйста, скажите мне, если что-то не ясно.

1 Ответ

2 голосов
/ 22 апреля 2020

Если у вас есть контроль над настройками, вы можете установить seqkit и запустить в файле FASTA следующее:

$ seqkit rmdup -s < in.fa > out.fa

Если у вас есть несколько файлов, вы можете объединить их и введите их как стандартный ввод:

$ seqkit rmdup -s < <(cat inA.fa ... inN.fa) > out.fa

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

Чтобы избежать сторонних зависимостей и понять, как удаляются дубликаты, можно использовать awk.

Идея состоит в том, чтобы прочитать все записи FASTA одну за другой в ассоциативный массив (или таблицу ha sh, также называемую «словарь» в Python), только если последовательность отсутствует в массиве.

Например, начиная с однострочного файла FASTA in.fa, который выглядит следующим образом:

>test1
ATAT
>test2
CGCG
>test3
ATAT
>test4
GCCT

Мы можем удалить дубликаты, сохранив первый заголовок Примерно так:

$ awk 'BEGIN {i = 1;} { if ($1 ~ /^>/) { tmp = h[i]; h[i] = $1; } else if (!a[$1]) { s[i] = $1; a[$1] = "1"; i++; } else { h[i] = tmp; } } END { for (j = 1; j < i; j++) { print h[j]; print s[j]; } }' < in.fa > out.fa
$ cat out.fa
>test1
ATAT
>test2
CGCG
>test4
GCCT

Требуется немного знаний о awk, если вам нужны модификации. Этот подход также зависит от того, как структурированы ваши файлы FASTA (записи с последовательностями в одну строку или несколько строк и т. Д. c.), Хотя обычно довольно легко преобразовать файлы FASTA в указанную выше структуру (по одной строке для заголовка и sequence).

Любой табличный подход ha sh также использует довольно много памяти (я думаю, что seqkit, вероятно, делает такой же компромисс для этой конкретной задачи, но я не смотрел на источник) , Это может быть проблемой для очень больших файлов FASTA.

Возможно, лучше использовать seqkit, если у вас есть локальная среда, в которой вы можете установить программное обеспечение. Если у вас настроена блокировка IT, то awk будет работать и для этой задачи, так как она поставляется с большинством Unix-систем из коробки.

...