diff --git a/scripts/test/README.md b/scripts/test/README.md index ab5643a92..2095ed418 100644 --- a/scripts/test/README.md +++ b/scripts/test/README.md @@ -11,23 +11,29 @@ From root: sh scripts/test/simple_tests.sh -#### Regression tests +#### pytest suite -They use the Python3 infrastructure in the `scripts/test` subdirectory, including `bt2face.py`, `dataface.py` and `btdata.py`. +The `test_*.py` modules in `scripts/test` exercise the `bowtie2`, `bowtie2-build` and `bowtie2-inspect` command lines, using small references generated from a fixed seed and the lambda phage example. Shared helpers are in `bt2_helpers.py` and fixtures in `conftest.py`. Requires Python 3 and `pytest`. Set `BOWTIE2_BIN_DIR` to test binaries built somewhere else. From root: - python3 scripts/test/regressions.py --verbose + python3 -m pytest scripts/test + +#### Regression tests + +Tests for specific GitHub issues are in `scripts/test/test_regressions.py` and run as part of the pytest suite. Val Antonescu originally set these up. + +From root: -Val Antonescu originally set these up. + python3 -m pytest scripts/test/test_regressions.py #### Big index test -Builds an index consisting of both human and mouse genomes, pushing the genome size above the 2^32 limit, and necessitating a "big" 64-bit index. This takes a lot of time and RAM. +`scripts/test/test_large_genome.py` aligns reads cut from beyond, and across, the 2^32 offset of a prebuilt large (`.bt2l`) index whose reference exceeds 2^32 bases, such as a combined human and mouse index. Building such an index takes hours and tens of GB of memory, so the tests are skipped unless `BOWTIE2_LARGE_INDEX` names its prefix. Running them needs enough memory to load that index. From root: - python3 scripts/test/large_idx.py --verbose + BOWTIE2_LARGE_INDEX=/path/to/hs_mm python3 -m pytest scripts/test/test_large_genome.py #### Randomized tests diff --git a/scripts/test/big_data/reads/human_reads.fa b/scripts/test/big_data/reads/human_reads.fa deleted file mode 100644 index 1c901bd9c..000000000 --- a/scripts/test/big_data/reads/human_reads.fa +++ /dev/null @@ -1,8 +0,0 @@ ->Fragment_1 -ACAAATGCACTGCTAGGCACCACCCCCAGTTCTAGAATCACACCAGCCAG -TTCACCCTCCAGATGGTTCACCCTCAACTTCATAAAAGTTCCCTACCTAA -TCTACTGACAGGCTCATCCCCGACCTAATTTTAAAGATTTCCTAGGAGCT -GCAGTGGGAATCCTGGACCTCAGCCTGGACAAAGAACAGCTGCAGGTCAT -TCTCATGTGTGGACACAGAAGCTCTGCCTGCCTTTGCTGGCCAGCTGGGC -TGAGCGGGCCTGGGAATTAAGGCTGCAGGGTTGGTCCCAGGCAGTCTTGC - diff --git a/scripts/test/big_data/reads/mouse_reads.fa b/scripts/test/big_data/reads/mouse_reads.fa deleted file mode 100644 index 1f104d597..000000000 --- a/scripts/test/big_data/reads/mouse_reads.fa +++ /dev/null @@ -1,8 +0,0 @@ ->Fragment_2 -GAATAGACAGGATGAGTATTGTCCTGTAATTTCATAAATTTGTTATTACA -CTCTTGTAGTTCATGTTTAAAACACTTTCTGTCACCTTTATCTTCTGAAA -ATTTACTATGTTATGAATTTCTTGAAGTTTTGATAGTTCGGGAAAAAGTT -TTAGTGTCCTCTTCTAAAGGCTTGCCTCTCCTCATTTCTTGAGAAATCTG -TGTTTTCCAGTTTTTGTAATCATTAATTTTGATtgtgtgtgtgtgtgtgt -gtgtgGTTTGTTTGTTGTTTTTTTTGCTACTTAAACGATATTGTTAACAC - diff --git a/scripts/test/bt2face.py b/scripts/test/bt2face.py deleted file mode 100755 index 8799491a6..000000000 --- a/scripts/test/bt2face.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python3 - -import os -import logging -import subprocess - - -class BowtieSuite(object): - - def __init__(self,bt2_path=''): - curr_path = os.path.realpath(bt2_path) - self.bowtie_bin = os.path.join(curr_path,'bowtie2') - self.bowtie_build = os.path.join(curr_path,'bowtie2-build') - self.bowtie_inspect = os.path.join(curr_path,'bowtie2-inspect') - logging.info('bowtie2 path: ' + self.bowtie_bin) - logging.info('bowtie2-build path: ' + self.bowtie_build) - logging.info('bowtie2-inspect path: ' + self.bowtie_inspect) - - - def run(self, *args): - cmd = self.bowtie_bin + " " + " ".join([i for i in args]) - logging.debug('align cmd: ' + cmd) - return(subprocess.call(cmd,shell=True)) - - - def silent_run(self, *args): - cmd = self.bowtie_bin + " " + " ".join([i for i in args]) - logging.debug('align cmd (silent): ' + cmd) - return(subprocess.call(cmd,shell=True,stderr=open(os.devnull, 'w'))) - - - def build(self, *args): - cmd = self.bowtie_build + " " + " ".join([i for i in args]) - curr_dir = os.getcwd() - os.chdir(os.path.dirname(cmd.split()[-2])) - logging.debug('build cmd: ' + cmd) - ret = subprocess.check_call(cmd,shell=True) - os.chdir(curr_dir) - return(ret) diff --git a/scripts/test/btdata.py b/scripts/test/btdata.py deleted file mode 100755 index 7ec9f3cba..000000000 --- a/scripts/test/btdata.py +++ /dev/null @@ -1,192 +0,0 @@ -#!/usr/bin/env python3 -""" -Note: This would look so much better replaced by XML or at least JSON. But - is not worth to do it for now. -""" - -import os -import gzip -import urllib.request, urllib.error, urllib.parse -import logging - - -class LargeTestsData(object): - """ - Large index tests use quite big datasets. Make sure these - are present before starting the time consuming tests. - """ - - def __init__(self,bt2_path=''): - self.data_dir = 'big_data' - curr_path = os.path.realpath(bt2_path) - - curr_path = os.path.join(curr_path,'scripts') - curr_path = os.path.join(curr_path,'test') - self.data_dir_path = os.path.join(curr_path,self.data_dir) - self.reads_dir_path = os.path.join(curr_path,'reads') - - try: - os.stat(self.data_dir_path) - except: - logging.error("Cannot find the working datadir %s!" % self.data_dir_path) - raise - - self.genomes = dict() - self.genomes['human'] = dict() - hm = self.genomes['human'] - hm['link'] = "ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/" - hm['ref_name'] = 'human.fa' - hm['chromosomes'] = [] - chromosomes = hm['chromosomes'] - for i in range(1,22): - chromosomes.append('chr%d' % i) - chromosomes.extend(['chrX', 'chrY', 'chrM']) - - self.genomes['mouse'] = dict() - ms = self.genomes['mouse'] - ms['link'] = "ftp://hgdownload.cse.ucsc.edu/goldenPath/rn4/chromosomes" - ms['ref_name'] = 'mouse.fa' - ms['chromosomes'] = [] - chromosomes = ms['chromosomes'] - for i in range(1,21): - chromosomes.append('chr%d' % i) - chromosomes.extend(['chrX', 'chrM']) - - self.joint_genomes = dict() - self.joint_genomes['ms_hum'] = dict() - mh = self.joint_genomes['ms_hum'] - mh['link'] = None - mh['ref_name'] = 'ms_hum.fa' - mh['genomes'] = ['human','mouse'] - - self.init_data() - - - - def init_data(self): - """ Try and init the data we need. - """ - for genome,gdata in list(self.genomes.items()): - gn_path = os.path.join(self.data_dir_path,genome) - gn_fasta = os.path.join(gn_path,gdata['ref_name']) - if not os.path.exists(gn_fasta): - self._get_genome(genome) - self._build_genome(genome) - - for genome,gdata in list(self.joint_genomes.items()): - gn_path = os.path.join(self.data_dir_path,genome) - gn_fasta = os.path.join(gn_path,gdata['ref_name']) - if not os.path.exists(gn_fasta): - self._build_joint_genome(genome) - - - - def _get_genome(self,genome): - g = self.genomes[genome] - gn_path = os.path.join(self.data_dir_path,genome) - - if not os.path.exists(gn_path): - os.mkdir(gn_path) - - logging.info("Downloading genome: %s " % genome) - - for chrs in g['chromosomes']: - chr_file = chrs + ".fa.gz" - fname = os.path.join(gn_path,chr_file) - - if os.path.exists(fname): - logging.info("Skip %s (already present)" % chr_file) - continue - - uri = g['link'] + r"/" + chr_file - logging.info("file: %s" % chr_file) - - try: - f = open(fname,'wb') - u = urllib.request.urlopen(uri) - f.write(u.read()) - except: - f.close() - os.remove(fname) - os.close(u.fileno()) - raise - else: - os.close(u.fileno()) - u.close() - f.close() - - - - def _build_genome(self,genome): - g = self.genomes[genome] - gn_path = os.path.join(self.data_dir_path,genome) - gn_fasta = os.path.join(gn_path,g['ref_name']) - - logging.info("Building fasta file for genome: %s" % genome) - - f_gn = open(gn_fasta,'wb') - - for chrs in g['chromosomes']: - chr_file = chrs + ".fa.gz" - fname = os.path.join(gn_path,chr_file) - - try: - f_chr = gzip.open(fname,'rb') - f_gn.write(f_chr.read()) - except: - f_chr_close() - f_gn.close() - os.remove(gn_fasta) - raise - else: - f_chr.close() - - f_gn.close() - - - - def _build_joint_genome(self,genome): - jg = self.joint_genomes[genome] - jgn_path = os.path.join(self.data_dir_path,genome) - jgn_fasta = os.path.join(jgn_path,jg['ref_name']) - - if not os.path.exists(jgn_path): - os.mkdir(jgn_path) - - logging.info("Building fasta file for genome: %s" % genome) - - f_jg = open(jgn_fasta,'wb') - for g in jg['genomes']: - gn_path = os.path.join(self.data_dir_path,g) - fasta_file = os.path.join(gn_path,self.genomes[g]['ref_name']) - try: - fin = open(fasta_file,'rb') - f_jg.write(fin.read()) - except: - fin.close() - f_jg.close() - os.remove(jgn_fasta) - raise - else: - fin.close() - - f_jg.close() - - - -class ExampleData(object): - """ The example data. - """ - - def __init__(self,bt2_path=''): - curr_path = os.path.realpath(bt2_path) - curr_path = os.path.join(curr_path,'example') - self.index_dir_path = os.path.join(curr_path,'index') - self.reads_dir_path = os.path.join(curr_path,'reads') - self.ref_dir_path = os.path.join(curr_path,'reference') - - try: - os.stat(curr_path) - except: - logging.error("Cannot find the example datadir %s!" % curr_path) - raise diff --git a/scripts/test/dataface.py b/scripts/test/dataface.py deleted file mode 100755 index 0a52e938c..000000000 --- a/scripts/test/dataface.py +++ /dev/null @@ -1,81 +0,0 @@ -#!/usr/bin/env python3 - -import os -import logging - - -class DataFace(object): - """ Some data IFace bowtie can work with. - """ - - def size(self): - raise NotImplementedError("size() needs to be implemented!") - - - -class SamFile(DataFace): - - def __init__(self,sam_desc): - self.file_desc = sam_desc - - - def size(self): - if hasattr(self,'no_frags'): - return self.no_frags - - count = 0 - try: - fh = open(self.file_desc,"r") - for line in fh: - if line[0] != "@": - count += 1 - except: - logging.error("Exception reading sam file!") - fh.close() - raise - else: - fh.close() - - self.no_frags = count - return count - - - -class FastaQFile(DataFace): - - def __init__(self,fastq_desc): - self.file_desc = fastq_desc - - - def size(self): - if hasattr(self,'no_seqs'): - return self.no_seqs - - count = 0 - try: - fh = open(self.file_desc,"r") - lno = 0 - for line in fh: - if lno == 0: - if line[0] != '@': - raise TypeError("This does not look like a fastq file!") - count += 1 - if lno == 2: - if line[0] != '+': - raise TypeError("This does not look like a fastq file!") - lno = (lno + 1) % 4 - except: - logging.error("Exception reading fastq file!") - fh.close() - raise - else: - fh.close() - - self.no_seqs = count - return self.no_seqs - - - - - - \ No newline at end of file diff --git a/scripts/test/large_idx.py b/scripts/test/large_idx.py deleted file mode 100755 index 055cc0628..000000000 --- a/scripts/test/large_idx.py +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env python3 - -import os -import gzip -import urllib.request, urllib.error, urllib.parse -import inspect -import unittest -import logging -import bt2face -import btdata -import subprocess -from optparse import OptionParser - - - -class TestLargeIndex(unittest.TestCase): - """ - Main fixture for large index testing. - """ - - def test_human(self): - wdir = g_bdata.data_dir_path - wdir = os.path.join(wdir,'human') - rdir = g_bdata.reads_dir_path - genome = g_bdata.genomes['human'] - genome_fasta = os.path.join(wdir,genome['ref_name']) - genome_index = os.path.join(wdir,'human') - reads = os.path.join(rdir,'human_reads.fa') - ret = g_bt.build("%s human" % genome_fasta) - self.assertEqual(ret,0) - args = "-x %s -f -U %s" % (genome_index,reads) - ret = g_bt.run(args) - self.assertEqual(ret,0) - - - def test_mouse(self): - wdir = g_bdata.data_dir_path - wdir = os.path.join(wdir,'mouse') - rdir = g_bdata.reads_dir_path - genome = g_bdata.genomes['mouse'] - genome_fasta = os.path.join(wdir,genome['ref_name']) - genome_index = os.path.join(wdir,'mouse') - reads = os.path.join(rdir,'mouse_reads.fa') - ret = g_bt.build("%s mouse" % genome_fasta) - self.assertEqual(ret,0) - args = "-x %s -f -U %s" % (genome_index,reads) - ret = g_bt.run(args) - self.assertEqual(ret,0) - - - def test_large_index(self): - wdir = g_bdata.data_dir_path - wdir = os.path.join(wdir,'ms_hum') - rdir = g_bdata.reads_dir_path - genome = g_bdata.joint_genomes['ms_hum'] - genome_fasta = os.path.join(wdir,genome['ref_name']) - genome_index = os.path.join(wdir,'ms_hum') - reads_human = os.path.join(rdir,'human_reads.fa') - reads_mouse = os.path.join(rdir,'mouse_reads.fa') - ret = g_bt.build("%s ms_hum" % genome_fasta) - self.assertEqual(ret,0) - args = "-x %s -f -U %s" % (genome_index,reads_human) - ret = g_bt.run(args) - args = "-x %s -f -U %s" % (genome_index,reads_mouse) - ret = g_bt.run(args) - self.assertEqual(ret,0) - - - -def get_suite(): - tests = ['test_human','test_mouse','test_large_index'] - return unittest.TestSuite(list(map(TestLargeIndex,tests))) - - - -def parse_args(): - usage = " %prog [options] \n\n" - usage += "Warning, this test runs some VERY resource consuming tests.\n" - parser = OptionParser(usage=usage) - parser.add_option("-v", "--verbose", - action="store_true",dest="verbose", default=False, - help="Print more info about each test.") - - (options, args) = parser.parse_args() - return options - - -g_bdata = None -g_bt = None - -if __name__ == "__main__": - logging.basicConfig(format='%(levelname)s:%(message)s',level=logging.ERROR) - options = parse_args() - - runner = unittest.TextTestRunner() - if options.verbose: - logging.getLogger().setLevel(level=logging.INFO) - runner = unittest.TextTestRunner(verbosity=2) - - src_file_path = os.path.realpath(inspect.getsourcefile(parse_args)) - curr_path = os.path.dirname(src_file_path) - bw2_subdir = 'bowtie2' - - i = curr_path.find(bw2_subdir) - bt2_path = curr_path[:i+len(bw2_subdir)] - - g_bdata = btdata.LargeTestsData(bt2_path) - g_bt = bt2face.BowtieSuite(bt2_path) - runner.run(get_suite()) - diff --git a/scripts/test/regressions.py b/scripts/test/regressions.py deleted file mode 100755 index 78da1ab61..000000000 --- a/scripts/test/regressions.py +++ /dev/null @@ -1,207 +0,0 @@ -#!/usr/bin/env python3 - -import os -import inspect -import unittest -import logging -import shutil -import bt2face -import dataface -import btdata -from optparse import OptionParser - - -class TestRegressions(unittest.TestCase): - """ - Main regression fixture. - """ - - def test_288(self): - """ Check if --un option works correctly when used with --no-unal - """ - # run --un and record unaligned file size and orig file size - # run --no-unal and record file size - # check if --no-unali + --un files size equal the orig file size - out_sam = 'test288_out.sam' - not_algn = 'test288_nal.fastq' - - ref_index = os.path.join(g_bdata.index_dir_path,'lambda_virus') - reads = os.path.join(g_bdata.reads_dir_path,'longreads.fq') - - args = "--quiet -x %s -U %s -a --un %s -S %s" % (ref_index,reads,not_algn,out_sam) - ret = g_bt.run(args) - self.assertEqual(ret,0) - sam_size = dataface.SamFile(out_sam).size() - no_algn_size = dataface.FastaQFile(not_algn).size() - args = "--quiet -x %s -U %s -a --un %s --no-unal -S %s" % (ref_index,reads,not_algn,out_sam) - ret = g_bt.run(args) - self.assertEqual(ret,0) - no_al_sam_size = dataface.SamFile(out_sam).size() - self.assertEqual(no_al_sam_size + no_algn_size, sam_size) - os.remove(out_sam) - os.remove(not_algn) - - - def test_279(self): - """ Check if --un-conc option works correctly when used with --no-unal - """ - out_no_conc = 'test279_nconc.fq' - out_conc = 'test279_conc.fq' - out_sam = 'test279.sam' - out_p1_nc = 'test279_nconc.1.fq' - out_p2_nc = 'test279_nconc.2.fq' - out_p1_conc = 'test279_conc.1.fq' - out_p2_conc = 'test279_conc.2.fq' - ref_index = os.path.join(g_bdata.index_dir_path,'lambda_virus') - pairs_1 = os.path.join(g_bdata.reads_dir_path,'reads_1.fq') - pairs_2 = os.path.join(g_bdata.reads_dir_path,'reads_2.fq') - - args = "--quiet -x %s -1 %s -2 %s --un-conc %s --dovetail --al-conc %s -S %s" % (ref_index,pairs_1,pairs_2,out_no_conc,out_conc,out_sam) - ret = g_bt.run(args) - self.assertEqual(ret,0) - conc_p1_size_step1 = dataface.FastaQFile(out_p1_conc).size() - conc_p2_size_step1 = dataface.FastaQFile(out_p2_conc).size() - nc_p1_size_step1 = dataface.FastaQFile(out_p1_nc).size() - nc_p2_size_step1 = dataface.FastaQFile(out_p2_nc).size() - self.assertEqual(conc_p1_size_step1, conc_p2_size_step1, "Number of concordant sequences in pass 1 should match.") - self.assertEqual(nc_p1_size_step1, nc_p2_size_step1, "Number of non-concordant sequences in pass 1 should match.") - - args = "--quiet -x %s -1 %s -2 %s --no-unal --un-conc %s --dovetail --al-conc %s -S %s" % (ref_index,pairs_1,pairs_2,out_no_conc,out_conc,out_sam) - ret = g_bt.run(args) - self.assertEqual(ret,0) - conc_p1_size_step2 = dataface.FastaQFile(out_p1_conc).size() - conc_p2_size_step2 = dataface.FastaQFile(out_p2_conc).size() - nc_p1_size_step2 = dataface.FastaQFile(out_p1_nc).size() - nc_p2_size_step2 = dataface.FastaQFile(out_p2_nc).size() - self.assertEqual(conc_p1_size_step2, conc_p2_size_step2, "Number of concordant sequences in pass 2 should match.") - self.assertEqual(nc_p1_size_step2, nc_p2_size_step2, "Number of non-concordant sequences in pass 2 should match.") - self.assertEqual(conc_p1_size_step1, conc_p1_size_step2, "Number of concordant sequences in both steps should match.") - self.assertEqual(conc_p2_size_step1, conc_p2_size_step2, "Number of concordant sequences in both steps should match.") - self.assertEqual(nc_p1_size_step1, nc_p1_size_step2, "Number of non-concordant sequences in both steps should match.") - self.assertEqual(nc_p2_size_step1, nc_p2_size_step2, "Number of non-concordant sequences in both steps should match.") - - self.assertTrue(1,"temp") - os.remove(out_p1_conc) - os.remove(out_p2_conc) - os.remove(out_p1_nc) - os.remove(out_p2_nc) - os.remove(out_sam) - - - def test_version(self): - """ Check if --version option works correctly. - """ - args = "--version" - ret = g_bt.silent_run(args) - self.assertEqual(ret,0) - - - def test_x_option(self): - """ Check if -x is indeed mandatory. - """ - ref_index = os.path.join(g_bdata.index_dir_path,'lambda_virus') - reads = os.path.join(g_bdata.reads_dir_path,'longreads.fq') - out_sam = 'test_x_option.sam' - - args = "%s -U %s -S %s " % (ref_index,reads,out_sam) - ret = g_bt.silent_run(args) - self.assertNotEqual(ret,0) - args = "-x %s -U %s -S %s " % (ref_index,reads,out_sam) - ret = g_bt.silent_run(args) - self.assertEqual(ret, 0) - os.remove(out_sam) - - - def test_313(self): - """ Check if --un-conc, --al-conc, --un, --al do not fail in case of a dot - character present in their path name. - """ - ref_index = os.path.join(g_bdata.index_dir_path,'lambda_virus') - pairs_1 = os.path.join(g_bdata.reads_dir_path,'reads_1.fq') - pairs_2 = os.path.join(g_bdata.reads_dir_path,'reads_2.fq') - no_dot_dir = 'nodotdir' - dot_dir = 'v.12_.test' - un_basename = 'bname' - out_spec = os.path.join(no_dot_dir,un_basename) - out_spec_dot= os.path.join(dot_dir,un_basename) - args = "--quiet -x %s -1 %s -2 %s --un-conc %s -S %s" % (ref_index,pairs_1,pairs_2,out_spec,os.devnull) - - try: - os.makedirs(no_dot_dir) - os.makedirs(dot_dir) - except: - pass - - ret = g_bt.silent_run(args) - self.assertEqual(ret, 0) - mate1_count = dataface.FastaQFile(out_spec + '.1').size() - args = "--quiet -x %s -1 %s -2 %s --un-conc %s -S %s" % (ref_index,pairs_1,pairs_2,out_spec_dot,os.devnull) - ret = g_bt.silent_run(args) - self.assertEqual(ret, 0, "--un-conc should work with dots in its path.") - mate1_count_dot = dataface.FastaQFile(out_spec_dot + '.1').size() - self.assertEqual(mate1_count, mate1_count_dot, "Number of seqs should be the same.") - - args = "--quiet -x %s -1 %s -2 %s --un-conc %s -S %s" % (ref_index,pairs_1,pairs_2,dot_dir,os.devnull) - ret = g_bt.silent_run(args) - self.assertEqual(ret, 0, "--un-conc should work with a path only.") - mate1_fs = os.path.join(dot_dir,"un-conc-mate.1") - self.assertTrue(os.path.isfile(mate1_fs), "--un-conc output is missing for path only case.") - mate1_count_dot_only = dataface.FastaQFile(mate1_fs).size() - self.assertEqual(mate1_count, mate1_count_dot_only, "Number of seqs should be the same.") - - args= "--quiet -x %s -1 %s -2 %s --un-conc %s --al-conc %s --un %s -S %s" % (ref_index,pairs_1,pairs_2,dot_dir,dot_dir,dot_dir,os.devnull) - ret = g_bt.run(args) - self.assertEqual(ret, 0) - mate1_al = os.path.join(dot_dir,"al-conc-mate.1") - mate1_un = os.path.join(dot_dir,"un-seqs") - self.assertTrue( - os.path.isfile(mate1_fs) and - os.path.isfile(mate1_al) and - os.path.isfile(mate1_un), - "--un-conc, --al-conc and --un should not clobber" - ) - shutil.rmtree(no_dot_dir) - shutil.rmtree(dot_dir) - - - -def get_suite(): - loader = unittest.TestLoader() - return loader.loadTestsFromTestCase(TestRegressions) - - -def parse_args(): - usage = " %prog [options] \n\n" - usage += "Some regression tests.\n" - parser = OptionParser(usage=usage) - parser.add_option("-v", "--verbose", - action="store_true",dest="verbose", default=False, - help="Print more info about each test.") - - (options, args) = parser.parse_args() - return options - - -g_bdata = None -g_bt = None - -if __name__ == "__main__": - logging.basicConfig(format='%(levelname)s:%(message)s',level=logging.ERROR) - options = parse_args() - - runner = unittest.TextTestRunner() - if options.verbose: - logging.getLogger().setLevel(level=logging.INFO) - runner = unittest.TextTestRunner(verbosity=2) - - src_file_path = os.path.realpath(inspect.getsourcefile(parse_args)) - curr_path = os.path.dirname(src_file_path) - bw2_subdir = 'bowtie2' - - i = curr_path.find(bw2_subdir) - bt2_path = curr_path[:i+len(bw2_subdir)] - - g_bdata = btdata.ExampleData(bt2_path) - g_bt = bt2face.BowtieSuite(bt2_path) - runner.run(get_suite()) - diff --git a/scripts/test/test_large_genome.py b/scripts/test/test_large_genome.py new file mode 100644 index 000000000..76db8e6ed --- /dev/null +++ b/scripts/test/test_large_genome.py @@ -0,0 +1,224 @@ +"""Alignment against a prebuilt large (.bt2l) index of more than 2^32 bases. + +The rest of the suite checks the large-index code path on small genomes +forced into the .bt2l format with --large-index. Only an index whose joined +reference really exceeds 2^32 bases exercises offsets that do not fit in 32 +bits, and building one takes hours and tens of GB of memory, so these tests +are skipped unless BOWTIE2_LARGE_INDEX names the prefix of such an index +(for example a human plus mouse index on a cluster). + +Reads are cut from the index's own reference, obtained by streaming +`bowtie2-inspect` one sequence at a time, at positions past the 2^32 offset +in the joined reference (which leaves out runs of N) and straddling it. Each +read must align perfectly, and every read whose alignment is unique must +land where it was cut. + + BOWTIE2_LARGE_INDEX=/path/to/hs_mm python3 -m pytest scripts/test/test_large_genome.py + +Optional settings: + + BOWTIE2_LARGE_INDEX_THREADS -p for the alignment run (default 4) + BOWTIE2_LARGE_INDEX_OFFSET joined offset the reads must reach (default + 2^32); lowering it lets the test logic be + checked against a small .bt2l index +""" + +import os +import subprocess + +import pytest + +import bt2_helpers as h + + +INDEX = os.environ.get('BOWTIE2_LARGE_INDEX') +THRESHOLD = int(os.environ.get('BOWTIE2_LARGE_INDEX_OFFSET', 2 ** 32)) +THREADS = os.environ.get('BOWTIE2_LARGE_INDEX_THREADS', '4') + +READ_LEN = 100 +FRAG_LEN = 300 +FAR_READS = 150 +CONTROL_READS = 5 + +pytestmark = pytest.mark.skipif(not INDEX, + reason='set BOWTIE2_LARGE_INDEX to a .bt2l index of more than 2^32 bases') + + +def inspect_sequences(bt2, prefix): + """Yield (name, sequence as bytes) for each reference sequence, streaming + bowtie2-inspect so that only one sequence is held in memory.""" + proc = subprocess.Popen([bt2.inspect, prefix], stdout=subprocess.PIPE) + name, seq = None, bytearray() + try: + for line in proc.stdout: + line = line.rstrip(b'\r\n') + if line.startswith(b'>'): + if name is not None: + yield name, bytes(seq) + name, seq = line[1:].decode().split()[0], bytearray() + else: + seq.extend(line) + if name is not None: + yield name, bytes(seq) + finally: + proc.stdout.close() + assert proc.wait() == 0, 'bowtie2-inspect failed on %s' % prefix + + +def joined_offset(seq, start, p): + """Offset in the joined reference of position p of seq, where start is the + joined offset of seq's first non-N base.""" + return start + p - seq.count(b'N', 0, p) + + +def position_at(seq, start, target): + """Smallest position p of seq whose joined offset is at least target.""" + lo, hi = 0, len(seq) + while lo < hi: + mid = (lo + hi) // 2 + if joined_offset(seq, start, mid) < target: + lo = mid + 1 + else: + hi = mid + return lo + + +def is_clean(seq, q, length): + return 0 <= q <= len(seq) - length and seq[q:q + length].translate(None, b'ACGT') == b'' + + +def clean_window(seq, p, length, limit=10000): + """First position at or after p where seq[p:p+length] is all A/C/G/T, or None.""" + for q in range(p, min(p + limit, len(seq) - length + 1)): + if is_clean(seq, q, length): + return q + return None + + +def straddling_window(seq, p, length): + """A position q with q < p < q+length and seq[q:q+length] all A/C/G/T, the + one with p nearest the middle, or None.""" + for d in range(length // 2): + for q in (p - length // 2 - d, p - length // 2 + d): + if q < p < q + length and is_clean(seq, q, length): + return q + return None + + +@pytest.fixture(scope='module') +def large_index(bt2): + assert os.path.exists(INDEX + '.1.bt2l'), \ + 'BOWTIE2_LARGE_INDEX=%s: no %s.1.bt2l' % (INDEX, INDEX) + return INDEX + + +@pytest.fixture(scope='module') +def cut_reads(bt2, large_index): + """Reads and pairs cut from the reference, with where they came from. + + Returns (unpaired, pairs, total joined length). unpaired holds + (name, seq, rname, pos, joined offset) and pairs holds + (name, mate1, mate2, rname, pos1, pos2, joined offset); positions are + 1-based as in SAM.""" + # Reads beyond the threshold are cut at evenly spaced joined offsets, so + # most come from long sequences rather than from the many short, often + # repetitive unplaced contigs an assembly may end with. The summary's + # lengths include runs of N, so the spacing is an overestimate at worst. + summary = bt2.run('bowtie2-inspect', ['-s', large_index]).stdout + total_len = sum(int(l.split('\t')[-1]) for l in summary.splitlines() + if l.startswith('Sequence-')) + spacing = max(FRAG_LEN, (total_len - THRESHOLD) // FAR_READS) + unpaired, pairs = [], [] + start, next_pick = 0, THRESHOLD + for rname, seq in inspect_sequences(bt2, large_index): + end = joined_offset(seq, start, len(seq)) + picks = [] + if not unpaired: + picks += [len(seq) * (i + 1) // (CONTROL_READS + 1) for i in range(CONTROL_READS)] + if start <= THRESHOLD < end: + # A read straddling the threshold. + q = straddling_window(seq, position_at(seq, start, THRESHOLD), READ_LEN) + if q is not None: + picks.append(q) + while next_pick < end: + picks.append(position_at(seq, start, next_pick)) + next_pick += spacing + for p in picks: + q = clean_window(seq, p, READ_LEN) + if q is None: + continue + off = joined_offset(seq, start, q) + name = 'r%d_%s_%d' % (len(unpaired), rname, q + 1) + unpaired.append((name, seq[q:q + READ_LEN].decode(), rname, q + 1, off)) + f = clean_window(seq, p, FRAG_LEN) + if f is not None and off >= THRESHOLD: + frag = seq[f:f + FRAG_LEN].decode() + pairs.append(('p%d_%s_%d' % (len(pairs), rname, f + 1), frag[:READ_LEN], + h.revcomp(frag[-READ_LEN:]), rname, f + 1, f + FRAG_LEN - READ_LEN + 1, + joined_offset(seq, start, f))) + start = end + return unpaired, pairs, start + + +def test_index_is_larger_than_threshold(cut_reads): + """The joined reference is longer than the offset the tests target, so + reads beyond it exist.""" + _, _, total = cut_reads + assert total > THRESHOLD, ('the index holds %d bases of non-N reference, not more ' + 'than %d' % (total, THRESHOLD)) + + +def test_large_offsets_align_where_cut(bt2, large_index, cut_reads, tmp_path): + """Unpaired reads cut beyond the threshold, and across it, align perfectly, + and pairs cut beyond it align concordantly. Uniquely aligned reads and + pairs land where they were cut, pairs with the right TLEN.""" + unpaired, pairs, _ = cut_reads + far = [r for r in unpaired if r[4] >= THRESHOLD] + assert far, 'no read could be cut beyond offset %d' % THRESHOLD + assert any(r[4] < THRESHOLD < r[4] + READ_LEN for r in unpaired), \ + 'no read straddles offset %d' % THRESHOLD + assert pairs, 'no pair could be cut beyond offset %d' % THRESHOLD + + u = h.write_fastq(str(tmp_path / 'u.fq'), [(r[0], r[1]) for r in unpaired]) + m1 = h.write_fastq(str(tmp_path / 'm1.fq'), [(p[0], p[1]) for p in pairs]) + m2 = h.write_fastq(str(tmp_path / 'm2.fq'), [(p[0], p[2]) for p in pairs]) + # Unpaired and paired reads go in separate runs: mixing them in one run + # with -p > 1 loses the pairs or hangs. + recs = {} + for name, reads in (('u', ['-U', u]), ('p', ['-1', m1, '-2', m2])): + sam = str(tmp_path / ('%s.sam' % name)) + bt2.run('bowtie2', ['-x', large_index, '-p', THREADS, '-S', sam] + reads, + timeout=7200) + with open(sam) as f: + recs.update(h.by_name(h.primary(h.parse_sam(f.read())[1]))) + + def unique(rec): + return 'XS' not in rec.tags or rec.tags['XS'] < rec.tags['AS'] + + # Repeats are common in large genomes, so only a quarter of the reads are + # required to align uniquely; every one that does must be in place. + located = 0 + for name, _, rname, pos, off in unpaired: + (rec,) = recs[name] + assert not rec.is_unmapped and rec.tags['AS'] == 0, rec + if unique(rec): + assert (rec.rname, rec.pos) == (rname, pos), rec + located += off >= THRESHOLD + assert located * 4 >= len(far), \ + 'only %d of %d reads beyond the threshold aligned uniquely' % (located, len(far)) + + # For a pair in a repeat, bowtie2 may report a concordant placement that is + # not the best one (one mate with mismatches and XS:i above AS:i), so a + # perfect score is required only of pairs whose mates both align uniquely. + located = 0 + for name, _, _, rname, pos1, pos2, _ in pairs: + r1, r2 = sorted(recs[name], key=lambda r: not r.flag & h.FIRST_IN_PAIR) + for r in (r1, r2): + assert not r.is_unmapped and r.tags['YT'] == 'CP', r + if unique(r1) and unique(r2): + assert r1.tags['AS'] == 0 and r2.tags['AS'] == 0, (r1, r2) + assert (r1.rname, r1.pos, r1.tlen) == (rname, pos1, FRAG_LEN), r1 + assert (r2.rname, r2.pos, r2.tlen) == (rname, pos2, -FRAG_LEN), r2 + located += 1 + assert located * 4 >= len(pairs), \ + 'only %d of %d pairs aligned uniquely' % (located, len(pairs)) diff --git a/scripts/test/test_regressions.py b/scripts/test/test_regressions.py new file mode 100644 index 000000000..a24f015ce --- /dev/null +++ b/scripts/test/test_regressions.py @@ -0,0 +1,118 @@ +"""Regression tests for specific bowtie2 issues. + +Each test names the GitHub issue (or behaviour) it guards and uses the +lambda phage example shipped in example/. New regression tests belong here +when they do not fit naturally in one of the topic modules. + +Run from the repository root after `make`: + + python3 -m pytest scripts/test/test_regressions.py +""" + +import os + +import bt2_helpers as h + + +def sam_records(path): + with open(path) as f: + return h.parse_sam(f.read())[1] + + +def read_names(path): + return [r[0] for r in h.read_fastq(path)] + + +def file_bytes(path): + with open(path, 'rb') as f: + return f.read() + + +def test_issue_288_un_with_no_unal(bt2, lambda_index, example_reads, tmp_path): + """--un writes the same reads with or without --no-unal (issue 288). + + With -a, every read that fails to align is one unaligned SAM record. Those + reads are exactly the ones --un writes, and --no-unal removes exactly + those records from the SAM and nothing else.""" + reads = example_reads['longreads'] + full_sam, full_un = str(tmp_path / 'full.sam'), str(tmp_path / 'full_un.fq') + nounal_sam, nounal_un = str(tmp_path / 'nounal.sam'), str(tmp_path / 'nounal_un.fq') + base = ['--quiet', '-x', lambda_index, '-U', reads, '-a'] + bt2.run('bowtie2', base + ['--un', full_un, '-S', full_sam]) + bt2.run('bowtie2', base + ['--un', nounal_un, '--no-unal', '-S', nounal_sam]) + full = sam_records(full_sam) + unaligned = [r.qname for r in full if r.is_unmapped] + assert unaligned, 'expected some unaligned long reads' + assert [n.split()[0] for n in read_names(full_un)] == unaligned + assert file_bytes(nounal_un) == file_bytes(full_un) + assert [r.line for r in sam_records(nounal_sam)] == \ + [r.line for r in full if not r.is_unmapped] + + +def test_issue_279_un_conc_with_no_unal(bt2, lambda_index, example_reads, tmp_path): + """--un-conc and --al-conc write the same pairs with or without --no-unal (issue 279). + + Both mate files of each output hold the same number of reads, and + --al-conc holds exactly the pairs the SAM reports as concordant.""" + outs = {} + for name, extra in (('plain', []), ('nounal', ['--no-unal'])): + d = tmp_path / name + d.mkdir() + sam = str(d / 'out.sam') + bt2.run('bowtie2', ['--quiet', '-x', lambda_index, '-1', example_reads['reads_1'], + '-2', example_reads['reads_2'], '--dovetail', '--un-conc', str(d / 'nconc.fq'), + '--al-conc', str(d / 'conc.fq'), '-S', sam] + extra) + outs[name] = d + for kind in ('conc', 'nconc'): + m1 = read_names(str(d / ('%s.1.fq' % kind))) + m2 = read_names(str(d / ('%s.2.fq' % kind))) + assert len(m1) == len(m2), (name, kind) + if name == 'plain': + concordant = sorted(set(r.qname for r in sam_records(sam) if r.tags.get('YT') == 'CP')) + conc = sorted(set(n.split()[0].rsplit('/', 1)[0] for n in read_names(str(d / 'conc.1.fq')))) + assert conc == concordant + assert concordant, 'expected some concordant pairs' + for f in ('conc.1.fq', 'conc.2.fq', 'nconc.1.fq', 'nconc.2.fq'): + assert file_bytes(str(outs['nounal'] / f)) == file_bytes(str(outs['plain'] / f)), f + + +def test_issue_313_dots_in_output_paths(bt2, lambda_index, example_reads, tmp_path): + """--un-conc, --al-conc and --un work when the output path has a dot in a + directory name, and when given a bare directory (issue 313). + + A bare directory gets the default names un-conc-mate.1/.2, + al-conc-mate.1/.2 and un-seqs, and the outputs do not overwrite each + other.""" + pairs = ['--quiet', '-x', lambda_index, '-1', example_reads['reads_1'], + '-2', example_reads['reads_2'], '-S', os.devnull] + plain_dir, dot_dir = tmp_path / 'nodotdir', tmp_path / 'v.12_.test' + plain_dir.mkdir() + dot_dir.mkdir() + + bt2.run('bowtie2', pairs + ['--un-conc', str(plain_dir / 'bname')]) + expected = file_bytes(str(plain_dir / 'bname.1')) + assert read_names(str(plain_dir / 'bname.1')), 'expected some non-concordant pairs' + + bt2.run('bowtie2', pairs + ['--un-conc', str(dot_dir / 'bname')]) + assert file_bytes(str(dot_dir / 'bname.1')) == expected + + bt2.run('bowtie2', pairs + ['--un-conc', str(dot_dir)]) + assert file_bytes(str(dot_dir / 'un-conc-mate.1')) == expected + + bt2.run('bowtie2', pairs + ['--un-conc', str(dot_dir), '--al-conc', str(dot_dir), + '--un', str(dot_dir)]) + for name in ('un-conc-mate.1', 'un-conc-mate.2', 'al-conc-mate.1', 'al-conc-mate.2', + 'un-seqs'): + assert (dot_dir / name).is_file(), name + assert file_bytes(str(dot_dir / 'un-conc-mate.1')) == expected + + +def test_index_requires_x(bt2, lambda_index, example_reads, tmp_path): + """An index given as a positional argument instead of with -x is an error, + and the same command with -x succeeds.""" + sam = str(tmp_path / 'out.sam') + reads = ['-U', example_reads['longreads'], '-S', sam] + proc = bt2.run('bowtie2', [lambda_index] + reads, check=False) + assert proc.returncode != 0 + bt2.run('bowtie2', ['-x', lambda_index] + reads) + assert sam_records(sam)