Обход файла Fatsa в C / Stack Overflow
Я ищу написать программу на C / C ++ для обхода файла Fasta, отформатированного как:
>ID and header information
SEQUENCE1
>ID and header information
SEQUENCE2
и так далее
чтобы найти все уникальные последовательности (проверьте, является ли подмножество любой другой последовательности)
и записать уникальные последовательности (и все заголовки) в выходной файл.
Мой подход был:
- Скопируйте все последовательности в массив / список в начале (более эффективный способ сделать это?)
- Возьмите заголовок, добавьте его в выходной файл, сравните последовательность для этого заголовка со всем в списке / массиве. Если уникален, напишите это под заголовком, если не удалите это.
Тем не менее, я немного не уверен, как правильно подходить к чтению строк. Мне нужно прочитать верхнюю строку для заголовка, а затем «вернуть?» на следующую строку, чтобы прочитать последовательность. Иногда последовательность занимает более двух строк, поэтому я бы использовал > (из примера выше) в качестве разделителя? Если я использую C ++, я думаю, я бы использовал iostreams для достижения этой цели?
Если бы кто-нибудь мог подтолкнуть меня в правильном направлении относительно того, как я хотел бы прочитать информацию, которой я должен манипулировать / как провести сравнение, это было бы очень ценно.
Решение
Во-первых, вместо того, чтобы писать свою собственную программу чтения FASTA, вы, вероятно, захотите использовать то, что существует, например, см.: http://lh3lh3.users.sourceforge.net/parsefastq.shtml
Внутренне вы будете иметь последовательности без перевода строки, и это, вероятно, полезно. Я думаю, что самый простой подход с высокого уровня
- зациклите на fasta и запишите последовательности в файл
- сортировать этот файл
- с отсортированным файлом становится легче выбирать подпоследовательности, поэтому напишите программу для поиска «уникальных идентификаторов»
- Используя уникальные идентификаторы, вернитесь к исходному файлу и получите любую дополнительную информацию, которая вам нужна.
Другие решения
Ваш подход пригоден для использования. Ниже приведена реализация этого.
Тем не менее, я немного не уверен, как подходить к чтению строк
в должном порядке. … Иногда последовательность занимает более двух строк, поэтому я бы использовал>(из примера выше) в качестве разделителя?
Вот так; кроме того, есть только EOF, который должен быть проверен.
Я написал функцию getd() для этого, который читает однострочное описание или объединенные строки данных последовательности и возвращает указатель на строку, которую он выделил.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
char *getd()
{
char *s = NULL, *eol;
size_t i = 0, size = 1; // 1 for null character
int c;
#define MAXLINE 80 // recommended max. line length; longer lines are okay
do // read single-line description or concatenated lines of sequence data
{
do // read a line (until '\n')
{
s = realloc(s, size += MAXLINE+1); // +1 for newline character
if (!s) puts("out of memory"), exit(1);
if (!fgets(s+i, size, stdin)) break;
eol = strchr(s+i, '\n');
i += MAXLINE+1;
} while (!eol);
if (!i) { free(s); return NULL; } // nothing read
if (*s == '>') return s; // single-line description
i = eol-s;
ungetc(c = getchar(), stdin); // peek at next character
} while (c != '>' && c != EOF);
return s;
}
int main()
{
char *s;
struct seq { char *head, *data; } *seq = NULL;
int n = 0, i, j;
while (s = getd())
if (*s == '>')
{ // new sequence: make room, store header
seq = realloc(seq, ++n * sizeof *seq);
if (!seq) puts("out of memory"), exit(1);
seq[n-1] = (struct seq){ s, "" };
}
else
if (n) // store sequence data if at least one header present
seq[n-1].data = s;
for (i = 0; i < n; ++i)
{
const int max = 70; // reformat output data to that line length max.
printf("%s", seq[i].head);
for (s = seq[i].data, j = 0; j < n; ++j)
if (j != i) // compare sequence to the others, delete if not unique
if (strstr(seq[j].data, s)) { s = seq[i].data = ""; break; }
for (; strlen(s) > max && s[max] != '\n'; s += max)
printf("%.*s\n", max, s);
printf("%s", s);
}
}