如何在 C++ 中滑动读取多行 fasta 文件的窗口

Question

我有一个 C++ 代码，它读取一个 fasta 文件（单行读取）并在读取时滑动一个 l 窗口以生成 kmers。 像这样：

 >NC_000913.3-2320800
 GACGAAGCTGCCGACCAGCGTTTTATGTCAGTGCGGACATCAGCACTGCGTGAACAATTTGCCTGGCTGCGCGAGAACGGTTATCAACCGGTCAGTATTG

现在我有另一个 fasta 文件，它读取了不止一行这样的内容：

  >LR2957 4054ae1b-13f5-89b1-e7b5-d64cabb8badb NC_000913.3,+strand,4451319-4510860 length=55312 error-free_length=59451 read_identity=74.36%
  TGTAGTCCGTCAAGTTACGTTATTGCTACGTCTATCAGGGAAGTCAACCTGCCTGCAATA
  TGGTAGATAAATCCTATTATGCCGCGAGACAACCCTTGGCTTCCTACACGCGCAGTGGAG

这个函数是：

 void Sliding_window_l (const char *ptr, size_t length) 
 {   
    size_t p=0;
   /*find start of a read*/
   for(; ptr[p]!='>' && p<length; p++) {/*noop*/}
   kmer_t kmer = 0;

   while(p<length) {
     //printf("hello second %c\n", ptr[p]);  
     assert(ptr[p]=='>'); /*this will be true*/

     /*skip till newline*/
     for(; p<length && ptr[p]!='\n'; p++) {/*noop*/ }
     p++; /*skip the newline*/

    if(p+LMER_LENGTH > length) break; /*too short a read*/
    kmer = 0;
    int i;
    for(i=0; ptr[p]!='\n' && i<LMER_LENGTH-1; i++) {
       kmer = lmer_shift(kmer, char_to_el(ptr[p++]));
      //kmer = kmer_cons(kmer, i, char_to_el(ptr[p++]));
    }

  while(p<length && ptr[p]!='\n') {
      kmer = lmer_shift(kmer, char_to_el(ptr[p++]));
      lmer_frequency[kmer]++;

  }
  p++; /*skip the newline*/
 }

}

该函数对于单行读取运行良好，但对于多行读取，我收到如下错误：

 void Sliding_window_l(const char*, size_t): Assertion `ptr[p]=='>'' failed.

任何人都可以指导我如何修改该函数，以便它也可以用于多行读取。非常感谢。

Answer 1

看看你的代码：

    while(p < length) {
        // The first time the assertion is true
        assert(ptr[p]=='>');

        for(; p<length && ptr[p]!='\n'; p++) {/*noop*/ }
        p++;
        // Now p points at the beginning of the first line of payload

        // more code...

        while(p < length && ptr[p] != '\n') {
            //...
        }
        // You've finished the first line
        p++; /*skip the newline*/
        // Now p may point to the beginning of another line of the payload
        // The next iteration would start with p pointing NOT to the '>'
    }

你的意思是断言应该在循环之前完成吗？

无论如何，如果我正确理解了您的代码，您将丢失与上一行结尾和下一行开头重叠的 kmers。

指南：

1. 简化您的生活。逐行读取，将有效负载连接到一大串核苷酸。每当您找到输入的结尾（或者您找到下一个以 '>' 开头的评论）时，您就可以开始处理连接的读取。
2. 避免const char*，更喜欢std::string：您使用的是C++，不是吗？
3. 重新格式化您的代码，使其更具可读性：糟糕的格式是万恶之源。