File: RefGenome.cpp

package info (click to toggle)
libseqlib 1.1.2+dfsg-1~bpo9+1
  • links: PTS, VCS
  • area: main
  • in suites: stretch-backports
  • size: 1,460 kB
  • sloc: cpp: 7,176; sh: 805; makefile: 60
file content (61 lines) | stat: -rw-r--r-- 1,440 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include "SeqLib/RefGenome.h"

#include <stdexcept>
#include "SeqLib/SeqLibUtils.h"

namespace SeqLib {

  bool RefGenome::LoadIndex(const std::string& file) {

    // clear the old one
    if (index)  
      fai_destroy(index);
    
    index = NULL;
    
    // check that its readable
    if (!read_access_test(file)) {
      return false;
      //throw std::invalid_argument("RefGenome: file not found - " + file);
    }
    
    // load it in
    index = fai_load(file.c_str());

    if (!index)
      return false;

    return true;

  }
  
  std::string RefGenome::QueryRegion(const std::string& chr_name, int32_t p1, int32_t p2) const {
    
    // check that we have a loaded index
    if (!index) 
      throw std::invalid_argument("RefGenome::queryRegion index not loaded");

    // check input is OK
    if (p1 > p2)
      throw std::invalid_argument("RefGenome::queryRegion p1 must be <= p2");
    if (p1 < 0)
      throw std::invalid_argument("RefGenome::queryRegion p1 must be >= 0");

    int len;
    char * f = faidx_fetch_seq(index, const_cast<char*>(chr_name.c_str()), p1, p2, &len);

    if (!f)
      throw std::invalid_argument("RefGenome::queryRegion - Could not find valid sequence");

    std::string out(f);

    free(f);

    if (out.empty())
      throw std::invalid_argument("RefGenome::queryRegion - Returning empty query on " + chr_name + ":" + tostring(p1) + "-" + tostring(p2));

    return (out);

  }

}