Дублируй FASTA, сохраняй seq id

Мне нужно отформатировать файлы для инструмента идентификации miRNA (miREAP).

У меня есть файл fasta в следующем формате:

>seqID_1
CCCGGCCGTCGAGGC
>seqID_2
AGGGCACGCCTGCCTGGGCGTCACGC
>seqID_3
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA
>seqID_4
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA
>seqID_5
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA
>seqID_6
AGGGCACGCCTGCCTGGGCGTCACGC

Я хочу посчитать, сколько раз каждая последовательность встречается, и добавить это число в строку seqID. Счетчик для каждой последовательности и исходный идентификатор, относящийся к последовательности, должны появляться в файле только один раз, например:

>seqID_1 1
CCCGGCCGTCGAGGC
>seqID_2 2
AGGGCACGCCTGCCTGGGCGTCACGC
>seqID_3 3
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA

Fastx_collapser делает то, что хотел бы я ( http://hannonlab.cshl.edu/fastx_toolkit/index.html). Однако вместо того, чтобы поддерживать seqID, он возвращает:

>1 1
CCCGGCCGTCGAGGC
>2 2
AGGGCACGCCTGCCTGGGCGTCACGC
>3 3
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA

Это означает, что связь между моей последовательностью, seqID и местоположением картирования генома потеряна. (Каждый seqID соответствует последовательности в моем файле fasta и месте отображения генома в отдельном сгенерированном Bowtie2 файле.sam)

Есть ли простой способ сделать нужную дедупликацию в командной строке?

Спасибо!

1 ответ

Решение

Линеаризовать и отсортировать /uniq -c

 awk '/^>/ {if(N>0) printf("\n"); ++N; printf("%s ",$0);next;} {printf("%s",$0);} END { printf("\n");}'  input.fa  | \
  sort -t ' ' -k2,2 | uniq -f 1 -c |\
  awk '{printf("%s_%s\n%s\n",$2,$1,$3);}'

>seqID_2_2
AGGGCACGCCTGCCTGGGCGTCACGC
>seqID_1_1
CCCGGCCGTCGAGGC
>seqID_3_3
CCGCATCAGGTCTCCAAGGTGAACAGCCTCTGGTCGA
Другие вопросы по тегам