Add a max_size parameter to the load command. Fixes #23

Author: Ben Fulton

cafe/cafe_commands.cpp

@@ -667,6 +667,7 @@ struct load_args get_load_arguments(vector<Argument> pargs)
 	args.num_threads = 0;
 	args.num_random_samples = 0;
 	args.pvalue = -1;
+	args.max_size = -1;
 
 	for (size_t i = 0; i < pargs.size(); i++)
 	{
@@ -678,6 +679,9 @@ struct load_args get_load_arguments(vector<Argument> pargs)
 		if (!strcmp(parg->opt, "-r"))
 			sscanf(parg->argv[0], "%d", &args.num_random_samples);
 
+		if (!strcmp(parg->opt, "-max_size"))
+			sscanf(parg->argv[0], "%d", &args.max_size);
+
 		if (!strcmp(parg->opt, "-p"))
 			sscanf(parg->argv[0], "%lf", &args.pvalue);
 
@@ -733,6 +737,8 @@ bool endsWith(std::string const &fullString, std::string const &ending) {
 \brief Loads families from a family file with a defined format
 *
 * Takes six arguments: -t, -r, -p, -l, -i, and -filter
+*
+* -filter ensures that there is at least one copy at the root (using parsimony) for each family. 
 */
 int cafe_cmd_load(Globals& globals, std::vector<std::string> tokens)
 {
@@ -764,7 +770,7 @@ int cafe_cmd_load(Globals& globals, std::vector<std::string> tokens)
 
 	param->str_fdata = string_new_with_string(args.family_file_name.c_str());
   ifstream ifst(args.family_file_name);
-	param->pfamily = load_gene_families(ifst, 1, separator);
+	param->pfamily = load_gene_families(ifst, separator, args.max_size);
 	if (param->pfamily == NULL) 
 		throw runtime_error("Failed to load file\n");
 
@@ -1948,7 +1954,7 @@ int cafe_cmd_accuracy(Globals& globals, std::vector<std::string> tokens)
 	{
 		// read in truth data
     ifstream ifst(truthfile);
-		pCafeFamily truthfamily = load_gene_families(ifst, 1, '\t');
+		pCafeFamily truthfamily = load_gene_families(ifst, '\t', -1);
 		if (truthfamily == NULL) {
 			fprintf(stderr, "failed to read in true values %s\n", truthfile.c_str());
 			return -1;
@@ -2100,7 +2106,7 @@ int cafe_cmd_cvfamily(Globals& globals, std::vector<std::string> tokens)
 
 		// read in training data
     ifstream ifst(trainfile);
-		pCafeFamily tmpfamily = load_gene_families(ifst, 1, '\t');
+		pCafeFamily tmpfamily = load_gene_families(ifst, '\t', -1);
 		if (tmpfamily == NULL) {
 			fprintf(stderr, "failed to read in training data %s\n", trainfile);
 			return -1;
@@ -2185,7 +2191,7 @@ int cafe_cmd_cvspecies(Globals& globals, std::vector<std::string> tokens)
 
 			// read in training data
       std::ifstream ifst(trainfile);
-			pCafeFamily tmpfamily = load_gene_families(ifst, 1, '\t');
+			pCafeFamily tmpfamily = load_gene_families(ifst, '\t', -1);
 			if (tmpfamily == NULL) {
 				fprintf(stderr, "failed to read in training data %s\n", trainfile);
 				fprintf(stderr, "did you load the family data with the cross-validation option (load -i <familyfile> -cv)?\n");

cafe/cafe_commands.h

@@ -72,6 +72,7 @@ void tree_set_branch_lengths(pCafeTree pcafe, std::vector<int> lengths);
 struct load_args {
 	int num_threads;
 	int num_random_samples;
+	int max_size;
 	double pvalue;
 	bool filter;
 	std::string log_file_name;

cafe/cross_validator.cpp

@@ -259,7 +259,7 @@ double cross_validator::validate_by_family(pCafeParam param, const char* queryfi
 
   // read in validation data
   std::ifstream ifst(truthfile);
-  pCafeFamily truthfamily = load_gene_families(ifst, 1, '\t');
+  pCafeFamily truthfamily = load_gene_families(ifst, '\t', -1);
   if (truthfamily == NULL) {
     fprintf(stderr, "failed to read in true values %s\n", truthfile);
     return -1;

cafe/gene_family.cpp

@@ -99,7 +99,12 @@ void cafe_family_free(pCafeFamily pcf)
   pcf = NULL;
 }
 
-pCafeFamily load_gene_families(std::istream& ist, int bpatcheck, char separator)
+int to_int(string s)
+{
+	return stoi(s);
+}
+
+pCafeFamily load_gene_families(std::istream& ist, char separator, int max_size)
 {
   char buf[STRING_BUF_SIZE];
   if (!ist)
@@ -122,18 +127,21 @@ pCafeFamily load_gene_families(std::istream& ist, int bpatcheck, char separator)
     if (!ist)
       break;
 
-    cafe_family_add_item(pcf, string_split(buf, separator));
+	vector<string> values = string_split(buf, separator);
+	vector<int> sizes(values.size() - 2);
+	transform(values.begin() + 2, values.end(), sizes.begin(), to_int);
+	if (max_size < 0 || *max_element(sizes.begin(), sizes.end()) <= max_size)
+		cafe_family_add_item(pcf, values[1], values[0], sizes);
   }
-  if (bpatcheck)
-  {
+
     __cafe_famliy_check_the_pattern(pcf);
-  }
+
   return pcf;
 }
 
 /// Data array is assumes to contain a description, an identifier, and a set of integers
 /// giving the family size in species order
-void cafe_family_add_item(pCafeFamily pcf, const vector<string>& data)
+void cafe_family_add_item(pCafeFamily pcf, std::string id, std::string desc, std::vector<int> sizes)
 {
   pCafeFamilyItem pitem = (pCafeFamilyItem)memory_new(1, sizeof(CafeFamilyItem));
   pitem->count = (int*)calloc(pcf->num_species, sizeof(int));
@@ -143,27 +151,18 @@ void cafe_family_add_item(pCafeFamily pcf, const vector<string>& data)
   pitem->mu = NULL;
   pitem->holder = 1;
 
-  if (data.size() != size_t(pcf->num_species + 2))
+  if (sizes.size() != size_t(pcf->num_species))
   {
-    std::cerr << "Inconsistency in column count: expected " << pcf->num_species + 2 << ", but found " << data.size();
+    std::cerr << "Inconsistency in column count: expected " << pcf->num_species + 2 << ", but found " << sizes.size() + 2;
   }
 
-  pitem->desc = new char[data[0].size()+1];
-  strcpy(pitem->desc, data[0].c_str());
-  pitem->id = new char[data[1].size() + 1];
-  strcpy(pitem->id, data[1].c_str());
+  pitem->desc = new char[desc.size()+1];
+  strcpy(pitem->desc, desc.c_str());
+  pitem->id = new char[id.size() + 1];
+  strcpy(pitem->id, id.c_str());
+
+  copy(sizes.begin(), sizes.end(), pitem->count);
 
-  for (int j = 0; j < pcf->num_species; j++)
-  {
-    if (data[j + 2].empty())
-    {
-      pitem->count[j] = 0;
-    }
-    else
-    {
-      pitem->count[j] = atoi(data[j + 2].c_str());
-    }
-  }
   pcf->max_size = max(pcf->max_size, *std::max_element(pitem->count, pitem->count + pcf->num_species));
   arraylist_add(pcf->flist, pitem);
 }

cafe/gene_family.h

@@ -16,10 +16,10 @@ const int COMMA_AS_WHITESPACE = 1;
 std::vector<std::string> tokenize(std::string s, int flags);
 
 pCafeFamily cafe_family_init(const std::vector<std::string>& species_list);
-pCafeFamily load_gene_families(std::istream& ist, int bpatcheck, char separator);
+pCafeFamily load_gene_families(std::istream& ist, char separator, int max_size);
 void cafe_family_free(pCafeFamily pcf);
 
-void cafe_family_add_item(pCafeFamily pcf, const std::vector<std::string>& data);
+void cafe_family_add_item(pCafeFamily pcf, std::string id, std::string desc, std::vector<int> sizes);
 void cafe_family_item_free(pCafeFamilyItem pitem);
 
 int log_cluster_membership(pCafeFamily pcf, int k_value, double **p_z_membership, std::ostream& log);

tests/command_tests.cpp

@@ -189,12 +189,21 @@ TEST(CommandTests, get_load_arguments)
 	struct load_args args = get_load_arguments(build_argument_list(command));
 	LONGS_EQUAL(1, args.num_threads);
 	LONGS_EQUAL(2, args.num_random_samples);
+	LONGS_EQUAL(-1, args.max_size);
 	DOUBLES_EQUAL(.05, args.pvalue, .000001);
 	STRCMP_EQUAL("log.txt", args.log_file_name.c_str());
 	STRCMP_EQUAL("fam.txt", args.family_file_name.c_str());
 	CHECK(!args.filter);
 }
 
+TEST(CommandTests, get_load_arguments_filter_and_max_size)
+{
+	vector<string> command = tokenize("load -filter -max_size 100", REGULAR_WHITESPACE);
+	struct load_args args = get_load_arguments(build_argument_list(command));
+	LONGS_EQUAL(100, args.max_size);
+	CHECK(args.filter);
+}
+
 TEST(CommandTests, cafe_cmd_load)
 {
 	try
@@ -261,7 +270,7 @@ TEST(CommandTests, cafe_cmd_tree_syncs_family_if_loaded)
 
 	std::vector<std::string> species = { "chimp", "human", "mouse", "rat", "dog" };
 	globals.param.pfamily = cafe_family_init(species);
-	cafe_family_add_item(globals.param.pfamily, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(globals.param.pfamily, "id", "description", {3, 5, 7, 11, 13 });
 
 	LONGS_EQUAL(-1, globals.param.pfamily->index[0]);
 	cafe_cmd_tree(globals, tokens);
@@ -297,7 +306,7 @@ void prepare_viterbi(CafeParam& param)
 	param.lambda = lambdas;
 
 	param.pfamily = cafe_family_init({ "chimp", "human", "mouse", "rat", "dog" });
-	cafe_family_add_item(param.pfamily, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(param.pfamily, "id", "description", { 3, 5, 7, 11, 13 });
 
 	cafe_family_set_species_index(param.pfamily, param.pcafe);
 

tests/family_tests.cpp

@@ -40,11 +40,11 @@ TEST_GROUP(FamilyTests)
 TEST(FamilyTests, TestCafeFamily)
 {
 	std::ifstream ifst("Nonexistent.tab");
-	POINTERS_EQUAL(NULL, load_gene_families(ifst, 1, '\t'));
+	POINTERS_EQUAL(NULL, load_gene_families(ifst, '\t', -1));
 
   pCafeFamily fam = cafe_family_init({ "Dog", "Chimp", "Human", "Mouse", "Rat" });
 
-  cafe_family_add_item(fam, { "OTOPETRIN", "ENSF1", "3", "5", "7", "11", "13" });
+  cafe_family_add_item(fam, "ENSF1", "OTOPETRIN", { 3, 5, 7, 11, 13 });
 
 	LONGS_EQUAL(5, fam->num_species);
 	LONGS_EQUAL(1, fam->flist->size);
@@ -65,10 +65,10 @@ TEST(FamilyTests, TestCafeFamily)
 TEST(FamilyTests, load_gene_families_tab_separated)
 {
   std::ifstream ifst("Nonexistent.tab");
-  POINTERS_EQUAL(NULL, load_gene_families(ifst, 1, '\t'));
+  POINTERS_EQUAL(NULL, load_gene_families(ifst, '\t', -1));
 
   std::istringstream ist("FAMILYDESC\tFAMILY\tA\tB\nOTOPETRIN\tENSF2\t2\t10\n");
-  pCafeFamily fam = load_gene_families(ist, 0, '\t');
+  pCafeFamily fam = load_gene_families(ist, '\t', -1);
   LONGS_EQUAL(2, fam->num_species);
   STRCMP_EQUAL("A", fam->species[0]);
   STRCMP_EQUAL("B", fam->species[1]);
@@ -83,10 +83,25 @@ TEST(FamilyTests, load_gene_families_tab_separated)
   cafe_family_free(fam);
 }
 
+TEST(FamilyTests, filter_families_greater_than_max_size)
+{
+	std::istringstream ist("FAMILYDESC\tFAMILY\tA\tB\ndesc\tid1\t2\t50\ndesc\tid2\t20\t10\n");
+	pCafeFamily fam = load_gene_families(ist, '\t', 25);
+	LONGS_EQUAL(2, fam->num_species);
+	STRCMP_EQUAL("A", fam->species[0]);
+	STRCMP_EQUAL("B", fam->species[1]);
+
+	LONGS_EQUAL(1, fam->flist->size);
+
+	CHECK(get_family_item(fam, "id1") == NULL);
+	CHECK(get_family_item(fam, "id2") != NULL);
+	cafe_family_free(fam);
+}
+
 TEST(FamilyTests, load_gene_families_comma_separated)
 {
   std::istringstream ist("FAMILY DESC,FAMILY,A,B\nOXYSTEROL BINDING RELATED,ENSF2,2,10\n");
-  pCafeFamily fam = load_gene_families(ist, 0, ',');
+  pCafeFamily fam = load_gene_families(ist, ',', -1);
   LONGS_EQUAL(2, fam->num_species);
   STRCMP_EQUAL("A", fam->species[0]);
   STRCMP_EQUAL("B", fam->species[1]);
@@ -106,7 +121,7 @@ TEST(FamilyTests, load_gene_families_is_not_whitespace_separated)
   std::istringstream ist("FAMILYDESC FAMILY A B C\n");
   try
   {
-    pCafeFamily fam = load_gene_families(ist, 0, '\t');
+    pCafeFamily fam = load_gene_families(ist, '\t', -1);
     FAIL("Exception should have been thrown");
     cafe_family_free(fam);
   }
@@ -119,7 +134,7 @@ TEST(FamilyTests, load_gene_families_is_not_whitespace_separated)
 TEST(FamilyTests, cafe_family_set_size)
 {
 	pCafeFamily pcf = cafe_family_init({"chimp", "human", "mouse", "rat", "dog" });
-	cafe_family_add_item(pcf, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
 
 	pCafeTree ptree = create_tree();
 	cafe_family_set_species_index(pcf, ptree);
@@ -155,7 +170,7 @@ TEST(FamilyTests, cafe_family_add_item)
 {
 	pCafeFamily pcf = cafe_family_init({"chimp", "human", "mouse", "rat", "dog" });
 
-	cafe_family_add_item(pcf, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
 	LONGS_EQUAL(1, pcf->flist->size);
 	pCafeFamilyItem pitem = (pCafeFamilyItem)arraylist_get(pcf->flist, 0);
 	STRCMP_EQUAL("description", pitem->desc);
@@ -203,7 +218,7 @@ TEST(FamilyTests, cafe_family_set_size_with_family_forced)
 	pCafeFamily pcf = cafe_family_init({"chimp", "human", "mouse", "rat", "dog" });
 	pCafeTree cafe_tree = create_tree();
 	cafe_family_set_species_index(pcf, cafe_tree);
-	cafe_family_add_item(pcf, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
 
 	// SUT
 	cafe_family_set_size_with_family_forced(pcf, 0, cafe_tree);
@@ -218,7 +233,7 @@ TEST(FamilyTests, write_family_gainloss)
 {
 	std::ostringstream ost;
 	pCafeFamily pcf = cafe_family_init({"chimp", "human", "mouse", "rat", "dog" });
-	cafe_family_add_item(pcf, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
 	pCafeTree tree1 = create_tree();
 	cafe_family_set_species_index(pcf, tree1);
 	pCafeTree tree2 = cafe_tree_copy(tree1);
@@ -237,11 +252,11 @@ TEST(FamilyTests, write_family)
 {
 	std::ostringstream ost;
 	pCafeFamily pcf = cafe_family_init({"chimp", "human", "mouse", "rat", "dog" });
-	cafe_family_add_item(pcf, { "my_description", "my_id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
 
 	write_family(ost, pcf);
 	STRCMP_CONTAINS("Desc\tFamily ID\tchimp\thuman\tmouse\trat\tdog\n", ost.str().c_str())
-	STRCMP_CONTAINS("my_description\tmy_id\t3\t5\t7\t11\t13\n", ost.str().c_str())
+	STRCMP_CONTAINS("description\tid\t3\t5\t7\t11\t13\n", ost.str().c_str())
   cafe_family_free(pcf);
 
 }
@@ -250,13 +265,13 @@ TEST(FamilyTests, log_cluster_membership)
 {
   std::ostringstream ost;
   pCafeFamily pcf = cafe_family_init({ "chimp", "human", "mouse", "rat", "dog" });
-  cafe_family_add_item(pcf, { "my_description", "my_id", "3", "5", "7", "11", "13" });
+  cafe_family_add_item(pcf, "id", "description", { 3, 5, 7, 11, 13 });
   double **z = (double **)memory_new_2dim(2, 2, sizeof(double));
   z[0][0] = 7.7;
   z[0][1] = 11.13;
   log_cluster_membership(pcf, 2, z, ost);
   STRCMP_CONTAINS("The Number of families : 1\n", ost.str().c_str());
-  STRCMP_CONTAINS("family my_id:  7.7 11.13", ost.str().c_str());
+  STRCMP_CONTAINS("family id:  7.7 11.13", ost.str().c_str());
 
   cafe_family_free(pcf);
 }

tests/lambda_tests.cpp

@@ -103,7 +103,7 @@ TEST(LambdaTests, TestCmdLambda)
 
   globals.param.pfamily = cafe_family_init({"Dog", "Chimp", "Human", "Mouse", "Rat" });
 
-  cafe_family_add_item(globals.param.pfamily, { "OTOPETRIN", "ENSF00000002390", "7", "7", "13", "18", "7" });
+  cafe_family_add_item(globals.param.pfamily, "ENSF00000002390", "OTOPETRIN", { 7, 7, 13, 18, 7 });
 
 	init_family_size(&globals.param.family_size, globals.param.pfamily->max_size);
 	cafe_tree_set_parameters(globals.param.pcafe, &globals.param.family_size, 0);
@@ -573,10 +573,10 @@ TEST(LambdaTests, best_lambda_by_fminsearch)
 	globals.param.pcafe = cafe_tree_new(tree, &range, 0, 0);
 	globals.param.pfamily = cafe_family_init({ "A", "B", "C", "D" });
 
-	cafe_family_add_item(globals.param.pfamily, { "1", "ENS01", "5", "10", "2", "6" });
-	cafe_family_add_item(globals.param.pfamily, { "2", "ENS02", "5", "10", "2", "6" });
-	cafe_family_add_item(globals.param.pfamily, { "3", "ENS03", "5", "10", "2", "6" });
-	cafe_family_add_item(globals.param.pfamily, { "4", "ENS04", "5", "10", "2", "6" });
+	cafe_family_add_item(globals.param.pfamily, "ENS01", "1", { 5, 10, 2, 6 });
+	cafe_family_add_item(globals.param.pfamily, "ENS02", "2", { 5, 10, 2, 6 });
+	cafe_family_add_item(globals.param.pfamily, "ENS03", "3", { 5, 10, 2, 6 });
+	cafe_family_add_item(globals.param.pfamily, "ENS04", "4", { 5, 10, 2, 6 });
 
 	init_family_size(&globals.param.family_size, globals.param.pfamily->max_size);
 	cafe_tree_set_parameters(globals.param.pcafe, &globals.param.family_size, 0);
@@ -638,7 +638,7 @@ TEST(LambdaTests, get_posterior)
 
 	param.pfamily = cafe_family_init({ "chimp", "human", "mouse", "rat", "dog" });
 	cafe_family_set_species_index(param.pfamily, param.pcafe);
-	cafe_family_add_item(param.pfamily, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(param.pfamily, "id", "description", { 3, 5, 7, 11, 13 });
 
 	param.ML = (double*)memory_new(15, sizeof(double));
 	param.MAP = (double*)memory_new(15, sizeof(double));

tests/report_tests.cpp

@@ -206,7 +206,7 @@ TEST(ReportTests, write_families_line)
   globals.param.pcafe = tree;
   globals.param.pfamily = cafe_family_init({ "chimp", "human", "mouse", "rat", "dog" });
   cafe_family_set_species_index(globals.param.pfamily, tree);
-  cafe_family_add_item(globals.param.pfamily, { "description", "id", "3", "5", "7", "11", "13" });
+  cafe_family_add_item(globals.param.pfamily, "id", "description", { 3, 5, 7, 11, 13 });
 
   viterbi_parameters* v = globals.viterbi;
   v->num_nodes = 6;
@@ -250,7 +250,7 @@ TEST(ReportTests, family_line_item_sets_pvalues)
 	pCafeTree tree = create_tree(range);
 	globals.param.pcafe = tree;
 	cafe_family_set_species_index(globals.param.pfamily, tree);
-	cafe_family_add_item(globals.param.pfamily, { "description", "id", "3", "5", "7", "11", "13" });
+	cafe_family_add_item(globals.param.pfamily, "id", "description", { 3, 5, 7, 11, 13 });
 
 	viterbi_parameters* v = globals.viterbi;
 	v->num_nodes = 6;