kmermap_maxbin2/Profiler.cpp at main · catkeyy/kmermap_maxbin2 · 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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#include "Profiler.h"

Profiler::Profiler(int input_kmerlen, const char *seq, kmerMap *input_kmap)
{
	kmerlen = input_kmerlen;
	kmap = input_kmap;
	tempKmer = (char*)malloc(sizeof(char) * (kmerlen + 1));
	memset(tempKmer, '\0', sizeof(char) * (kmerlen + 1));
	profile = (double*)malloc(sizeof(double) * kmap->getEntryNum());//entryNum es el tamano del vector de la secuencia
	memset(profile, '\0', sizeof(double) * kmap->getEntryNum());
	computeProfile(seq);
	total_weight = 1;
}

Profiler::~Profiler()
{
	free(tempKmer);
	free(profile);
}

double* Profiler::getProfile()
{
	return(profile);
}

float Profiler::getPercentN()
{
	return(percent_N);
}

float Profiler::getTargetProfile(char *ch)
{
	return((float)(profile[kmap->getMapping(ch)]));
}

void Profiler::reset()
{
	total_weight = 0;
	memset(profile, '\0', sizeof(double) * kmap->getEntryNum());
}

//parte de get profile
void Profiler::addProfile(Profiler *pf, long double weight)
{
	int i;
	for (i = 0; i < kmap->getEntryNum(); i++)
	{
		profile[i] = profile[i] + (pf->profile[i] * (double)weight);
	}
	total_weight = total_weight + weight;;
}

void Profiler::calcProfile()
{
	int i;
	if (total_weight == 0)
	{
		printf("Total weight is 0! Set all profiles to 0. Nan encountered?\n");
	}
	for (i = 0; i < kmap->getEntryNum(); i++)
	{
		if (total_weight == 0)
		{
			profile[i] = 0;
		}
		else
		{
			profile[i] = profile[i] / (double)total_weight;
		}
	}
	total_weight = 1;
}


//aqui calcula el perfil (el vector para cada read)
/*
profile[i] =
	(número de veces que aparece el k-mer i en la secuencia) / (total de k-mers en la secuencia)
*/
void Profiler::computeProfile(const char *seq)
{
	int len, i, j, N, total;
	N = 0;
	total = 0;
	len = strlen(seq);
	for (i = 0; i < len - kmerlen + 1; i++)
	{
		//nota: N es cuando es una base desconocida (no se pasó bien al archivo)
		if (seq[i] == 'N')
		{
			N++;
		}
		memcpy(tempKmer, seq + i, kmerlen);
		j = kmap->getMapping(tempKmer);
		if (j != -1)
		{
			total++;
			profile[j]++;
			//j = kmap->getReverseMapping(tempKmer);
			//total++;
			//profile[j]++;
		}
	}
	j = kmap->getEntryNum();
	percent_N = (float)N / (float)(len - kmerlen + 1);
	if (percent_N == 1)
	{
		for (i = 0; i < j; i++)
		{
			profile[i] = 0;
		}
	}
	else if (total == 0)
	{
		for (i = 0; i < j; i++)
		{
			profile[i] = 0;
		}
	}
	else
	{
		for (i = 0; i < j; i++)
		{
			profile[i] = profile[i] / total;
		}
	}
}