#!/usr/bin/env python #file test_split_libraries.py __author__ = "Rob Knight" __copyright__ = "Copyright 2011, The QIIME Project" #consider project name __credits__ = ["Rob Knight", "William Walters"] #remember to add yourself __license__ = "GPL" __version__ = "1.4.0" __maintainer__ = "William Walters" __email__ = "william.a.walters@colorado.edu" __status__ = "Release" from os.path import exists from StringIO import StringIO from numpy import array from shutil import rmtree from cogent import DNA from cogent.util.unit_test import TestCase, main from qiime.util import get_tmp_filename from cogent.util.misc import remove_files from qiime.split_libraries import ( expand_degeneracies, get_infile, count_mismatches, ok_mm_primer, check_map, fasta_ids, count_ambig, split_seq, primer_exceeds_mismatches, check_barcode, make_histograms, SeqQualBad, seq_exceeds_homopolymers, check_window_qual_scores, check_seqs, local_align_primer_seq, preprocess ) from qiime.parse import parse_qual_score from qiime.util import create_dir class FakeOutFile(object): def __init__(self): self.data = "" def write(self,s): self.data += s class TopLevelTests(TestCase): """Tests of top-level functions""" def setUp(self): """ """ self.in_seqs_variable_len_bc1 = in_seqs_variable_len_bc1 self.bc_map_variable_len_bc1 = bc_map_variable_len_bc1 self.primer_seq_lens_variable_len_bc1 = primer_seq_lens_variable_len_bc1 self.all_primers_variable_len_bc1 = all_primers_variable_len_bc1 self.expected_fasta_variable_len_bc1 = expected_fasta_variable_len_bc1 self.in_seqs_variable_len_bc2 = in_seqs_variable_len_bc1 self.bc_map_variable_len_bc2 = bc_map_variable_len_bc2 self.primer_seq_lens_variable_len_bc2 = primer_seq_lens_variable_len_bc2 self.all_primers_variable_len_bc2 = all_primers_variable_len_bc2 self.expected_fasta_variable_len_bc2 = expected_fasta_variable_len_bc1 self.in_seqs_fixed_len_bc1 = in_seqs_fixed_len_bc1 self.bc_map_fixed_len_bc1 = bc_map_fixed_len_bc1 self.primer_seq_lens_fixed_len_bc1 = primer_seq_lens_fixed_len_bc1 self.all_primers_fixed_len_bc1 = all_primers_fixed_len_bc1 self.expected_fasta_fixed_len_bc1 = expected_fasta_fixed_len_bc1 self.in_seqs_fixed_len_bc2 = in_seqs_fixed_len_bc1 self.bc_map_fixed_len_bc2 = bc_map_fixed_len_bc1 self.primer_seq_lens_fixed_len_bc2 = primer_seq_lens_fixed_len_bc1 self.all_primers_fixed_len_bc2 = all_primers_fixed_len_bc1 self.expected_fasta_fixed_len_bc2 = expected_fasta_fixed_len_bc2 self.expected_fasta_fixed_len_bc1_no_primers = \ expected_fasta_fixed_len_bc1_no_primers self.reverse_primers_fixed_len_bc1 =\ reverse_primers_fixed_len_bc1 self.in_seqs_reverse_primers =\ in_seqs_reverse_primers self.expected_in_seqs_reverse_primers =\ expected_in_seqs_reverse_primers self.in_seqs_reverse_primers_mismatch =\ in_seqs_reverse_primers_mismatch self.expected_in_seqs_reverse_primers_mismatch =\ expected_in_seqs_reverse_primers_mismatch self.expected_in_seqs_reverse_primers_full_remove =\ expected_in_seqs_reverse_primers_full_remove self.expected_in_seqs_reverse_primers_mismatch_allowed =\ expected_in_seqs_reverse_primers_mismatch_allowed self.expected_fasta_fixed_len_bc1_sliding_window =\ expected_fasta_fixed_len_bc1_sliding_window self.in_seqs_fixed_len_extra_bc = in_seqs_fixed_len_extra_bc self.expected_fasta_extra_bc = expected_fasta_extra_bc self.expected_fasta_mad = expected_fasta_mad self.sample_seqs_fna_file = sample_seqs_fna_file self.expected_fasta_fixed_added_demultiplex =\ expected_fasta_fixed_added_demultiplex self.in_seqs_added_demultiplex = in_seqs_added_demultiplex self.bc_map_added_demultiplex = bc_map_added_demultiplex self.bc_map_added_demultiplex_group = bc_map_added_demultiplex_group self.expected_fasta_added_demultiplex_group =\ expected_fasta_added_demultiplex_group self.in_seqs_fixed_len_bc1_qual_scores =\ in_seqs_fixed_len_bc1_qual_scores self.sample_mapping = sample_mapping self.sample_fasta_file = get_tmp_filename(prefix = "sample_seqs_", suffix = ".fasta") seq_file = open(self.sample_fasta_file, 'w') seq_file.write("\n".join(self.in_seqs_fixed_len_extra_bc)) seq_file.close() self.sample_qual_file = get_tmp_filename(prefix = "sample_qual_", suffix = ".qual") qual_file = open(self.sample_qual_file, "w") qual_file.write("\n".join(self.in_seqs_fixed_len_bc1_qual_scores)) qual_file.close() self.sample_mapping_file = get_tmp_filename(prefix = "sample_mapping_", suffix = ".txt") map_file = open(self.sample_mapping_file, "w") map_file.write(self.sample_mapping) map_file.close() self.output_dir = get_tmp_filename(prefix = "split_libraries_", suffix = "/") create_dir(self.output_dir) self._files_to_remove = \ [self.sample_fasta_file, self.sample_qual_file, self.sample_mapping_file] def tearDown(self): if self._files_to_remove: remove_files(self._files_to_remove) if exists(self.output_dir): rmtree(self.output_dir) def test_check_window_qual_scores(self): """check_window_qual_scores returns False, index if window below qual threshold.""" scores1 = [8,8,8,8,8,8,8,2,2,2,2,2] self.assertEqual(check_window_qual_scores(scores1, 5, 5), (False, 5)) self.assertEqual(check_window_qual_scores(scores1, 10, 5), (True, 2)) # windowsize larger than qual score list works self.assertEqual(check_window_qual_scores(scores1, 100, 5), (True, 0)) self.assertEqual(check_window_qual_scores([], 5, 1), True) #check each base in its own window self.assertEqual(check_window_qual_scores(scores1, 1, 2), (True, 11)) self.assertEqual(check_window_qual_scores(scores1, 1, 5), (False, 7)) def test_expand_degeneracies(self): """generate_possibilities should make possible strings""" self.assertEqual(expand_degeneracies(['ACG']), ['ACG']) self.assertEqual(expand_degeneracies(['RGY']), ['AGT', 'AGC', 'GGT', 'GGC']) def test_get_infile(self): """get_infile should return filehandle""" pass #not practically testable, but obvious file I/O def test_count_mismatches(self): """count_mismatches should count mismatches correctly""" self.assertEqual(count_mismatches('GG','GG',10), 0) self.assertEqual(count_mismatches('GGG','AAA',10), 3) self.assertEqual(count_mismatches('GGG','AAA',1), 2) def test_ok_mm_primer(self): """ok_mm_primer should test that primer is within max mismatches""" primers = ['AAAA', 'GGGG'] self.assertEqual(ok_mm_primer('AAAA', primers, 0), True) self.assertEqual(ok_mm_primer('AAAA', primers, 3), True) self.assertEqual(ok_mm_primer('CCCC', primers, 0), False) self.assertEqual(ok_mm_primer('CCCA', primers, 3), True) self.assertEqual(ok_mm_primer('CCCA', primers, 2), False) self.assertEqual(ok_mm_primer('CCGG', primers, 2), True) self.assertEqual(ok_mm_primer('CCGA', primers, 2), False) def test_check_map(self): """check_map should return valid barcodes as expected""" s = """#SampleID\tBarcodeSequence\tLinkerPrimerSequence\tX\tDescription #fake data x\tAA\tAC\t3\tsample_x y\t"AC"\tAC\t4\t"sample_y" z\tGG\tGC\t5\tsample_z""" f = StringIO(s) f.name='test.xls' headers, id_map, barcode_to_sample_id, warnings, errors, \ primer_seqs_lens, all_primers = check_map(f, disable_primer_check=False) self.assertEqual(barcode_to_sample_id, {'AA':'x','AC':'y','GG':'z'}) self.assertEqual(errors, []) self.assertEqual(warnings, []) def test_check_map_primer_pool(self): """check_map should handle primer pools as expected""" s = """#SampleID\tBarcodeSequence\tLinkerPrimerSequence\tX\tDescription #fake data x\tAA\tAC\t3\tsample_x y\t"AC"\tAT,DC\t4\t"sample_y" z\tGG\tGC\t5\tsample_z""" f = StringIO(s) f.name='test.xls' headers, id_map, barcode_to_sample_id, warnings, errors, \ primer_seqs_lens, all_primers = check_map(f, disable_primer_check=False) self.assertEqual(barcode_to_sample_id, {'AA':'x','AC':'y','GG':'z'}) self.assertEqual(errors, []) self.assertEqual(warnings, []) # Returns all possible primers with lengths associated. expected_all_primers = {'AC': 2, 'GC': 2, 'AT': 2, 'TC': 2} self.assertEqual(all_primers, expected_all_primers) # Returns all primers associated with each barcode. expected_primer_seqs_lens = {'AA': {'AC': 2}, 'GG': {'GC': 2}, 'AC': {'AC': 2, 'GC': 2, 'AT': 2, 'TC': 2}} self.assertEqual(primer_seqs_lens, expected_primer_seqs_lens) def test_fasta_ids(self): """fasta_ids should return list of ids in fasta files, no dups""" first = StringIO('>x\nACT\n>y\nAAA') first_copy = StringIO('>x\nACT\n>y\nAAA') second = StringIO('>a\nGGG\n>b\nCCC') self.assertEqual(fasta_ids([first, second]), set(['x','y','a','b'])) first.seek(0) #need to reset so we can read it again self.assertRaises(ValueError, fasta_ids, [first,first_copy]) def test_count_ambig(self): """count_ambig should count ambiguous bases in seq""" s = 'ACC' s2 = 'RNY' s3 = 'NA' self.assertEqual(count_ambig(s), 0) self.assertEqual(count_ambig(s2), 3) self.assertEqual(count_ambig(s3), 1) self.assertEqual(count_ambig(''), 0) def test_split_seq(self): """split_seq should split seq into pieces""" seq = 'AAAACCCCCGTGTGTGT' barcode, primer, remainder = split_seq(seq, 4, 5) self.assertEqual(barcode, 'AAAA') self.assertEqual(primer, 'CCCCC') self.assertEqual(remainder, 'GTGTGTGT') def test_primer_exceeds_mismatches(self): """primer_exceeds_mismatches returns True if too many mismatches""" primers = ['AAAA', 'TTTT'] exact = 'AAAA' mismatch_ok = 'AAAT' mismatch_bad = 'GGGG' self.assertEqual(primer_exceeds_mismatches(exact, primers, 0), False) self.assertEqual(primer_exceeds_mismatches(mismatch_ok, primers, 1), False) self.assertEqual(primer_exceeds_mismatches(mismatch_bad, primers, 2), True) # Tests for local alignment functions def test_local_align_primer_seq_fwd_rev_match(self): "local_align function can handle fwd/rev primers with no mismatches" # forward primer primer = DNA.makeSequence('TAGC',Name='F5') seq = 'TAGC' # mismatch_count, hit_start expected = (0,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) primer = DNA.makeSequence('TAGC',Name='F5') seq = 'TAGCCCCC' # mismatch_count, hit_start expected = (0,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) primer = DNA.makeSequence('TAGC',Name='F5') seq = 'CCCTAGCCCCC' # mismatch_count, hit_start expected = (0,3) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) # different length primer primer = DNA.makeSequence('GTTTAGC',Name='F5') seq = 'GTTTAGC' # mismatch_count, hit_start expected = (0,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) primer = DNA.makeSequence('GCTC',Name='R5') seq = 'TAGCCCCC' # mismatch_count, hit_start expected = (1,2) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) primer = DNA.makeSequence('GCTA',Name='R5') seq = 'CCCTAGCCCCC' # mismatch_count, hit_start expected = (1,1) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) def test_local_align_primer_seq_fwd_rev_match_ambig(self): "local_align function can handle fwd/rev primers with ambigs" primer = DNA.makeSequence('TASC',Name='F5') seq = 'TAGC' # primer_hit, target, mismatch_count, hit_start expected = (0,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) def test_local_align_primer_seq_mm(self): "local_align function can handle fwd/rev primers with mismatches" # forward primer primer = DNA.makeSequence('AAAAACTTTTT',Name='F5') seq = 'AAAAAGTTTTT' # mismatch_count, hit_start expected = (1,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) # forward primer primer = DNA.makeSequence('AAAACCTTTTT',Name='F5') seq = 'AAAAAGTTTTT' # mismatch_count, hit_start expected = (2,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) def test_local_align_primer_seq_indels_middle(self): "local_align function can handle fwd/rev primers with indels in middle of seq" # Insertion in target sequence primer = DNA.makeSequence('CGAATCGCTATCG',Name='F5') seq = 'CGAATCTGCTATCG' # mismatch count, hit_start expected = (1,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) # Deletion in target sequence primer = DNA.makeSequence('CGAATCGCTATCG',Name='F5') seq = 'CGAATGCTATCG' # mismatch_count, hit_start expected = (1,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) def test_local_align_primer_seq_multiple_mismatch_indel(self): "local_align function can handle fwd/rev primers with indels and mismatches" # multiple insertions primer = DNA.makeSequence('ATCGGGCGATCATT',Name='F5') seq = 'ATCGGGTTCGATCATT' # mismatch_count, hit_start expected = (2,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) # two deletions primer = DNA.makeSequence('ACGGTACAGTGG',Name='F5') seq = 'ACGGCAGTGG' # mismatch_count, hit_start expected = (2,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) # deletion and mismatch primer = DNA.makeSequence('CATCGTCGATCA',Name='F5') seq = 'CCTCGTGATCA' # mismatch_count, hit_start expected = (2,0) actual = local_align_primer_seq(primer,seq) self.assertEqual(actual,expected) def test_seq_exceeds_homopolymers(self): """seq_exceeds_homopolymers returns True if too many homopolymers""" self.assertEqual(seq_exceeds_homopolymers('AAACGA',3), False) self.assertEqual(seq_exceeds_homopolymers('AAACGA',2), True) self.assertEqual(seq_exceeds_homopolymers('AAACGA',1), True) self.assertEqual(seq_exceeds_homopolymers('AAACGATTTT',3), True) def test_check_barcode(self): """check_barcode should return False if barcode ok, True otherwise""" self.assertEqual(check_barcode('AA', None, ['AA']), (False, 'AA', \ False)) self.assertEqual(check_barcode('GCATCGTCCACA', 'golay_12', ['GCATCGTCAACA']), (2, 'GCATCGTCAACA', True)) # num errors for golay code is currently in bits self.assertEqual(check_barcode('GGTT', 4, ['TTTT']), (2, 'TTTT', True)) def test_make_histograms(self): """make_histograms should make correct histograms""" pre_lengths = [100, 110, 105, 130, 135] post_lengths = [130, 135] pre_hist, post_hist, bin_edges = \ make_histograms(pre_lengths, post_lengths) self.assertEqual(pre_hist, array([2,1,0,2])) self.assertEqual(post_hist, array([0,0,0,2])) self.assertEqual(bin_edges, array([100,110,120,130,140])) def test_check_seqs_sliding_window(self): """check_seqs handles sliding window truncations/removal """ in_seqs = self.in_seqs_fixed_len_bc1 bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_fixed_len_bc1_sliding_window out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings=parse_qual_score(self.in_seqs_fixed_len_bc1_qual_scores), filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False, qual_score_window=5, discard_bad_windows=False, min_qual_score=25, min_seq_len=200) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_variable_len_bc(self): """check_seqs handles variable length barcodes """ # Simple test with variable length primers in_seqs = self.in_seqs_variable_len_bc1 bc_map = self.bc_map_variable_len_bc1 primer_seq_lens = self.primer_seq_lens_variable_len_bc1 all_primers = self.all_primers_variable_len_bc1 expected = self.expected_fasta_variable_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=None, keep_primer=False, keep_barcode=False, barcode_type="variable_length", max_bc_errors=0, retain_unassigned_reads=False, attempt_bc_correction=False, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Second test, includes truncated form of one of the barcodes-the # longest barcode should be found first in_seqs = self.in_seqs_variable_len_bc2 bc_map = self.bc_map_variable_len_bc2 primer_seq_lens = self.primer_seq_lens_variable_len_bc2 all_primers = self.all_primers_variable_len_bc2 expected = self.expected_fasta_variable_len_bc2 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=None, keep_primer=False, keep_barcode=False, barcode_type="variable_length", max_bc_errors=0, retain_unassigned_reads=False, attempt_bc_correction=False, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_fixed_len_bc(self): """check_seqs handles fixed length barcodes """ # Third test, fixed length barcodes, fixed length primers + one # degenerate test. Should correct one of the passed barcodes in_seqs = self.in_seqs_fixed_len_bc1 bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Fourth test-set max_bc_errors to 0, and allow some primer mismatches in_seqs = self.in_seqs_fixed_len_bc2 bc_map = self.bc_map_fixed_len_bc2 primer_seq_lens = self.primer_seq_lens_fixed_len_bc2 all_primers = self.all_primers_fixed_len_bc2 expected = self.expected_fasta_fixed_len_bc2 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=0.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=1, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_no_primers(self): """check_seqs handles disabled primers """ # Fifth test, no primers, fixed length barcodes # Should correct one of the passed barcodes in_seqs = self.in_seqs_fixed_len_bc1 bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = {} all_primers = {} expected = self.expected_fasta_fixed_len_bc1_no_primers out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=True, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_reverse_primers(self): """check_seqs handles truncating reverse primers """ # Initial test, should truncate all seqs in_seqs = self.in_seqs_reverse_primers bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_in_seqs_reverse_primers rev_primers_test = self.reverse_primers_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'truncate_only', rev_primers = rev_primers_test, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Second test with a mismatch in seq a, should not find reverse primer # and will write out entire sequence. in_seqs = self.in_seqs_reverse_primers_mismatch bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_in_seqs_reverse_primers_mismatch rev_primers_test = self.reverse_primers_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'truncate_only', rev_primers = rev_primers_test, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # With reverse_primer_mismatches allowed set to 1, # should restore truncation. in_seqs = self.in_seqs_reverse_primers_mismatch bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_in_seqs_reverse_primers_mismatch_allowed rev_primers_test = self.reverse_primers_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'truncate_only', rev_primers = rev_primers_test, qual_out = False, reverse_primer_mismatches=1) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Testing truncate_remove, which should not write sequences where # the reverse primer is not found in_seqs = self.in_seqs_reverse_primers bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_in_seqs_reverse_primers_full_remove rev_primers_test = self.reverse_primers_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'truncate_remove', rev_primers = rev_primers_test, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Testing truncate_remove, with reverse_primer_mismatches set to 1 # should allow all 4 sequences to be written, truncated in_seqs = self.in_seqs_reverse_primers_mismatch bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_in_seqs_reverse_primers_mismatch_allowed rev_primers_test = self.reverse_primers_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=1, disable_primer_check=False, reverse_primers = 'truncate_remove', rev_primers = rev_primers_test, qual_out = False, reverse_primer_mismatches=1) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_qual_out(self): """ check_seqs handles optional quality output file """ in_seqs = self.in_seqs_fixed_len_bc1 bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = expected_qual_fixed_len_bc1 out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) qual_out_f = FakeOutFile() actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings=parse_qual_score(self.in_seqs_fixed_len_bc1_qual_scores), filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = qual_out_f) self.assertEqual(qual_out_f.data,expected) def test_check_seqs_retain_unassigned_reads(self): """ check_seqs handles retaining Unassigned reads """ in_seqs = self.in_seqs_fixed_len_extra_bc bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_extra_bc out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=True, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_median_abs_dev(self): """ check_seqs handles median absolute deviation calculations """ in_seqs = self.in_seqs_fixed_len_bc1 bc_map = self.bc_map_fixed_len_bc1 primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_mad out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False, median_length_filtering=2.0) out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_check_seqs_added_demultiplex(self): """check_seqs handles added demultiplex field""" # Test added demultiplex for the run_prefix in_seqs = self.in_seqs_added_demultiplex bc_map = self.bc_map_added_demultiplex primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_fixed_added_demultiplex out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False, added_demultiplex_field='run_prefix') out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) # Demultiplex by the 'group' in the fasta label in_seqs = self.in_seqs_added_demultiplex bc_map = self.bc_map_added_demultiplex_group primer_seq_lens = self.primer_seq_lens_fixed_len_bc1 all_primers = self.all_primers_fixed_len_bc1 expected = self.expected_fasta_added_demultiplex_group out_f = open(get_tmp_filename(prefix = "sample_seqs_", suffix = ".fna.tmp"), "w") self._files_to_remove.append(out_f.name.replace('.tmp','')) actual = check_seqs( fasta_out=out_f, fasta_files = [in_seqs], starting_ix=0, valid_map = bc_map, qual_mappings={}, filters=[], barcode_len=12, keep_primer=False, keep_barcode=False, barcode_type="golay_12", max_bc_errors=1.5, retain_unassigned_reads=False, attempt_bc_correction=True, primer_seqs_lens=primer_seq_lens, all_primers=all_primers, max_primer_mm=0, disable_primer_check=False, reverse_primers = 'disable', rev_primers = {}, qual_out = False, added_demultiplex_field='group') out_f = open(out_f.name.replace('.tmp',''), "U") actual_results = '\n'.join([line.strip() for line in out_f]) self.assertEqual(actual_results, expected) def test_preprocess(self): """ Overall module functionality test """ # Should discard all reads due to sequence length being too short fasta_files = [self.sample_fasta_file] qual_files = [self.sample_qual_file] mapping_file = self.sample_mapping_file barcode_type="golay_12" min_seq_len=200 max_seq_len=1000 min_qual_score=25 starting_ix=1 keep_primer=False max_ambig=0 max_primer_mm=1 trim_seq_len=True dir_prefix=self.output_dir max_bc_errors=2 max_homopolymer=4 retain_unassigned_reads=False keep_barcode=False attempt_bc_correction=True qual_score_window=0 disable_primer_check=False reverse_primers='disable' record_qual_scores=False discard_bad_windows=False median_length_filtering=None added_demultiplex_field=None preprocess(fasta_files, qual_files, mapping_file, barcode_type, min_seq_len, max_seq_len, min_qual_score, starting_ix, keep_primer, max_ambig, max_primer_mm, trim_seq_len, dir_prefix, max_bc_errors, max_homopolymer, retain_unassigned_reads, keep_barcode, attempt_bc_correction, qual_score_window, disable_primer_check, reverse_primers, record_qual_scores, discard_bad_windows, median_length_filtering, added_demultiplex_field) output_seqs = open(dir_prefix + "seqs.fna", "U") output_log = open(dir_prefix + "split_library_log.txt", "U") output_histograms = open(dir_prefix + "histograms.txt", "U") actual_seqs = [line for line in output_seqs] actual_log = [line for line in output_log] actual_histograms = [line for line in output_histograms] expected_seqs = [] expected_log = ['Number raw input seqs\t6\n', '\n', 'Length outside bounds of 200 and 1000\t6\n', 'Num ambiguous bases exceeds limit of 0\t0\n', 'Missing Qual Score\t0\n', 'Mean qual score below minimum of 25\t0\n', 'Max homopolymer run exceeds limit of 4\t0\n', 'Num mismatches in primer exceeds limit of 1: 0\n', '\n', 'Raw len min/max/avg\t25.0/35.0/31.2\n', '\n', 'Barcodes corrected/not\t0/0\n', 'Uncorrected barcodes will not be written to the output fasta file.\n', 'Corrected barcodes will be written with the appropriate barcode category.\n', 'Corrected but unassigned sequences will not be written unless --retain_unassigned_reads is enabled.\n', '\n', 'Total valid barcodes that are not in mapping file\t0\n', 'Sequences associated with valid barcodes that are not in the mapping file will not be written.\n', '\n', 'Barcodes in mapping file\n', 'Sample\tSequence Count\tBarcode\n', 's2\t0\tAGAGTCCTGAGC\n', 's1\t0\tACACATGTCTAC\n', 's3\t0\tAACTGTGCGTAC\n', '\n', 'Total number seqs written\t0'] expected_histograms = ['Length\tBefore\tAfter\n', '20\t2\t0\n', '30\t4\t0'] self.assertEqual(actual_seqs, expected_seqs) self.assertEqual(actual_log, expected_log) self.assertEqual(actual_histograms, expected_histograms) # With minimal length at 5, should retain 4 sequences fasta_files = [self.sample_fasta_file] qual_files = [self.sample_qual_file] mapping_file = self.sample_mapping_file barcode_type="golay_12" min_seq_len=5 max_seq_len=1000 min_qual_score=25 starting_ix=1 keep_primer=False max_ambig=0 max_primer_mm=0 trim_seq_len=False dir_prefix=self.output_dir max_bc_errors=2 max_homopolymer=4 retain_unassigned_reads=False keep_barcode=False attempt_bc_correction=True qual_score_window=0 disable_primer_check=False reverse_primers='disable' record_qual_scores=False discard_bad_windows=False median_length_filtering=None added_demultiplex_field=None preprocess(fasta_files, qual_files, mapping_file, barcode_type, min_seq_len, max_seq_len, min_qual_score, starting_ix, keep_primer, max_ambig, max_primer_mm, trim_seq_len, dir_prefix, max_bc_errors, max_homopolymer, retain_unassigned_reads, keep_barcode, attempt_bc_correction, qual_score_window, disable_primer_check, reverse_primers, record_qual_scores, discard_bad_windows, median_length_filtering, added_demultiplex_field) output_seqs = open(dir_prefix + "seqs.fna", "U") output_log = open(dir_prefix + "split_library_log.txt", "U") output_histograms = open(dir_prefix + "histograms.txt", "U") actual_seqs = [line for line in output_seqs] actual_log = [line for line in output_log] actual_histograms = [line for line in output_histograms] expected_seqs = ['>s1_1 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\n', 'CCCTTATATATATAT\n', '>s2_2 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\n', 'CCCTTTCCA\n', '>s3_3 c orig_bc=AACTGTGCGTAC new_bc=AACTGTGCGTAC bc_diffs=0\n', 'AACCGGCCGGTT\n', '>s1_4 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\n', 'CCCTTACTATATAT\n'] expected_log = ['Number raw input seqs\t6\n', '\n', 'Length outside bounds of 5 and 1000\t0\n', 'Num ambiguous bases exceeds limit of 0\t0\n', 'Missing Qual Score\t0\n', 'Mean qual score below minimum of 25\t0\n', 'Max homopolymer run exceeds limit of 4\t0\n', 'Num mismatches in primer exceeds limit of 0: 2\n', '\n', 'Raw len min/max/avg\t25.0/35.0/31.2\n', 'Wrote len min/max/avg\t29.0/35.0/32.5\n', '\n', 'Barcodes corrected/not\t1/0\n', 'Uncorrected barcodes will not be written to the output fasta file.\n', 'Corrected barcodes will be written with the appropriate barcode category.\n', 'Corrected but unassigned sequences will not be written unless --retain_unassigned_reads is enabled.\n', '\n', 'Total valid barcodes that are not in mapping file\t0\n', 'Sequences associated with valid barcodes that are not in the mapping file will not be written.\n', '\n', 'Barcodes in mapping file\n', 'Num Samples\t3\n', 'Sample ct min/max/mean: 1 / 2 / 1.33\n', 'Sample\tSequence Count\tBarcode\n', 's1\t2\tACACATGTCTAC\n', 's2\t1\tAGAGTCCTGAGC\n', 's3\t1\tAACTGTGCGTAC\n', '\n', 'Total number seqs written\t4'] expected_histograms = ['Length\tBefore\tAfter\n', '20\t2\t1\n', '30\t4\t3'] self.assertEqual(actual_seqs, expected_seqs) self.assertEqual(actual_log, expected_log) self.assertEqual(actual_histograms, expected_histograms) # Added sliding window should discard read "b" fasta_files = [self.sample_fasta_file] qual_files = [self.sample_qual_file] mapping_file = self.sample_mapping_file barcode_type="golay_12" min_seq_len=5 max_seq_len=1000 min_qual_score=22 starting_ix=1 keep_primer=False max_ambig=0 max_primer_mm=0 trim_seq_len=False dir_prefix=self.output_dir max_bc_errors=2 max_homopolymer=4 retain_unassigned_reads=False keep_barcode=False attempt_bc_correction=True qual_score_window=3 disable_primer_check=False reverse_primers='disable' record_qual_scores=False discard_bad_windows=True median_length_filtering=None added_demultiplex_field=None reverse_primer_mismatches=0 preprocess(fasta_files, qual_files, mapping_file, barcode_type, min_seq_len, max_seq_len, min_qual_score, starting_ix, keep_primer, max_ambig, max_primer_mm, trim_seq_len, dir_prefix, max_bc_errors, max_homopolymer, retain_unassigned_reads, keep_barcode, attempt_bc_correction, qual_score_window, disable_primer_check, reverse_primers, reverse_primer_mismatches, record_qual_scores, discard_bad_windows, median_length_filtering, added_demultiplex_field) output_seqs = open(dir_prefix + "seqs.fna", "U") output_log = open(dir_prefix + "split_library_log.txt", "U") output_histograms = open(dir_prefix + "histograms.txt", "U") actual_seqs = [line for line in output_seqs] actual_log = [line for line in output_log] actual_histograms = [line for line in output_histograms] expected_seqs = ['>s1_1 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\n', 'CCCTTATATATATAT\n', '>s3_2 c orig_bc=AACTGTGCGTAC new_bc=AACTGTGCGTAC bc_diffs=0\n', 'AACCGGCCGGTT\n', '>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\n', 'CCCTTACTATATAT\n'] expected_log = ['Number raw input seqs\t6\n', '\n', 'Length outside bounds of 5 and 1000\t0\n', 'Num ambiguous bases exceeds limit of 0\t0\n', 'Missing Qual Score\t0\n', 'Mean qual score below minimum of 22\t0\n', 'Max homopolymer run exceeds limit of 4\t0\n', 'Num mismatches in primer exceeds limit of 0: 2\n', '\n', 'Size of quality score window, in base pairs: 3\n', 'Number of sequences where a low quality score window was detected: 1\n', 'Sequences with a low quality score were not written, -g option enabled.\n', '\n', 'Raw len min/max/avg\t25.0/35.0/31.2\n', 'Wrote len min/max/avg\t32.0/35.0/33.7\n', '\n', 'Barcodes corrected/not\t1/0\n', 'Uncorrected barcodes will not be written to the output fasta file.\n', 'Corrected barcodes will be written with the appropriate barcode category.\n', 'Corrected but unassigned sequences will not be written unless --retain_unassigned_reads is enabled.\n', '\n', 'Total valid barcodes that are not in mapping file\t0\n', 'Sequences associated with valid barcodes that are not in the mapping file will not be written.\n', '\n', 'Barcodes in mapping file\n', 'Num Samples\t2\n', 'Sample ct min/max/mean: 1 / 2 / 1.50\n', 'Sample\tSequence Count\tBarcode\n', 's1\t2\tACACATGTCTAC\n', 's3\t1\tAACTGTGCGTAC\n', 's2\t0\tAGAGTCCTGAGC\n', '\n', 'Total number seqs written\t3'] expected_histograms = ['Length\tBefore\tAfter\n', '20\t2\t0\n', '30\t4\t3'] self.assertEqual(actual_seqs, expected_seqs) self.assertEqual(actual_log, expected_log) self.assertEqual(actual_histograms, expected_histograms) in_seqs_variable_len_bc1 = """>a ACCGGTCCGGACCCTTATATATATAT >b AGGAGTCCGGACCCTTTCCA >c ATTAACCCGGAAACCGGCCGGTT >d ACCGGTCCGGACCCTTACTATATAT >e TTTTGTCCGGACCCTTACTATATAT >d_primer_error ACCTGGTCCGGACCCTTACTATATAT """.split() bc_map_variable_len_bc1 = {'ACC':'s1','AGGA':'s2','ATTA':'s3'} primer_seq_lens_variable_len_bc1 = {'ACC':{'GGTCCGGA':8}, 'AGGA':{'GTCCGGA':7}, 'ATTA':{'ACCCGGA':7}} all_primers_variable_len_bc1 = {'GGTCCGGA':8,'GTCCGGA':7,'ACCCGGA':7} expected_fasta_variable_len_bc1 = """>s1_0 a orig_bc=ACC new_bc=ACC bc_diffs=0 CCCTTATATATATAT >s2_1 b orig_bc=AGGA new_bc=AGGA bc_diffs=0 CCCTTTCCA >s3_2 c orig_bc=ATTA new_bc=ATTA bc_diffs=0 AACCGGCCGGTT >s1_3 d orig_bc=ACC new_bc=ACC bc_diffs=0 CCCTTACTATATAT""" bc_map_variable_len_bc2 = {'ACC':'s1','AGGA':'s2','ATTA':'s3','AGG':'s4'} primer_seq_lens_variable_len_bc2 = {'ACC':{'GGTCCGGA':8}, 'AGGA':{'GTCCGGA':7}, 'ATTA':{'ACCCGGA':7}, 'AGG':{'AGTCCGGA':8}} all_primers_variable_len_bc2 = {'GGTCCGGA':8,'GTCCGGA':7,'ACCCGGA':7, 'AGTCCGGA':8} # Fixed barcode test data in_seqs_fixed_len_bc1 = """>a ACACATGTCTACGGTCCGGACCCTTATATATATAT >b AGAGTCCTGAGCGGTCCGGACCCTTTCCA >c AATCGTGACTCGGGTCTGGAAACCGGCCGGTT >d ACTCATGTCTACGGTCCGGACCCTTACTATATAT >e_no_barcode_match TTTTGTCCGGACCCTTACTATATAT >d_primer_error AGAGTCCTGAGCGGTCCGGTACGTTTACTGGA """.split('\n') sample_mapping = """#SampleID\tBarcodeSequence\tLinkerPrimerSequence\tTreatment\tDescription # Test mapping for split_libraries.py unit tests s1\tACACATGTCTAC\tGGTCCGGA\tControl\ts1_mouse s2\tAGAGTCCTGAGC\tGGTCCGGA\tControl\ts2_mouse s3\tAACTGTGCGTAC\tGGTCYGGA\tFasted\ts3_mouse """ # Fixed barcode test data, with one valid BC not in mapping data in_seqs_fixed_len_extra_bc = """>a ACACATGTCTACGGTCCGGACCCTTATATATATAT >b AGAGTCCTGAGCGGTCCGGACCCTTTCCA >c AACTGTGCGTACGGTCTGGAAACCGGCCGGTT >d ACTCATGTCTACGGTCCGGACCCTTACTATATAT >e_no_barcode_match TTTTGTCCGGACCCTTACTATATAT >d_primer_error AGAGTCCTGAGCGGTCCGGTACGTTTACTGGA """.split('\n') expected_fasta_extra_bc = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nCCCTTATATATATAT\n>s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\nCCCTTTCCA\n>Unassigned_2 c orig_bc=AACTGTGCGTAC new_bc=AACTGTGCGTAC bc_diffs=0\nAACCGGCCGGTT\n>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\nCCCTTACTATATAT""" in_seqs_fixed_len_bc1_qual_scores = """>a 37 37 37 37 37 37 37 37 37 37 37 37 37 37 37 37 37 36 36 33 33 33 36 37 37 37 37 37 37 40 40 40 39 39 38 >b 35 31 31 23 23 23 31 21 21 21 35 35 37 37 37 36 36 36 36 36 36 37 37 37 37 37 37 37 37 >c 37 37 37 37 37 37 37 37 37 37 37 37 37 37 36 32 32 32 36 37 35 32 32 32 32 32 32 32 32 36 37 37 >d 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37 37 37 37 37 >e_no_barcode_match 23 23 23 31 21 21 21 35 35 37 37 37 36 36 36 36 36 36 37 37 37 37 37 37 37 >d_primer_error 31 31 23 23 23 31 21 21 21 35 35 37 37 37 36 36 36 36 36 36 37 37 37 37 37 37 37 37 37 37 37 37""".split('\n') # These test data have equal length barcodes, one degenerate primer bc_map_fixed_len_bc1 = {'ACACATGTCTAC':'s1','AGAGTCCTGAGC':'s2', 'AATCGTGACTCG':'s3'} primer_seq_lens_fixed_len_bc1 = {'ACACATGTCTAC':{'GGTCCGGA':8}, 'AGAGTCCTGAGC':{'GGTCCGGA':8}, 'AATCGTGACTCG':{'GGTCCGGA':8,'GGTCTGGA':8}} all_primers_fixed_len_bc1 = {'GGTCCGGA':8,'GGTCTGGA':8} # added_demultiplex barcode test data in_seqs_added_demultiplex = """>mintberrycrunch group=coonandfriends AATCGTGACTCGGGTCCGGACCCTTATATATATAT >mysterion group=coonandfriends AGAGTCCTGAGCGGTCCGGACCCTTTCCA >mansquito group=professorchaos AATCGTGACTCGGGTCTGGAAACCGGCCGGTT >mintberrycrunch group=coonandfriends AATCGTGACTTGGGTCCGGACCCTTACTATATAT >e_no_barcode_match group=professorchaos TTTTGTCCGGACCCTTACTATATAT >d_primer_error group=coonandfriends AGAGTCCTGAGCGGTCCGGTACGTTTACTGGA """.split('\n') # Modified barcode map to test added_demultiplex option bc_map_added_demultiplex = {'AATCGTGACTCG,mintberry':'s1', 'AGAGTCCTGAGC,myst':'s2', 'AATCGTGACTCG,mansq':'s3'} # Modified barcode map to test added_demultiplex from 'group' in label bc_map_added_demultiplex_group = {'AATCGTGACTCG,coonandfriends':'s1', 'AGAGTCCTGAGC,coonandfriends':'s2', 'AATCGTGACTCG,coonandfriends':'s3'} expected_fasta_fixed_len_bc1 = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0 CCCTTATATATATAT >s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 CCCTTTCCA >s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0 AACCGGCCGGTT >s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1 CCCTTACTATATAT""" # Added demultiplex options expected_fasta_fixed_added_demultiplex = """>s1_0 mintberrycrunch orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG,mintberry bc_diffs=0\nCCCTTATATATATAT\n>s2_1 mysterion orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC,myst bc_diffs=0\nCCCTTTCCA\n>s3_2 mansquito orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG,mansq bc_diffs=0\nAACCGGCCGGTT\n>s1_3 mintberrycrunch orig_bc=AATCGTGACTTG new_bc=AATCGTGACTCG,mintberry bc_diffs=1\nCCCTTACTATATAT""" expected_fasta_added_demultiplex_group = """>s3_0 mintberrycrunch orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG,coonandfriends bc_diffs=0\nCCCTTATATATATAT\n>s2_1 mysterion orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC,coonandfriends bc_diffs=0\nCCCTTTCCA\n>s3_2 mintberrycrunch orig_bc=AATCGTGACTTG new_bc=AATCGTGACTCG,coonandfriends bc_diffs=1\nCCCTTACTATATAT""" # Poor quality window results in second sequence being removed expected_fasta_fixed_len_bc1_sliding_window = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0 CCCTTATATATATAT >s3_1 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0 AACCGGCCGGTT >s1_2 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1 CCCTTACTATATAT""" expected_qual_fixed_len_bc1 = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0 33 33 36 37 37 37 37 37 37 40 40 40 39 39 38 >s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 36 37 37 37 37 37 37 37 37 >s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0 35 32 32 32 32 32 32 32 32 36 37 37 >s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1 37 37 37 37 37 37 37 37 37 37 37 37 37 37 """ # Will be longer because primers are no longer sliced off or checked for # mismatches expected_fasta_fixed_len_bc1_no_primers = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0 GGTCCGGACCCTTATATATATAT >s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 GGTCCGGACCCTTTCCA >s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0 GGTCTGGAAACCGGCCGGTT >s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1 GGTCCGGACCCTTACTATATAT >s2_4 d_primer_error orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 GGTCCGGTACGTTTACTGGA""" # Equal length barcodes, primers, should give different results # Due to parameter changes regarding barcode changes, primer mismatches expected_fasta_fixed_len_bc2 = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0 CCCTTATATATATAT >s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 CCCTTTCCA >s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0 AACCGGCCGGTT >s2_3 d_primer_error orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0 ACGTTTACTGGA""" # These test data have equal length barcodes reverse_primers_fixed_len_bc1 = {'ACACATGTCTAC':'CTTATAT', 'AGAGTCCTGAGC':'GCCCTTT', 'AATCGTGACTCG':'AGTACC'} # Fixed barcode, reverse primers test data in_seqs_reverse_primers = """>a ACACATGTCTACGGTCCGGAGTACCATGATCGGCCCTTATATATATAT >b AGAGTCCTGAGCGGTCCGGAACCGTCGGATCAGCCCTTTCCA >c AATCGTGACTCGGGTCTGGAAACCGATCGACCATAGTACCGGCCGGTT >d ACTCATGTCTACGGTCCGGAATACGTTACGTCCCTTACTATATAT >e_no_barcode_match TTTTGTCCGGACCCTTACTATATAT >d_primer_error AGAGTCCTGAGCGGGGAGGTACGTTTACTGGA """.split() # expected to find the reverse rprimers for seqs a, b, and c, but does not # find the reverse primer for d so writes out whole sequence following forward # primer. expected_in_seqs_reverse_primers = '>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nGTACCATGATCGGCC\n>s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\nACCGTCGGATCA\n>s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0\nAACCGATCGACCAT\n>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\nATACGTTACGTCCCTTACTATATAT' # Will truncate sequence d properly if mismatch in primer allowed expected_in_seqs_reverse_primers_mismatch_allowed = '>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nGTACCATGATCGGCC\n>s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\nACCGTCGGATCA\n>s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0\nAACCGATCGACCAT\n>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\nATACGTTACGTCC' # Fixed barcode, reverse primers test data in_seqs_reverse_primers_mismatch = """>a ACACATGTCTACGGTCCGGAGTACCATGATCGGCCCCTATATATATAT >b AGAGTCCTGAGCGGTCCGGAACCGTCGGATCAGCCCTTTCCA >c AATCGTGACTCGGGTCTGGAAACCGATCGACCATAGTACCGGCCGGTT >d ACTCATGTCTACGGTCCGGAATACGTTACGTCCCTTACATATCCAT >e_no_barcode_match TTTTGTCCGGACCCTTACTATATAT >d_primer_error AGAGTCCTGAGCGGTCCTTTACGCCCACTGGA """.split() # expected to find the reverse rprimers for seqs b, and c, will not find # seq a's reverse primer due to mismatch and will write the whole sequence expected_in_seqs_reverse_primers_mismatch = '>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nGTACCATGATCGGCCCCTATATATATAT\n>s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\nACCGTCGGATCA\n>s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0\nAACCGATCGACCAT\n>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\nATACGTTACGTCCCTTACATATCCAT' # Will not write the d sequence, as the primer mismatches. expected_in_seqs_reverse_primers_full_remove = '>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nGTACCATGATCGGCC\n>s2_1 b orig_bc=AGAGTCCTGAGC new_bc=AGAGTCCTGAGC bc_diffs=0\nACCGTCGGATCA\n>s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0\nAACCGATCGACCAT' # Sample seqs.fna output file sample_seqs_fna_file = """>a1 testing ACACATGTCTACGGTCCGGAACGACGACGAGCGAGGGTAGC >b2 testing more AGAGTCCTGAGCGGTCCGGAACAGACAGGGAGAGACAGAA >a3 testing ACAACAGACGAGTTAGACCAA >d4 AATCGTGACTCGGGTCTGGACAGACGAGAACGAGTTACAGACCAGA""" expected_fasta_mad = """>s1_0 a orig_bc=ACACATGTCTAC new_bc=ACACATGTCTAC bc_diffs=0\nCCCTTATATATATAT\n>s3_2 c orig_bc=AATCGTGACTCG new_bc=AATCGTGACTCG bc_diffs=0\nAACCGGCCGGTT\n>s1_3 d orig_bc=ACTCATGTCTAC new_bc=ACACATGTCTAC bc_diffs=1\nCCCTTACTATATAT""" class SeqQualBadTests(TestCase): """Tests of the SeqQualBad class""" def test_init(self): """SeqQualBad should init OK""" sq = SeqQualBad('Q', None) self.assertEqual(sq.Name, 'Q') self.assertEqual(sq.F, None) self.assertEqual(sq.FailedIds, []) def test_call(self): """SeqQualBad should apply correct function""" f = lambda id_, seq, qual: len(seq) > 3 s1 = 'aa' s2 = 'aaaa' sq = SeqQualBad('Q', f) self.assertEqual(sq('x',s1, [1,2,3]), False) self.assertEqual(sq('y',s2, [1,2,3]), True) self.assertEqual(sq.FailedIds, ['y']) def test_str(self): """SeqQualBad should apply correct function""" f = lambda id_, seq, qual: len(seq) > 3 s1 = 'aa' s2 = 'aaaa' sq = SeqQualBad('Q', f) self.assertEqual(sq('x',s1, [1,2,3]), False) self.assertEqual(str(sq), 'Q\t0') self.assertEqual(sq('y',s2, [1,2,3]), True) self.assertEqual(str(sq), 'Q\t1') if __name__ =='__main__': main()