diff --git a/CMakeLists.txt b/CMakeLists.txt index 3efe616..496eaf6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,7 +54,7 @@ if("${CMAKE_CXX_COMPILER_ID}" MATCHES "GNU") list(APPEND TGT_WARN_FLAGS "-Wno-int-in-bool-context") endif() - if(GCC_VERSION VERSION_GREATER_EQUAL "9.1") + if(GCC_VERSION VERSION_GREATER_EQUAL "9.2") list(APPEND TGT_WARN_FLAGS "-Wno-deprecated-copy") endif() diff --git a/include/FASTAParser.hpp b/include/FASTAParser.hpp index ddcfc6a..1292cb5 100644 --- a/include/FASTAParser.hpp +++ b/include/FASTAParser.hpp @@ -12,11 +12,12 @@ class FASTAParser { public: FASTAParser(); FASTAParser(const std::string& fname); - void populateTargets(std::vector& transcripts); + void populateTargets(std::vector& transcripts, uint32_t readLength); void populateIntronTargets( std::vector& refs, std::string& intronFileName, - std::unordered_map& transcriptNameMap + std::unordered_map& transcriptNameMap, + uint32_t readLength ) ; void updateTranscriptLevelIntron(std::vector& transcripts, std::unordered_map& transcriptNameMap diff --git a/include/MatrixParser.hpp b/include/MatrixParser.hpp index 7c1adec..ed5a581 100644 --- a/include/MatrixParser.hpp +++ b/include/MatrixParser.hpp @@ -123,13 +123,13 @@ public : int ind = rand() % trVec.size() ; auto& fixTrInfo = trVec[ind] ; - if((fixTrInfo.end - fixTrInfo.start) < READ_LEN){ + if((fixTrInfo.end - fixTrInfo.start) < refInfo.readLength){ std::cerr << "Should not happend REPORT IT \n" ; std::exit(1) ; } int mid = fixTrInfo.start + std::round((fixTrInfo.end - fixTrInfo.start)/2) ; - if((refInfo.transcripts[fixTrInfo.tid].RefLength - mid) < READ_LEN){ + if((refInfo.transcripts[fixTrInfo.tid].RefLength - mid) < refInfo.readLength){ checked++; continue ; } diff --git a/include/ProgOpts.hpp b/include/ProgOpts.hpp index e88e9b6..3f0926b 100644 --- a/include/ProgOpts.hpp +++ b/include/ProgOpts.hpp @@ -135,6 +135,7 @@ class EstimateOptions { std::string gene2txpFile{""} ; std::string eqClassFolder{""} ; std::string refFile ; + uint32_t ReadLength{100}; std::string bfhFile{""} ; std::string outDir{""} ; diff --git a/include/ReferenceInfo.hpp b/include/ReferenceInfo.hpp index 52ccffa..5377b3c 100644 --- a/include/ReferenceInfo.hpp +++ b/include/ReferenceInfo.hpp @@ -17,15 +17,17 @@ class Reference{ Reference( std::string& fastaFileIn, std::string& gene2txpFileIn, + uint32_t readLengthIn, std::shared_ptr& consoleLogIn ){ consoleLog = consoleLogIn ; gene2txpFile = gene2txpFileIn ; fastaFile = fastaFileIn ; FASTAParser fastaParser(fastaFile) ; + readLength = readLengthIn; { - fastaParser.populateTargets(transcripts) ; + fastaParser.populateTargets(transcripts, readLength) ; size_t trId{0} ; for(auto& tr : transcripts){ transcriptNameMap[tr.RefName] = trId++ ; @@ -43,7 +45,8 @@ class Reference{ fastaParser.populateIntronTargets( transcripts, intronFastaFile, - transcriptNameMap + transcriptNameMap, + readLength ) ; consoleLog->info("Intron file {} is read", intronFastaFile); }else{ @@ -210,6 +213,7 @@ class Reference{ std::unordered_map transcript2geneMap ; std::unordered_map transcriptNameMap ; size_t numOfTranscripts ; + uint32_t readLength; std::shared_ptr consoleLog; } ; diff --git a/include/macros.hpp b/include/macros.hpp index cb906db..495922a 100644 --- a/include/macros.hpp +++ b/include/macros.hpp @@ -5,12 +5,12 @@ #define MAX_FRAGLENGTH 1000 #define MAX_QUEUE_SIZE 20 -#define READ_LEN 100 +//#define READ_LEN 100 #define MAXNUM 100000 -#define FRAGMENT_END_DIST 404 + READ_LEN -#define FRAGMENT_START_DIST 53 + READ_LEN -#define FRAGMENT_RANGE (FRAGMENT_END_DIST - FRAGMENT_START_DIST) +// #define FRAGMENT_END_DIST 404 + READ_LEN +// #define FRAGMENT_START_DIST 53 + READ_LEN +// #define FRAGMENT_RANGE (FRAGMENT_END_DIST - FRAGMENT_START_DIST) //#define CB_LENGTH 16 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 228a2c8..a7d11c8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,12 +2,12 @@ include_directories( ${TBB_INCLUDE_DIRS} ) -link_directories( - ${GAT_SOURCE_DIR}/external/install/lib - DESTINATION $(TBB_LIBRARY_DIRS) - FILES_MATCHING PATTERN "libtbb*.${SHARED_LIB_EXTENSION}*" +# link_directories( +# ${GAT_SOURCE_DIR}/external/install/lib +# DESTINATION $(TBB_LIBRARY_DIRS) +# FILES_MATCHING PATTERN "libtbb*.${SHARED_LIB_EXTENSION}*" -) +# ) set( GRAPH_LIB_SRC diff --git a/src/FASTAParser.cpp b/src/FASTAParser.cpp index 412b639..f612dcc 100644 --- a/src/FASTAParser.cpp +++ b/src/FASTAParser.cpp @@ -145,7 +145,8 @@ void FASTAParser::updateTranscriptLevelIntron( void FASTAParser::populateTargets( - std::vector& refs + std::vector& refs, + uint32_t readLength ) { using single_parser = fastx_parser::FastxParser; @@ -192,7 +193,7 @@ void FASTAParser::populateTargets( // std::string& seq = j->data[i].seq; std::string& seq = read.seq; size_t readLen = seq.length(); - if(readLen < READ_LEN) + if(readLen < readLength) continue ; // Replace non-ACGT bases @@ -233,7 +234,8 @@ void FASTAParser::populateTargets( void FASTAParser::populateIntronTargets( std::vector& refs, std::string& intronFileName, - std::unordered_map& transcriptNameMap + std::unordered_map& transcriptNameMap, + uint32_t readLength ) { using single_parser = fastx_parser::FastxParser; @@ -279,7 +281,7 @@ void FASTAParser::populateIntronTargets( // std::string& seq = j->data[i].seq; std::string& seq = read.seq; size_t readLen = seq.length(); - if(readLen < READ_LEN) + if(readLen < readLength) continue ; // Replace non-ACGT bases diff --git a/src/GFAReader.cpp b/src/GFAReader.cpp index 2ca17fd..d04db7c 100644 --- a/src/GFAReader.cpp +++ b/src/GFAReader.cpp @@ -89,7 +89,7 @@ void GFAReader::updateEqClass( if(unitigMap.find(contigId) != unitigMap.end()){ auto sizeOfUnitig = unitigMap[contigInfo.first].size() ; - if(sizeOfUnitig < READ_LEN){ + if(sizeOfUnitig < refInfo.readLength){ std::cerr << sizeOfUnitig << " --- possible twopaco bug !!!!\n" ; std::exit(1) ; } @@ -166,7 +166,7 @@ void GFAReader::parseFile( bool foundOverlapSize{false}; // size_t overlapsize = READ_LEN-1; - size_t overlapsize = READ_LEN + 1; + size_t overlapsize = refInfo.readLength + 1; consoleLog->info("Predicted overlap size: {}", overlapsize); file.reset(new std::ifstream(gfaFileName_)) ; diff --git a/src/MatrixParser.cpp b/src/MatrixParser.cpp index ea77cd4..ab39c1a 100644 --- a/src/MatrixParser.cpp +++ b/src/MatrixParser.cpp @@ -866,7 +866,7 @@ void DataMatrix::loadAlevinData( auto tcInfoVec = dbgPtr->eqClassMap[seg][tid] ; for(auto tInfo : tcInfoVec){ if(refInfo.transcripts[tid].RefLength - tInfo.eposInContig <= MAX_FRAGLENGTH){ - if(tInfo.eposInContig - tInfo.sposInContig < READ_LEN){ + if(tInfo.eposInContig - tInfo.sposInContig < refInfo.readLength){ consoleLog->error("encountered a contig shorter than read length", "this is not permitted currently" ); diff --git a/src/Minnow.cpp b/src/Minnow.cpp index 6edd5c8..40139f9 100644 --- a/src/Minnow.cpp +++ b/src/Minnow.cpp @@ -107,6 +107,9 @@ int main(int argc, char* argv[]) { value(ensure_file_exists,"reference", estimateOpt.refFile)) % "transcript fasta file", + (option("--ReadLength") & + value("Read length", estimateOpt.ReadLength)) % "read length by default is 100", + (required("--g2t") & value(ensure_file_exists,"gene_tr", estimateOpt.gene2txpFile)) % "tab separated list of Gene to Transcirpt mapping", @@ -163,7 +166,7 @@ int main(int argc, char* argv[]) { (option("--CBLength") & value("Cell barcode length", simulateOpt.CBLength)) % "Cell barcode length by default is 16", (option("--UMILength") & value("UMI length length", simulateOpt.UMILength)) % "Cell barcode length by default is 10", - (option("--ReadLength") & value("Read length", simulateOpt.UMILength)) % "read length by default is 100", + (option("--ReadLength") & value("Read length", simulateOpt.ReadLength)) % "read length by default is 100", // (option("--alevin-mode").set(simulateOpt.alevinMode, true)) % // "The program would assume that the input matrix is obtained from Alevin", diff --git a/src/MinnowEstimate.cpp b/src/MinnowEstimate.cpp index 07719d5..cdf6164 100644 --- a/src/MinnowEstimate.cpp +++ b/src/MinnowEstimate.cpp @@ -27,12 +27,14 @@ int minnowEstimate(EstimateOptions& eopts){ auto gene2txpFile = eopts.gene2txpFile; auto outDir = eopts.outDir; auto bfhFile = eopts.bfhFile; + auto readLength = eopts.ReadLength; consoleLog->info("Reading reference sequences ...") ; Reference refInfo( refFileName, gene2txpFile, + readLength, consoleLog ) ; refInfo.updateGene2TxpMap(); diff --git a/src/MinnowSimulate.cpp b/src/MinnowSimulate.cpp index d655540..4a2c66b 100644 --- a/src/MinnowSimulate.cpp +++ b/src/MinnowSimulate.cpp @@ -43,7 +43,7 @@ using SpinLockT = std::mutex; #define MAX_FRAGLENGTH 1000 #define MAX_QUEUE_SIZE 20 -#define READ_LEN 100 +// #define READ_LEN 100 #define FRAGMENT_END_DIST 404 + READ_LEN // this is from the empirical e^6 #define FRAGMENT_START_DIST 53 + READ_LEN // this is from the empirical limit e^4 @@ -52,6 +52,7 @@ using SpinLockT = std::mutex; uint32_t CB_LENGTH ; uint32_t UMI_LENGTH ; uint32_t POOL_SIZE ; +uint32_t READ_LEN ; #define _verbose(fmt, args...) fprintf(stderr, fmt, ##args) @@ -1913,6 +1914,8 @@ void minnowSimulate(SimulateOptions& simOpts){ CB_LENGTH = expConfig.barcodeLength ; UMI_LENGTH = expConfig.umiLength ; POOL_SIZE = expConfig.maxValue ; + READ_LEN = simOpts.ReadLength ; + uint32_t readLength = simOpts.ReadLength; //auto& numOfSampleCells = simOpts.sampleCells ; @@ -1939,6 +1942,7 @@ void minnowSimulate(SimulateOptions& simOpts){ Reference refInfo( refFileName, gene2txpFile, + readLength, consoleLog ) ; consoleLog->info("Reference sequence is loaded ...") ; diff --git a/src/MinnowValidator.cpp b/src/MinnowValidator.cpp index f20b82a..0282dc5 100644 --- a/src/MinnowValidator.cpp +++ b/src/MinnowValidator.cpp @@ -29,7 +29,8 @@ #define _verbose(fmt, args...) fprintf(stderr, fmt, ##args) -#define READ_LEN 100 +//#define READ_LEN 100 +uint32_t READ_LEN ; struct ValidateOpt{ std::string gfaFile{""} ; @@ -38,6 +39,7 @@ struct ValidateOpt{ std::string outFile{""} ; int numThreads{0} ; int edit_max_lim{0} ; + uint32_t read_length{100} ; } ; @@ -45,6 +47,7 @@ struct DumpOpt{ std::string fastqFile{""} ; std::string t2gFile{""} ; std::string outFile{""} ; + uint32_t read_length{100}; } ; @@ -53,6 +56,7 @@ struct ExtractOpt{ std::string rightFastqFile{""} ; std::string queryFileName{""} ; std::string outFile{""} ; + uint32_t read_length{100}; } ; @@ -108,6 +112,7 @@ void queryCellName( std::string& leftFastq, std::string& rightFastq, std::string& queryFileName, + uint32_t read_length, std::string& outFile ){ @@ -117,6 +122,8 @@ void queryCellName( std::string right_out_file = outFile + "_2.fastq.gz" ; std::string left_out_file = outFile + "_1.fastq.gz" ; + READ_LEN = read_length; + zstr::ofstream leftFileStream(left_out_file.c_str(), std::ios::out ) ; zstr::ofstream rightFileStream(right_out_file.c_str(), std::ios::out ) ; @@ -300,12 +307,15 @@ void refValidate( std::string& fastqFile, std::string& referenceFile, int& edit_max_lim, + uint32_t read_length, std::string& outFile ){ + READ_LEN = read_length; + std::vector transcripts; std::cout << "Reading reference file\n" ; FASTAParser fastaParser(referenceFile) ; - fastaParser.populateTargets(transcripts) ; + fastaParser.populateTargets(transcripts, READ_LEN) ; std::unordered_map trMap ; trMap.reserve(transcripts.size()) ; @@ -416,12 +426,14 @@ void gfaValidate( std::string& gfaFile, std::string& fastqFile, int& edit_max_lim, + uint32_t read_length, std::string& outFile, std::shared_ptr& consoleLog ){ GFAReader gfaObj(gfaFile, consoleLog) ; + READ_LEN = read_length; gfaObj.readUnitigs() ; std::map editDistanceMap ; @@ -517,6 +529,7 @@ void validate(ValidateOpt& valOpts){ valOpts.fastqFile, valOpts.referenceFile, valOpts.edit_max_lim, + valOpts.read_length, valOpts.outFile ) ; }else{ @@ -524,6 +537,7 @@ void validate(ValidateOpt& valOpts){ valOpts.gfaFile, valOpts.fastqFile, valOpts.edit_max_lim, + valOpts.read_length, valOpts.outFile, consoleLog ) ; @@ -536,6 +550,7 @@ void extract(ExtractOpt& valOpts){ valOpts.leftFastqFile, valOpts.rightFastqFile, valOpts.queryFileName, + valOpts.read_length, valOpts.outFile ) ; @@ -568,6 +583,10 @@ int main(int argc, char* argv[]) { value("cell-list", extOpt.queryFileName)) % "list of cell files", + (option("-l", "--read-length") & + value("read-length", extOpt.read_length)) % + "read length in bp", + (option("-o", "--output") & value("output", extOpt.outFile)) % "output fastq file" @@ -585,6 +604,10 @@ int main(int argc, char* argv[]) { (option("-t", "--reference") & value("transcript", dumpOpt.t2gFile)) % "tsv file with transcript to gene mapping", + + (option("-l", "--read-length") & + value("read-length", extOpt.read_length)) % + "read length in bp", (option("-o", "--output") & value("output", dumpOpt.outFile)) % @@ -614,6 +637,10 @@ int main(int argc, char* argv[]) { (option("-p", "--num-threads") & value("number-threads", valOpts.numThreads)) % "number of Threads needed for the parsing", + + (option("-l", "--read-length") & + value("read-length", extOpt.read_length)) % + "read length in bp", (option("-o", "--output") & value("output", valOpts.outFile)) %