SHOGUN
v1.1.0
|
00001 00002 #include <stdio.h> 00003 #include <string.h> 00004 00005 #include <shogun/mathematics/Math.h> 00006 #include <shogun/lib/config.h> 00007 #include <shogun/io/SGIO.h> 00008 #include <shogun/structure/IntronList.h> 00009 00010 using namespace shogun; 00011 00012 CIntronList::CIntronList() 00013 :CSGObject() 00014 { 00015 m_length = 0; 00016 m_all_pos = NULL; 00017 m_intron_list = NULL; 00018 m_quality_list = NULL; 00019 } 00020 CIntronList::~CIntronList() 00021 { 00022 for (int i=0; i<m_length; i++) 00023 { 00024 SG_FREE(m_intron_list[i]); 00025 SG_FREE(m_quality_list[i]); 00026 } 00027 SG_FREE(m_intron_list); 00028 SG_FREE(m_quality_list); 00029 SG_FREE(m_all_pos); 00030 } 00031 void CIntronList::init_list(int32_t* all_pos, int32_t len) 00032 { 00033 m_length = len; 00034 m_all_pos = SG_MALLOC(int32_t, len); 00035 memcpy(m_all_pos, all_pos, len*sizeof(int32_t)); 00036 m_intron_list = SG_MALLOC(int32_t*, len); 00037 m_quality_list = SG_MALLOC(int32_t*, len); 00038 if (m_intron_list==NULL||m_quality_list==NULL) 00039 SG_ERROR("IntronList: Out of mem 1"); 00040 //initialize all elements with an array of length one 00041 int32_t* one; 00042 for (int i=0;i<m_length;i++) 00043 { 00044 one = SG_MALLOC(int32_t, 1);//use malloc here because mem can be increased efficiently with realloc later 00045 m_intron_list[i] = one; 00046 m_intron_list[i][0] = 1; 00047 one = SG_MALLOC(int32_t, 1); 00048 m_quality_list[i] = one; 00049 m_quality_list[i][0] = 1; 00050 } 00051 } 00052 void CIntronList::read_introns(int32_t* start_pos, int32_t* end_pos, int32_t* quality, int32_t len) 00053 { 00054 int k=0; 00055 for(int i=0;i<m_length;i++)//iterate over candidate positions 00056 { 00057 while (k<len) 00058 { 00059 //SG_PRINT("i:%i, m_all_pos[i]:%i, k:%i, end_pos[k]: %i\n", i, m_all_pos[i], k, end_pos[k]) ; 00060 if (k>0) 00061 if (end_pos[k]<end_pos[k-1]) 00062 SG_ERROR("end pos array is not sorted: end_pos[k-1]:%i end_pos[k]:%i\n", end_pos[k-1], end_pos[k]); 00063 if (end_pos[k]>=m_all_pos[i]) 00064 break; 00065 else 00066 k++; 00067 00068 } 00069 while (k<len && end_pos[k]==m_all_pos[i]) 00070 { 00071 //SG_PRINT("\tk:%i, end_pos[k]: %i, start_pos[k]:%i\n", k, end_pos[k], start_pos[k]) ; 00072 ASSERT(start_pos[k]<end_pos[k]); 00073 ASSERT(end_pos[k]<=m_all_pos[m_length-1]); 00074 // intron list 00075 //------------ 00076 int32_t from_list_len = m_intron_list[i][0]; 00077 int32_t* new_list = SG_REALLOC(int32_t, m_intron_list[i], (from_list_len+1)); 00078 if (new_list == NULL) 00079 SG_ERROR("IntronList: Out of mem 4"); 00080 new_list[from_list_len]= start_pos[k]; 00081 new_list[0]++; 00082 m_intron_list[i] = new_list; 00083 // quality list 00084 //-------------- 00085 int32_t q_list_len = m_quality_list[i][0]; 00086 //SG_PRINT("\t q_list_len:%i, from_list_len:%i \n",q_list_len, from_list_len); 00087 ASSERT(q_list_len==from_list_len); 00088 new_list = SG_REALLOC(int32_t, m_quality_list[i], (q_list_len+1)); 00089 if (new_list == NULL) 00090 SG_ERROR("IntronList: Out of mem 5"); 00091 new_list[q_list_len]= quality[k]; 00092 new_list[0]++; 00093 m_quality_list[i] = new_list; 00094 00095 k++; 00096 } 00097 } 00098 } 00103 void CIntronList::get_intron_support(int32_t* values, int32_t from_pos, int32_t to_pos) 00104 { 00105 if (from_pos>=m_length) 00106 SG_ERROR("from_pos (%i) is not < m_length (%i)\n",to_pos, m_length); 00107 if (to_pos>=m_length) 00108 SG_ERROR("to_pos (%i) is not < m_length (%i)\n",to_pos, m_length); 00109 int32_t* from_list = m_intron_list[to_pos]; 00110 int32_t* q_list = m_quality_list[to_pos]; 00111 00112 //SG_PRINT("from_list[0]: %i\n", from_list[0]); 00113 00114 int32_t coverage = 0; 00115 int32_t quality = 0; 00116 for (int i=1;i<from_list[0]; i++) 00117 { 00118 //SG_PRINT("from_list[%i]: %i, m_all_pos[from_pos]:%i\n", i, from_list[i], m_all_pos[from_pos]); 00119 if (from_list[i]==m_all_pos[from_pos]) 00120 { 00121 //SG_PRINT("found intron: %i->%i\n", from_pos, to_pos ); 00122 coverage = coverage+1; 00123 quality = CMath::max(quality, q_list[i]); 00124 } 00125 } 00126 values[0] = coverage; 00127 values[1] = quality; 00128 }