Дублируй 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