Repository: acorg/dark-matter

Issue:
Add various options to some scripts
- Add `--sampleName` and `--readGroup` options to `Bowtie2.align` and corresponding options to `run-bowtie2.py`.
- Add `--id` and `--idLambda` options to `make-consensus.py`.
- Add `--aligner` and `--alignerOptions` to `compare-sequences.py`.

And add support for aligning using MAFFT to `compare-sequences.py`.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c6b2654..e78ba44 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,13 @@
+## 3.1.48 Feb 21, 2020
+
+
+
+## 3.1.47 Feb 19, 2020
+
+Added `--sampleName` and `--readGroup` options to `Bowtie2.align`. Added
+`--id` and `--idLambda` options to `make-consensus.py` to make it possible
+to set (or adjust) the name of the generated consensus.
+
 ## 3.1.46 Feb 19, 2020
 
 Trivial change to setting of `tempdir` in `run-bowtie2.py`.
diff --git a/bin/compare-sequences.py b/bin/compare-sequences.py
index d92d09c..47e4a70 100755
--- a/bin/compare-sequences.py
+++ b/bin/compare-sequences.py
@@ -7,6 +7,9 @@ import argparse
 from os.path import join
 from math import log10
 from subprocess import CalledProcessError
+import multiprocessing
+from tempfile import mkdtemp
+from shutil import rmtree
 
 from dark.dna import compareDNAReads, matchToString
 from dark.fasta import FastaReads
@@ -15,29 +18,61 @@ from dark.reads import (Reads, addFASTACommandLineOptions,
 from dark.process import Executor
 from dark.utils import parseRangeString
 
+MAFFT_DEFAULT_ARGS = '--globalpair --maxiterate 1000 --preservecase'
+MAFFT_ALGORITHMS_URL = (
+    'https://mafft.cbrc.jp/alignment/software/algorithms/algorithms.html')
+NEEDLE_DEFAULT_ARGS = 'auto'
 
-def needle(reads, verbose=False):
+
+def mafft(reads, verbose=False, options=None, threads=None):
     """
-    Run a Needleman-Wunsch alignment and return the two sequences.
+    Run a MAFT alignment and return the two sequences.
 
     @param reads: An iterable of two reads.
     @param verbose: If C{True} print progress info to sys.stderr.
+    @param options: A C{str} of options to pass to mafft.
     @return: A C{Reads} instance with the two aligned sequences.
     """
-    from tempfile import mkdtemp
-    from shutil import rmtree
+    tempdir = mkdtemp()
+
+    infile = join(tempdir, 'input.fasta')
+    out = join(tempdir, 'result.fasta')
+
+    reads.save(infile)
+
+    if verbose:
+        print('Running mafft.', file=sys.stderr)
+
+    Executor().execute("mafft --thread %d %s '%s' > '%s'" % (
+        threads, options or '', infile, out))
 
-    dir = mkdtemp()
+    # Use 'list' in the following to force reading the FASTA from disk.
+    result = Reads(list(FastaReads(out)))
+    rmtree(tempdir)
+
+    return result
 
-    file1 = join(dir, 'file1.fasta')
+
+def needle(reads, verbose=False, options=None):
+    """
+    Run a Needleman-Wunsch alignment and return the two sequences.
+
+    @param reads: An iterable of two reads.
+    @param verbose: If C{True} print progress info to sys.stderr.
+    @param options: Additional options to pass to needle.
+    @return: A C{Reads} instance with the two aligned sequences.
+    """
+    tempdir = mkdtemp()
+
+    file1 = join(tempdir, 'file1.fasta')
     with open(file1, 'w') as fp:
         print(reads[0].toString('fasta'), end='', file=fp)
 
-    file2 = join(dir, 'file2.fasta')
+    file2 = join(tempdir, 'file2.fasta')
     with open(file2, 'w') as fp:
         print(reads[1].toString('fasta'), end='', file=fp)
 
-    out = join(dir, 'result.fasta')
+    out = join(tempdir, 'result.fasta')
 
     def useStderr(e):
         return "Sequences too big. Try 'stretcher'" not in e.stderr
@@ -45,9 +80,10 @@ def needle(reads, verbose=False):
     if verbose:
         print('Running needle.', file=sys.stderr)
     try:
-        Executor().execute("needle -asequence '%s' -bsequence '%s' -auto "
-                           "-outfile '%s' -aformat fasta" % (
-                               file1, file2, out), useStderr=useStderr)
+        Executor().execute(
+            "needle -asequence '%s' -bsequence '%s' %s "
+            "-outfile '%s' -aformat fasta" % (
+                file1, file2, options or '', out), useStderr=useStderr)
     except CalledProcessError as e:
         if useStderr(e):
             raise
@@ -62,7 +98,7 @@ def needle(reads, verbose=False):
 
     # Use 'list' in the following to force reading the FASTA from disk.
     result = Reads(list(FastaReads(out)))
-    rmtree(dir)
+    rmtree(tempdir)
 
     return result
 
@@ -81,7 +117,27 @@ parser.add_argument(
 
 parser.add_argument(
     '--align', default=False, action='store_true',
-    help='If given, use needle to do a global alignment of the two sequences.')
+    help=('If given, use mafft (the default) or needle (according to the '
+          'algorithm selected by --aligner) to align the two sequences.'))
+
+parser.add_argument(
+    '--aligner', default='mafft', choices=('mafft', 'needle'),
+    help='The alignment algorithm to use.')
+
+parser.add_argument(
+    '--alignerOptions',
+    help=('Optional arguments to pass to the alignment algorithm. If the '
+          'aligner is mafft, the default options are "%s". If needle, "%s". '
+          'Do not try to set the number of threads here - use the --threads '
+          'argument instead. If you are using mafft, see %s for some possible '
+          'option combinations.' %
+          (MAFFT_DEFAULT_ARGS, NEEDLE_DEFAULT_ARGS, MAFFT_ALGORITHMS_URL)))
+
+parser.add_argument(
+    '--threads', type=int, default=multiprocessing.cpu_count(),
+    help=('The number of threads to use when running the aligner (if --align '
+          'is used and the alignment algorithm can make use of multiple '
+          'threads (mafft can, needle cannot)).'))
 
 parser.add_argument(
     '--alignmentFile',
@@ -136,13 +192,27 @@ if args.alignmentFile:
 if args.align:
     len1, len2 = map(len, reads)
     if len1 == len2:
-        print('Pre-alignment, sequence lengths were identical: %s' % len1)
+        print('Pre-alignment, sequence lengths were identical: %d' % len1)
     else:
         print('Pre-alignment, sequence lengths: %d, %d (difference %d)' % (
             len1, len2, abs(len1 - len2)))
 
-    # Align.
-    reads = needle(reads, args.verbose)
+    if args.aligner == 'mafft':
+        # Be careful in examining args.alignerOptions because we want the
+        # user to be able to pass an empty string (so check against None
+        # before deciding to use the default.)
+        options = (MAFFT_DEFAULT_ARGS if args.alignerOptions is None
+                   else args.alignerOptions)
+        reads = mafft(reads, args.verbose, options=options,
+                      threads=args.threads)
+    else:
+        assert args.aligner == 'needle'
+        # Be careful in examining args.alignerOptions because we want the
+        # user to be able to pass an empty string (so check against None
+        # before deciding to use the default.)
+        options = (NEEDLE_DEFAULT_ARGS if args.alignerOptions is None
+                   else args.alignerOptions)
+        reads = needle(reads, args.verbose, options=options)
 
     if args.alignmentFile:
         assert reads.save(args.alignmentFile) == 2
diff --git a/bin/make-consensus.py b/bin/make-consensus.py
index 74e56a9..4f3dc28 100755
--- a/bin/make-consensus.py
+++ b/bin/make-consensus.py
@@ -5,6 +5,7 @@ import argparse
 from tempfile import mkdtemp
 from os.path import join
 
+from dark.fasta import FastaReads
 from dark.process import Executor
 
 
@@ -17,22 +18,40 @@ def main():
         '--reference', required=True,
         help='The reference FASTA file.')
 
-    parser.add_argument(
+    group = parser.add_mutually_exclusive_group(required=True)
+
+    group.add_argument(
         '--bam',
         help=('The BAM file from which the consensus should be made. Only '
               'used if a VCF file is not provided.'))
 
+    group.add_argument(
+        '--vcfFile',
+        help=('The VCF file. If omitted, bcftools will be used to make a VCF '
+              'file from the BAM file.'))
+
+    group = parser.add_mutually_exclusive_group()
+
+    group.add_argument(
+        '--id',
+        help=('The id to use in the consensus sequence in the output FASTA. '
+              'If not given, the reference sequence id will be used.'))
+
+    group.add_argument(
+        '--idLambda', metavar='LAMBDA-FUNCTION',
+        help=('A one-argument function taking and returning a read id. '
+              'This can be used to set the id of the reference sequence based '
+              'on the id of the reference sequence (the function will be '
+              'called with the id of the reference sequence). E.g., '
+              '--idLambda "lambda id: id.split(\'_\')[0]" or '
+              '--idLambda "lambda id: id[:10] + \'-consensus\'".'))
+
     parser.add_argument(
         '--sample',
         help=('The name of the sample (from the @RG SM tag in the original '
-              'alignment BAM file for which a consensus should be made. '
-              'If not given, the first sample name in the VCF file will be '
-              'used.'))
-
-    parser.add_argument(
-        '--vcfFile',
-        help=('The VCF file. If ommitted, bcftools will be used to make a VCF '
-              'file.'))
+              'alignment BAM file) for which a consensus should be made. '
+              'If not given, the first sample name (from the #CHROM header) '
+              'in the VCF file will be used.'))
 
     parser.add_argument(
         '--dryRun', default=False, action='store_true',
@@ -75,16 +94,25 @@ def main():
     if args.sample:
         sample = args.sample
     else:
-        result = e.execute("gunzip -c '%s' | egrep '^#CHROM' | cut -f10" %
-                           vcfFile)
-
-        sample = result.stdout.split('\n')[0]
+        result = e.execute(
+            "gunzip -c '%s' | egrep -m 1 '^#CHROM' | cut -f10" % vcfFile)
+        sample = result.stdout.strip()
 
+    consensusFile = join(tempdir, 'consensus.fasta')
     result = e.execute(
-        "bcftools consensus --sample '%s' --fasta-ref '%s' '%s'" %
-        (sample, args.reference, vcfFile))
+        "bcftools consensus --sample '%s' --fasta-ref '%s' '%s' > '%s'" %
+        (sample, args.reference, vcfFile, consensusFile))
+
+    consensus = list(FastaReads(consensusFile))[0]
+
+    if args.id is not None:
+        consensus.id = args.id
+    elif args.idLambda is not None:
+        idLambda = eval(args.idLambda)
+        consensus.id = idLambda(args.id)
+
+    print(consensus.toString('fasta'), end='')
 
-    print(result.stdout, end='')
     if result.stderr:
         print(result.stderr, end='', file=sys.stderr)
 
diff --git a/dark/__init__.py b/dark/__init__.py
index 0826037..6688c56 100644
--- a/dark/__init__.py
+++ b/dark/__init__.py
@@ -7,4 +7,4 @@ if sys.version_info < (2, 7):
 # will not be found by the version() function in ../setup.py
 #
 # Remember to update ../CHANGELOG.md describing what's new in each version.
-__version__ = '3.1.46'
+__version__ = '3.1.48'
diff --git a/dark/bowtie2.py b/dark/bowtie2.py
index 0f9d987..2da2bc6 100644
--- a/dark/bowtie2.py
+++ b/dark/bowtie2.py
@@ -59,7 +59,8 @@ class Bowtie2(object):
         self._indexCalled = True
 
     def align(self, bowtie2Args='--no-unal', fastq1=None, fastq2=None,
-              threads=None, discardSAM=False):
+              threads=None, discardSAM=False, readGroup='orig',
+              sampleName='orig'):
         """
         Run Bowtie2.
         """
@@ -73,15 +74,16 @@ class Bowtie2(object):
 
         if fastq1 and fastq2:
             self._executor.execute(
-                "bowtie2 %s --threads %d --rg-id orig --rg SM:orig -x '%s' "
+                "bowtie2 %s --threads %d --rg-id '%s' --rg 'SM:%s' -x '%s' "
                 "-1 '%s' -2 '%s' > '%s'" % (
-                    bowtie2Args, nThreads, self._indexFile, fastq1, fastq2,
-                    samFile))
+                    bowtie2Args, nThreads, readGroup, sampleName,
+                    self._indexFile, fastq1, fastq2, samFile))
         elif fastq1:
             self._executor.execute(
-                "bowtie2 %s --threads %d --rg-id orig --rg SM:orig -x '%s' "
+                "bowtie2 %s --threads %d --rg-id '%s' --rg 'SM:%s' -x '%s' "
                 "-U '%s' > '%s'" % (
-                    bowtie2Args, nThreads, self._indexFile, fastq1, samFile))
+                    bowtie2Args, nThreads, readGroup, sampleName,
+                    self._indexFile, fastq1, samFile))
         else:
             raise ValueError('At least fastq1 must be passed.')
 
@@ -200,15 +202,14 @@ class Bowtie2(object):
         inFile = self._bamFile if which == 'BAM' else self._samFile
         vcfFile = vcfFile or join(self.tempdir, 'output.g.vcf.gz')
 
-        if not referenceFasta:
+        if referenceFasta is None:
             if self._reference:
                 self._report('Using %s as a reference.' % self._reference)
+                referenceFasta = self._reference
             else:
                 raise ValueError('No reference was passed, given to the '
                                  'Bowtie2 __init__, or used in buildIndex.')
 
-        referenceFasta = referenceFasta or self._reference
-
         indexFile = referenceFasta + '.fai'
         if os.path.exists(indexFile):
             removeIndex = False
diff --git a/dark/filter.py b/dark/filter.py
index 7c97870..749d341 100644
--- a/dark/filter.py
+++ b/dark/filter.py
@@ -180,49 +180,50 @@ def addFASTAFilteringCommandLineOptions(parser):
     @param parser: An C{argparse.ArgumentParser} instance.
     """
     parser.add_argument(
-        '--minLength', type=int,
+        '--minLength', type=int, metavar='N',
         help='The minimum sequence length')
 
     parser.add_argument(
-        '--maxLength', type=int,
+        '--maxLength', type=int, metavar='N',
         help='The maximum sequence length')
 
     parser.add_argument(
-        '--whitelist', action='append',
+        '--whitelist', action='append', metavar='SEQUENCE-ID',
         help='Sequence titles (ids) that should be whitelisted')
 
     parser.add_argument(
-        '--blacklist', action='append',
+        '--blacklist', action='append', metavar='SEQUENCE-ID',
         help='Sequence titles (ids) that should be blacklisted')
 
     parser.add_argument(
-        '--whitelistFile',
+        '--whitelistFile', metavar='SEQUENCE-ID-FILE',
         help=('The name of a file that contains sequence titles (ids) that '
               'should be whitelisted, one per line'))
 
     parser.add_argument(
-        '--blacklistFile',
+        '--blacklistFile', metavar='SEQUENCE-ID-FILE',
         help=('The name of a file that contains sequence titles (ids) that '
               'should be blacklisted, one per line'))
 
     parser.add_argument(
-        '--titleRegex', help='A regex that sequence titles (ids) must match.')
+        '--titleRegex', metavar='REGEX',
+        help='A regex that sequence titles (ids) must match.')
 
     parser.add_argument(
-        '--negativeTitleRegex',
+        '--negativeTitleRegex', metavar='REGEX',
         help='A regex that sequence titles (ids) must not match.')
 
     # A mutually exclusive group for --keepSequences and --removeSequences.
     group = parser.add_mutually_exclusive_group()
 
     group.add_argument(
-        '--keepSequences',
+        '--keepSequences', metavar='NUMBER,RANGE,...',
         help=('Specify (1-based) ranges of sequence numbers that should be '
               'kept. E.g., --keepSequences 1-3,5 will output just the 1st, '
               '2nd, 3rd, and 5th sequences. All others will be omitted.'))
 
     group.add_argument(
-        '--removeSequences',
+        '--removeSequences', metavar='NUMBER,RANGE,...',
         help=('Specify (1-based) ranges of sequence numbers that should be '
               'removed. E.g., --removeSequences 1-3,5 will output all but the '
               '1st, 2nd, 3rd, and 5th sequences. All others will be ouput.'))
@@ -252,26 +253,26 @@ def addFASTAFilteringCommandLineOptions(parser):
     # See the docstring for dark.reads.Reads.filter for more detail on
     # randomSubset.
     parser.add_argument(
-        '--randomSubset', type=int,
+        '--randomSubset', type=int, metavar='N',
         help=('An integer giving the number of sequences that should be kept. '
               'These will be selected at random.'))
 
     # See the docstring for dark.reads.Reads.filter for more detail on
     # trueLength.
     parser.add_argument(
-        '--trueLength', type=int,
+        '--trueLength', type=int, metavar='N',
         help=('The number of reads in the FASTA input. Only to be used with '
               'randomSubset'))
 
     parser.add_argument(
-        '--sampleFraction', type=float,
+        '--sampleFraction', type=float, metavar='FRACTION',
         help=('A [0.0, 1.0] C{float} indicating a fraction of the reads that '
               'should be allowed to pass through the filter. The sample size '
               'will only be approximately the product of the sample fraction '
               'and the number of reads. The sample is taken at random.'))
 
     parser.add_argument(
-        '--sequenceNumbersFile',
+        '--sequenceNumbersFile', metavar='FILENAME',
         help=('A file of (1-based) sequence numbers to retain. Numbers must '
               'be one per line.'))
 
diff --git a/dark/sam.py b/dark/sam.py
index 5511b18..bbfdd31 100644
--- a/dark/sam.py
+++ b/dark/sam.py
@@ -303,19 +303,19 @@ class SAMFilter(object):
                   'failures will be included in the output.'))
 
         parser.add_argument(
-            '--minScore', type=float,
+            '--minScore', type=float, metavar='FLOAT',
             help=('If given, alignments with --scoreTag (default AS) values '
                   'less than this value will not be output. If given, '
                   'alignments that do not have a score will not be output.'))
 
         parser.add_argument(
-            '--maxScore', type=float,
+            '--maxScore', type=float, metavar='FLOAT',
             help=('If given, alignments with --scoreTag (default AS) values '
                   'greater than this value will not be output. If given, '
                   'alignments that do not have a score will not be output.'))
 
         parser.add_argument(
-            '--scoreTag', default='AS',
+            '--scoreTag', default='AS', metavar='TAG',
             help=('The alignment tag to extract for --minScore and --maxScore '
                   'comparisons.'))
 


Test patch:
diff --git a/test/test_bowtie2.py b/test/test_bowtie2.py
index 04d4176..1b0a61d 100644
--- a/test/test_bowtie2.py
+++ b/test/test_bowtie2.py
@@ -110,8 +110,48 @@ class TestBowtie2(TestCase):
         bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
         bt.buildIndex('file.fasta')
         bt.align(fastq1='file1.fastq', threads=4)
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
+        self.assertEqual(
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'orig' --rg 'SM:orig' "
+            "-x '/tmp/xxx/index' "
+            "-U 'file1.fastq' > '/tmp/xxx/result.sam'", e.log[-1])
+
+    @patch('os.path.exists')
+    def testAlignOneFASTQWithReadGroup(self, existsMock):
+        """
+        Making a Bowtie index from a file name and running an alignment with
+        one FASTQ file and specifying a read group must result in the expected
+        commands being run and output being produced.
+        """
+        e = Executor(dryRun=True)
+        fp = StringIO()
+        bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
+        bt.buildIndex('file.fasta')
+        bt.align(fastq1='file1.fastq', threads=4, readGroup='xxx')
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
         self.assertEqual(
-            "$ bowtie2 --no-unal --threads 4 --rg-id orig --rg SM:orig "
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'xxx' --rg 'SM:orig' "
+            "-x '/tmp/xxx/index' "
+            "-U 'file1.fastq' > '/tmp/xxx/result.sam'", e.log[-1])
+
+    @patch('os.path.exists')
+    def testAlignOneFASTQWithSampleName(self, existsMock):
+        """
+        Making a Bowtie index from a file name and running an alignment with
+        one FASTQ file and specifying a sample name must result in the expected
+        commands being run and output being produced.
+        """
+        e = Executor(dryRun=True)
+        fp = StringIO()
+        bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
+        bt.buildIndex('file.fasta')
+        bt.align(fastq1='file1.fastq', threads=4, sampleName='xxx')
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
+        self.assertEqual(
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'orig' --rg 'SM:xxx' "
             "-x '/tmp/xxx/index' "
             "-U 'file1.fastq' > '/tmp/xxx/result.sam'", e.log[-1])
 
@@ -126,9 +166,11 @@ class TestBowtie2(TestCase):
         fp = StringIO()
         bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
         bt.buildIndex('file.fasta')
-        bt.align(fastq1='file1.fastq', threads=4, bowtie2Args='--up --down')
+        bt.align(fastq1='file1.fastq', threads=4, bowtie2Args='--up --dn')
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
         self.assertEqual(
-            "$ bowtie2 --up --down --threads 4 --rg-id orig --rg SM:orig -x "
+            "$ bowtie2 --up --dn --threads 4 --rg-id 'orig' --rg 'SM:orig' -x "
             "'/tmp/xxx/index' -U 'file1.fastq' > '/tmp/xxx/result.sam'",
             e.log[-1])
 
@@ -144,8 +186,52 @@ class TestBowtie2(TestCase):
         bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
         bt.buildIndex('file.fasta')
         bt.align(fastq1='file1.fastq', fastq2='file2.fastq', threads=4)
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
+        self.assertEqual(
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'orig' --rg 'SM:orig' -x "
+            "'/tmp/xxx/index' "
+            "-1 'file1.fastq' -2 'file2.fastq' > '/tmp/xxx/result.sam'",
+            e.log[-1])
+
+    @patch('os.path.exists')
+    def testAlignTwoFASTQsWithReadGroup(self, existsMock):
+        """
+        Making a Bowtie index from a file name and running an alignment with
+        two FASTQ files and a read group must result in the expected commands
+        being run and output being produced.
+        """
+        e = Executor(dryRun=True)
+        fp = StringIO()
+        bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
+        bt.buildIndex('file.fasta')
+        bt.align(fastq1='file1.fastq', fastq2='file2.fastq', threads=4,
+                 readGroup='xxx')
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
+        self.assertEqual(
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'xxx' --rg 'SM:orig' -x "
+            "'/tmp/xxx/index' "
+            "-1 'file1.fastq' -2 'file2.fastq' > '/tmp/xxx/result.sam'",
+            e.log[-1])
+
+    @patch('os.path.exists')
+    def testAlignTwoFASTQsWithSampleName(self, existsMock):
+        """
+        Making a Bowtie index from a file name and running an alignment with
+        two FASTQ files and a sample name must result in the expected commands
+        being run and output being produced.
+        """
+        e = Executor(dryRun=True)
+        fp = StringIO()
+        bt = Bowtie2(executor=e, dryRun=True, verboseFp=fp)
+        bt.buildIndex('file.fasta')
+        bt.align(fastq1='file1.fastq', fastq2='file2.fastq', threads=4,
+                 sampleName='xxx')
+        log = fp.getvalue()
+        self.assertTrue(log.endswith('\nAligning with Bowtie2.\n'))
         self.assertEqual(
-            "$ bowtie2 --no-unal --threads 4 --rg-id orig --rg SM:orig -x "
+            "$ bowtie2 --no-unal --threads 4 --rg-id 'orig' --rg 'SM:xxx' -x "
             "'/tmp/xxx/index' "
             "-1 'file1.fastq' -2 'file2.fastq' > '/tmp/xxx/result.sam'",
             e.log[-1])