libStatGen Software
1
|
00001 #include <iostream> 00002 #include <time.h> 00003 #include "GenomeSequence.h" 00004 #include "InputFile.h" 00005 00006 void readDbsnp(mmapArrayBool_t& dbSNP, const char* fileName, GenomeSequence& ref); 00007 00008 int main(int argc, char ** argv) 00009 { 00010 // time_t startTime; 00011 // time_t endTime; 00012 // startTime = time(NULL); 00013 GenomeSequence* refPtr = new GenomeSequence("testFiles/chr1_partial.fa"); 00014 // endTime = time(NULL); 00015 00016 // std::cerr << "Time to read reference: " << endTime - startTime << std::endl; 00017 if(refPtr == NULL) 00018 { 00019 std::cerr << "Failed to read the reference\n"; 00020 return(-1); 00021 } 00022 std::cerr << "\nStandard VCF DBSNP test\n"; 00023 mmapArrayBool_t dbsnpArray1; 00024 const char* dbsnpFileName = "testFiles/dbsnp.vcf"; 00025 00026 // startTime = time(NULL); 00027 refPtr->loadDBSNP(dbsnpArray1, dbsnpFileName); 00028 // endTime = time(NULL); 00029 // std::cerr << "Time to read dbsnp through reference: " << endTime - startTime << std::endl; 00030 00031 00032 genomeIndex_t mapPos = 00033 refPtr->getGenomePosition("1", 10233); 00034 std::cerr << "dbsnp " << mapPos << ": " 00035 << dbsnpArray1[mapPos] << std::endl; 00036 std::cerr << "dbsnp " << mapPos+1 << ": " 00037 << dbsnpArray1[mapPos+1] << std::endl; 00038 std::cerr << "dbsnp " << mapPos+2 << ": " 00039 << dbsnpArray1[mapPos+2] << std::endl; 00040 00041 00042 std::cerr << "\nGZIP VCF DBSNP test\n"; 00043 00044 mmapArrayBool_t dbsnpArray2; 00045 dbsnpFileName = "testFiles/dbsnp.vcf.gz"; 00046 00047 // startTime = time(NULL); 00048 refPtr->loadDBSNP(dbsnpArray2, dbsnpFileName); 00049 // endTime = time(NULL); 00050 // std::cerr << "Time to read dbsnp through reference: " << endTime - startTime << std::endl; 00051 00052 00053 mapPos = refPtr->getGenomePosition("1", 10233); 00054 std::cerr << "dbsnp " << mapPos << ": " 00055 << dbsnpArray2[mapPos] << std::endl; 00056 std::cerr << "dbsnp " << mapPos+1 << ": " 00057 << dbsnpArray2[mapPos+1] << std::endl; 00058 std::cerr << "dbsnp " << mapPos+2 << ": " 00059 << dbsnpArray2[mapPos+2] << std::endl; 00060 return(0); 00061 }