GDTools COMPARE TSV output headers can vary

[jvera888]

Hi!

Breseq is a great tool! thank you for all your work.

I'm not sure if this reported behavior is a bug or a feature, so this is a question with a potential bug report.

question/issue: When using using GDTools COMPARE on multiple samples to produce TSV file output, I will usually get the same (40) headers, but occasionally get a different number. Some headers go missing (e.g. 'new_seq'), while new ones appear (e.g. 'repeat_length'). This behavior is somewhat unpredictable (i.e. does not always correspond to types of mutations present; lots of other headers always appear even when column is empty) and thus perhaps is unintended? It definitely makes parsing the table problematic. Here are a couple of different header row examples:

aa_new_seq	aa_position	aa_ref_seq	clone	codon_new_seq	codon_number	codon_position	codon_position_is_indeterminate	codon_ref_seq	gene_name	gene_position	gene_product	gene_strand	genes_inactivated	genes_overlapping	genes_promoter	locus_tag	locus_tags_inactivated	locus_tags_overlapping	locus_tags_promoter	mutation_category	mutator_status	new_read_count	new_read_count_basis	population	position	position_end	position_start	ref_read_count	ref_read_count_basis	ref_seq	repeat_length	repeat_new_copies	repeat_ref_copies	repeat_seq	seq_id	size	snp_type	time	title	transl_table	treatment	type
									PM_10340/PM_10350	intergenic (+150/+5)	ester cyclase/hypothetical protein	>/<				PM_10340/PM_10350				small_indel		12	1		964623	964663	964623	5	2	41-bp					PM_chr1	41		-1	output			DEL
									PM39400_25310	coding (99-107/837 nt)	DNA-directed RNA polymerase beta subunit	>	PM39400_25310			PM39400_25310	PM_25310			small_indel		79	1		2361336	2361344	2361336	11	2	AGTAGCCCC	9	2	3	AGTAGCCCC	PM_chr1	9		-1	output			DEL

aa_new_seq	aa_position	aa_ref_seq	clone	codon_new_seq	codon_number	codon_position	codon_position_is_indeterminate	codon_ref_seq	gene_name	gene_position	gene_product	gene_strand	genes_inactivated	genes_overlapping	genes_promoter	locus_tag	locus_tags_inactivated	locus_tags_overlapping	locus_tags_promoter	mutation_category	mutator_status	new_read_count	new_read_count_basis	new_seq	population	position	position_end	position_start	ref_read_count	ref_read_count_basis	ref_seq	seq_id	size	snp_type	time	title	transl_table	treatment	type
T	263	T		ACG	263	3		ACA	PM39400_71000	789	IS200/IS605 family transposase ISBce3	<		PM39400_71000		PM39400_71000		PM39400_71000		snp_synonymous		9	1	C		382	382	382	0	1	T	PM39400_pla12		synonymous	-1	PM39400	1		SNP
									PM39400_10340/PM39400_10350	intergenic (+150/+5)	ester cyclase/hypothetical protein	>/<				PM39400_10340/PM39400_10350				small_indel		10	1			964623	964663	964623	1	2	41-bp	PM39400_chr1	41		-1	PE39400-G4			DEL

Thanks for your help!

[jeffreybarrick]

Fields appearing or not is due to the presence of certain mutations or their annotations with certain named fields. For example, indels in simple sequence repeats add these fields: repeat_length, repeat_new_copies, repeat_ref_copies, repeat_seq. When an INS is present, then new_seq is there. ETC.

The TSV/CSV output is meant as an overall dump of all information. It would be hard to keep this up-to-date as other fields potentially get added by "fixing" the header to always be the same and include all possibilities.

If using R or Python, it should be easy to load and access data by the column names rather than by column position so that this isn't a problem. The "core" columns, which are the only ones you likely care about, will always be there. Is there a reason that won't work for you?

[dannagifford]

I think this causes issues if one were trying to merge TSV/CSV outputs via base R rbind(), which will complain about mismatching column names, but probably not if using dplyr/tidyverse bind_rows(), which put NA values for any missing columns. FWIW, it's also problematic for joining files with UNIX cat.

Though there can also be a problem with bind_rows() too due to mismatching data types--e.g. depending on what's in each column output by gdtools, a particular column might get typed as either chr, dbl, or lgl (logical, the default type if every entry is NA or empty iirc) when read in by readr's read_csv() (NB: if you're hitting this problem, sometimes setting guess_max=Inf in read_csv() can help! This makes R check every row before deciding what data type a column is.)

It would be helpful somewhere to have these as a workaround to incorporate into R scripts.

1. a list of all the possible column names that gdtools can possibly output
2. a list of possible data types for each column name so that you can specify col_types rather than relying on readr's guessing function.

Then we could do something like recreate the missing columns and specify a col_types string for read_csv().

I just remembered that this is what inspired my bodge_rows() hack...

[jvera888]

Hi,

Thanks for the responses!

I think @dannagifford 's suggestion of specifying potential output columns would be sufficient to design workarounds for the variable header columns. Perhaps listed somewhere in the documentation or even a json/yaml file?

[jeffreybarrick]

From my side, there are maintainability issues with specifying a set number of columns that will always be in the CSV/TSV output. One is that GenomeDiff format allows the user or another program to add custom fields to mutations/evidence items. They will be preserved and output currently by gdtools COMPARE/ANNOTATE.

If you pass all the GenomeDiff files you want to a single gdtools COMPARE/ANNOTATE command, then it will do the full table merge for you and output a large TSV file. Then, there will be no problem loading that and looking for or filtering to the columns (by name) that your code expects.

Is there something in your use case that is not covered by doing things this way?

[dannagifford]

Maintainability is very important, so I totally get where you're coming from Jeff (and we're really grateful that you do maintain it so well!). As users, could we parse the potential column names from the source code? (which files(s)?)

These are some use cases for being able to merge CSVs directly:

1. Someone may want to combine output from two different reference genomes. e.g. We're currently doing experiments with 8 different Pseudomonas species. Cannot run a single gdtools COMPARE for this, I think?
2. Sometimes one no longer has the output.gd evidence files any more, so can't rerun gdtools COMPARE (usually because of our cluster's /scratch file system deletion policy).
3. Someone might want to include gdtools COMPARE output as part of a pipeline. Inconsistent headers make it challenging.

Best wishes,
Danna

[jeffreybarrick]

There are lots of potential column names and they exist several places in the source code.

The "core" ones can be found at the top of genome_diff_entry.h. Anything here is lower-cased to get the corresponding output column/field names.

  const map<gd_entry_type, vector<string> > line_specification = make_map<gd_entry_type, vector<string> >
  //! seq_id and positions are already parameters in cDiffEntry
  //## mutations
  (SNP,make_vector<string> (SEQ_ID)(POSITION)(NEW_SEQ))
  (SUB,make_vector<string> (SEQ_ID)(POSITION)(SIZE)(NEW_SEQ))
  (DEL,make_vector<string> (SEQ_ID)(POSITION)(SIZE))
  (INS,make_vector<string> (SEQ_ID)(POSITION)(NEW_SEQ))
  (MOB,make_vector<string> (SEQ_ID)(POSITION)(REPEAT_NAME)(STRAND)(DUPLICATION_SIZE))
  (AMP,make_vector<string> (SEQ_ID)(POSITION)(SIZE)(NEW_COPY_NUMBER))
  (INV,make_vector<string> (SEQ_ID)(POSITION)(SIZE))
  (CON,make_vector<string> (SEQ_ID)(POSITION)(SIZE)(REGION))
  (INT,make_vector<string> (SEQ_ID)(POSITION)(SIZE)(REGION))
  
  //## evidence
  (RA,make_vector<string> (SEQ_ID)(POSITION)(INSERT_POSITION)(REF_BASE)(NEW_BASE))
  (MC,make_vector<string> (SEQ_ID)(START)(END)(START_RANGE)(END_RANGE))
  (JC,make_vector<string> (SIDE_1_SEQ_ID)(SIDE_1_POSITION)(SIDE_1_STRAND)(SIDE_2_SEQ_ID)(SIDE_2_POSITION)(SIDE_2_STRAND)(OVERLAP))
  (CN,make_vector<string> (SEQ_ID)(START)(END)(COPY_NUMBER))
  (UN,make_vector<string> (SEQ_ID)(START)(END))

...

  ; // end line specifications
  
  // These specs include addition fields used when determined equal mutations and sorting!
  // IMPORTANT: They include fields that may not always be defined.
  // NOTE: 'unique' gets added to ALL specs for determining uniqueness
  const map<gd_entry_type, vector<string> > extended_line_specification = make_map<gd_entry_type, vector<string> >
  (INS,make_vector<string> (SEQ_ID)(POSITION)(INSERT_POSITION)(NEW_SEQ))
  (MOB,make_vector<string> (SEQ_ID)(POSITION)(REPEAT_NAME)(STRAND)(DUPLICATION_SIZE)(INS_START)(INS_END)(DEL_START)(DEL_END)(MOB_REGION))
  (AMP,make_vector<string> (SEQ_ID)(POSITION)(SIZE)(NEW_COPY_NUMBER)(MEDIATED)(MEDIATED_STRAND)(MOB_REGION))
  (RA,make_vector<string>  (SEQ_ID)(POSITION)(INSERT_POSITION)(REF_BASE)(NEW_BASE))
  (JC,make_vector<string>  (SIDE_1_SEQ_ID)(SIDE_1_POSITION)(SIDE_1_STRAND)(SIDE_2_SEQ_ID)(SIDE_2_POSITION)(SIDE_2_STRAND)(OVERLAP)(UNIQUE_READ_SEQUENCE))
  ;

This block defines which are non-text fields for validation and sorting. Everything else is treated as free text internally by breseq.

  enum diff_entry_field_variable_t {
    kDiffEntryFieldVariableType_BaseSequence,
    kDiffEntryFieldVariableType_PositiveInteger,
    kDiffEntryFieldVariableType_NonNegativeInteger,
    kDiffEntryFieldVariableType_Integer,
    kDiffEntryFieldVariableType_Strand, // must be -1 or +1
    kDiffEntryFieldVariableType_PositiveInteger_ReverseSort,
  };
  
  map<string, diff_entry_field_variable_t > diff_entry_field_variable_types = make_map<string, diff_entry_field_variable_t>
  (POSITION, kDiffEntryFieldVariableType_PositiveInteger)
  (START, kDiffEntryFieldVariableType_PositiveInteger)
  (END, kDiffEntryFieldVariableType_PositiveInteger)
  (SIZE, kDiffEntryFieldVariableType_PositiveInteger_ReverseSort)
  (STRAND, kDiffEntryFieldVariableType_Strand)
  (NEW_SEQ, kDiffEntryFieldVariableType_BaseSequence)
  (NEW_COPY_NUMBER, kDiffEntryFieldVariableType_PositiveInteger)
  (DUPLICATION_SIZE, kDiffEntryFieldVariableType_Integer)
  (SIZE, kDiffEntryFieldVariableType_NonNegativeInteger)
  (DEL_START, kDiffEntryFieldVariableType_PositiveInteger)
  (DEL_END, kDiffEntryFieldVariableType_PositiveInteger)
  (INS_START, kDiffEntryFieldVariableType_BaseSequence)
  (INS_END, kDiffEntryFieldVariableType_BaseSequence)
  (INSERT_POSITION, kDiffEntryFieldVariableType_Integer)
  (MEDIATED_STRAND, kDiffEntryFieldVariableType_Strand)
  //JC item
  (SIDE_1_POSITION, kDiffEntryFieldVariableType_PositiveInteger)
  (SIDE_1_STRAND, kDiffEntryFieldVariableType_Strand)
  (SIDE_2_POSITION, kDiffEntryFieldVariableType_PositiveInteger)
  (SIDE_2_STRAND, kDiffEntryFieldVariableType_Strand)
  (OVERLAP, kDiffEntryFieldVariableType_Integer)
  (UNIQUE_READ_SEQUENCE, kDiffEntryFieldVariableType_BaseSequence)
  ;

[jeffreybarrick]

The solution I would propose is still to code something in R/Python that takes a set of CSV files to merge fills in the columns that a downstream pipeline expects with blanks or NAs when they are not present. I'd rather help you code a robust version of this and maintain that than have gdtools COMPARE print out extra columns that will often be empty given how many of these potential fields that there are (there are many defined for the annotations) and how some will only be used for specific purposes or by specific users.

It is possible this could be a utility script distributed as part of the genomediff Python package:

https://github.com/barricklab/genomediff-python
https://pypi.org/project/genomediff/

[abrozzi]

Hi, @jeffreybarrick. I am on it. I have started a simple R package that reads and parse "evidence.gd" files

Coming soon.

HTH
-A