Align nucleotide sequences with mrtrans and the protein alignment.
195672
def run_mrtrans(align_fasta, recs, work_dir, outfmt="paml"):
    
    align_file = op.join(work_dir, "prot-align.fasta")
    nuc_file = op.join(work_dir, "nuc.fasta")
    output_file = op.join(work_dir, "nuc-align.mrtrans")

    # make the prot_align file and nucleotide file
    align_h0 = open(align_file + "0", "w")
    align_h0.write(str(align_fasta))
    align_h0.close()
    prot_seqs = {}
    i = 0
    for rec in SeqIO.parse(align_h0.name, "fasta"):
        prot_seqs[i] = rec.seq
        i += 1
    align_h = open(align_file, "w")
    for i, rec in enumerate(recs):
        if len(rec.id) > 30:
            rec.id = rec.id[:28] + "_" + str(i)
            rec.description = ""
        print(">{0}\n{1}".format(rec.id, prot_seqs[i]), file=align_h)
    align_h.close()
    SeqIO.write(recs, file(nuc_file, "w"), "fasta")

    # run the program
    cl = MrTransCommandline(align_file, nuc_file, output_file, outfmt=outfmt)
    r, e = cl.run()
    if e is None:
        print("\tpal2nal:", cl, file=sys.stderr)
        return output_file
    elif e.read().find("could not translate") >= 0:
        print("***pal2nal could not translate", file=sys.stderr)
        return None
678433
def aligned_read_to_read(read, revcomp=True, qual=None, ignore_quality=False):
    ''''''
    if read.qual is None or ignore_quality:
        if qual is None or ignore_quality:
            seq = pyfastaq.sequences.Fasta(read.qname, common.decode(read.seq))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)
    else:
        if qual is None:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), common.decode(read.qual))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)

    if read.is_reverse and revcomp:
        seq.revcomp()

    return seq
1041065
def read_alignments(self, input_dir, file_format, header_grep=None):
        

        compression = ['', 'gz', 'bz2']

        if file_format == 'fasta':
            extensions = ['fa', 'fas', 'fasta']

        elif file_format == 'phylip':
            extensions = ['phy']

        else:
            extensions = []

        extensions = list('.'.join([x,y]) if y else x for x,y in itertools.product(extensions, compression))

        files = fileIO.glob_by_extensions(input_dir, extensions)
        files.sort(key=SORT_KEY)
        self._input_files = files
        records = []

        if self.show_progress:
            pbar = setup_progressbar("Loading files", len(files), simple_progress=True)
            pbar.start()

        for i, f in enumerate(files):
            if f.endswith('.gz') or f.endswith('.bz2'):
                fd, tmpfile = tempfile.mkstemp()
                with fileIO.freader(f, f.endswith('.gz'), f.endswith('.bz2')) as reader,\
                     fileIO.fwriter(tmpfile) as writer:
                    for line in reader:
                        if ISPY3:
                            line = line.decode()
                        writer.write(line)
                try:
                    record = Alignment(tmpfile, file_format, True)
                except ValueError:
                    record = Alignment(tmpfile, file_format, False)
                finally:
                    os.close(fd)
                    os.unlink(tmpfile)

            else:
                try:
                    record = Alignment(f, file_format, True)
                except RuntimeError:
                    record = Alignment(f, file_format, False)

            if header_grep:
                try:
                    datatype = 'dna' if record.is_dna() else 'protein'

                    record = Alignment([(header_grep(x), y) for (x, y) in record.get_sequences()], datatype)

                except TypeError:
                    raise TypeError("Couldn't apply header_grep to header\n"
                                    "alignment number={}, name={}\n"
                                    "header_grep={}".format(i, fileIO.strip_extensions(f), header_grep))
                except RuntimeError:
                    print('RuntimeError occurred processing alignment number={}, name={}'
                          .format(i, fileIO.strip_extensions(f)))
                    raise

            record.name = (fileIO.strip_extensions(f))
            records.append(record)
            if self.show_progress:
                pbar.update(i)
        if self.show_progress:
            pbar.finish()
        return records
1071843
def subalignment(alnfle, subtype, alntype="fasta"):
    
    aln = AlignIO.read(alnfle, alntype)
    alnlen = aln.get_alignment_length()
    nseq = len(aln)
    subaln = None
    subalnfile = alnfle.rsplit(".", 1)[0] + "_{0}.{1}".format(subtype, alntype)

    if subtype == "synonymous":
        for j in range( 0, alnlen, 3 ):
            aa = None
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in CODON_TRANSLATION:
                    break
                if aa and CODON_TRANSLATION[codon] != aa:
                    break
                else:
                    aa = CODON_TRANSLATION[codon]
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subtype == "fourfold":
        for j in range( 0, alnlen, 3 ):
            for i in range(nseq):
                codon = str(aln[i, j: j + 3].seq)
                if codon not in FOURFOLD:
                    break
            else:
                if subaln is None:
                    subaln = aln[:, j: j + 3]
                else:
                    subaln += aln[:, j: j + 3]

    if subaln:
        AlignIO.write(subaln, subalnfile, alntype)
        return subalnfile
    else:
        print("No sites {0} selected.".format(subtype), file=sys.stderr)
        return None
1179580
def protein_statistics(self):
        """"""

        # TODO: can i use get_dict here instead
        d = {}

        d['id'] = self.id
        d['sequences'] = [x.id for x in self.sequences]
        d['num_sequences'] = self.num_sequences
        if self.representative_sequence:
            d['representative_sequence'] = self.representative_sequence.id
            d['repseq_gene_name'] = self.representative_sequence.gene_name
            d['repseq_uniprot'] = self.representative_sequence.uniprot
            d['repseq_description'] = self.representative_sequence.description
        d['num_structures'] = self.num_structures
        d['experimental_structures'] = [x.id for x in self.get_experimental_structures()]
        d['num_experimental_structures'] = self.num_structures_experimental
        d['homology_models'] = [x.id for x in self.get_homology_models()]
        d['num_homology_models'] = self.num_structures_homology
        if self.representative_structure:
            d['representative_structure'] = self.representative_structure.id
            d['representative_chain'] = self.representative_chain
            d['representative_chain_seq_coverage'] = self.representative_chain_seq_coverage
            d['repstruct_description'] = self.description
            if self.representative_structure.is_experimental:
                d['repstruct_resolution'] = self.representative_structure.resolution
        d['num_sequence_alignments'] = len(self.sequence_alignments)
        d['num_structure_alignments'] = len(self.structure_alignments)

        return d
290402
def run_align(*data):
    
    work_dir = dd.get_work_dir(data[0][0])
    out_dir = op.join(work_dir, "seqcluster", "prepare")
    seq_out = op.join(out_dir, "seqs.fastq")
    bam_dir = op.join(work_dir, "align")
    new_bam_file = op.join(bam_dir, "seqs.bam")
    tools = dd.get_expression_caller(data[0][0])
    if not file_exists(new_bam_file):
        sample = process_alignment(data[0][0], [seq_out, None])
        bam_file = dd.get_work_bam(sample[0][0])
        shutil.move(bam_file, new_bam_file)
        shutil.move(bam_file + ".bai", new_bam_file + ".bai")
        shutil.rmtree(op.join(bam_dir, sample[0][0]["rgnames"]['sample']))
    for sample in data:
        # sample[0]["align_bam"] = sample[0]["clean_fastq"]
        sample[0]["cluster_bam"] = new_bam_file

    if "mirdeep2" in tools:
        novel_db = mirdeep.run(data)
    return data
705163
def _retranslate(seq):
        '''
        
        '''
        if len(seq['vdj_nt']) % 3 != 0:
            trunc = len(seq['vdj_nt']) % 3
            seq['vdj_nt'] = seq['vdj_nt'][:-trunc]
        seq['vdj_aa'] = Seq(seq['vdj_nt'], generic_dna).translate()
186795
def get_query_align(hit, contig):
      

      # Getting data needed to extract sequences
      query_seq = hit['query_string']
      homo_seq = hit['homo_string']
      sbjct_start = int(hit['sbjct_start'])
      sbjct_end = int(hit['sbjct_end'])
      query_start = int(hit['query_start'])
      query_end = int(hit['query_end'])
      length = int(hit['sbjct_length'])

      # If the alignment doesn't start at the first position data is
      # added to the begnning
      if sbjct_start != 1:
         missing = sbjct_start - 1

         if(query_start >= missing and hit['strand'] != 1
            or hit['strand'] == 1 and missing <= (len(contig) - query_end)):

            # Getting the query sequence.
            # If the the hit is on the other strand the characters
            # are reversed.
            if hit['strand'] == 1:
               start_pos = query_end
               end_pos = query_end + missing
               chars = contig[start_pos:end_pos]
               chars = Blaster.reversecomplement(chars)
            else:
               start_pos = query_start - missing - 1
               end_pos = query_start - 1
               chars = contig[start_pos:end_pos]

            query_seq = chars + str(query_seq)
         else:
            # Getting the query sequence.
            # If the the hit is on the other strand the characters
            # are reversed.
            if hit['strand'] == 1:
               if query_end == len(contig):
                  query_seq = "-" * missing + str(query_seq)
               else:
                  start_pos = query_end
                  chars = contig[start_pos:]
                  chars = Blaster.reversecomplement(chars)

                  query_seq = ("-" * (missing - len(chars))
                               + chars + str(query_seq))
            elif query_start < 3:
               query_seq = "-" * missing + str(query_seq)
            else:
               end_pos = query_start - 2
               chars = contig[0:end_pos]

               query_seq = ("-" * (missing - len(chars))
                            + chars + str(query_seq))

         # Adding to the homo sequence
         spaces = " " * missing
         homo_seq = str(spaces) + str(homo_seq)

      # If the alignment dosen't end and the last position data is
      # added to the end
      if sbjct_end < length:
         missing = length - sbjct_end

         if(missing <= (len(contig) - query_end) and hit['strand'] != 1
            or hit['strand'] == 1 and query_start >= missing):

            # Getting the query sequence.
            # If the the hit is on the other strand the characters
            # are reversed.
            if hit['strand'] == 1:
               start_pos = query_start - missing - 1
               end_pos = query_start - 1
               chars = contig[start_pos:end_pos]
               chars = Blaster.reversecomplement(chars)
            else:
               start_pos = query_end
               end_pos = query_end + missing
               chars = contig[start_pos:end_pos]

            query_seq = query_seq + chars
         else:
            # If the hit is on the other strand the characters are reversed
            if hit['strand'] == 1:
               if query_start < 3:
                  query_seq = query_seq + "-" * missing
               else:
                  end_pos = query_start - 2
                  chars = contig[0:end_pos]
                  chars = Blaster.reversecomplement(chars)

                  query_seq = (query_seq
                               + chars + "-" * (missing - len(chars)))
            elif query_end == len(contig):
               query_seq = query_seq + "-" * missing
            else:
               start_pos = query_end
               chars = contig[start_pos:]

               query_seq = query_seq + chars + "-" * (missing - len(chars))

         # Adding to the homo sequence
         spaces = " " * int(missing)
         homo_seq = str(homo_seq) + str(spaces)

      return query_seq, homo_seq
157097
private static synchronized Map<String, String> getReverseNucleotideTemplateMap(
      String notationSource) {
    Map<String, String> map = new TreeMap<String, String>(String.CASE_INSENSITIVE_ORDER);
    Map<String, String> normalMap = nucleotideTemplates.get(notationSource);
    if (normalMap != null) {
      Set<String> set = normalMap.keySet();
      for (Iterator i = set.iterator(); i.hasNext();) {
        String tmp = (String) i.next();
        map.put(normalMap.get(tmp), tmp);
      }
    }
    return map;
  }
921357
def align(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(align.__doc__)
    p.add_option("--rnaseq", default=False, action="store_true",
                 help="Input is RNA-seq reads, turn splicing on")
    p.add_option("--native", default=False, action="store_true",
                 help="Convert GSNAP output to NATIVE format")
    p.set_home("eddyyeh")
    p.set_outdir()
    p.set_cpus()
    opts, args = p.parse_args(args)

    if len(args) == 2:
        logging.debug("Single-end alignment")
    elif len(args) == 3:
        logging.debug("Paired-end alignment")
    else:
        sys.exit(not p.print_help())

    dbfile, readfile = args[:2]
    outdir = opts.outdir
    assert op.exists(dbfile) and op.exists(readfile)
    prefix = get_prefix(readfile, dbfile)
    logfile = op.join(outdir, prefix + ".log")
    gsnapfile = op.join(outdir, prefix + ".gsnap")
    nativefile = gsnapfile.rsplit(".", 1)[0] + ".unique.native"

    if not need_update((dbfile, readfile), gsnapfile):
        logging.error("`{0}` exists. `gsnap` already run.".format(gsnapfile))
    else:
        dbdir, dbname = check_index(dbfile)
        cmd = "gsnap -D {0} -d {1}".format(dbdir, dbname)
        cmd += " -B 5 -m 0.1 -i 2 -n 3"  # memory, mismatch, indel penalty, nhits
        if opts.rnaseq:
            cmd += " -N 1"
        cmd += " -t {0}".format(opts.cpus)
        cmd += " --gmap-mode none --nofails"
        if readfile.endswith(".gz"):
            cmd += " --gunzip"
        try:
            offset = "sanger" if guessoffset([readfile]) == 33 else "illumina"
            cmd += " --quality-protocol {0}".format(offset)
        except AssertionError:
            pass
        cmd += " " + " ".join(args[1:])
        sh(cmd, outfile=gsnapfile, errfile=logfile)

    if opts.native:
        EYHOME = opts.eddyyeh_home
        if need_update(gsnapfile, nativefile):
            cmd = op.join(EYHOME, "convert2native.pl")
            cmd += " --gsnap {0} -o {1}".format(gsnapfile, nativefile)
            cmd += " -proc {0}".format(opts.cpus)
            sh(cmd)

    return gsnapfile, logfile
756788
def _misalign_split(self,alns):
    
    total = []
    z = 0
    for x in alns:
      z += 1
      exon_num = z
      if self._alignment.strand == '-':
        exon_num = (len(alns)-z)+1
      buffer = {'query':x['query'][0],'target':x['target'][0],'query_quality':x['query_quality'][0],'exon':exon_num}
      if buffer['query'] == '-': buffer['nt'] = buffer['target']
      elif buffer['target'] == '-': buffer['nt'] = buffer['query']
      elif buffer['query'] == buffer['target']: buffer['nt'] = buffer['query']
      elif buffer['query'] != buffer['target']: buffer['nt'] = '*'
      else:
	sys.stderr.write("WARNING unkonwn case\n")
      for i in range(1,len(x['query'])):
        qchar = x['query'][i]
        tchar = x['target'][i]
        qualchar = x['query_quality'][i]
        if qchar != tchar and (qchar != '-' and tchar != '-'):
          #classic mismatch
          #print 'mismatch'
          #print buffer
          total.append(buffer)
          buffer = {'query':qchar,'target':tchar,'query_quality':qualchar,'exon':exon_num}
          buffer['nt'] = '*'
        elif qchar == buffer['nt'] or tchar == buffer['nt']:
          # its a homopolymer match
          buffer['query'] += qchar

          buffer['target'] += tchar
          buffer['query_quality'] += qualchar
          #print 'homopoly'
        else:
          #print 'new thing'
          #print buffer
          total.append(buffer)
          buffer = {'query':qchar,'target':tchar,'query_quality':qualchar,'exon':exon_num}
          if qchar == '-': buffer['nt'] = tchar
          else:  buffer['nt'] = qchar
      total.append(buffer)
    result = [AlignmentErrors.HPAGroup(self,y) for y in total]
    return result
571696
def protein_header_split_generator(elements, headers, ns):
    """"""
    for el in elements:
        header_not_matching = False
        for protein in el.findall('{%s}protein_id' % ns['xmlns']):
            if not any((re.search(h, protein.text) for h in headers)):
                header_not_matching = True
                break
        if header_not_matching:
            formatting.clear_el(el)
        else:
            yield formatting.string_and_clear(el, ns)
1187703
def main():
    
    parser = argparse.ArgumentParser(
        description='monoseq: pretty-printing DNA and protein sequences',
        epilog='If INPUT is in FASTA format, each record is pretty-printed '
        'after printing its name and ANNOTATION (if supplied) is used by '
        'matching chromosome/record name. If INPUT contains a raw sequence, '
        'only the first chromosome in ANNOTATION is used.')
    parser.add_argument(
        'sequence_file', metavar='INPUT', nargs='?', default=sys.stdin,
        type=argparse.FileType('r'), help='file to read sequence(s) from, '
        'can be in FASTA format (default: standard input)')
    parser.add_argument(
        '-b', '--block-length', metavar='LENGTH', dest='block_length',
        type=int, default=10, help='block length in letters (default: 10)')
    parser.add_argument(
        '-l', '--blocks-per-line', metavar='BLOCKS', dest='blocks_per_line',
        type=int, default=6, help='blocks per line (default: 6)')
    parser.add_argument(
        '-a', '--annotation', metavar='POS', dest='annotation', nargs=2,
        action='append', type=int, help='first and last positions of '
        'subsequence to annotate (allowed more than once)')
    parser.add_argument(
        '-e', '--bed', metavar='ANNOTATION', dest='annotation_file',
        type=argparse.FileType('r'), help='file to read annotation from in '
        'BED format')

    args = parser.parse_args()
    pprint(_until_eof(args.sequence_file), annotation=args.annotation,
           annotation_file=args.annotation_file,
           block_length=args.block_length,
           blocks_per_line=args.blocks_per_line)
837484
def translate(dnaseq,host='human',fmtout=str,tax_id=None):
    
    if isinstance(dnaseq,str): 
        dnaseq=Seq.Seq(dnaseq,Alphabet.generic_dna)
    if tax_id is None:
        tax_id=1 # stanndard codon table. ref http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc25
    prtseq=dnaseq.translate(table=tax_id)
    if fmtout is str:
        return str(prtseq)
    else:
        return prtseq
959650
public static <S extends Sequence<S>> Alignment<S> alignGlobalLinear(LinearGapAlignmentScoring scoring,
                                                                         S seq1, S seq2) {
        if (seq1.getAlphabet() != seq2.getAlphabet() ||
                seq1.getAlphabet() != scoring.getAlphabet())
            throw new IllegalArgumentException("Different alphabets.");

        int size1 = seq1.size() + 1,
                size2 = seq2.size() + 1;
        int[] matrix = new int[size1 * (seq2.size() + 1)];

        for (int i = 0; i < size2; ++i)
            matrix[i] = scoring.getGapPenalty() * i;

        for (int j = 1; j < size1; ++j)
            matrix[size2 * j] = scoring.getGapPenalty() * j;

        int i1, i2,
                match, delete, insert;

        for (i1 = 0; i1 < seq1.size(); ++i1)
            for (i2 = 0; i2 < seq2.size(); ++i2) {
                match = matrix[i1 * size2 + i2] +
                        scoring.getScore(seq1.codeAt(i1), seq2.codeAt(i2));
                delete = matrix[i1 * size2 + i2 + 1] + scoring.getGapPenalty();
                insert = matrix[(i1 + 1) * size2 + i2] + scoring.getGapPenalty();
                matrix[(i1 + 1) * size2 + i2 + 1] = Math.max(match, Math.max(delete, insert));
            }

        MutationsBuilder<S> builder = new MutationsBuilder<>(seq1.getAlphabet(), true);

        i1 = seq1.size() - 1;
        i2 = seq2.size() - 1;
        int score = matrix[(i1 + 1) * size2 + i2 + 1];

        while (i1 >= 0 || i2 >= 0) {
            if (i1 >= 0 && i2 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] == matrix[i1 * size2 + i2] +
                            scoring.getScore(seq1.codeAt(i1), seq2.codeAt(i2))) {
                if (seq1.codeAt(i1) != seq2.codeAt(i2))
                    builder.appendSubstitution(i1, seq1.codeAt(i1), seq2.codeAt(i2));
                --i1;
                --i2;
            } else if (i1 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] ==
                            matrix[i1 * size2 + i2 + 1] + scoring.getGapPenalty()) {
                builder.appendDeletion(i1, seq1.codeAt(i1));
                i1--;
            } else if (i2 >= 0 &&
                    matrix[(i1 + 1) * size2 + i2 + 1] ==
                            matrix[(i1 + 1) * size2 + i2] + scoring.getGapPenalty()) {
                builder.appendInsertion(i1 + 1, seq2.codeAt(i2));
                i2--;
            } else
                throw new RuntimeException();
        }

        return new Alignment<>(seq1, builder.createAndDestroy(), new Range(0, seq1.size()), new Range(0, seq2.size()),
                score);
    }
693187
def find_seq_homologues(self, return_raw=False):
        
        
        
        # Get sequence from structure/chain
        # We could use Bio.PDB.Polypeptide?
        
        from Bio.SCOP.Raf import to_one_letter_code
        
        s = self
        seq_iter = s.get_residues()
        seq_str = ''

        for aa in seq_iter:
            if aa.resname in to_one_letter_code:
                seq_str += to_one_letter_code[aa.resname]
        
        # Use BLAST to find homologous sequences with associated
        # structures in the PDB database.
        # Perhaps include local BLAST?
        
        from Bio.Blast.NCBIWWW import qblast

        # Adapt for short query sequences if needed
        # From http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#shortp
        
        if len(seq_str) < 15:
            word_size = 2
            expect = 20000
            matrix_name = 'PAM30'
            filter = None
        else:
            word_size = 3
            expect = 10.0
            matrix_name = 'BLOSUM62'
            filter = 'SEG'
             
        query_result = qblast(  "blastp", "pdb", seq_str, 
                                word_size, expect, matrix_name, filter)
        
        if return_raw:
           return query_result 
            
        # Parse BLAST result to yield results
        # PDBID : (E-Value, Identity, Positives, Gaps, Alignment)
        
        from Bio.Blast import NCBIXML
        
        blast_records = NCBIXML.read(query_result)
        
        results = []
        
        for alignment in blast_records.alignments:
            for hsps in alignment.hsps:
                id_perc = "%s/%s" %(hsps.identities, alignment.length)
                pos_perc = "%s/%s" %(hsps.positives, alignment.length)
                gaps_perc = "%s/%s" %(hsps.gaps, alignment.length)

                pdb_id = alignment.title.split('|')[3]
                e_value = hsps.expect
                
                results.append( ( pdb_id,
                                "%2.5e" %e_value,
                                id_perc, 
                                pos_perc, 
                                gaps_perc,
                                '\n'.join([hsps.query, hsps.match, hsps.sbjct]) ) )

        return results
710082
def align(fastq_file, pair_file, ref_file, names, align_dir, data,
          extra_args=None):
    
    num_hits = 1
    if data["analysis"].lower().startswith("smallrna-seq"):
        num_hits = 1000
    config = data['config']
    out_file = os.path.join(align_dir, "{0}-sort.bam".format(dd.get_sample_name(data)))
    if data.get("align_split"):
        final_file = out_file
        out_file, data = alignprep.setup_combine(final_file, data)
        fastq_file, pair_file = alignprep.split_namedpipe_cls(fastq_file, pair_file, data)
    else:
        final_file = None
        if fastq_file.endswith(".gz"):
            fastq_file = "<(gunzip -c %s)" % fastq_file
            if pair_file:
                pair_file = "<(gunzip -c %s)" % pair_file

    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with postalign.tobam_cl(data, out_file, pair_file is not None) as (tobam_cl, tx_out_file):
            cl = [config_utils.get_program("bowtie", config)]
            cl += _bowtie_args_from_config(data)
            cl += extra_args if extra_args is not None else []
            cl += ["-q",
                   "-v", 2,
                   "-k", num_hits,
                   "-X", 2000, # default is too selective for most data
                   "--best",
                   "--strata",
                   "--sam",
                   ref_file]
            if pair_file:
                cl += ["-1", fastq_file, "-2", pair_file]
            else:
                cl += [fastq_file]
            cl = [str(i) for i in cl]
            fix_rg_cmd = r"samtools addreplacerg -r '%s' -" % novoalign.get_rg_info(data["rgnames"])
            if fix_rg_cmd:
                cmd = " ".join(cl) + " | " + fix_rg_cmd + " | " + tobam_cl
            else:
                cmd = " ".join(cl) + " | " + tobam_cl
            do.run(cmd, "Running Bowtie on %s and %s." % (fastq_file, pair_file), data)
    return out_file