Repository: ComparativeGenomicsToolkit/taffy (language: python)

Issue:
taffy norm seems to miss straightforward case
input maf
```
##maf version=1

a
s	x	0	9	+	50	CAAATAAGG
s	y	0	8	+	50	CAAATAAG-
s	z	0	9	+	50	CAAATAAGG

a
s	y	8	1	+	50	A

a
s	x	9	41	+	50	CTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
s	y	9	41	+	50	CTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
s	z	9	41	+	50	CTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
```

running this command doesn't change anything
```
taffy view -i input.maf | taffy norm -k
```




Discussion / hints:
This is the "raw" output of `vg2maf`.  So as it stands, normalization is able to connect consecutive nodes, but does not handle the variant. 

```
##maf version=1

a
s	x	0	8	+	50	CAAATAAG
s	y	0	8	+	50	CAAATAAG
s	z	0	8	+	50	CAAATAAG

a
s	x	8	1	+	50	G
s	z	8	1	+	50	G

a
s	y	8	1	+	50	A

a
s	x	9	1	+	50	C
s	y	9	1	+	50	C
s	z	9	1	+	50	C

a
s	x	10	3	+	50	TTG
s	y	10	3	+	50	TTG
s	z	10	3	+	50	TTG
```

![tiny1](https://github.com/ComparativeGenomicsToolkit/taffy/assets/901102/0d7eb13e-6736-4d00-bb10-8ce597be9a38)

Taffy norm doesn't merge these blocks because of the -q parameter, which specifies a minimum fraction of shared sequences between two blocks to merge. By default this is set of 0.6, setting it to 0, i.e.:

./bin/taffy view -i temp.maf | ./bin/taffy norm -q 0 | ./bin/taffy view -m

I get:

##maf version=1

a
s       x       9       41      +       50      ----------CTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
s       x       0       9       +       50      CAAATAAGG------------------------------------------
s       y       0       50      +       50      CAAATAAG-ACTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
s       z       9       41      +       50      ----------CTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG
s       z       0       9       +       50      CAAATAAGG------------------------------------------


which does the merge, but now I need to figure out why it doesn't actually properly merge together the x and z rows


Fix patch:
diff --git a/Makefile b/Makefile
index d9e8603..0c0e6c9 100644
--- a/Makefile
+++ b/Makefile
@@ -24,8 +24,8 @@ ${sonLibDir}/sonLib.a : sonLib
 
 ${sonLibDir}/cuTest.a : sonLib
 
-${LIBDIR}/libstTaf.a : ${libTests} ${libHeaders} ${srcDir}/alignment_block.o ${srcDir}/line_iterator.o ${srcDir}/maf.o ${srcDir}/paf.o ${srcDir}/ond.o ${srcDir}/taf.o ${srcDir}/merge_adjacent_alignments.o ${srcDir}/prefix_sort.o ${srcDir}/tai.o ${libHeaders} ${stTafDependencies}
-	${AR} rc libstTaf.a ${srcDir}/alignment_block.o ${srcDir}/line_iterator.o ${srcDir}/maf.o ${srcDir}/paf.o ${srcDir}/ond.o ${srcDir}/taf.o ${srcDir}/merge_adjacent_alignments.o ${srcDir}/prefix_sort.o ${srcDir}/tai.o
+${LIBDIR}/libstTaf.a : ${libTests} ${libHeaders} ${srcDir}/alignment_block.o ${srcDir}/line_iterator.o ${srcDir}/maf.o ${srcDir}/paf.o ${srcDir}/ond.o ${srcDir}/taf.o ${srcDir}/add_gap_bases.o ${srcDir}/merge_adjacent_alignments.o ${srcDir}/prefix_sort.o ${srcDir}/tai.o ${libHeaders} ${stTafDependencies}
+	${AR} rc libstTaf.a ${srcDir}/alignment_block.o ${srcDir}/line_iterator.o ${srcDir}/maf.o ${srcDir}/paf.o ${srcDir}/ond.o ${srcDir}/taf.o ${srcDir}/add_gap_bases.o ${srcDir}/merge_adjacent_alignments.o ${srcDir}/prefix_sort.o ${srcDir}/tai.o
 	mv libstTaf.a ${LIBDIR}/
 
 ${srcDir}/alignment_block.o : ${srcDir}/alignment_block.c ${libHeaders}
@@ -43,6 +43,9 @@ ${srcDir}/paf.o : ${srcDir}/paf.c ${libHeaders}
 ${srcDir}/ond.o : ${srcDir}/ond.c ${libHeaders}
 	${CC} ${CFLAGS} ${LDFLAGS} -o ${srcDir}/ond.o -c ${srcDir}/ond.c
 
+${srcDir}/add_gap_bases.o : ${srcDir}/add_gap_bases.cpp ${libHeaders}
+	${CXX} ${CPPFLAGS} ${CXXFLAGS} -o ${srcDir}/add_gap_bases.o -c ${srcDir}/add_gap_bases.cpp
+
 ${srcDir}/merge_adjacent_alignments.o : ${srcDir}/merge_adjacent_alignments.c ${libHeaders}
 	${CC} ${CFLAGS} ${LDFLAGS} -o ${srcDir}/merge_adjacent_alignments.o -c ${srcDir}/merge_adjacent_alignments.c
 
diff --git a/README.md b/README.md
index c97ac41..6c84309 100644
--- a/README.md
+++ b/README.md
@@ -269,11 +269,12 @@ There is also a utility to merge together short alignment blocks to create a mor
 
 For example, to normalize a maf file do the following:
 
-    taffy view -i MAF_FILE | taffy add-gap-bases SEQUENCE_FILES | taffy norm -k > out.maf
+    taffy view -i MAF_FILE | taffy norm -k -b SEQUENCE_FILES out.maf
 
-`taffy view` converts MAF_FILE into taf, `taffy add-gap-bases` adds in missing
-unaligned sequences between maf blocks and `taffy norm` then merges together the blocks. The 
-`-k` option causes the output to be in maf format.
+`taffy view` converts MAF_FILE into taf, `taffy norm` then merges together blocks. The 
+`-k` option causes the output to be in maf format. The `-b`
+ option reads in underlying sequence files and is used to 
+retrieve any sequences that are unaligned between two blocks that is necessary to include in stitching together adjacent blocks. This uses the same method as `taffy add-gap-bases` to add these unaligned sequences.
 
 ## Taffy Sort
 
@@ -281,7 +282,11 @@ It can be useful to sort the rows of an alignment. For this we have `taffy sort`
 
     taffy view -i MAF_FILE | taffy sort -n SORT_FILE | taffy view -m
 
-Taffy view first converts the input maf to TAF, taffy sort then sorts the rows of the taf according to a given file and finally the last taffy view pipes the output back to maf. The sort file is a sequence of sequence name prefixes, so that each sequence name in the alignment has, uniquely, one of the strings in the sort file as a prefix. The order of these prefixes is then used to sort the rows.
+Taffy view first converts the input maf to TAF, taffy sort then sorts the rows of the taf according to a given file and finally the last taffy view pipes the output back to maf. The sort file is a sequence of sequence name prefixes, so that each sequence name in the alignment has, uniquely, one of the strings in the sort file as a prefix. The order of these prefixes is then used to sort the rows. Taffy sort also supports filtering, to remove selected rows of an alignment. For example:
+
+    taffy view -i MAF_FILE | taffy sort -f FILTER_FILE -n SORT_FILE | taffy view -m
+
+Where here we additionally remove any rows with sequence names with a prefix contained in the given FILTER_FILE.
 
 # Referenced-based MAF/TAF and Indexing
 
diff --git a/taf_add_gap_bases.cpp b/taf_add_gap_bases.cpp
index b68709a..1d19763 100644
--- a/taf_add_gap_bases.cpp
+++ b/taf_add_gap_bases.cpp
@@ -15,7 +15,6 @@ extern "C" {
 #endif
 
 static int64_t maximum_gap_string_length = 50;
-static int64_t repeat_coordinates_every_n_columns = 1000;
 
 static void usage() {
     fprintf(stderr, "taffy add_gap_bases SEQ_FILExN [options]\n");    
@@ -26,81 +25,9 @@ static void usage() {
     fprintf(stderr, "-m --maximumGapStringLength : The maximum size of a gap string to add, be default: %" PRIi64 "\n",
             maximum_gap_string_length);
     fprintf(stderr, "-l --logLevel : Set the log level\n");
-    fprintf(stderr, "-s --repeatCoordinatesEveryNColumns : Repeat coordinates of each sequence at least every n columns. By default: %" PRIi64 "\n", repeat_coordinates_every_n_columns);
-    fprintf(stderr, "-c --useCompression : Write the output using bgzip compression.\n");
     fprintf(stderr, "-h --help : Print this help message\n");
 }
 
-static void add_to_hash(void *fastas, const char *fasta_header, const char *sequence, int64_t length) {
-    if(stHash_search((stHash *)fastas, (void *)fasta_header) != NULL) {
-        // c++ gives an angry warning if we try to send our string literal directly to st_errAbort, so we do this
-        char msg[8192];
-        snprintf(msg, 8192, "Found duplicate sequence header: %s\n", fasta_header);
-        st_errAbort(msg);
-    }
-    stHash_insert((stHash *)fastas, stString_copy(fasta_header), stString_copy(sequence));
-}
-
-// get a dna interval either from the fastas hash file or from the hal_handle
-// note the string returned needs to be freed
-static char *get_sequence_fragment(const char* sequence_name, int64_t start, int64_t length, stHash *fastas, int hal_handle, stSet *hal_species) {
-    char *fragment = NULL;
-    if (fastas) {
-        assert(hal_handle == -1);
-        char *seq = (char*)stHash_search(fastas, (void*)sequence_name);
-        if(seq != NULL) {
-            fragment = stString_getSubString(seq, start, length);
-        }
-    } else {
-#ifdef USE_HAL
-        assert(fastas == NULL);
-        char* species_name = extract_genome_name(sequence_name, hal_species, NULL);
-        char* chrom_name = (char*)sequence_name + strlen(species_name) + 1;
-        fragment = halGetDna(hal_handle, species_name, chrom_name, start, start + length, NULL);
-        free(species_name);
-#else
-        assert(false);
-#endif
-    }
-
-    return fragment;
-}
-
-
-static void add_gap_strings(Alignment *p_alignment, Alignment *alignment, stHash *fastas, int hal_handle, stSet *hal_species) {
-    Alignment_Row *row = alignment->row;
-    while(row != NULL) {
-        if(row->l_row != NULL && alignment_row_is_predecessor(row->l_row, row)) {
-            int64_t gap_length = row->start - (row->l_row->start + row->l_row->length);
-            if(gap_length <= maximum_gap_string_length && row->left_gap_sequence == NULL) {
-                char* seq_interval = NULL;
-                int64_t i = row->l_row->start + row->l_row->length;
-                assert(i >= 0 && i < row->sequence_length);
-                if(row->strand) {
-                    seq_interval = get_sequence_fragment(row->sequence_name, i, gap_length, fastas, hal_handle, hal_species);
-                }
-                else { // Case sequence is on the negative strand
-                    assert(row->sequence_length - i - gap_length >= 0);
-                    char *s = get_sequence_fragment(row->sequence_name,
-                                                    row->sequence_length - i - gap_length, gap_length, fastas, hal_handle, hal_species);
-                    if(s != NULL) {
-                        seq_interval = stString_reverseComplementString(s);
-                        free(s);
-                    }
-                }
-                if(seq_interval == NULL) {
-                    st_logDebug("[taf] Missing sequence for gap, seq name: %s, skipping!\n", row->sequence_name);
-                }
-                else {
-                    row->left_gap_sequence = seq_interval;
-                    assert(strlen(row->left_gap_sequence) == (size_t)gap_length);
-                }
-            }
-        }
-        row = row->n_row;
-    }
-}
-
 int taf_add_gap_bases_main(int argc, char *argv[]) {
     time_t startTime = time(NULL);
 
@@ -112,7 +39,6 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
     char *outputFile = NULL;
     char *hal_file = NULL;
     bool run_length_encode_bases = 0;
-    bool use_compression = 0;
 
     ///////////////////////////////////////////////////////////////////////////
     // Parse the inputs
@@ -125,12 +51,10 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
                                                 { "halFile", required_argument, 0, 'a' },
                                                 { "help", no_argument, 0, 'h' },
                                                 { "maximumGapStringLength", required_argument, 0, 'm' },
-                                                { "repeatCoordinatesEveryNColumns", required_argument, 0, 's' },
-                                                { "useCompression", no_argument, 0, 'c' },
                                                 { 0, 0, 0, 0 } };
 
         int option_index = 0;
-        int64_t key = getopt_long(argc, argv, "l:i:o:a:hm:s:c", long_options, &option_index);
+        int64_t key = getopt_long(argc, argv, "l:i:o:a:hm:", long_options, &option_index);
         if (key == -1) {
             break;
         }
@@ -151,15 +75,9 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
             case 'h':
                 usage();
                 return 0;
-            case 's':
-                repeat_coordinates_every_n_columns = atol(optarg);
-                break;
             case 'm':
                 maximum_gap_string_length = atol(optarg);
                 break;
-            case 'c':
-                use_compression = 1;
-                break;
             default:
                 usage();
                 return 1;
@@ -190,8 +108,6 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
         st_logInfo("Number of input FASTA files : %ld\n", argc - optind);
     }            
     st_logInfo("Maximum maximum gap string length : %" PRIi64 "\n", maximum_gap_string_length);
-    st_logInfo("Repeat coordinates every n bases : %" PRIi64 "\n", repeat_coordinates_every_n_columns);
-    st_logInfo("Write compressed output : %s\n", use_compression ? "true" : "false");
 
     //////////////////////////////////////////////
     // Read in the sequence files
@@ -201,27 +117,10 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
     stSet *hal_species = NULL;
     int hal_handle = -1;
     if (optind < argc) {
-        fastas = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, free, free);
-        while(optind < argc) {
-            char *seq_file = argv[optind++];
-            st_logInfo("Parsing sequence file : %s\n", seq_file);
-            FILE *fh = fopen(seq_file, "r");
-            fastaReadToFunction(fh, fastas, add_to_hash);
-            fclose(fh);
-        }
-        st_logInfo("Finished parsing sequence files\n");
-    } else {
-#ifdef USE_HAL
-        hal_handle = halOpen(hal_file, NULL);
-        // make a table of every species in the hal. this will be used to try to parse
-        // contig names with lots of dots in them where we don't know which dot separates genome from chrom...
-        hal_species = stSet_construct3(stHash_stringKey, stHash_stringEqualKey, free);
-        struct hal_species_t *species_list = halGetSpecies(hal_handle, NULL);
-        for (struct hal_species_t *species = species_list; species != NULL; species = species->next) {
-            stSet_insert(hal_species, stString_copy(species->name));
-        }
-        halFreeSpeciesList(species_list);
-#endif
+        fastas = load_sequences_from_fasta_files(&(argv[optind]), argc - optind);
+    }
+    else {
+        hal_species = load_sequences_from_hal_file(hal_file, &hal_handle);
     }
 
     //////////////////////////////////////////////
@@ -229,15 +128,11 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
     //////////////////////////////////////////////
 
     FILE *input = inputFile == NULL ? stdin : fopen(inputFile, "r");
-    LW *output = LW_construct(outputFile == NULL ? stdout : fopen(outputFile, "w"), use_compression);
+    LW *output = LW_construct(outputFile == NULL ? stdout : fopen(outputFile, "w"), 0);
     LI *li = LI_construct(input);
 
     // Pass the header line to determine parameters and write the updated taf header
-    Tag *tag = taf_read_header(li);
-    Tag *t = tag_find(tag, (char*)"run_length_encode_bases");
-    if(t != NULL && strcmp(t->value, "1") == 0) {
-        run_length_encode_bases = 1;
-    }
+    Tag *tag = taf_read_header_2(li, &run_length_encode_bases);
     taf_write_header(tag, output);
     tag_destruct(tag);
 
@@ -245,11 +140,11 @@ int taf_add_gap_bases_main(int argc, char *argv[]) {
     while((alignment = taf_read_block(p_alignment, run_length_encode_bases, li)) != NULL) {
         // Add in the gap strings if there is a previous block
         if(p_alignment != NULL) {
-            add_gap_strings(p_alignment, alignment, fastas, hal_handle, hal_species);
+            alignment_add_gap_strings(p_alignment, alignment, fastas, hal_handle, hal_species, maximum_gap_string_length);
         }
 
         // Write the block
-        taf_write_block(p_alignment, alignment, run_length_encode_bases, repeat_coordinates_every_n_columns, output);
+        taf_write_block(p_alignment, alignment, run_length_encode_bases, -1, output);
 
         // Clean up the previous alignment
         if(p_alignment != NULL) {
diff --git a/taf_norm.c b/taf_norm.c
index 33ab919..125e875 100644
--- a/taf_norm.c
+++ b/taf_norm.c
@@ -11,7 +11,8 @@
 
 int64_t maximum_block_length_to_merge = 200;
 int64_t maximum_gap_length = 30;
-float fraction_shared_rows = 0.6;
+int64_t minimum_shared_rows = 1;
+float fraction_shared_rows = 0.0;
 int64_t repeat_coordinates_every_n_columns = 1000;
 
 static void usage(void) {
@@ -23,10 +24,13 @@ static void usage(void) {
     fprintf(stderr, "-k --maf : Print maf output instead of taf\n");
     fprintf(stderr, "-m --maximumBlockLengthToMerge : Only merge together any two adjacent blocks if one or both is less than this many bases long, by default: %" PRIi64 "\n", maximum_block_length_to_merge);
     fprintf(stderr, "-n --maximumGapLength : Only merge together two adjacent blocks if the total number of unaligned bases between the blocks is less than this many bases, by default: %" PRIi64 "\n", maximum_gap_length);
+    fprintf(stderr, "-Q --minimumSharedRows : The minimum number of rows between two blocks that need to be shared for a merge, default: %" PRIi64 "\n", minimum_shared_rows);
     fprintf(stderr, "-q --fractionSharedRows : The fraction of rows between two blocks that need to be shared for a merge, default: %f\n", fraction_shared_rows);
     fprintf(stderr, "-d --filterGapCausingDupes : Reduce the number of MAF blocks by filtering out rows that induce gaps > maximumGapLength. Rows are only filtered out if they are duplications (contig of same name appears elsewhere in block, or contig with same prefix up to \".\" appears in the same block).\n");
     fprintf(stderr, "-s --repeatCoordinatesEveryNColumns : Repeat coordinates of each sequence at least every n columns. By default: %" PRIi64 "\n", repeat_coordinates_every_n_columns);
     fprintf(stderr, "-c --useCompression : Write the output using bgzip compression.\n");
+    fprintf(stderr, "-a --halFile : HAL file for extracting gap sequence (MAF must be created with hal2maf *without* --onlySequenceNames)\n");
+    fprintf(stderr, "-b --seqFiles : Fasta files for extracting gap sequence. Do not specify both this option and --halFile\n");
     fprintf(stderr, "-h --help : Print this help message\n");
 }
 
@@ -139,14 +143,6 @@ static bool greedy_prune_by_gap(Alignment *alignment, int64_t maximum_gap_length
             if (row == row_to_prune) {
                 assert(p_row != NULL);
                 p_row->n_row = row->n_row;
-                if (row->l_row) {
-                    row->l_row->r_row = NULL;
-                } 
-                if (row->r_row) {
-                    row->r_row->l_row = NULL;
-                    free(row->r_row->left_gap_sequence);
-                    row->r_row->left_gap_sequence = NULL;                    
-                }
                 alignment_row_destruct(row);
                 ++to_prune_idx;
                 row_to_prune = to_prune_idx < stList_length(to_prune) ? stList_get(to_prune, to_prune_idx) : NULL;
@@ -168,6 +164,7 @@ static bool greedy_prune_by_gap(Alignment *alignment, int64_t maximum_gap_length
     return pruned;
 }
 
+
 int taf_norm_main(int argc, char *argv[]) {
     time_t startTime = time(NULL);
 
@@ -181,6 +178,8 @@ int taf_norm_main(int argc, char *argv[]) {
     bool output_maf = 0;
     bool use_compression = 0;
     bool filter_gap_causing_dupes = 0;
+    stList *fasta_files = stList_construct();
+    char *hal_file = NULL;
 
     ///////////////////////////////////////////////////////////////////////////
     // Parse the inputs
@@ -195,13 +194,16 @@ int taf_norm_main(int argc, char *argv[]) {
                                                 { "maximumBlockLengthToMerge", required_argument, 0, 'm' },
                                                 { "maximumGapLength", required_argument, 0, 'n' },
                                                 { "fractionSharedRows", required_argument, 0, 'q' },
+                                                { "minimumSharedRows", required_argument, 0, 'Q' },
                                                 { "filterGapCausingDupes", no_argument, 0, 'd' },
                                                 { "repeatCoordinatesEveryNColumns", required_argument, 0, 's' },
                                                 { "useCompression", no_argument, 0, 'c' },
+                                                { "halFile", required_argument, 0, 'a' },
+                                                { "seqFiles", required_argument, 0, 'b' },
                                                 { 0, 0, 0, 0 } };
 
         int option_index = 0;
-        int64_t key = getopt_long(argc, argv, "l:i:o:hcm:n:dkq:s:", long_options, &option_index);
+        int64_t key = getopt_long(argc, argv, "l:i:o:hcm:n:dkQ:q:s:a:b:", long_options, &option_index);
         if (key == -1) {
             break;
         }
@@ -231,6 +233,9 @@ int taf_norm_main(int argc, char *argv[]) {
             case 'd':
                 filter_gap_causing_dupes = 1;
                 break;
+            case 'Q':
+                minimum_shared_rows = atol(optarg);
+                break;
             case 'q':
                 fraction_shared_rows = atof(optarg);
                 break;
@@ -240,6 +245,16 @@ int taf_norm_main(int argc, char *argv[]) {
             case 's':
                 repeat_coordinates_every_n_columns = atol(optarg);
                 break;
+            case 'a':
+                hal_file = optarg;
+                break;
+            case 'b':
+                // Parse the set of sequence files (this is a bit fragile - files can not start with a '-' character)
+                optind--;
+                for( ;optind < argc && *argv[optind] != '-'; optind++){
+                    stList_append(fasta_files, argv[optind]);
+                }
+                break;
             default:
                 usage();
                 return 1;
@@ -260,6 +275,26 @@ int taf_norm_main(int argc, char *argv[]) {
     st_logInfo("Repeat coordinates every n bases : %" PRIi64 "\n", repeat_coordinates_every_n_columns);
     st_logInfo("Fraction shared rows to merge adjacent blocks : %f\n", fraction_shared_rows);
     st_logInfo("Write compressed output : %s\n", use_compression ? "true" : "false");
+    if (hal_file) {
+        st_logInfo("HAL file string : %s\n", hal_file);
+    } else {
+        st_logInfo("Number of input FASTA files : %ld\n", argc - optind);
+    }
+
+    //////////////////////////////////////////////
+    // Read in the sequences if joining over unaligned gaps
+    //////////////////////////////////////////////
+
+    stHash *fastas_map = NULL;
+    stSet *hal_species = NULL;
+    int hal_handle = -1;
+    if (hal_file) {
+        hal_species = load_sequences_from_hal_file(hal_file, &hal_handle);
+    }
+    else if(stList_length(fasta_files) > 0) {
+        fastas_map = load_sequences_from_fasta_files(stList_getBackingArray(fasta_files), stList_length(fasta_files));
+        stList_destruct(fasta_files);
+    }
 
     //////////////////////////////////////////////
     // Read in the taf blocks and merge blocks that are sufficiently small
@@ -270,13 +305,9 @@ int taf_norm_main(int argc, char *argv[]) {
     LI *li = LI_construct(input);
 
     // Pass the header line to determine parameters and write the updated taf header
-    Tag *tag = taf_read_header(li);
-    Tag *t = tag_find(tag, "run_length_encode_bases");
-    if(t != NULL && strcmp(t->value, "1") == 0) {
-        run_length_encode_bases = 1;
-        if(output_maf) { // Remove this tag from the maf output as not relevant
-            tag = tag_remove(tag, "run_length_encode_bases");
-        }
+    Tag *tag = taf_read_header_2(li, &run_length_encode_bases);
+    if(output_maf && run_length_encode_bases) { // Remove this tag from the maf output as not relevant
+        tag = tag_remove(tag, "run_length_encode_bases");
     }
     output_maf ? maf_write_header(tag, output) : taf_write_header(tag, output);
     tag_destruct(tag);
@@ -284,10 +315,15 @@ int taf_norm_main(int argc, char *argv[]) {
     Alignment *alignment, *p_alignment = NULL, *p_p_alignment = NULL;
     while((alignment = get_next_taf_block(li, run_length_encode_bases)) != NULL) {
         if(p_alignment != NULL) {
+            // First realign the rows in case we in the process of merging prior blocks we have
+            // identified rows that can be merged
+            alignment_link_adjacent(p_alignment, alignment, 1);
+
             bool merged = false;
             int64_t common_rows = alignment_number_of_common_rows(p_alignment, alignment);
             int64_t total_rows = alignment->row_number + p_alignment->row_number - common_rows;
-            if (common_rows >= total_rows * fraction_shared_rows &&
+            if (common_rows >= minimum_shared_rows &&
+                common_rows >= total_rows * fraction_shared_rows &&
                 (alignment_length(p_alignment) <= maximum_block_length_to_merge ||
                  alignment_length(alignment) <= maximum_block_length_to_merge)) {
                 int64_t total_gap = alignment_total_gap_length(p_alignment);
@@ -297,7 +333,10 @@ int taf_norm_main(int argc, char *argv[]) {
                     total_gap = alignment_total_gap_length(p_alignment);
                     assert(was_pruned == (total_gap <= maximum_gap_length));
                 }
-                if (total_gap <= maximum_gap_length) {                    
+                if (total_gap <= maximum_gap_length) {
+                    if(hal_species || fastas_map) { // Now add in any gap bases if sequences are provided
+                        alignment_add_gap_strings(p_alignment, alignment, fastas_map, hal_handle, hal_species, -1);
+                    }
                     p_alignment = alignment_merge_adjacent(p_alignment, alignment);
                     merged = true;
                 }
@@ -332,6 +371,13 @@ int taf_norm_main(int argc, char *argv[]) {
     }
     LW_destruct(output, outputFile != NULL);
 
+    if (fastas_map) {
+        stHash_destruct(fastas_map);
+    }
+    if (hal_species) {
+        stSet_destruct(hal_species);
+    }
+
     st_logInfo("taffy norm is done, %" PRIi64 " seconds have elapsed\n", time(NULL) - startTime);
 
     //while(1);
diff --git a/taf_sort.c b/taf_sort.c
index 3041020..71cc106 100644
--- a/taf_sort.c
+++ b/taf_sort.c
@@ -16,11 +16,44 @@ static void usage(void) {
     fprintf(stderr, "-i --inputFile : Input TAF or MAF file. If not specified reads from stdin\n");
     fprintf(stderr, "-o --outputFile : Output file. If not specified outputs to stdout\n");
     fprintf(stderr, "-n --sortFile : File in which each line is a prefix of a sequence name. Rows are sorted accordingly, \n"
-                    "with any ties broken by lexicographic sort of the suffixes. Can not be None\n");
+                    "with any ties broken by lexicographic sort of the suffixes.\n");
+    fprintf(stderr, "-f --filterFile : Remove any rows with sequences matching a prefix in this file\n");
     fprintf(stderr, "-l --logLevel : Set the log level\n");
     fprintf(stderr, "-h --help : Print this help message\n");
 }
 
+stList *load_sort_file(char *sort_file) {
+    if (sort_file == NULL) {
+        return NULL;
+    }
+    FILE *sort_fh = fopen(sort_file, "r");
+    if (sort_fh == NULL) {
+        fprintf(stderr, "Unable to open sort/filter file: %s\n", sort_file);
+        return 1;
+    }
+    stList *prefixes_to_sort_by = sequence_prefix_load(sort_fh);
+    st_logInfo("Loaded the sort/filter file, got %i rows\n", (int) stList_length(prefixes_to_sort_by));
+    return prefixes_to_sort_by;
+}
+
+void process_alignment_block(Alignment *pp_alignment, Alignment *p_alignment, stList *prefixes_to_filter_by,
+                             stList *prefixes_to_sort_by, bool run_length_encode_bases, LW *output) {
+    if(p_alignment) {
+        if(prefixes_to_filter_by) { //Remove rows matching a prefix
+            alignment_filter_the_rows(p_alignment, prefixes_to_filter_by);
+        }
+        if(prefixes_to_sort_by) { // Sort the alignment block rows
+            alignment_sort_the_rows(pp_alignment, p_alignment, prefixes_to_sort_by);
+        }
+        // Write the block
+        taf_write_block(pp_alignment, p_alignment,
+                        run_length_encode_bases, -1, output); // Write the block
+    }
+    if(pp_alignment != NULL) {
+        alignment_destruct(pp_alignment, 1); // Delete the left most block
+    }
+}
+
 int taf_sort_main(int argc, char *argv[]) {
     time_t startTime = time(NULL);
 
@@ -31,21 +64,23 @@ int taf_sort_main(int argc, char *argv[]) {
     char *input_file = NULL;
     char *output_file = NULL;
     char *sort_file = NULL;
+    char *filter_file = NULL;
 
     ///////////////////////////////////////////////////////////////////////////
     // Parse the inputs
     ///////////////////////////////////////////////////////////////////////////
 
     while (1) {
-        static struct option long_options[] = { { "logLevel", required_argument, 0, 'l' },
-                                                { "inputFile", required_argument, 0, 'i' },
-                                                { "outputFile", required_argument, 0, 'o' },
-                                                { "sortFile", required_argument, 0, 'n' },
-                                                { "help", no_argument, 0, 'h' },
-                                                { 0, 0, 0, 0 } };
+        static struct option long_options[] = {{"logLevel",   required_argument, 0, 'l'},
+                                               {"inputFile",  required_argument, 0, 'i'},
+                                               {"outputFile", required_argument, 0, 'o'},
+                                               {"sortFile",   required_argument, 0, 'n'},
+                                               {"filterFile", required_argument, 0, 'f'},
+                                               {"help",       no_argument,       0, 'h'},
+                                               {0, 0,                            0, 0}};
 
         int option_index = 0;
-        int64_t key = getopt_long(argc, argv, "l:i:o:n:h", long_options, &option_index);
+        int64_t key = getopt_long(argc, argv, "l:i:o:n:hf:", long_options, &option_index);
         if (key == -1) {
             break;
         }
@@ -63,6 +98,9 @@ int taf_sort_main(int argc, char *argv[]) {
             case 'n':
                 sort_file = optarg;
                 break;
+            case 'f':
+                filter_file = optarg;
+                break;
             case 'h':
                 usage();
                 return 0;
@@ -80,6 +118,7 @@ int taf_sort_main(int argc, char *argv[]) {
     st_logInfo("Input file string : %s\n", input_file);
     st_logInfo("Output file string : %s\n", output_file);
     st_logInfo("Sort file string : %s\n", sort_file);
+    st_logInfo("Filter file string : %s\n", filter_file);
 
     //////////////////////////////////////////////
     // Read in the taf/maf blocks and sort order file
@@ -102,22 +141,12 @@ int taf_sort_main(int argc, char *argv[]) {
     LW *output = LW_construct(output_fh, 0);
 
     // Sort file
-    if (sort_file == NULL) {
-        fprintf(stderr, "No sort file specified!\n");
-        return 1;
-    }
-    FILE *sort_fh = fopen(sort_file, "r");
-    if (sort_fh == NULL) {
-        fprintf(stderr, "Unable to open sort file: %s\n", sort_file);
-        return 1;
-    }
-    stList *prefixes_to_sort_by = sequence_prefix_load(sort_fh);
-    st_logInfo("Loaded the sort file, got %i rows\n", (int)stList_length(prefixes_to_sort_by));
+    stList *prefixes_to_sort_by = load_sort_file(sort_file);
+    stList *prefixes_to_filter_by = load_sort_file(filter_file);
 
     // Parse the header
-    Tag *tag = taf_read_header(li);
-    Tag *t = tag_find(tag, "run_length_encode_bases");
-    bool run_length_encode_bases = t != NULL && strcmp(t->value, "1") == 0;
+    bool run_length_encode_bases;
+    Tag *tag = taf_read_header_2(li, &run_length_encode_bases);
 
     // Write the header
     taf_write_header(tag, output);
@@ -126,24 +155,22 @@ int taf_sort_main(int argc, char *argv[]) {
     // Write the alignment blocks
     Alignment *alignment, *p_alignment = NULL, *pp_alignment = NULL;
     while((alignment = taf_read_block(p_alignment, run_length_encode_bases, li)) != NULL) {
-        if(p_alignment) {
-            // Sort the alignment block rows
-            alignment_sort_the_rows(pp_alignment, p_alignment, prefixes_to_sort_by);
-            // Write the block
-            taf_write_block(pp_alignment, p_alignment, run_length_encode_bases, -1, output); // Write the block
-        }
+        process_alignment_block(pp_alignment, p_alignment, prefixes_to_filter_by,
+                                prefixes_to_sort_by, run_length_encode_bases, output);
         pp_alignment = p_alignment;
         p_alignment = alignment;
     }
     if(p_alignment) { // Write the final block
-        alignment_sort_the_rows(pp_alignment, p_alignment, prefixes_to_sort_by);
-        taf_write_block(pp_alignment, p_alignment, run_length_encode_bases, -1, output); // Write the block
+        process_alignment_block(pp_alignment, p_alignment, prefixes_to_filter_by,
+                                prefixes_to_sort_by, run_length_encode_bases, output);
+        alignment_destruct(p_alignment, 1);
     }
 
     //////////////////////////////////////////////
     // Cleanup
     //////////////////////////////////////////////
 
+
     LI_destruct(li);
     if(input_file != NULL) {
         fclose(input);
diff --git a/taffy/impl/add_gap_bases.cpp b/taffy/impl/add_gap_bases.cpp
new file mode 100644
index 0000000..de26d7d
--- /dev/null
+++ b/taffy/impl/add_gap_bases.cpp
@@ -0,0 +1,111 @@
+extern "C" {
+#include "taf.h"
+#include "sonLib.h"
+}
+#include "bioioC.h"
+#include <getopt.h>
+#include <time.h>
+#ifdef USE_HAL
+#include "halBlockViz.h"
+#endif
+
+static void add_to_hash(void *fastas, const char *fasta_header, const char *sequence, int64_t length) {
+    if(stHash_search((stHash *)fastas, (void *)fasta_header) != NULL) {
+        // c++ gives an angry warning if we try to send our string literal directly to st_errAbort, so we do this
+        char msg[8192];
+        snprintf(msg, 8192, "Found duplicate sequence header: %s\n", fasta_header);
+        st_errAbort(msg);
+    }
+    stHash_insert((stHash *)fastas, stString_copy(fasta_header), stString_copy(sequence));
+}
+
+// get a dna interval either from the fastas hash file or from the hal_handle
+// note the string returned needs to be freed
+static char *get_sequence_fragment(const char* sequence_name, int64_t start, int64_t length, stHash *fastas, int hal_handle, stSet *hal_species) {
+    char *fragment = NULL;
+    if (fastas) {
+        assert(hal_handle == -1);
+        char *seq = (char*)stHash_search(fastas, (void*)sequence_name);
+        if(seq != NULL) {
+            fragment = stString_getSubString(seq, start, length);
+        }
+    } else {
+#ifdef USE_HAL
+        assert(fastas == NULL);
+        char* species_name = extract_genome_name(sequence_name, hal_species, NULL);
+        char* chrom_name = (char*)sequence_name + strlen(species_name) + 1;
+        fragment = halGetDna(hal_handle, species_name, chrom_name, start, start + length, NULL);
+        free(species_name);
+#else
+        assert(false);
+#endif
+    }
+
+    return fragment;
+}
+
+void alignment_add_gap_strings(Alignment *p_alignment, Alignment *alignment, stHash *fastas, int hal_handle, stSet *hal_species,
+                     int64_t maximum_gap_string_length) {
+    Alignment_Row *row = alignment->row;
+    while(row != NULL) {
+        if(row->l_row != NULL && alignment_row_is_predecessor(row->l_row, row)) {
+            int64_t gap_length = row->start - (row->l_row->start + row->l_row->length);
+            if((maximum_gap_string_length < 0 || gap_length <= maximum_gap_string_length) && row->left_gap_sequence == NULL) {
+                char* seq_interval = NULL;
+                int64_t i = row->l_row->start + row->l_row->length;
+                assert(i >= 0 && i < row->sequence_length);
+                if(row->strand) {
+                    seq_interval = get_sequence_fragment(row->sequence_name, i, gap_length, fastas, hal_handle, hal_species);
+                }
+                else { // Case sequence is on the negative strand
+                    assert(row->sequence_length - i - gap_length >= 0);
+                    char *s = get_sequence_fragment(row->sequence_name,
+                                                    row->sequence_length - i - gap_length, gap_length, fastas, hal_handle, hal_species);
+                    if(s != NULL) {
+                        seq_interval = stString_reverseComplementString(s);
+                        free(s);
+                    }
+                }
+                if(seq_interval == NULL) {
+                    st_logDebug("[taf] Missing sequence for gap, seq name: %s, skipping!\n", row->sequence_name);
+                }
+                else {
+                    row->left_gap_sequence = seq_interval;
+                    assert(strlen(row->left_gap_sequence) == (size_t)gap_length);
+                }
+            }
+        }
+        row = row->n_row;
+    }
+}
+
+stHash *load_sequences_from_fasta_files(char **seq_file_names, int64_t seq_file_number) {
+    stHash *fastas = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, free, free);
+    for (int64_t i = 0; i < seq_file_number; i++) { //optind < argc) {
+        char *seq_file = seq_file_names[i];
+        st_logInfo("Parsing sequence file : %s\n", seq_file);
+        FILE *fh = fopen(seq_file, "r");
+        fastaReadToFunction(fh, fastas, add_to_hash);
+        fclose(fh);
+    }
+    st_logInfo("Finished parsing sequence fasta files\n");
+    return fastas;
+}
+
+stSet *load_sequences_from_hal_file(char *hal_file, int *hal_handle) {
+    stSet *hal_species = NULL;
+    st_logInfo("Parsing hal file : %s\n", hal_file);
+#ifdef USE_HAL
+    *hal_handle = halOpen(hal_file, NULL);
+    // make a table of every species in the hal. this will be used to try to parse
+    // contig names with lots of dots in them where we don't know which dot separates genome from chrom...
+    hal_species = stSet_construct3(stHash_stringKey, stHash_stringEqualKey, free);
+    struct hal_species_t *species_list = halGetSpecies(*hal_handle, NULL);
+    for (struct hal_species_t *species = species_list; species != NULL; species = species->next) {
+        stSet_insert(hal_species, stString_copy(species->name));
+    }
+    halFreeSpeciesList(species_list);
+#endif
+    st_logInfo("Finished parsing hal file\n");
+    return hal_species;
+}
diff --git a/taffy/impl/alignment_block.c b/taffy/impl/alignment_block.c
index 601f3b1..2c07f2c 100644
--- a/taffy/impl/alignment_block.c
+++ b/taffy/impl/alignment_block.c
@@ -104,9 +104,15 @@ Tag *tag_parse(char *tag_string, char *delimiter, Tag *p_tag) {
     return tag;
 }
 
-
-
 void alignment_row_destruct(Alignment_Row *row) {
+    if(row->l_row != NULL) { // If there is a preceding, left row then unlink it
+        assert(row->l_row->r_row == row);
+        row->l_row->r_row = NULL;
+    }
+    if(row->r_row != NULL) { // If there is a proceeding, right row then unlink it
+        assert(row->r_row->l_row == row);
+        row->r_row->l_row = NULL;
+    }
     if(row->bases != NULL) {
         free(row->bases);
     }
@@ -191,6 +197,11 @@ void alignment_link_adjacent(Alignment *left_alignment, Alignment *right_alignme
     Alignment_Row *row = left_alignment->row;
     while(row != NULL) {
         row->r_row = NULL;
+        // Clean up any gap sequences because once unlinked we can not guarantee continuity
+        if(row->left_gap_sequence != NULL) {
+            free(row->left_gap_sequence);
+            row->left_gap_sequence = NULL;
+        }
         row = row->n_row;
     }
     row = right_alignment->row;
diff --git a/taffy/impl/merge_adjacent_alignments.c b/taffy/impl/merge_adjacent_alignments.c
index 9713ae2..eb2467f 100644
--- a/taffy/impl/merge_adjacent_alignments.c
+++ b/taffy/impl/merge_adjacent_alignments.c
@@ -254,29 +254,34 @@ Alignment *alignment_merge_adjacent(Alignment *left_alignment, Alignment *right_
             l_row->bases = bases;
         }
         else {
+            Alignment_Row *r_row = l_row->r_row;
+
             // Check the rows agree coordinate wise
-            assert(strcmp(l_row->sequence_name, l_row->r_row->sequence_name) == 0);
-            assert(l_row->strand == l_row->r_row->strand);
-            assert(l_row->start + l_row->length <= l_row->r_row->start);
+            assert(strcmp(l_row->sequence_name, r_row->sequence_name) == 0);
+            assert(l_row->strand == r_row->strand);
+            assert(l_row->start + l_row->length <= r_row->start);
 
             // Is not a deletion, so merge together two adjacent rows
-            assert(l_row->r_row->left_gap_sequence != NULL);
-            assert(strlen(l_row->r_row->left_gap_sequence) == interstitial_alignment_length);
-            char *bases = stString_print("%s%s%s", l_row->bases, l_row->r_row->left_gap_sequence, l_row->r_row->bases);
+            assert(r_row->left_gap_sequence != NULL);
+            assert(strlen(r_row->left_gap_sequence) == interstitial_alignment_length);
+            char *bases = stString_print("%s%s%s", l_row->bases, r_row->left_gap_sequence, r_row->bases);
             free(l_row->bases); // clean up
             l_row->bases = bases;
 
             // Update the left row's length coordinate
-            int64_t interstitial_bases = l_row->r_row->start - (l_row->start + l_row->length);
-            l_row->length += interstitial_bases + l_row->r_row->length;
+            int64_t interstitial_bases = r_row->start - (l_row->start + l_row->length);
+            l_row->length += interstitial_bases + r_row->length;
             // Update the l_row's r_row pointer...
-            if(l_row->r_row->r_row != NULL) { // Check pointers are correct
-                assert(l_row->r_row->r_row->l_row == l_row->r_row);
+            if(r_row->r_row != NULL) { // Check pointers are correct
+                assert(r_row->r_row->l_row == r_row);
             }
-            l_row->r_row = l_row->r_row->r_row;
+            l_row->r_row = r_row->r_row;
             if(l_row->r_row != NULL) {
                 l_row->r_row->l_row = l_row;
             }
+            // Null r_row's left and right pointers
+            r_row->l_row = NULL;
+            r_row->r_row = NULL;
         }
 
         l_row = l_row->n_row; // Move to the next left alignment row
diff --git a/taffy/impl/prefix_sort.c b/taffy/impl/prefix_sort.c
index 2fca853..66afc89 100644
--- a/taffy/impl/prefix_sort.c
+++ b/taffy/impl/prefix_sort.c
@@ -93,6 +93,23 @@ void alignment_sort_the_rows(Alignment *p_alignment, Alignment *alignment, stLis
     }
 }
 
+void alignment_filter_the_rows(Alignment *alignment, stList *prefixes_to_filter_by) {
+    Alignment_Row *row = alignment->row, **p_row = &(alignment->row);
+    while(row != NULL) {
+        if(alignment_row_get_closest_sequence_prefix(row, prefixes_to_filter_by) != -1) { // Filter the row
+            alignment->row_number--; // Reduce the number of rows
+            assert(alignment->row_number >= 0);
+            *p_row = row->n_row; // Update the previous link to point at the row after the current row
+            alignment_row_destruct(row); // Clean up the row
+            row = (*p_row); // Set row to point at the next row
+        }
+        else { // Keep the row, do nothing
+            p_row = &(row->n_row);  // Update pointers
+            row = row->n_row;
+        }
+    }
+}
+
 void alignment_show_only_lineage_differences(Alignment *alignment, char mask_char, stList *sequence_prefixes, stList *tree_nodes) {
     // First create map of tree nodes to bases
     stHash *tree_nodes_to_bases = stHash_construct2(NULL, (void (*)(void *))stList_destruct);
diff --git a/taffy/impl/taf.c b/taffy/impl/taf.c
index ad620b9..2d7ec77 100644
--- a/taffy/impl/taf.c
+++ b/taffy/impl/taf.c
@@ -84,11 +84,7 @@ static Alignment *parse_coordinates_and_establish_block(Alignment *p_block, stLi
             alignment->row_number--;
             Alignment_Row *r = *row;
             *row = r->n_row;
-            // Fix left pointer of row in previous block
-            assert(r->l_row != NULL);
-            r->l_row->r_row = NULL;
             // Now delete the row
-            r->n_row = NULL;
             alignment_row_destruct(r);
         } else if(op_type[0] == 'g') { // Is making a gap without the sequence specified
             int64_t gap_length = atol(stList_get(tokens, j++)); // Get the index of the affected row
@@ -276,6 +272,14 @@ Tag *taf_read_header(LI *li) {
     return tag;
 }
 
+Tag *taf_read_header_2(LI *li, bool *run_length_encode_bases) {
+    Tag *tag = taf_read_header(li);
+    Tag *t = tag_find(tag, (char *) "run_length_encode_bases");
+    *run_length_encode_bases = t != NULL && strcmp(t->value, "1") == 0;
+    return tag;
+}
+
+
 static void write_base(char base, int64_t base_count, LW *lw, bool run_length_encode_bases, bool color_bases) {
     if(base != '\0') {
         if(run_length_encode_bases) {
diff --git a/taffy/inc/taf.h b/taffy/inc/taf.h
index b7a0c6f..6afdec5 100644
--- a/taffy/inc/taf.h
+++ b/taffy/inc/taf.h
@@ -203,6 +203,11 @@ void paf_write_block(Alignment *alignment, LW *lw, bool all_to_all, bool cs_ciga
  */
 Tag *taf_read_header(LI *li);
 
+/*
+ * Read a taf header line, check if the run_length_encode_bases flag is set.
+ */
+Tag *taf_read_header_2(LI *li, bool *run_length_encode_bases);
+
 /*
  * Read a taf block - that is a column with column coordinates and all subsequent coordinate-less columns that
  * are considered to be part of the block.
@@ -333,5 +338,26 @@ int64_t alignment_row_get_closest_sequence_prefix(Alignment_Row *row, stList *pr
  */
 void alignment_sort_the_rows(Alignment *p_alignment, Alignment *alignment, stList *prefixes_to_sort_by);
 
+/*
+ * Removes any rows from the alignment whose sequence name prefix matches a string in the prefixes_to_filtet_by list
+ */
+void alignment_filter_the_rows(Alignment *alignment, stList *prefixes_to_filter_by);
+
+/*
+ * Load sequences in fasta files into a hash from sequence names to sequences
+ */
+stHash *load_sequences_from_fasta_files(char **seq_file_names, int64_t seq_file_number);
+
+/*
+ * Load sequences in hal file into memory.
+ */
+stSet *load_sequences_from_hal_file(char *hal_file, int *hal_handle);
+
+/*
+ * Add any gap strings between representing unaligned sequences between rows of alignment and p_alignment.
+ */
+void alignment_add_gap_strings(Alignment *p_alignment, Alignment *alignment, stHash *fastas, int hal_handle, stSet *hal_species,
+                               int64_t maximum_gap_string_length);
+
 #endif /* STTAF_H_ */
 


Test patch:
diff --git a/tests/evolverMammals.maf.mini.filtered b/tests/evolverMammals.maf.mini.filtered
new file mode 100644
index 0000000..ddebac0
--- /dev/null
+++ b/tests/evolverMammals.maf.mini.filtered
@@ -0,0 +1,42 @@
+##maf version=1 scoring=N/A
+
+a
+s	simCow_chr6.simCow.chr6	5045	50	-	602619	GTGAAGCTCAGTTGATGCTGGATTGGGAACTCATGAGTTAAGCTGTAAGC
+s	simDog_chr6.simDog.chr6	589129	50	+	593897	GTCAAGCTCAGTTGGTGCTGGATTAAGAATTCATGAGTTAGGCTGCAGTC
+s	simHuman_chr6.simHuman.chr6	597375	50	+	601863	GTCAAGCTCAGTAGATATTGGATTAGGAATTCATAAGTTAACCTGTAGCC
+s	simMouse_chr6.simMouse.chr6	630640	50	+	636262	GTCAAGCATTGTACATACTAGATTGGACATTCATGGATGACAATGTGACT
+s	simRat_chr6.simRat.chr6	642153	50	+	647215	GTCAAGCTCTGTAAATAGTAGATTGGACATTCATGGATGAAACTGTGCCT
+
+a
+s	simCow_chr6.simCow.chr6	5095	1	-	602619	T
+s	simDog_chr6.simDog.chr6	589179	1	+	593897	A
+s	simHuman_chr6.simHuman.chr6	597425	1	+	601863	A
+s	simMouse_chr6.simMouse.chr6	630691	1	+	636262	A
+s	simRat_chr6.simRat.chr6	642204	1	+	647215	A
+
+a
+s	simCow_chr6.simCow.chr6	5097	7	-	602619	ATATCCG
+s	simDog_chr6.simDog.chr6	589180	7	+	593897	ACATCTG
+s	simHuman_chr6.simHuman.chr6	597426	7	+	601863	ACATCGG
+s	simMouse_chr6.simMouse.chr6	630692	7	+	636262	ACAAAAA
+s	simRat_chr6.simRat.chr6	642205	7	+	647215	ACAAAAA
+
+a
+s	simCow_chr6.simCow.chr6	5104	29	-	602619	ATGTATTCTCTATAATAGTTTCTAATCTT
+s	simDog_chr6.simDog.chr6	589187	5	+	593897	ATGTA------------------------
+s	simHuman_chr6.simHuman.chr6	597433	29	+	601863	GTATGCTCTATCTAATATTTAGTACTCTT
+s	simMouse_chr6.simMouse.chr6	630701	19	+	636262	AAATCTACTTCCTATTATC----------
+s	simRat_chr6.simRat.chr6	642214	19	+	647215	AAATCTTCCTCCTATTATC----------
+
+a
+s	simCow_chr6.simCow.chr6	5133	13	-	602619	AGTCGATCTGGGA
+s	simDog_chr6.simDog.chr6	589192	3	+	593897	----------GGA
+s	simHuman_chr6.simHuman.chr6	597463	13	+	601863	AGTCGATGTGGGA
+
+a
+s	simCow_chr6.simCow.chr6	5146	18	-	602619	GTACAGTCTGAGTATGAG
+s	simDog_chr6.simDog.chr6	589195	18	+	593897	ATGCAGTTTGAGTATGAC
+s	simHuman_chr6.simHuman.chr6	597482	18	+	601863	GTAAAGTTTAATTATGTC
+s	simMouse_chr6.simMouse.chr6	630720	13	+	636262	----AATTTGATTAT-TC
+s	simRat_chr6.simRat.chr6	642233	8	+	647215	----AGTTTGAT------
+
diff --git a/tests/evolverMammals.maf.mini.filtered.normalized b/tests/evolverMammals.maf.mini.filtered.normalized
new file mode 100644
index 0000000..a18fecd
--- /dev/null
+++ b/tests/evolverMammals.maf.mini.filtered.normalized
@@ -0,0 +1,9 @@
+##maf version=1 scoring=N/A
+
+a
+s	simCow_chr6.simCow.chr6	5045	119	-	602619	GTGAAGCTCAGTTGATGCTGGATTGGGAACTCATGAGTTAAGCTGTAAGC-TGATATCCG--ATGTATTCTCTATAATAGTTTCTAATCTT-AGTCGATCTGGGA------GTACAGTCTGAGTATGAG
+s	simDog_chr6.simDog.chr6	589129	84	+	593897	GTCAAGCTCAGTTGGTGCTGGATTAAGAATTCATGAGTTAGGCTGCAGTC-A-ACATCTG--ATGTA-----------------------------------GGA------ATGCAGTTTGAGTATGAC
+s	simHuman_chr6.simHuman.chr6	597375	125	+	601863	GTCAAGCTCAGTAGATATTGGATTAGGAATTCATAAGTTAACCTGTAGCC-A-ACATCGG--GTATGCTCTATCTAATATTTAGTACTCTTTAGTCGATGTGGGAACAGGTGTAAAGTTTAATTATGTC
+s	simMouse_chr6.simMouse.chr6	630640	93	+	636262	GTCAAGCATTGTACATACTAGATTGGACATTCATGGATGACAATGTGACTGA-ACAAAAAAAAAATCTACTTCCTATTATC----------------------------------AATTTGATTAT-TC
+s	simRat_chr6.simRat.chr6	642153	88	+	647215	GTCAAGCTCTGTAAATAGTAGATTGGACATTCATGGATGAAACTGTGCCTGA-ACAAAAAAAAAATCTTCCTCCTATTATC----------------------------------AGTTTGAT------
+
diff --git a/tests/filter_file.txt b/tests/filter_file.txt
new file mode 100644
index 0000000..42b64b5
--- /dev/null
+++ b/tests/filter_file.txt
@@ -0,0 +1,2 @@
+Anc
+mr
diff --git a/tests/normalize_test.c b/tests/normalize_test.c
index c461f68..c55093b 100644
--- a/tests/normalize_test.c
+++ b/tests/normalize_test.c
@@ -110,7 +110,7 @@ static void test_maf_norm_to_maf(CuTest *testCase) {
     // Example maf file
     char *example_file = "./tests/evolverMammals.maf";
     char *output_file = "./tests/evolverMammals.maf.norm";
-    int i = st_system("./bin/taffy view -i %s | ./bin/taffy add-gap-bases ./tests/seqs/* | ./bin/taffy norm -k > %s",
+    int i = st_system("./bin/taffy view -i %s | ./bin/taffy norm -k -b ./tests/seqs/* > %s",
               example_file, output_file);
     CuAssertIntEquals(testCase, 0, i); // return value should be zero
 
@@ -188,11 +188,27 @@ static void test_dupe_filter(CuTest *testCase) {
     CuAssertIntEquals(testCase, 0, diff_ret); // return value should be zero if files sames
 }
 
+static void test_norm_pipeline(CuTest *testCase) {
+    /*
+     * Run taf_norm -d on really simple test maf to make sure that it drops a gap-causing dupe
+     */
+    char *input_file = "./tests/evolverMammals.maf.mini";
+    char *filter_file = "./tests/filter_file.txt";
+    char *output_file = "./tests/evolverMammals.maf.norm_pipeline";
+    int i = st_system("./bin/taffy view -i %s | ./bin/taffy sort -f %s | ./bin/taffy norm -b ./tests/seqs/* -k > %s",
+                      input_file, filter_file, output_file);
+    CuAssertIntEquals(testCase, 0, i); // return value should be zero
+    char *truth_file = "./tests/evolverMammals.maf.mini.filtered.normalized";
+    int diff_ret = st_system("diff %s %s", output_file, truth_file);
+    CuAssertIntEquals(testCase, 0, diff_ret); // return value should be zero if files sames
+}
+
 CuSuite* normalize_test_suite(void) {
     CuSuite* suite = CuSuiteNew();
     SUITE_ADD_TEST(suite, test_normalize);
     SUITE_ADD_TEST(suite, test_maf_norm);
     SUITE_ADD_TEST(suite, test_maf_norm_to_maf);
     SUITE_ADD_TEST(suite, test_dupe_filter);
+    SUITE_ADD_TEST(suite, test_norm_pipeline);
     return suite;
 }
diff --git a/tests/sort_test.c b/tests/sort_test.c
index b999c75..bd89d7c 100644
--- a/tests/sort_test.c
+++ b/tests/sort_test.c
@@ -18,8 +18,24 @@ static void test_sort(CuTest *testCase) {
     }
 }
 
+static void test_filter(CuTest *testCase) {
+    {
+        char *example_file = "./tests/evolverMammals.maf.mini";
+        char *output_file = "./tests/filter_test.maf.out";
+        char *truth_file = "./tests/evolverMammals.maf.mini.filtered";
+        char *filter_file = "./tests/filter_file.txt";
+        int i = st_system("./bin/taffy view -i %s | ./bin/taffy sort -f %s | ./bin/taffy view -m > %s",
+                          example_file, filter_file, output_file);
+        CuAssertIntEquals(testCase, 0, i); // return value should be zero
+        int diff_ret = st_system("diff %s %s", output_file, truth_file);
+        CuAssertIntEquals(testCase, 0, diff_ret); // return value should be zero if files sames
+        st_system("rm -f %s", output_file);
+    }
+}
+
 CuSuite* sort_test_suite(void) {
     CuSuite* suite = CuSuiteNew();
     SUITE_ADD_TEST(suite, test_sort);
+    SUITE_ADD_TEST(suite, test_filter);
     return suite;
 }