diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 043afdb..d365123 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -47,7 +47,7 @@ jobs: if [ `uname` == Darwin ]; then meson setup --buildtype=release --unity on builddir else - env CC=clang CXX=clang meson setup --buildtype=release --unity on builddir + env CC=clang CXX=clang++ meson setup --buildtype=release --unity on builddir fi - name: meson_compile diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..36726cb --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,38 @@ +cmake_minimum_required(VERSION 3.0) + +project(PretextGraph LANGUAGES CXX C) + +set(CMAKE_CXX_STANDARD 11) + +set(version "0.0.6") + +if (NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Debug) + message(STATUS "No build type selected, default to Debug") +endif() +message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") + +link_directories( + ${CMAKE_SOURCE_DIR}/libdeflate +) +add_executable(PretextGraph PretextGraph.cpp) + +target_include_directories( + PretextGraph + PRIVATE + wrapper + libdeflate +) + +target_link_libraries( + PretextGraph + PRIVATE + pthread + deflate +) + +if (CMAKE_BUILD_TYPE STREQUAL "Debug") + target_compile_definitions(PretextGraph PRIVATE DEBUG) +endif() + +target_compile_definitions(PretextGraph PRIVATE PV="${version}") diff --git a/Header.h b/Header.h index a437880..c0db3bb 100644 --- a/Header.h +++ b/Header.h @@ -31,6 +31,15 @@ SOFTWARE. #include #endif + +// cpp includes +#include +#include +#include +#include +#include // used to store the data_type_dic + + #include #include #include @@ -112,22 +121,23 @@ typedef size_t memptr; #define EndMain return(0) #ifndef _WIN32 -#include -#include -#include + #include + #include + #include #else -#define __atomic_fetch_add(x, y, z) _InterlockedExchangeAdd(x, y) -#define __atomic_add_fetch(x, y, z) (y + _InterlockedExchangeAdd(x, y)) -#define __sync_fetch_and_add(x, y) _InterlockedExchangeAdd(x, y) -#define __sync_fetch_and_sub(x, y) _InterlockedExchangeAdd(x, -y) -#define __atomic_store(x, y, z) _InterlockedCompareExchange(x, *y, *x) + #define __atomic_fetch_add(x, y, z) _InterlockedExchangeAdd(x, y) + #define __atomic_add_fetch(x, y, z) (y + _InterlockedExchangeAdd(x, y)) + #define __sync_fetch_and_add(x, y) _InterlockedExchangeAdd(x, y) + #define __sync_fetch_and_sub(x, y) _InterlockedExchangeAdd(x, -y) + #define __atomic_store(x, y, z) _InterlockedCompareExchange(x, *y, *x) #endif #ifndef _WIN32 -#define ThreadFence __asm__ volatile("" ::: "memory") + #define ThreadFence __asm__ volatile("" ::: "memory") #else -#define ThreadFence _mm_mfence() + #define ThreadFence _mm_mfence() #endif + #define FenceIn(x) ThreadFence; \ x; \ ThreadFence @@ -268,10 +278,10 @@ thread_context }; #ifdef DEBUG -#include -#define Assert(x) assert(x) + #include + #define Assert(x) assert(x) #else -#define Assert(x) + #define Assert(x) #endif #define KiloByte(x) 1024*x @@ -763,9 +773,8 @@ ThreadPoolWait(thread_pool *threadPool) { LockMutex(threadPool->threadCountLock); - while (threadPool->jobQueue.len || threadPool->numThreadsWorking) - { - WaitOnCond(threadPool->threadsAllIdle, threadPool->threadCountLock); + while (threadPool->jobQueue.len || threadPool->numThreadsWorking) { + WaitOnCond(threadPool->threadsAllIdle, threadPool->threadCountLock); } UnlockMutex(threadPool->threadCountLock); @@ -934,9 +943,9 @@ StringToInt_Check(u08 *string, u32 *result, u08 stringTerminator = '\0') while(--strLen > 0 && goodResult) { - *result += (u32)(*--string - '0') * pow; - goodResult = (*string >= '0' && *string <= '9'); - pow *= 10; + *result += (u32)(*--string - '0') * pow; + goodResult = (*string >= '0' && *string <= '9'); + pow *= 10; } *result += (u32)(*--string - '0') * pow; @@ -1037,18 +1046,21 @@ PushStringIntoIntArray(u32 *intArray, u32 arrayLength, u08 *string, u08 stringTe u08 *stringToInt = (u08 *)intArray; u32 stringLength = 0; + // Copy the string into the integer array until the string terminator is reached or the integer array is full while (*string != stringTerminator && stringLength < (arrayLength << 2)) { *(stringToInt++) = *(string++); ++stringLength; } - + + // Pad the integer array with zeros until it is 4 byte aligned while (stringLength & 3) { *(stringToInt++) = 0; ++stringLength; } + // make sure the integer array is zeroed out for ( u32 index = (stringLength >> 2); index < arrayLength; ++index ) @@ -1056,7 +1068,7 @@ PushStringIntoIntArray(u32 *intArray, u32 arrayLength, u08 *string, u08 stringTe intArray[index] = 0; } - return(string); + return(string); // return the line buffer pointer } global_function diff --git a/PretextGraph.cpp b/PretextGraph.cpp index 933065d..9f0883c 100644 --- a/PretextGraph.cpp +++ b/PretextGraph.cpp @@ -98,6 +98,8 @@ Third-party software and resources used in this project, along with their respec WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. )thirdparty"; + + global_variable u08 Status_Marco_Expression_Sponge = 0; @@ -138,6 +140,7 @@ global_variable thread_pool * Thread_Pool; +/* #pragma clang diagnostic push #pragma GCC diagnostic ignored "-Wreserved-id-macro" #pragma GCC diagnostic ignored "-Wsign-conversion" @@ -151,9 +154,60 @@ Thread_Pool; #define STB_SPRINTF_IMPLEMENTATION #include "stb_sprintf.h" #pragma clang diagnostic pop +*/ + +#ifdef __clang__ // Clang pragma + #pragma clang diagnostic push + #pragma clang diagnostic ignored "-Wreserved-id-macro" + #pragma clang diagnostic ignored "-Wsign-conversion" + #pragma clang diagnostic ignored "-Wcast-align" + #pragma clang diagnostic ignored "-Wextra-semi-stmt" + #pragma clang diagnostic ignored "-Wunused-parameter" + #pragma clang diagnostic ignored "-Wconditional-uninitialized" + #pragma clang diagnostic ignored "-Wdouble-promotion" + #pragma clang diagnostic ignored "-Wpadded" + #pragma clang diagnostic ignored "-Wimplicit-fallthrough" +#elif defined(__GNUC__) // GCC pragma + #pragma GCC diagnostic push + // #pragma GCC diagnostic ignored "-Wreserved-id-macro" + #pragma GCC diagnostic ignored "-Wsign-conversion" + #pragma GCC diagnostic ignored "-Wcast-align" + #pragma GCC diagnostic ignored "-Wunused-parameter" + #pragma GCC diagnostic ignored "-Wdouble-promotion" + #pragma GCC diagnostic ignored "-Wpadded" + #pragma GCC diagnostic ignored "-Wimplicit-fallthrough" +#endif + +// Include the header +#define STB_SPRINTF_IMPLEMENTATION +#include "stb_sprintf.h" + +#ifdef __clang__ + #pragma clang diagnostic pop +#elif defined(__GNUC__) + #pragma GCC diagnostic pop +#endif #include "LineBufferQueue.cpp" + +u32 NUM_THREADS; +void set_num_threads() { // function to set the number of threads + // something interesting happened + // if set the num_threads as 1, the thread will be blocked while reading the file into line buffer + // TODO solve the thread block problem + printf("\n\n"); + #ifdef DEBUG + NUM_THREADS = 4; // after testing, if using 4 threads, the thread is blocked, haven't found the reason yet + PrintStatus("Debug mode, number of thread: %d\n", NUM_THREADS); + #else + NUM_THREADS = 4; // define the thread used in RELEASE mode + PrintStatus("Release mode, number of thread: %d\n", NUM_THREADS); + #endif // DEBUG + return ; +} + + struct contig { @@ -279,10 +333,41 @@ graph volatile s32 *values; }; +// used to store the graph values in f32 format +// actually, don't need the graph_f struct, just use the f32 array to store the values +// but leave this here will help the model to be easier to understand and +// easier to further extend the code +struct +graph_f +{ + f32 *values; +}; + +// the gap extension is different with the coverage and repeat density extensions +// thus we have to add this flag to distinguish the gap extension with the other extensions +// if there is "gap" in the extension name, then the flag will be set as true +// and during the ProcessLine function, the gap value will not be weighted according to the length of the bin and pixel + +// add a dictionary to store the data type +unsigned int data_type(0); +std::unordered_map data_type_dic{ // use this data_type + {"default", 0 }, + {"repeat_density", 1}, // as this is counted in every single bin, so we need to normalise this within the bin + {"gap", 2}, // +}; + global_variable graph * Graph; +global_variable +graph_f * +Graph_tmp; + +// add the mutex here to protect the graph->values while updating the values in multi-thread mode +// before updating the graph->value calculation into f32, Ed used the atomic operation to update the values +std::mutex mtx_global; // define the glabal mutex to protext graph_f->values while updating the values + global_variable volatile u08 Data_Added = 0; @@ -298,19 +383,24 @@ void ProcessLine(void *in) { u32 sizeGraphArray = Map_Properties->textureResolution * Map_Properties->numberOfTextures1D; - u64 binSize_genome = (u64)((f64)Map_Properties->totalGenomeLength / (f64)sizeGraphArray); + u64 bp_per_pixel = (u64)((f64)Map_Properties->totalGenomeLength / (f64)sizeGraphArray); // number of bps per pixel - line_buffer *buffer = (line_buffer *)in; + line_buffer *buffer = (line_buffer *)in; u08 *line = buffer->line; ForLoop(buffer->nLines) { u32 nameBuff[16]; - line = PushStringIntoIntArray(nameBuff, ArrayCount(nameBuff), line, '\t'); + + // std::string line_str((char *)line); + // std::string name_chr = line_str.substr(0, line_str.find('\t')); + // line_str.erase(0, name_chr.size() + 1); + + line = PushStringIntoIntArray(nameBuff, ArrayCount(nameBuff), line, '\t'); // the buffer of line moved to '\t' contig *cont = 0; - if (ContigHashTableLookup(nameBuff, ArrayCount(nameBuff), &cont) && cont) + if (ContigHashTableLookup(nameBuff, ArrayCount(nameBuff), &cont) && cont) // find the contig and get the contig pointer { - u64 prevLength_genome = (u64)((f64)cont->previousCumulativeLength * (f64)Map_Properties->totalGenomeLength); + u64 prevLength_genome = (u64)((f64)cont->previousCumulativeLength * (f64)Map_Properties->totalGenomeLength); // get the previous cumulative length of the contig and transfer that into bp u32 len = 0; while (*++line != '\t') ++len; @@ -320,35 +410,74 @@ ProcessLine(void *in) while (*++line != '\t') ++len; u64 to_genome = (u64)StringToInt(line, len) + prevLength_genome; - u64 sectionLength_genome = to_genome - from_genome; + u64 bp_left_in_this_bin = to_genome - from_genome; u32 value; - if (StringToInt_Check(line + 1, &value, '\n')) + if (StringToInt_Check(line + 1, &value, '\n')) // get the value of the bedgraph bin { Data_Added = 1; + u64 bin_size = bp_left_in_this_bin; // used to normalise the repeat density data + u32 from_pixel = (u32)(((f64)from_genome / (f64)Map_Properties->totalGenomeLength) * (f64)sizeGraphArray); // coordinate with unit of pixel + u32 to_pixel = (u32)(((f64)to_genome / (f64)Map_Properties->totalGenomeLength) * (f64)sizeGraphArray); - u32 from = (u32)(((f64)from_genome / (f64)Map_Properties->totalGenomeLength) * (f64)sizeGraphArray); - u32 to = (u32)(((f64)to_genome / (f64)Map_Properties->totalGenomeLength) * (f64)sizeGraphArray); + u64 bp_covered_in_current_pixel; + if ( (u64) (from_pixel + 1) * bp_per_pixel <= from_genome) { + bp_covered_in_current_pixel = 0; // because bp_per_pixel is a little bit smaller than the original value, so it can be a negtiave value, as this is an unsigned value, so it will a very large value if not set as 0 + } + else { + bp_covered_in_current_pixel = ((u64)(from_pixel + 1) * bp_per_pixel) - from_genome; // this is the bp covered in the current pixel + } - u64 binSize_ind = ((u64)(from + 1) * binSize_genome) - from_genome; + // if (from_pixel != to_pixel) { + // printf("from_pixel: %d, to_pixel: %d\n", from_pixel, to_pixel); + // } - for ( u32 index = from; - index <= to && index < sizeGraphArray; + for ( u32 index = from_pixel; + index <= to_pixel && index < sizeGraphArray; ++index ) { - u32 nThisBin = (u32)(Min(binSize_ind, sectionLength_genome)); - s32 valueToAdd = (s32)(value * nThisBin); - if (valueToAdd < 0) valueToAdd = s32_max; - - s32 oldValue = __atomic_fetch_add(Graph->values + index, valueToAdd, 0); - if ((s32_max - oldValue) < valueToAdd) - { - s32 cap = s32_max; - __atomic_store(Graph->values + index, &cap, 0); + /* + Iterate over all the pixels covered by the current bin. + + First calculate the number of bp of the first pixel that the current bin can cover, then add the value to graph->values. + Then update the number of bp's left in the current bin, and the number of bp's the current bin can cover. + */ + u32 nThisBin = (u32)(Min(bp_covered_in_current_pixel, bp_left_in_this_bin)); + // s32 valueToAdd = (s32)(value * nThisBin); + if (data_type == data_type_dic["gap"]) { // + // add the value to graph->values + std::unique_lock lock(mtx_global); + Graph_tmp->values[index] = Min(Max(0.f, (f32)value + Graph_tmp->values[index]), 1.0f); // if set the value vector as f32 array, then we can not use the atomic operation. If mutil-thread is used, then we have to use the mutex to protect the values + lock.unlock(); } - - sectionLength_genome -= (u64)nThisBin; - binSize_ind = binSize_genome; + else if (data_type == data_type_dic["repeat_density"]){ // normalise the repeat density data by the length of the bin, and scale that into 0 - 100 + f32 valueToAdd_f = (f32)value / (f32)bin_size * (f32)nThisBin / (f32)bp_per_pixel * 100.0f; + if (valueToAdd_f > 100.f) { + printf("Warning: valueToAdd_f is larger than 100: %f\n", valueToAdd_f); + } + std::unique_lock lock(mtx_global); + Graph_tmp->values[index] += valueToAdd_f; + lock.unlock(); + } + else { + f32 valueToAdd_f = (f32)value * (f32)nThisBin / (f32)bp_per_pixel; + // add the value to graph->values + std::unique_lock lock(mtx_global); + Graph_tmp->values[index] += valueToAdd_f; // if set the value vector as f32 array, then we can not use the atomic operation. If mutil-thread is used, then we have to use the mutex to protect the values + lock.unlock(); + } + // if (valueToAdd < 0) valueToAdd = s32_max; + + // s32 oldValue = __atomic_fetch_add(Graph->values + index, valueToAdd, 0); + + // if ((s32_max - oldValue) < valueToAdd) // if the value is overflow, set this as maximum value (s32_max) + // { + // s32 cap = s32_max; + // __atomic_store(Graph->values + index, &cap, 0); + // } + + bp_left_in_this_bin -= (u64)nThisBin; + bp_covered_in_current_pixel = bp_per_pixel; } if (!(__atomic_add_fetch(&Number_of_Lines_Read, 1, 0) & ((Pow2(Number_of_Lines_Print_Status_Every_Log2)) - 1))) @@ -395,13 +524,21 @@ global_function void NormaliseGraph_Thread(void *in) { + + /* + Normalise the graph data. + + Graph->values = Graph->values * (mapResolution / totalGenomeLength) + + */ normalise_graph_thread_data *data = (normalise_graph_thread_data *)in; u32 start = data->start; u32 nLanes = data->nLanes; u32 mapResolution = data->mapResolution; - f32 texelsPerBp = (f32)((f64)mapResolution / (f64)Map_Properties->totalGenomeLength); - s32 *dataPtr = ((s32 *)Graph->values) + start; + f32 texelsPerBp = (f32)((f64)mapResolution / (f64)Map_Properties->totalGenomeLength); // fraction of pixel per bp + s32 *dataPtr = ((s32 *)Graph->values) + start; // Graph->values with num_pixels s32 values which save the graph data + #ifdef UsingAVX __m256 normFactor = _mm256_set1_ps (texelsPerBp); ForLoop(nLanes) @@ -421,9 +558,9 @@ NormaliseGraph_Thread(void *in) dataPtr += 8; } #else - nLanes *= 2; + nLanes *= 2; // 2048 __m128 normFactor = _mm_set1_ps (texelsPerBp); - ForLoop(nLanes) + ForLoop(nLanes) // 4 values will be processed at one iteration { #pragma clang diagnostic push #pragma GCC diagnostic ignored "-Wcast-align" @@ -437,7 +574,7 @@ NormaliseGraph_Thread(void *in) #pragma GCC diagnostic ignored "-Wcast-align" *(__m128i *)dataPtr = intLane; #pragma clang diagnostic pop - dataPtr += 4; + dataPtr += 4; // process 4 * s32 } #endif } @@ -509,12 +646,16 @@ GrabStdIn() u32 bufferPtr = 0; read_pool *readPool = CreateReadPool(&Working_Set); - readPool->handle = -#ifdef DEBUG - open("test_in", O_RDONLY); -#else - STDIN_FILENO; -#endif + readPool->handle = STDIN_FILENO; // read from stdin + +// old version used to debug +// readPool->handle = +// // #ifdef DEBUG +// // open("data_for_test/repeat_density.bedgraph", O_RDONLY); +// // #else +// // STDIN_FILENO; +// // #endif +// STDIN_FILENO; // read from stdin while running u08 line[KiloByte(16)]; u32 linePtr = 0; @@ -542,9 +683,9 @@ GrabStdIn() line[128] = c; } - if ((u64)linePtr > (Line_Buffer_Size - bufferPtr)) + if ((u64)linePtr > (Line_Buffer_Size - bufferPtr)) // { - buffer->nLines = numLines; + buffer->nLines = numLines; // number of lines in the buffer ThreadPoolAddTask(Thread_Pool, ProcessLine, buffer); buffer = TakeLineBufferFromQueue_Wait(Line_Buffer_Queue); @@ -604,7 +745,9 @@ CopyFile(void *in) } MainArgs -{ +{ + set_num_threads(); + if (ArgCount == 1) { printf("\n%s\n\n", PretextGraph_Version); @@ -619,14 +762,15 @@ MainArgs } if (ArgCount == 2) - { - if (AreNullTerminatedStringsEqual((u08 *)ArgBuffer[1], (u08 *)"--licence")) + { + std::string arg1(ArgBuffer[1]); + if (arg1.find("licen") != std::string::npos) { printf("%s\n", Licence); exit(EXIT_SUCCESS); } - if (AreNullTerminatedStringsEqual((u08 *)ArgBuffer[1], (u08 *)"--thirdparty")) + if (arg1.find("third") != std::string::npos) { printf("%s\n", ThirdParty); exit(EXIT_SUCCESS); @@ -643,30 +787,42 @@ MainArgs u32 copyBufferSize = KiloByte(256); copy_file_data copyFileData; - for ( u32 index = 1; - index < (u32)ArgCount; - ++index ) - { - if (AreNullTerminatedStringsEqual((u08 *)ArgBuffer[index], (u08 *)"-i")) - { - ++index; - pretextFile = ArgBuffer[index]; + u32 index_arg = 1; + while (index_arg < (u32) ArgCount) { + std::string arg_tmp(ArgBuffer[index_arg]); + if (arg_tmp == "-i") { + ++index_arg; + if (index_arg >= (u32) ArgCount) { + PrintError("Input file required for key \'-i\'"); + returnCode = EXIT_FAILURE; + goto end; + } + pretextFile = ArgBuffer[index_arg]; } - else if (AreNullTerminatedStringsEqual((u08 *)ArgBuffer[index], (u08 *)"-n")) - { - ++index; - PushStringIntoIntArray(nameBuffer, ArrayCount(nameBuffer), (u08 *)ArgBuffer[index]); + else if (arg_tmp == "-n") { + ++index_arg; + if (index_arg >= (u32) ArgCount) { + PrintError("Extension name required for key \'-n\'"); + returnCode = EXIT_FAILURE; + goto end; + } + PushStringIntoIntArray(nameBuffer, ArrayCount(nameBuffer), (u08 *)ArgBuffer[index_arg]); } - else if (AreNullTerminatedStringsEqual((u08 *)ArgBuffer[index], (u08 *)"-o")) - { - ++index; - outputFile = ArgBuffer[index]; + else if ( arg_tmp == "-o") { + ++index_arg; + if (index_arg >= (u32) ArgCount) { + PrintError("Ouput file path required for key \'-o\'"); + returnCode = EXIT_FAILURE; + goto end; + } + outputFile = ArgBuffer[index_arg]; } + index_arg ++ ; } - if (pretextFile) - { - PrintStatus("Input file: \'%s\'", pretextFile); + if (pretextFile) { + fprintf(stdout, "[PretextGraph status] :: Pretext file: %s\n", pretextFile); + // PrintStatus("Input file: \'%s\' \n", pretextFile); // if the path of the file is too long, there can be segmentation fault } else { @@ -678,6 +834,21 @@ MainArgs if (nameBuffer[0]) { PrintStatus("Graph name: \'%s\'", (char *)nameBuffer); + + // update the data_type flag + std::string tmp_string((char *)nameBuffer); + if (tmp_string.find("gap") != std::string::npos) { + data_type = data_type_dic["gap"]; + PrintStatus("The extension of the graph name is gap, set the data_type to %d (%s).", data_type, "gap"); + } + else if (tmp_string.find("repeat_density") != std::string::npos) { + data_type = data_type_dic["repeat_density"]; + PrintStatus("The extension of the graph name is repeat_density, set the data_type to %d (%s).", data_type, "repeat_density"); + } + else { + data_type = data_type_dic["default"]; + PrintStatus("The extension of the graph name is default, set the data_type into %d (%s).", data_type, "default"); + } } else { @@ -687,14 +858,15 @@ MainArgs } CreateMemoryArena(Working_Set, MegaByte(128)); - Thread_Pool = ThreadPoolInit(&Working_Set, 4); + Thread_Pool = ThreadPoolInit(&Working_Set, NUM_THREADS); - if (outputFile) + if (outputFile) // if output file is specified, define the copy buffer. Copy the input file to output file. { copyBuffer = PushArray(Working_Set, u08, copyBufferSize); copyInFile = fopen((const char *)pretextFile, "rb"); copyOutFile = fopen((const char *)outputFile, "wb"); - PrintStatus("Output file: \'%s\'", outputFile); + // PrintStatus("Output file: \'%s\'\n", outputFile); + fprintf(stdout, "[PretextGraph status] :: Pretext file: %s\n", outputFile); } else { @@ -729,7 +901,14 @@ MainArgs file = 0; } - if (file) + if (!file) + { + PrintError("Invalid file format. \'%s\' is not a pretext file", pretextFile); + returnCode = EXIT_FAILURE; + goto end; + } + + /* Read file */ { u32 nBytesHeaderComp; u32 nBytesHeader; @@ -739,8 +918,15 @@ MainArgs u08 *compressionBuffer = PushArray(Working_Set, u08, nBytesHeaderComp); fread(compressionBuffer, 1, nBytesHeaderComp, file); - if (!libdeflate_deflate_decompress(decompressor, (const void *)compressionBuffer, nBytesHeaderComp, (void *)header, nBytesHeader, NULL)) - { + if (!libdeflate_deflate_decompress(decompressor, (const void *)compressionBuffer, nBytesHeaderComp, (void *)header, nBytesHeader, NULL)) // decompress the header + { + + if (outputFile && !(copyInFile && copyOutFile)) // outputfile exists, and not both of copyin and copyout files exist, then give the error + { + PrintError("Error copying input file"); + returnCode = EXIT_FAILURE; + } + if (!outputFile || (copyInFile && copyOutFile)) { if (outputFile) @@ -838,15 +1024,20 @@ MainArgs u32 mapResolution = Map_Properties->textureResolution * Map_Properties->numberOfTextures1D; Graph = PushStruct(Working_Set, graph); - u32 graphPlusNameBufferSize = (sizeof(s32) * mapResolution) + sizeof(nameBuffer); + u32 graphPlusNameBufferSize = (sizeof(f32) * mapResolution) + sizeof(nameBuffer); // add the f32 also into this buffer u08 *graphPlusNameBuffer = PushArray(Working_Set, u08, graphPlusNameBufferSize, 5); #pragma clang diagnostic push #pragma GCC diagnostic ignored "-Wcast-align" Graph->values = (volatile s32 *)(graphPlusNameBuffer + sizeof(nameBuffer)); //PushArray(Working_Set, volatile s32, mapResolution, 5); + // Graph->values_f = (volatile s32 *)(graphPlusNameBuffer + sizeof(nameBuffer) + (sizeof(s32) * mapResolution)); //PushArray(Working_Set, volatile s32, mapResolution, 5); #pragma clang diagnostic pop - memset((void *)Graph->values, 0, sizeof(s32) * mapResolution); + memset((void *)Graph->values, 0, sizeof(s32) * mapResolution); // initalise the graph values to 0 + // initialise the graph values_f and set values to 0. + Graph_tmp = (graph_f *) malloc(sizeof(graph_f)); + // Graph_tmp->values = (f32 *) malloc(sizeof(f32) * mapResolution); + Graph_tmp->values = (f32 *) calloc(mapResolution, sizeof(f32)); // make sure the values are set to 0 - ThreadPoolAddTask(Thread_Pool, GrabStdIn, 0); + ThreadPoolAddTask(Thread_Pool, GrabStdIn, 0); // reading the graph data ThreadPoolWait(Thread_Pool); fprintf(stdout, "\n"); @@ -858,6 +1049,7 @@ MainArgs exit(EXIT_FAILURE); } + /* PrintStatus("Normalising graph data..."); { normalise_graph_thread_data data[4]; @@ -866,12 +1058,12 @@ MainArgs data[2].mapResolution = mapResolution; data[3].mapResolution = mapResolution; - u32 nLanes = (mapResolution + 7) >> 3; - u32 halfLanes = nLanes >> 1; - u32 quaterLanes = halfLanes >> 1; + u32 nLanes = (mapResolution + 7) >> 3; // 4096， total number of units to process, 1 unit = 8 pixels + u32 halfLanes = nLanes >> 1; // 2048 half of the total units + u32 quaterLanes = halfLanes >> 1; // 1024 quater of the total units - data[0].nLanes = quaterLanes; - data[0].start = 0; + data[0].nLanes = quaterLanes; // 1024 + data[0].start = 0; data[1].nLanes = halfLanes - data[0].nLanes; data[1].start = data[0].nLanes << 3; @@ -888,6 +1080,19 @@ MainArgs ThreadPoolAddTask(Thread_Pool, NormaliseGraph_Thread, (data + 3)); ThreadPoolWait(Thread_Pool); } + */ + + PrintStatus("Transfer f32 to s32..."); + { + // as only the s32 values can be accepted by PretextView + ForLoop(mapResolution) + { + Graph->values[index] = (s32)(Graph_tmp->values[index]); + } + free(Graph_tmp->values); + Graph_tmp->values = NULL; + free(Graph_tmp); + } PrintStatus("Saving graph..."); { @@ -909,7 +1114,8 @@ MainArgs { *namePtr++ = nameBuffer[index]; } - + + // compress file by libdeflate, and write the compressed data into the file if (compressor && (compSize = (u32)libdeflate_deflate_compress(compressor, (const void *)graphPlusNameBuffer, graphPlusNameBufferSize, (void *)compBuffer, compBufferSize))) { fwrite(graphMagic, 1, sizeof(graphMagic), graphOutputFile); @@ -924,11 +1130,6 @@ MainArgs } PrintStatus("Done"); } - else - { - PrintError("Error copying input file"); - returnCode = EXIT_FAILURE; - } } else { @@ -936,11 +1137,6 @@ MainArgs returnCode = EXIT_FAILURE; } } - else - { - PrintError("\'%s\' is not a pretext file", pretextFile); - returnCode = EXIT_FAILURE; - } } else { diff --git a/README.md b/README.md index 4043bde..6fa4963 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,10 @@ # PretextGraph Converts bedgraph formatted data and embeds inside a Pretext contact map. + # Bioconda +**Note: currently seems not able to install this tool via Bioconda. (will be fixed soon I hope...)** + All commandline Pretext tools for Unix (Linux and Mac) are available on [bioconda](https://bioconda.github.io/).
The full suite of Pretext tools can be installed with @@ -20,15 +23,52 @@ Or, just PretextGraph can be installed with PretextGraph reads bedgraph formatted data from `stdin`, e.g:
```bash > zcat bedgraph.file.gz | PretextGraph -i input.pretext -n "graph name" +> PretextGraph -i input.pretext -n "graph name" < /path/to/bedgraph.file > bigWigToBedGraph bigwig.file /dev/stdout | PretextGraph -i input.pretext -n "graph name" ``` -Important: only non-negative integer data is supported. -# Options --i input Pretext file, required. Sequence names in the Pretext file must match sequence names in the bedgraph data; although relative sort order is unimportant.
--n graph name, required. A name for the graph.
+## Inject extensions in to an existing `.pretext` file +Usage: +```bash +PretextGraph -i /path/to/input/.pretext/file -n "name_of_extension" [-o /path/to/the/output/.pretext/file] < /path/to/extension/data/file +``` +NOTE: while using the `bedgraph` as the input for extension file, the newline character must be `\n`, and the separator must be `\t`. There is an example `repeat_density.bedgraph` [file](data_for_test). The file is as follows: +``` +#1_usercol 2_usercol 3_usercol N_density +chr1 0 10000 5107 +chr1 10000 20000 3579 +chr1 20000 30000 4850 +chr1 30000 40000 2643 +chr1 40000 50000 2309 +chr1 50000 60000 3308 +chr1 60000 70000 3605 +chr1 70000 80000 3908 +``` +***Important: only non-negative integer data is supported. All values in 4th column should be `int`.*** + +## NOTE: Special process different extensions + +There are different extensions types, such as the coverage, repeat_density, gap... + +Currently, there are 3 ways to transform the data into extension graphs: +```cpp +std::unordered_map data_type_dic{ // use this data_type + {"default", 0, }, + {"repeat_density", 1}, // as this is counted in every single bin, so we need to normalise this within the bin + {"gap", 2}, // +}; +``` +- `0`: default, just add the weighted value of every bin to the `graph->values[index]`; +- `1`: `repeat_density`, before add the value of every bin to the `graph->values[index]`, the value is normlised by the `bin_size` as the `repeat.bedgraph` counts number of repeat bps in one single bin; +- `2`: `gap`, if there is gap in the related pixel, then `graph->values[index]` is set to `1` (no matter how many gaps), if no gaps within the range related with the pixel, the value is `0`. + + --o output Pretext file, optional. If no output is specified the graph data will be appended to the input file.
+## Options +`-i` input Pretext file, required. Sequence names in the Pretext file must match sequence names in the bedgraph data; although relative sort order is unimportant.
+`-i` graph name, required. A name for the graph.
+ +`-o` output Pretext file, optional. If no output is specified the graph data will be appended to the input file.
# Requirments, running 4 cpu cores
@@ -43,6 +83,8 @@ PretextGraph uses the following third-party libraries:
* [stb_sprintf.h](https://github.com/nothings/stb/blob/master/stb_sprintf.h) # Installation + +## meson Requires: * clang >= 11.0.0 * meson >= 0.57.1 @@ -53,3 +95,15 @@ cd builddir meson compile meson test meson install +``` + +## cmake +```bash +git submodule update --init --recursive +cd libdeflate +make || exit 1 +cd .. +cmake -DCMAKE_BUILD_TYPE=Release -B build_cmake || exit 1 +cd build_cmake +make +``` \ No newline at end of file diff --git a/data_for_test/repeat_density.bedgraph b/data_for_test/repeat_density.bedgraph new file mode 100644 index 0000000..1759dd8 --- /dev/null +++ b/data_for_test/repeat_density.bedgraph @@ -0,0 +1,100 @@ +#1_usercol 2_usercol 3_usercol N_density +chr1 0 10000 5107 +chr1 10000 20000 3579 +chr1 20000 30000 4850 +chr1 30000 40000 2643 +chr1 40000 50000 2309 +chr1 50000 60000 3308 +chr1 60000 70000 3605 +chr1 70000 80000 3908 +chr1 80000 90000 6822 +chr1 90000 100000 8680 +chr1 100000 110000 7409 +chr1 110000 120000 2288 +chr1 120000 130000 4497 +chr1 130000 140000 3476 +chr1 140000 150000 4639 +chr1 150000 160000 4201 +chr1 160000 170000 4839 +chr1 170000 180000 6350 +chr1 180000 190000 4564 +chr1 190000 200000 4115 +chr1 200000 210000 4844 +chr1 210000 220000 7143 +chr1 220000 230000 4005 +chr1 230000 240000 2339 +chr1 240000 250000 2470 +chr1 250000 260000 1753 +chr1 260000 270000 2220 +chr1 270000 280000 1649 +chr1 280000 290000 2145 +chr1 290000 300000 1569 +chr1 300000 310000 2912 +chr1 310000 320000 5878 +chr1 320000 330000 6025 +chr1 330000 340000 2718 +chr1 340000 350000 854 +chr1 350000 360000 1429 +chr1 360000 370000 685 +chr1 370000 380000 2004 +chr1 380000 390000 1957 +chr1 390000 400000 1205 +chr1 400000 410000 3517 +chr1 410000 420000 3924 +chr1 420000 430000 4427 +chr1 430000 440000 3216 +chr1 440000 450000 6135 +chr1 450000 460000 1408 +chr1 460000 470000 2316 +chr1 470000 480000 1367 +chr1 480000 490000 3073 +chr1 490000 500000 2916 +chr1 500000 510000 3323 +chr1 510000 520000 1031 +chr1 520000 530000 4350 +chr1 530000 540000 4855 +chr1 540000 550000 3632 +chr1 550000 560000 2339 +chr1 560000 570000 2669 +chr1 570000 580000 1409 +chr1 580000 590000 2537 +chr1 590000 600000 885 +chr1 600000 610000 2455 +chr1 610000 620000 4682 +chr1 620000 630000 3525 +chr1 630000 640000 2582 +chr1 640000 650000 1614 +chr1 650000 660000 2512 +chr1 660000 670000 1017 +chr1 670000 680000 2319 +chr1 680000 690000 3164 +chr1 690000 700000 3724 +chr1 700000 710000 3952 +chr1 710000 720000 4422 +chr1 720000 730000 3271 +chr1 730000 740000 2052 +chr1 740000 750000 2391 +chr1 750000 760000 1667 +chr1 760000 770000 1425 +chr1 770000 780000 1533 +chr1 780000 790000 2451 +chr1 790000 800000 6692 +chr1 800000 810000 3113 +chr1 810000 820000 2741 +chr1 820000 830000 4815 +chr1 830000 840000 3300 +chr1 840000 850000 6067 +chr1 850000 860000 4311 +chr1 860000 870000 4285 +chr1 870000 880000 1951 +chr1 880000 890000 5826 +chr1 890000 900000 4019 +chr1 900000 910000 4603 +chr1 910000 920000 4332 +chr1 920000 930000 3480 +chr1 930000 940000 6847 +chr1 940000 950000 3758 +chr1 950000 960000 3227 +chr1 960000 970000 3491 +chr1 970000 980000 1452 +chr1 980000 990000 3280 diff --git a/meson.build b/meson.build index 5b475cb..21dcac0 100644 --- a/meson.build +++ b/meson.build @@ -4,7 +4,10 @@ project('PretextGraph', ['c', 'cpp'], version: '0.0.6' ) -deps = [dependency('threads')] +deps = [ + dependency('threads'), + # dependency('pybind11') + ] base_flags = ['-Ofast'] if get_option('buildtype') == 'debug' base_flags = ['-O0', '-g', '-DDEBUG'] @@ -12,7 +15,13 @@ if get_option('buildtype') == 'debug' endif base_flags += ['-DPV=' + meson.project_version()] -extensions = {'.avx2': ['-mavx2'], '.avx': ['-mavx'], '.sse42': ['-msse4.2'], '.sse41': ['-msse4.1']} +extensions = { + # '': [], + '.avx2': ['-mavx2'], + '.avx': ['-mavx'], + '.sse42': ['-msse4.2'], + '.sse41': ['-msse4.1'] + } foreach ext, flag : extensions flags = base_flags + flag