Changeset 3889


Ignore:
Timestamp:
Mar 25, 2020, 11:40:58 AM (2 months ago)
Author:
Peter
Message:

closes #941

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/bam.cc

    r3883 r3889  
    225225  }
    226226
     227
    227228  // just for consistency
    228229  utility::Aligner::Cigar cig;
     
    244245  seq.resize(120, 'A');
    245246  qual.resize(120, 0);
     247  bam.sequence(seq, qual);
     248  if (bam.sequence()!=seq) {
     249    suite.add(false);
     250    suite.err() << "incorrect sequence:" << bam.sequence() << "\n";
     251    suite.err() << "expected: " << seq << "\n";
     252  }
     253
     254
     255  // set arbitrary sequence
     256  seq = "CCCCGGTGCGAAAAAAATATATATATATT";
     257  qual.resize(seq.size(), 60);
    246258  bam.sequence(seq, qual);
    247259  if (bam.sequence()!=seq) {
  • trunk/yat/omic/BamRead.cc

    r3883 r3889  
    347347            static_cast<const void*>(aux()), aux_size());
    348348    // copy sequence
    349     for (size_t i=0; i<seq.size(); ++i)
    350       sequence(i, nt16_table(seq[i]));
    351     core().l_qseq = seq.size();
     349    const size_t N = seq.size();
     350    core().l_qseq = N;
     351    std::string::const_iterator first = seq.begin();
     352    std::string::const_iterator second = seq.begin()+1;
     353    for (size_t i=0; second<seq.end(); ++i) {
     354      sequence(i, nt16_table(*first), nt16_table(*second));
     355      first+=2;
     356      second+=2;
     357    }
     358    if (N % 2) {
     359      assert(seq.end() - first == 1);
     360      sequence(N-1, nt16_table(*first));
     361    }
    352362    // copy quality
    353363    memcpy(static_cast<void*>(bam_get_qual(bam_)),
     
    363373    assert(seq < 16);
    364374    bam1_seq_seti(bam_get_seq(bam_), i, seq);
     375  }
     376
     377
     378  void BamRead::sequence(size_t k, uint8_t first, uint8_t second)
     379  {
     380    assert(bam_);
     381    assert(first < 16);
     382    assert(second < 16);
     383    assert(2*k < sequence_length());
     384
     385    /*
     386    sequence(2*k, c1);
     387    sequence(2*k+1, c2);
     388
     389    is equivalent to
     390    s[(2*k)>>1] =
     391      (s[(2*k)>>1] & 0xf<<(((2*k)&1)<<2)) | (c1)<<((~(2*k)&1)<<2)
     392    s[(2*k+1)>>1] =
     393      (s[(2*k+1)>>1] & 0xf<<(((2*k+1)&1)<<2)) | (c2)<<((~(2*k+1)&1)<<2)
     394
     395    s[k] = (s[k] & 0xf<<(0<<2)) | (c1)<<(1<<2)
     396    s[k] = (s[k] & 0xf<<(1<<2)) | (c2)<<(0<<2)
     397
     398    s[k] = (s[k] & 0xf<<(0x0)) | (c1)<<(0x4)
     399    s[k] = (s[k] & 0xf<<(0x4)) | (c2)<<(0x0)
     400
     401    s[k] = (s[k] & 0xf ) | (c1<<0x4)
     402    s[k] = (s[k] & 0xf0) | c2
     403
     404    s[k] = (c1<<0x4) | c2
     405    */
     406
     407    bam_get_seq(bam_)[k] = (first << 4) | second;
    365408  }
    366409
  • trunk/yat/omic/BamRead.h

    r3883 r3889  
    385385    // ensure capacity of data pointer is (at least) n
    386386    void reserve(uint32_t n);
     387
     388    /**
     389       \param k which byte in underlying sequence container to set
     390       \param x 4-bit integer
     391       \param y 4-bit integer
     392
     393       equivalent to calling
     394       sequence(2*k, x);
     395       sequence(2*k+1, y);
     396    */
     397    void sequence(size_t k, uint8_t x, uint8_t y);
    387398
    388399    bam1_t* bam_;
Note: See TracChangeset for help on using the changeset viewer.