Skip to content

Instantly share code, notes, and snippets.

@holtgrewe
Created October 4, 2014 11:52
Show Gist options
  • Save holtgrewe/eb22e3a6bab6efd5c184 to your computer and use it in GitHub Desktop.
Save holtgrewe/eb22e3a6bab6efd5c184 to your computer and use it in GitHub Desktop.
Uploaded files from SeqAn Trac Ticket.
// ==========================================================================
// mini_bowtie
// ==========================================================================
// Copyright (c) 2006-2012, Knut Reinert, FU Berlin
// All rights reserved.
#include <iostream>
#include <time.h>
#include <seqan/basic.h>
#include <seqan/sequence.h>
#include <seqan/file.h>
#include <seqan/index.h>
#include <seqan/index/index_fm.h>
#include <seqan/store.h>
using namespace seqan;
int main(int argc, char *argv[])
{
clock_t indexing_time = clock();
// type definitions
typedef Dna5String TString;
typedef StringSet<TString> TStringSet;
typedef Index<TStringSet, FMIndex<> > TIndex;
typedef Iterator<TIndex, TopDown<ParentLinks<> > >::Type TIter;
// reading the command line arguments
if (argc < 2) {
std::cerr << "Invalid number of arguments." << std::endl
<< "USAGE: fmIndex GENOME.fasta" << std::endl;
return 1;
}
// declaration and initialization of the fragment store
FragmentStore<> fragStore;
// loading the contigs
if (!loadContigs(fragStore, argv[1])) {
std::cerr << "Loading GENOME failed" << std::endl;
return 1;
}
// copy the contigs to the string set "mySeq"
TStringSet mySeq;
for (unsigned i = 0; i < length(fragStore.contigStore); ++i) {
appendValue(mySeq, fragStore.contigStore[i].seq);
std::cout << "Sequence " << i << " has " << length(mySeq[i]) << " nucleotides" << std::endl;
}
//std::cout << mySeq[0] << std::endl;
// clear contigs to free memory
clearContigs(fragStore);
////////////////////////////////////////////////////////////////////////
// REVERSE the ref genome !!! REVERSE !!! REVERSE !!!
reverse(mySeq);
////////////////////////////////////////////////////////////////////////
// FM INDEXING
TIndex fmIndex(mySeq); // the index is not created yet
TIter it(fmIndex); // the index is created after this line
char indexFile[100];
strcpy(indexFile, argv[1]);
strcat(indexFile, ".index");
if (save(fmIndex, indexFile)) {
std::cout << "Index saved!" << std::endl;
}
else {
std::cout << "Failed to save index" << std::endl;
return 1;
}
indexing_time = clock() - indexing_time;
std::cout << "Total indexing time: " << ((double)indexing_time)/CLOCKS_PER_SEC << " seconds" << std::endl;
return 0;
}
// ==========================================================================
// mini_bowtie
// ==========================================================================
// Copyright (c) 2006-2012, Knut Reinert, FU Berlin
// All rights reserved.
#include <iostream>
#include <time.h>
#include <seqan/basic.h>
#include <seqan/sequence.h>
#include <seqan/file.h>
#include <seqan/index.h>
#include <seqan/index/index_fm.h>
#include <seqan/store.h>
using namespace seqan;
int main(int argc, char *argv[])
{
clock_t indexing_time = clock();
// type definitions
typedef Dna5String TString;
typedef StringSet<TString> TStringSet;
typedef Index<TStringSet, FMIndex<> > TIndex;
typedef Iterator<TIndex, TopDown<ParentLinks<> > >::Type TIter;
// reading the command line arguments
if (argc < 2) {
std::cerr << "Invalid number of arguments." << std::endl
<< "USAGE: fmIndex GENOME.fasta" << std::endl;
return 1;
}
// declaration and initialization of the fragment store
FragmentStore<> fragStore;
// loading the contigs
if (!loadContigs(fragStore, argv[1])) {
std::cerr << "Loading GENOME failed" << std::endl;
return 1;
}
// copy the contigs to the string set "mySeq"
TStringSet mySeq;
for (unsigned i = 0; i < length(fragStore.contigStore); ++i) {
appendValue(mySeq, fragStore.contigStore[i].seq);
std::cout << "Sequence " << i << " has " << length(mySeq[i]) << " nucleotides" << std::endl;
}
//std::cout << mySeq[0] << std::endl;
// clear contigs to free memory
clearContigs(fragStore);
////////////////////////////////////////////////////////////////////////
// REVERSE the ref genome !!! REVERSE !!! REVERSE !!!
reverse(mySeq);
////////////////////////////////////////////////////////////////////////
// FM INDEXING
TIndex fmIndex(mySeq); // the index is not created yet
TIter it(fmIndex); // the index is created after this line
char indexFile[100];
strcpy(indexFile, argv[1]);
strcat(indexFile, ".index");
if (save(fmIndex, indexFile)) {
std::cout << "Index saved!" << std::endl;
}
else {
std::cout << "Failed to save index" << std::endl;
return 1;
}
indexing_time = clock() - indexing_time;
std::cout << "Total indexing time: " << ((double)indexing_time)/CLOCKS_PER_SEC << " seconds" << std::endl;
return 0;
}
// ==========================================================================
// mini_bowtie
// ==========================================================================
// Copyright (c) 2006-2012, Knut Reinert, FU Berlin
// All rights reserved.
#include <iostream>
#include <time.h>
#include <seqan/basic.h>
#include <seqan/sequence.h>
#include <seqan/file.h>
#include <seqan/index.h>
#include <seqan/index/index_fm.h>
#include <seqan/store.h>
using namespace seqan;
int main(int argc, char *argv[])
{
clock_t indexing_time = clock();
// type definitions
typedef Dna5String TString;
typedef StringSet<TString> TStringSet;
typedef Index<TStringSet, FMIndex<> > TIndex;
typedef Iterator<TIndex, TopDown<ParentLinks<> > >::Type TIter;
// reading the command line arguments
if (argc < 2) {
std::cerr << "Invalid number of arguments." << std::endl
<< "USAGE: fmIndex GENOME.fasta" << std::endl;
return 1;
}
// declaration and initialization of the fragment store
FragmentStore<> fragStore;
// loading the contigs
if (!loadContigs(fragStore, argv[1])) {
std::cerr << "Loading GENOME failed" << std::endl;
return 1;
}
// copy the contigs to the string set "mySeq"
TStringSet mySeq;
for (unsigned i = 0; i < length(fragStore.contigStore); ++i) {
appendValue(mySeq, fragStore.contigStore[i].seq);
std::cout << "Sequence " << i << " has " << length(mySeq[i]) << " nucleotides" << std::endl;
}
//std::cout << mySeq[0] << std::endl;
// clear contigs to free memory
clearContigs(fragStore);
////////////////////////////////////////////////////////////////////////
// REVERSE the ref genome !!! REVERSE !!! REVERSE !!!
reverse(mySeq);
////////////////////////////////////////////////////////////////////////
// FM INDEXING
TIndex fmIndex(mySeq); // the index is not created yet
TIter it(fmIndex); // the index is created after this line
char indexFile[100];
strcpy(indexFile, argv[1]);
strcat(indexFile, ".index");
if (save(fmIndex, indexFile)) {
std::cout << "Index saved!" << std::endl;
}
else {
std::cout << "Failed to save index" << std::endl;
return 1;
}
indexing_time = clock() - indexing_time;
std::cout << "Total indexing time: " << ((double)indexing_time)/CLOCKS_PER_SEC << " seconds" << std::endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment