thermonucleotideBLAST/sequence_data_annot.cpp at master · jgans/thermonucleotideBLAST · GitHub

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#include "sequence_data.h"
#include "throw.h"

// Needed for sort()
#include <algorithm>

using namespace std;

void sequence_data::load_gbk(const std::string &m_filename)
{
	format = GBK;

	// Clear any existing annotation data
	mol.clear();
	size_t pos = 0;

	gzFile fin = gzopen(m_filename.c_str(), "r");

	// Increase the size of the internal zlib buffer used for decompression
	gzbuffer(fin, 32768);

	if(fin == NULL){
		THROW(__FILE__ ":sequence_data::load_gbk: Unable to open input file");
	}

	while(true){

		mol.push_back( DNAMol() );

		if(mol.back().load(fin, DNAMol::GBK, pos) == false){
			break;
		}

		if(mol.back().empty() == true){
			mol.pop_back();
		}
	}

	gzclose(fin);

	list<DNAMol>::const_iterator iter = mol.begin();
	const size_t num_seq = size();

	seq_length.resize(num_seq);

	for(size_t i = 0;i < num_seq;i++, iter++){

		// Don't use readdb_get_sequence_length -- it's too slow on large databases
		const unsigned int seq_len = iter->num_bases();

		seq_length[i] = make_pair(seq_len, i);
	}
}

void sequence_data::load_embl(const std::string &m_filename)
{
	format = EMBL;

	// Clear any existing annotation data
	mol.clear();
	size_t pos = 0;

	gzFile fin = gzopen(m_filename.c_str(), "r");

	// Increase the size of the internal zlib buffer used for decompression
	gzbuffer(fin, 32768);

	if(fin == NULL){
		THROW(__FILE__ ":sequence_data::load_embl: Unable to open input file");
	}

	while(true){

		mol.push_back( DNAMol() );

		if(mol.back().load(fin, DNAMol::EMBL, pos) == false){
			break;
		}
	}

	gzclose(fin);

	if(mol.back().empty() == true){
		mol.pop_back();
	}

	list<DNAMol>::const_iterator iter = mol.begin();
	const size_t num_seq = size();

	seq_length.resize(num_seq);

	for(size_t i = 0;i < num_seq;i++, iter++){

		// Don't use readdb_get_sequence_length -- it's too slow on large databases
		const unsigned int seq_len = iter->num_bases();

		seq_length[i] = make_pair(seq_len, i);
	}
}