Goncalo: Created page with " int main(int argc, char ** argv) { // Record start time for this run time_t start = time(NULL); // This function summarizes available options and proc..."

2013-10-28T17:29:35Z

Created page with "<source lang="cpp"> int main(int argc, char ** argv) { // Record start time for this run time_t start = time(NULL); // This function summarizes available options and proc..."

New page

<source lang="cpp">

int main(int argc, char ** argv)
{
// Record start time for this run
time_t start = time(NULL);

// This function summarizes available options and processes command line parameters
GreetUser(argc, argv);

// This function sets up the reference string, with bases ACGT encoded as 1234.
LoadReferenceGenome();

// This function prints a basic summary of the reference genome
DescribeReferenceGenome();

// Build reference genome indices
printf("Building reference indices...\n");

String mask1('1', indexComplexity), mask0('0', indexComplexity), mask;

MaqIndex index1(reference, mask = mask1 + mask1); printf(" Built index with mask %s\n", (const char *) mask);
MaqIndex index2(reference, mask = mask1 + mask0 + mask1); printf(" Built index with mask %s\n", (const char *) mask);
MaqIndex index3(reference, mask = mask1 + mask0 + mask0 + mask1); printf(" Built index with mask %s\n", (const char *) mask);
printf(" Done!\n\n");

// Map reads, one at a time
FastQ reads;
BestAlignment best;
String buffer;

int readCount = 0;
int mapCount = 0;
int mapMismatchesEq0 = 0;
int mapMismatchesEq1 = 0;
int progress = 0x100000;

reads.OpenArchive(shortReadFile);
IFILE sam = ifopen(samfile.IsEmpty() ? shortReadFile + ".sam" : samfile, "wb");

printf("Mapping reads\n");
while (!reads.EndOfFile())
{
reads.Load();
best.MaqAlignment(index1, index2, index3, reads);

if (reads.Length() == 0) continue;

if (++readCount == progress)
printf(" %d reads mapped so far ...\n", readCount), progress += 0x100000;

ReportAlignment(sam, reads, best);

if (best.mapQ > 30) mapCount++;
if (best.mismatches == 0) mapMismatchesEq0++;
if (best.mismatches == 1) mapMismatchesEq1++;
}

ifclose(sam);
reads.CloseArchive();

printf(" %d reads processed\n", readCount);
printf(" %d reads with mapping quality > 30\n", mapCount);
printf(" %d reads with zero mismatches\n", mapMismatchesEq0);
printf(" %d reads with one mismatches\n", mapMismatchesEq1);

// Log runtime and run date
time_t stop = time(NULL);
int seconds = (int) (stop - start);

printf("Run completed in %d hours, %d mins and %d seconds on %s\n",
seconds / 3600, (seconds % 3600) / 60, seconds % 60, ctime(&stop));
}

</source>

830 - BasicMapper::main() - Revision history

Goncalo: Created page with " int main(int argc, char ** argv) { // Record start time for this run time_t start = time(NULL); // This function summarizes available options and proc..."