Pick and extract different particle types separately

scripts/trace_amyloids.py

@@ -759,7 +759,7 @@ def run_compiled_model(data):
                 df = pd.DataFrame({
                     'rlnCoordinateX': sampled_points[:, 0],
                     'rlnCoordinateY': sampled_points[:, 1],
-                    'rlnParticleSelectionType':  sampled_points[:, 2]
+                    'rlnHelicalLineId':  sampled_points[:, 2]
                 })
                 starfile.write(df, outfile, overwrite=True)
             else:

src/displayer.cpp

@@ -636,6 +636,7 @@ int basisViewerWindow::fillPickerViewerCanvas(MultidimArray<RFLOAT> image, Multi
 											  RFLOAT _minimum_pick_fom, RFLOAT _min_fom, RFLOAT _max_fom)
 {
 	current_selection_type = 2; // Green
+	current_line_id = 1;
 
 	// Scroll bars
 	Fl_Scroll scroll(0, 0, w(), h());
@@ -2080,6 +2081,8 @@ void pickerViewerCanvas::draw()
 	long int icoord = 0;
     int my_prev_type = 0;
 	int xcoori_start, ycoori_start;
+	int current_line_id;
+	int prev_line_id = 0;
 	FOR_ALL_OBJECTS_IN_METADATA_TABLE(MDcoords)
 	{
 		icoord++;
@@ -2097,6 +2100,9 @@ void pickerViewerCanvas::draw()
 			if (fom < minimum_pick_fom) continue;
 		}
 
+		if (MDcoords.containsLabel(EMDL_PARTICLE_HELICAL_LINE_ID))
+			MDcoords.getValue(EMDL_PARTICLE_HELICAL_LINE_ID, current_line_id);
+
         if (color_label != EMDL_UNDEFINED)
 		{
 			RFLOAT colval;
@@ -2160,7 +2166,7 @@ void pickerViewerCanvas::draw()
 
         if (do_lines)
         {
-            if (mytype == my_prev_type)
+            if (current_line_id == prev_line_id)
             {
                 fl_line(xcoori_start, ycoori_start, xcoori, ycoori);
             }
@@ -2184,7 +2190,8 @@ void pickerViewerCanvas::draw()
                 }
             }
         }
-        my_prev_type = mytype;
+        //my_prev_type = mytype;
+		prev_line_id = current_line_id;
     }
 }
 
@@ -2196,7 +2203,7 @@ int pickerViewerCanvas::handle(int ev)
 	const int key = Fl::event_key();
     if (do_lines && has_dragged && ev == FL_RELEASE)
     {
-        current_selection_type++;
+        current_line_id++;
     }
     has_dragged = false;
 
@@ -2251,10 +2258,13 @@ int pickerViewerCanvas::handle(int ev)
 			MDcoords.setValue(EMDL_IMAGE_COORD_X, xcoor);
 			MDcoords.setValue(EMDL_IMAGE_COORD_Y, ycoor);
 			// No autopicking, but still always fill in the parameters for autopicking with dummy values (to prevent problems in joining autopicked and manually picked coordinates)
-			MDcoords.setValue(EMDL_PARTICLE_SELECTION_TYPE, iaux);
+			MDcoords.setValue(EMDL_PARTICLE_SELECTION_TYPE, current_selection_type);
 			MDcoords.setValue(EMDL_ORIENT_PSI, aux);
 			MDcoords.setValue(EMDL_PARTICLE_AUTOPICK_FOM, zero);
 
+			if (do_lines)
+				MDcoords.setValue(EMDL_PARTICLE_HELICAL_LINE_ID, current_line_id);
+
 			redraw();
             return 1;
 		}
@@ -2275,13 +2285,13 @@ int pickerViewerCanvas::handle(int ev)
 					if (do_lines)
                     {
                         int delval;
-                        MDcoords.getValue(EMDL_PARTICLE_SELECTION_TYPE, delval);
+                        MDcoords.getValue(EMDL_PARTICLE_HELICAL_LINE_ID, delval);
                         MetaDataTable MDout;
                         MDout.setName(MDcoords.getName());
                         FOR_ALL_OBJECTS_IN_METADATA_TABLE(MDcoords)
                         {
                             int val;
-                            MDcoords.getValue(EMDL_PARTICLE_SELECTION_TYPE, val);
+                            MDcoords.getValue(EMDL_PARTICLE_HELICAL_LINE_ID, val);
 
                             if (val != delval) MDout.addObject(MDcoords.getObject(current_object));
                         }
@@ -2450,10 +2460,10 @@ void pickerViewerCanvas::loadCoordinates(bool ask_filename)
         FOR_ALL_OBJECTS_IN_METADATA_TABLE(MDcoords)
         {
             int ifil;
-            MDcoords.getValue(EMDL_PARTICLE_SELECTION_TYPE, ifil);
+            MDcoords.getValue(EMDL_PARTICLE_HELICAL_LINE_ID, ifil);
             if (ifil > ifil_max) ifil_max=ifil;
         }
-        current_selection_type = ifil_max+1;
+        current_line_id = ifil_max+1;
     }
 
 	if (fn_color != "")

src/displayer.h

@@ -73,6 +73,7 @@ static int predrag_yc;
 static bool has_shift;
 static int preshift_ipos;
 static int current_selection_type;
+static int current_line_id;
 static int colour_scheme;
 static int fom_is_grey_instead;
 

src/gui_jobwindow.cpp

@@ -1076,6 +1076,9 @@ Pixels values higher than this many times the image stddev will be replaced with
 	group7->end();
 	guientries["do_fom_threshold"].cb_menu_i();
 
+	current_y += STEPY/2;
+	place("selection_type", TOGGLE_DEACTIVATE);
+
 	tab2->end();
 	tab3->begin();
 	tab3->label("Helix");

src/helix.cpp

@@ -1828,6 +1828,7 @@ void convertHelicalTubeCoordsToMetaDataTable(
 	MetaDataTable MD_in;
 	std::vector<RFLOAT> x1_coord_list, y1_coord_list, x2_coord_list, y2_coord_list, pitch_list;
 	std::vector<int> tube_id_list;
+	std::vector<int> selection_type_list;
 
 	// Check parameters and open files
 	if ( (nr_asu < 1) || (rise_A < 0.001) || (pixel_size_A < 0.01) )
@@ -1863,20 +1864,29 @@ void convertHelicalTubeCoordsToMetaDataTable(
     // Sjors added MDin_has_id and MDin_has_pitch to allow manual calculation of different cross-over distances to be carried onto the extracted segments...
     bool MDin_has_id = MD_in.containsLabel(EMDL_PARTICLE_HELICAL_TUBE_ID);
     bool MDin_has_pitch =  MD_in.containsLabel(EMDL_PARTICLE_HELICAL_TUBE_PITCH);
+	bool MDin_has_selection_type = MD_in.containsLabel(EMDL_PARTICLE_SELECTION_TYPE);
     x1_coord_list.clear();
     y1_coord_list.clear();
     x2_coord_list.clear();
     y2_coord_list.clear();
     tube_id_list.clear();
     pitch_list.clear();
+	selection_type_list.clear();
     MDobj_id = 0;
+	int selection_type = 0;
     FOR_ALL_OBJECTS_IN_METADATA_TABLE(MD_in)
     {
     	MDobj_id++;
     	MD_in.getValue(EMDL_IMAGE_COORD_X, xp);
     	MD_in.getValue(EMDL_IMAGE_COORD_Y, yp);
     	if (MDobj_id % 2)
     	{
+    		if (MDin_has_selection_type)
+    		{
+    			MD_in.getValue(EMDL_PARTICLE_SELECTION_TYPE, selection_type);
+    			selection_type_list.push_back(selection_type);
+    		}
+
     		x1_coord_list.push_back(xp);
     		y1_coord_list.push_back(yp);
     		if (MDin_has_id)
@@ -1917,6 +1927,8 @@ void convertHelicalTubeCoordsToMetaDataTable(
     	MD_out.addLabel(EMDL_PARTICLE_HELICAL_TUBE_ID);
     if (MDin_has_pitch)
     	MD_out.addLabel(EMDL_PARTICLE_HELICAL_TUBE_PITCH);
+	if (MDin_has_selection_type)
+		MD_out.addLabel(EMDL_PARTICLE_SELECTION_TYPE);
 
 	// Calculate all coordinates for helical segments
 	nr_segments = 0;
@@ -1937,6 +1949,10 @@ void convertHelicalTubeCoordsToMetaDataTable(
 		dx = step_pix * cos(psi_rad);
 		dy = step_pix * sin(psi_rad);
 
+    	int current_selection_type = 0;
+    	if (MDin_has_selection_type)
+    		current_selection_type = selection_type_list[tube_id];
+
     	if (!cut_into_segments)
     	{
             MD_out.addObject();
@@ -1952,6 +1968,8 @@ void convertHelicalTubeCoordsToMetaDataTable(
 				MD_out.setValue(EMDL_PARTICLE_HELICAL_TUBE_ID, id);
 			if (MDin_has_pitch)
 				MD_out.setValue(EMDL_PARTICLE_HELICAL_TUBE_PITCH, pitch);
+    		if (MDin_has_selection_type)
+    			MD_out.setValue(EMDL_PARTICLE_SELECTION_TYPE, current_selection_type);
 
 	        nr_segments++;
     		continue;
@@ -1997,6 +2015,8 @@ void convertHelicalTubeCoordsToMetaDataTable(
     				MD_out.setValue(EMDL_PARTICLE_HELICAL_TUBE_ID, id);
     			if (MDin_has_pitch)
     				MD_out.setValue(EMDL_PARTICLE_HELICAL_TUBE_PITCH, pitch);
+    			if (MDin_has_selection_type)
+    				MD_out.setValue(EMDL_PARTICLE_SELECTION_TYPE, current_selection_type);
    			nr_segments++;
     		}
     	}

src/metadata_label.h

@@ -498,6 +498,7 @@ enum EMDLabel
 	EMDL_PARTICLE_HELICAL_TUBE_PITCH,
 	EMDL_PARTICLE_HELICAL_TRACK_LENGTH, //deprecated
 	EMDL_PARTICLE_HELICAL_TRACK_LENGTH_ANGSTROM,
+	EMDL_PARTICLE_HELICAL_LINE_ID,
 	EMDL_PARTICLE_SELECTION_TYPE,
 	EMDL_PARTICLE_CLASS,
 	EMDL_PARTICLE_DLL,
@@ -1237,6 +1238,7 @@ class StaticInitialization
 		EMDL::addLabel(EMDL_PARTICLE_HELICAL_TUBE_PITCH, EMDL_DOUBLE, "rlnHelicalTubePitch", "Cross-over distance for a helical segment (A)");
 		EMDL::addLabel(EMDL_PARTICLE_HELICAL_TRACK_LENGTH, EMDL_DOUBLE, "rlnHelicalTrackLength", "Distance (in pix) from the position of this helical segment to the starting point of the tube");
 		EMDL::addLabel(EMDL_PARTICLE_HELICAL_TRACK_LENGTH_ANGSTROM, EMDL_DOUBLE, "rlnHelicalTrackLengthAngst", "Distance (in A) from the position of this helical segment to the starting point of the tube");
+		EMDL::addLabel(EMDL_PARTICLE_HELICAL_LINE_ID, EMDL_INT, "rlnHelicalLineId", "Helical line ID for a manually selected helical line");
 		EMDL::addLabel(EMDL_PARTICLE_CLASS, EMDL_INT, "rlnClassNumber", "Class number for which a particle has its highest probability");
 		EMDL::addLabel(EMDL_PARTICLE_SELECTION_TYPE, EMDL_INT, "rlnParticleSelectionType", "Selection type for manually picked particles");
 		EMDL::addLabel(EMDL_PARTICLE_DLL, EMDL_DOUBLE, "rlnLogLikeliContribution", "Contribution of a particle to the log-likelihood target function");

src/pipeline_jobs.cpp

@@ -2555,6 +2555,7 @@ Pixels values higher than this many times the image stddev will be replaced with
 	joboptions["rescale"] = JobOption("Re-scaled size (pixels): ", 128, 64, 512, 8, "The re-scaled value needs to be an even number");
 	joboptions["do_fom_threshold"] = JobOption("Use autopick FOM threshold?", false, "If set to Yes, only particles with rlnAutopickFigureOfMerit values below the threshold below will be extracted.");
 	joboptions["minimum_pick_fom"] = JobOption("Minimum autopick FOM: ", 0, -5, 10, 0.1, "The minimum value for the rlnAutopickFigureOfMerit for particles to be extracted.");
+	joboptions["selection_type"] = JobOption("Selected particles: ", job_extract_selection_type_options, 0, "Extract particles only of this selection type. Default: extract all particles regardless of selection type.");
 
 	joboptions["do_extract_helix"] = JobOption("Extract helical segments?", false, "Set to Yes if you want to extract helical segments. RELION (.star), EMAN2 (.box) and XIMDISP (.coords) formats of tube or segment coordinates are supported.");
 	joboptions["helical_tube_outer_diameter"] = JobOption("Tube diameter (A): ", 200, 100, 1000, 10, "Outer diameter (in Angstroms) of helical tubes. \
@@ -2673,6 +2674,12 @@ bool RelionJob::getCommandsExtractJob(std::string &outputname, std::vector<std::
 		command += " --minimum_pick_fom " + joboptions["minimum_pick_fom"].getString();
 	}
 
+	std::string selection_type = joboptions["selection_type"].getString();
+	if (selection_type != "All particles")
+	{
+		command += " --selection_type " + std::string(1, selection_type[0]);
+	}
+
 	if (joboptions["do_float16"].getBoolean())
 	{
 		command += " --float16 ";

src/pipeline_jobs.h

@@ -170,6 +170,16 @@ static const std::vector<std::string> job_blush_version_options{
         "amy-v1.0"
 };
 
+static const std::vector<std::string> job_extract_selection_type_options{
+	     "All particles",
+		 "1 (red)",
+	     "2 (green)",
+	     "3 (blue)",
+		 "4 (cyan)",
+		 "5 (magenta)",
+		 "6 (yellow)"
+};
+
 // To have a line on the GUI to change the minimum number of dedicated in a job
 static bool do_allow_change_minimum_dedicated;
 

src/preprocessing.cpp

@@ -661,7 +661,7 @@ void Preprocessing::readCoordinates(FileName fn_coord, MetaDataTable &MD)
 }
 
 void Preprocessing::addOneHelicalSegment(MetaDataTable &MD, RFLOAT xcoord, RFLOAT ycoord, int tube_id,
-                                  RFLOAT psi_prior, RFLOAT helix_length, RFLOAT psi_prior_flip_ratio)
+                                  RFLOAT psi_prior, RFLOAT helix_length, RFLOAT psi_prior_flip_ratio, int selection_type)
 {
     MD.addObject();
     MD.setValue(EMDL_IMAGE_COORD_X, xcoord);
@@ -671,7 +671,8 @@ void Preprocessing::addOneHelicalSegment(MetaDataTable &MD, RFLOAT xcoord, RFLOA
     MD.setValue(EMDL_ORIENT_PSI_PRIOR, psi_prior);
     MD.setValue(EMDL_PARTICLE_HELICAL_TRACK_LENGTH_ANGSTROM, helix_length);
     MD.setValue(EMDL_ORIENT_PSI_PRIOR_FLIP_RATIO, psi_prior_flip_ratio);
-
+	if (selection_type > 0)
+		MD.setValue(EMDL_PARTICLE_SELECTION_TYPE, selection_type);
 }
 
 void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
@@ -695,8 +696,8 @@ void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
         return;
 
     if ( (!MD_in.containsLabel(EMDL_IMAGE_COORD_X)) || (!MD_in.containsLabel(EMDL_IMAGE_COORD_Y)) ||
-            (!MD_in.containsLabel(EMDL_PARTICLE_SELECTION_TYPE)) )
-        REPORT_ERROR("Preprocessing::convertHelicalLineCoordsToMetaDataTable: Input STAR file does not contain X, Y coordinates or particle selection type! Error(s) in " + fn_in);
+            (!MD_in.containsLabel(EMDL_PARTICLE_HELICAL_LINE_ID)) )
+        REPORT_ERROR("Preprocessing::convertHelicalLineCoordsToMetaDataTable: Input STAR file does not contain X, Y coordinates or helical line ID! Error(s) in " + fn_in);
 
     MD_out.clear();
     MD_out.addLabel(EMDL_IMAGE_COORD_X);
@@ -706,6 +707,10 @@ void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
     MD_out.addLabel(EMDL_ORIENT_PSI_PRIOR);
     MD_out.addLabel(EMDL_PARTICLE_HELICAL_TRACK_LENGTH_ANGSTROM);
     MD_out.addLabel(EMDL_ORIENT_PSI_PRIOR_FLIP_RATIO);
+
+	if (MD_in.containsLabel(EMDL_PARTICLE_SELECTION_TYPE))
+		MD_out.addLabel(EMDL_PARTICLE_SELECTION_TYPE);
+
     RFLOAT psi_prior_flip_ratio = (bimodal_angular_priors) ? BIMODAL_PSI_PRIOR_FLIP_RATIO : UNIMODAL_PSI_PRIOR_FLIP_RATIO;
 
     int filament_id=0;
@@ -715,15 +720,19 @@ void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
     RFLOAT step_pix = step_A / pixel_size_A;
 
     total_tubes = total_segments = 0;
-    int my_type=0, my_prev_type= 0;
+    int current_line_id=0, prev_line_id= 0;
+	int selection_type = 0;
     bool is_first = false;
     RFLOAT my_xcoord, my_ycoord, prev_xcoord, prev_ycoord;
     FOR_ALL_OBJECTS_IN_METADATA_TABLE(MD_in)
     {
         MD_in.getValue(EMDL_IMAGE_COORD_X, my_xcoord);
         MD_in.getValue(EMDL_IMAGE_COORD_Y, my_ycoord);
-        MD_in.getValue(EMDL_PARTICLE_SELECTION_TYPE, my_type);
-        if (my_type != my_prev_type)
+        MD_in.getValue(EMDL_PARTICLE_HELICAL_LINE_ID, current_line_id);
+    	if (MD_in.containsLabel(EMDL_PARTICLE_SELECTION_TYPE))
+    		MD_in.getValue(EMDL_PARTICLE_SELECTION_TYPE, selection_type);
+
+        if (current_line_id != prev_line_id)
         {
             if (is_first)
             {
@@ -740,7 +749,7 @@ void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
             total_segments++;
             //std::cerr << "is_first= "<<is_first << " x,y= " << my_xcoord << " , " << my_ycoord<< " id= " << filament_id<< " l= "<<filament_length<< std::endl;
             addOneHelicalSegment(MD_out, my_xcoord, my_ycoord,
-                                 filament_id, 0., filament_length, psi_prior_flip_ratio); // this is wrong tilt prior for tilted data!!! (we might never use this anyway...)
+                                 filament_id, 0., filament_length, psi_prior_flip_ratio, selection_type); // this is wrong tilt prior for tilted data!!! (we might never use this anyway...)
         }
         else
         {
@@ -759,15 +768,15 @@ void Preprocessing::convertHelicalLineCoordsToMetaDataTable(
                 if (is_first) MD_out.setValue(EMDL_ORIENT_PSI_PRIOR, psi_prior, MD_out.numberOfObjects()-1);
                 //std::cerr << "added="<<added << " newx,y= " << new_x << " , " << new_y<< " psi= " << psi_prior << " id= " << filament_id<< " l= "<<filament_length<< std::endl;
                 addOneHelicalSegment(MD_out, new_x, new_y,filament_id,
-                                     psi_prior, filament_length, psi_prior_flip_ratio);
+                                     psi_prior, filament_length, psi_prior_flip_ratio, selection_type);
             }
             remaining_step = added - dist;
             is_first = false;
         }
 
         prev_xcoord = my_xcoord;
         prev_ycoord = my_ycoord;
-        my_prev_type = my_type;
+        prev_line_id = current_line_id;
     }
 
 }

src/preprocessing.h

@@ -215,7 +215,7 @@ class Preprocessing
 	void readCoordinates(FileName fn_coord, MetaDataTable &MD);
 
     void addOneHelicalSegment(MetaDataTable &MD, RFLOAT xcoord, RFLOAT ycoord, int tube_id,
-                       RFLOAT psi_prior, RFLOAT helix_length, RFLOAT psi_prior_flip_ratio);
+                       RFLOAT psi_prior, RFLOAT helix_length, RFLOAT psi_prior_flip_ratio, int selection_type);
 
     void convertHelicalLineCoordsToMetaDataTable(
 		FileName& fn_in,