830 - BasicMapper::main()

From Genome Analysis Wiki
Jump to navigationJump to search
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));
}