830 - BasicMapper::main()
From Genome Analysis Wiki
Jump to navigationJump to searchint 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));
}