Perform the alignment of non-mapped reads from previous tier.
486473
def tiered_alignment(in_bam, tier_num, multi_mappers, extra_args,
                     genome_build, pair_stats,
                     work_dir, dirs, config):
    
    nomap_fq1, nomap_fq2 = select_unaligned_read_pairs(in_bam, "tier{}".format(tier_num),
                                                       work_dir, config)
    if nomap_fq1 is not None:
        base_name = "{}-tier{}out".format(os.path.splitext(os.path.basename(in_bam))[0],
                                          tier_num)
        config = copy.deepcopy(config)
        dirs = copy.deepcopy(dirs)
        config["algorithm"]["bam_sort"] = "queryname"
        config["algorithm"]["multiple_mappers"] = multi_mappers
        config["algorithm"]["extra_align_args"] = ["-i", int(pair_stats["mean"]),
                                               int(pair_stats["std"])] + extra_args
        out_bam, ref_file = align_to_sort_bam(nomap_fq1, nomap_fq2,
                                              lane.rg_names(base_name, base_name, config),
                                              genome_build, "novoalign",
                                              dirs, config,
                                              dir_ext=os.path.join("hydra", os.path.split(nomap_fq1)[0]))
        return out_bam
    else:
        return None
644924
def _AlignDecodedDataOffset(self, decoded_data_offset):
    
    self._file_object.seek(0, os.SEEK_SET)

    self._decoder = self._GetDecoder()
    self._decoded_data = b''

    encoded_data_offset = 0
    encoded_data_size = self._file_object.get_size()

    while encoded_data_offset < encoded_data_size:
      read_count = self._ReadEncodedData(self._ENCODED_DATA_BUFFER_SIZE)
      if read_count == 0:
        break

      encoded_data_offset += read_count

      if decoded_data_offset < self._decoded_data_size:
        self._decoded_data_offset = decoded_data_offset
        break

      decoded_data_offset -= self._decoded_data_size
4470
def alignment_ranges(self):
     """"""
     if not self.is_aligned(): raise ValueError("you can't get alignment ranges from something that didn't align")
     if self._alignment_ranges: return self._alignment_ranges
     self._alignment_ranges = self._get_alignment_ranges()
     return self._alignment_ranges
211318
def fastaAlignmentRead(fasta, mapFn=(lambda x : x), l=None):
    
    if l is None:
        l = _getMultiFastaOffsets(fasta)
    else:
        l = l[:]
    seqNo = len(l)
    for i in xrange(0, seqNo):
        j = open(fasta, 'r')
        j.seek(l[i])
        l[i] = j
    column = [sys.maxint]*seqNo
    if seqNo != 0:
        while True:
            for j in xrange(0, seqNo):
                i = l[j].read(1)
                while i == '\n':
                    i = l[j].read(1)
                column[j] = i
            if column[0] == '>' or column[0] == '':
                for j in xrange(1, seqNo):
                    assert column[j] == '>' or column[j] == ''
                break
            for j in xrange(1, seqNo):
                 assert column[j] != '>' and column[j] != ''
                 column[j] = mapFn(column[j])
            yield column[:]
    for i in l:
        i.close()
363101
@Override
  public ArrayList<Read> getUnmappedMates(Read read) {
    if (read.getNumberReads() < 2 ||
        (read.hasNextMatePosition() && read.getNextMatePosition() != null) ||
        !read.hasAlignment() || read.getAlignment() == null ||
        !read.getAlignment().hasPosition() || read.getAlignment().getPosition() == null ||
        read.getAlignment().getPosition().getReferenceName() == null || read.getAlignment().getPosition().getReferenceName().isEmpty() ||
        read.getFragmentName() == null || read.getFragmentName().isEmpty()) {
      return null;
    }
    final String reference = read.getAlignment().getPosition().getReferenceName();
    final String key = getReadKey(read);
    
    Map<String, ArrayList<Read>> reads = unmappedReads.get(reference);
    if (reads != null) {
      final ArrayList<Read> mates = reads.get(key);
      if (mates != null && mates.size() > 1) {
        Collections.sort(mates, matesComparator);
      }
      return mates;
    }
    return null;
  }
1239141
def merge_pairs(self):
        
        logging.info('Merging paired reads')
        for sample in self.metadata:
            # Can only merge paired-end
            if len(sample.general.fastqfiles) == 2:
                # Set the name of the merged, and unmerged files
                sample.general.mergedreads = \
                    os.path.join(sample.general.outputdirectory, '{}_paired.fastq.gz'.format(sample.name))
                sample.general.unmergedforward = \
                    os.path.join(sample.general.outputdirectory, '{}_unpaired_R1.fastq.gz'.format(sample.name))
                sample.general.unmergedreverse = \
                    os.path.join(sample.general.outputdirectory, '{}_unpaired_R2.fastq.gz'.format(sample.name))
                try:
                    # Run the merging command - forward_in=sample.general.normalisedreads[0],
                    out, err, cmd = bbtools.bbmerge(forward_in=sorted(sample.general.trimmedcorrectedfastqfiles)[0],
                                                    merged_reads=sample.general.mergedreads,
                                                    returncmd=True,
                                                    outu1=sample.general.unmergedforward,
                                                    outu2=sample.general.unmergedreverse,
                                                    threads=self.cpus)
                    sample[self.analysistype].bbmergecmd = cmd
                    write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                except CalledProcessError:
                    delattr(sample.general, 'mergedreads')
                    delattr(sample.general, 'unmergedforward')
                    delattr(sample.general, 'unmergedreverse')
                except IndexError:
                    delattr(sample.general, 'mergedreads')
                    delattr(sample.general, 'unmergedforward')
                    delattr(sample.general, 'unmergedreverse')
            else:
                sample.general.mergedreads = sorted(sample.general.trimmedcorrectedfastqfiles)[0]
965669
def _load_data_alignment(self, chain1, chain2):
        
        parser = PDB.PDBParser(QUIET=True)
        ppb = PDB.PPBuilder()

        structure1 = parser.get_structure(chain1, self.pdb1)
        structure2 = parser.get_structure(chain2, self.pdb2)

        seq1 = str(ppb.build_peptides(structure1)[0].get_sequence())
        seq2 = str(ppb.build_peptides(structure2)[0].get_sequence())

        # Alignment parameters taken from PconsFold renumbering script.
        align = pairwise2.align.globalms(seq1, seq2, 2, -1, -0.5, -0.1)[0]
        indexes = set(i for i, (s1, s2) in enumerate(zip(align[0], align[1]))
                      if s1 != '-' and s2 != '-')
        coord1 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure1.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)
        coord2 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure2.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)

        self.coord1 = coord1
        self.coord2 = coord2
        self.N = len(seq1)
229202
def fix_nonref_positions(in_file, ref_file):
    
    ignore_chrs = ["."]
    ref2bit = twobit.TwoBitFile(open(ref_file))
    out_file = in_file.replace("-raw.vcf", ".vcf")

    with open(in_file) as in_handle:
        with open(out_file, "w") as out_handle:
            for line in in_handle:
                if line.startswith("#"):
                    out_handle.write(line)
                else:
                    parts = line.rstrip("\r\n").split("\t")
                    pos = int(parts[1])
                    # handle chr/non-chr naming
                    if parts[0] not in ref2bit.keys() and parts[0].replace("chr", "") in ref2bit.keys():
                        parts[0] = parts[0].replace("chr", "")
                    # handle X chromosome
                    elif parts[0] not in ref2bit.keys() and parts[0] == "23":
                        for test in ["X", "chrX"]:
                            if test in ref2bit.keys():
                                parts[0] == test
                    ref_base = None
                    if parts[0] not in ignore_chrs:
                        try:
                            ref_base = ref2bit[parts[0]].get(pos-1, pos).upper()
                        except Exception as msg:
                            print "Skipping line. Failed to retrieve reference base for %s\n%s" % (str(parts), msg)
                    parts = fix_vcf_line(parts, ref_base)
                    if parts is not None:
                        out_handle.write("\t".join(parts) + "\n")
        return out_file
100854
def align
      alignment_1 = []
      alignment_2 = []

      traceback do |i, j, flags|
        seq1_val = seq2_val = @skip_obj
        case flags
        when :align
          seq1_val = @seq1[i]
          seq2_val = @seq2[j]
        when :insert
          seq2_val = @seq2[j]
        when :delete
          seq1_val = @seq1[i]
        end
        alignment_1.unshift seq1_val
        alignment_2.unshift seq2_val
      end

      [alignment_1, alignment_2]
    end
527735
def align(self, out_path=None):
        """"""
        if out_path is None: out_path = self.prefix_path + '.aln'
        sh.muscle38("-in", self.path, "-out", out_path)
        return AlignedFASTA(out_path)
1256216
def DeserializeUnsigned(self, reader):
        
        self.Version = reader.ReadUInt32()
        self.PrevHash = reader.ReadUInt256()
        self.MerkleRoot = reader.ReadUInt256()
        self.Timestamp = reader.ReadUInt32()
        self.Index = reader.ReadUInt32()
        self.ConsensusData = reader.ReadUInt64()
        self.NextConsensus = reader.ReadUInt160()
959604
def align(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(align.__doc__)
    p.set_firstN(firstN=0)
    p.add_option("--full", default=False, action="store_true",
                 help="Enforce end-to-end alignment [default: local]")
    p.add_option("--reorder", default=False, action="store_true",
                 help="Keep the input read order [default: %default]")
    p.add_option("--null", default=False, action="store_true",
                 help="Do not write to SAM/BAM output")
    p.add_option("--fasta", default=False, action="store_true",
                 help="Query reads are FASTA")
    p.set_cutoff(cutoff=800)
    p.set_mateorientation(mateorientation="+-")
    p.set_sam_options(bowtie=True)

    opts, args = p.parse_args(args)
    extra = opts.extra
    mo = opts.mateorientation
    if mo == '+-':
        extra += ""
    elif mo == '-+':
        extra += "--rf"
    else:
        extra += "--ff"

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

    firstN = opts.firstN
    mapped = opts.mapped
    unmapped = opts.unmapped
    fasta = opts.fasta
    gl = "--end-to-end" if opts.full else "--local"

    dbfile, readfile = args[0:2]
    dbfile = check_index(dbfile)
    prefix = get_prefix(readfile, dbfile)
    samfile, mapped, unmapped = get_samfile(readfile, dbfile, bowtie=True,
                                            mapped=mapped, unmapped=unmapped,
                                            bam=opts.bam)
    logfile = prefix + ".log"
    if not fasta:
        offset = guessoffset([readfile])

    if not need_update(dbfile, samfile):
        logging.error("`{0}` exists. `bowtie2` already run.".format(samfile))
        return samfile, logfile

    cmd = "bowtie2 -x {0}".format(dbfile)
    if PE:
        r1, r2 = args[1:3]
        cmd += " -1 {0} -2 {1}".format(r1, r2)
        cmd += " --maxins {0}".format(opts.cutoff)
        mtag, utag = "--al-conc", "--un-conc"
    else:
        cmd += " -U {0}".format(readfile)
        mtag, utag = "--al", "--un"

    if mapped:
        cmd += " {0} {1}".format(mtag, mapped)
    if unmapped:
        cmd += " {0} {1}".format(utag, unmapped)

    if firstN:
        cmd += " --upto {0}".format(firstN)
    cmd += " -p {0}".format(opts.cpus)
    if fasta:
        cmd += " -f"
    else:
        cmd += " --phred{0}".format(offset)
    cmd += " {0}".format(gl)
    if opts.reorder:
        cmd += " --reorder"

    cmd += " {0}".format(extra)
    # Finally the log
    cmd += " 2> {0}".format(logfile)

    if opts.null:
        samfile = "/dev/null"

    cmd = output_bam(cmd, samfile)
    sh(cmd)
    print(open(logfile).read(), file=sys.stderr)

    return samfile, logfile
299692
def write_alignment(self, filename, file_format, interleaved=None):
        
        if file_format == 'phylip':
            file_format = 'phylip-relaxed'
        AlignIO.write(self._msa, filename, file_format)
48365
public Alignment<S> invert(S sequence2) {
        return new Alignment<>(sequence2, getRelativeMutations().invert().move(sequence2Range.getFrom()),
                sequence2Range, sequence1Range, score);
    }
103937
def _getReadAlignments(
            self, reference, start, end, readGroupSet, readGroup):
        
        # TODO If reference is None, return against all references,
        # including unmapped reads.
        samFile = self.getFileHandle(self._dataUrl)
        referenceName = reference.getLocalId().encode()
        # TODO deal with errors from htslib
        start, end = self.sanitizeAlignmentFileFetch(start, end)
        readAlignments = samFile.fetch(referenceName, start, end)
        for readAlignment in readAlignments:
            tags = dict(readAlignment.tags)
            if readGroup is None:
                if 'RG' in tags:
                    alignmentReadGroupLocalId = tags['RG']
                    readGroupCompoundId = datamodel.ReadGroupCompoundId(
                        readGroupSet.getCompoundId(),
                        str(alignmentReadGroupLocalId))
                yield self.convertReadAlignment(
                    readAlignment, readGroupSet, str(readGroupCompoundId))
            else:
                if self._filterReads:
                    if 'RG' in tags and tags['RG'] == self._localId:
                        yield self.convertReadAlignment(
                            readAlignment, readGroupSet,
                            str(readGroup.getCompoundId()))
                else:
                    yield self.convertReadAlignment(
                        readAlignment, readGroupSet,
                        str(readGroup.getCompoundId()))
958031
def alignments(self):
        
        referenceIds = self.referenceIds
        dropUnmapped = self.dropUnmapped
        dropSecondary = self.dropSecondary
        dropSupplementary = self.dropSupplementary
        dropDuplicates = self.dropDuplicates
        keepQCFailures = self.keepQCFailures
        storeQueryIds = self.storeQueryIds
        filterRead = self.filterRead
        minScore = self.minScore
        maxScore = self.maxScore
        scoreTag = self.scoreTag

        if storeQueryIds:
            self.queryIds = queryIds = set()

        lastAlignment = None
        count = 0
        with samfile(self.filename) as samAlignment:
            for count, alignment in enumerate(samAlignment.fetch(), start=1):
                if storeQueryIds:
                    queryIds.add(alignment.query_name)

                if minScore is not None or maxScore is not None:
                    try:
                        score = alignment.get_tag(scoreTag)
                    except KeyError:
                        continue
                    else:
                        if ((minScore is not None and score < minScore) or
                                (maxScore is not None and score > maxScore)):
                            continue

                # Secondary and supplementary alignments may have a '*'
                # (pysam returns this as None) SEQ field, indicating that
                # the previous sequence should be used. This is best
                # practice according to section 2.5.2 of
                # https://samtools.github.io/hts-specs/SAMv1.pdf So we use
                # the last alignment query and quality strings if we get
                # None as a query sequence.
                if alignment.query_sequence is None:
                    if lastAlignment is None:
                        raise InvalidSAM(
                            'pysam produced an alignment (number %d) with no '
                            'query sequence without previously giving an '
                            'alignment with a sequence.' % count)
                    # Use the previous query sequence and quality. I'm not
                    # making the call to _hardClip dependent on
                    # alignment.cigartuples (as in the else clause below)
                    # because I don't think it's possible for
                    # alignment.cigartuples to be None in this case. If we
                    # have a second match on a query, then it must be
                    # aligned to something (i.e., it cannot be unmapped
                    # with no CIGAR string). The assertion will tell us if
                    # this is ever not the case.
                    assert alignment.cigartuples
                    (alignment.query_sequence,
                     alignment.query_qualities, _) = _hardClip(
                         lastAlignment.query_sequence,
                         lastAlignment.query_qualities,
                         alignment.cigartuples)
                else:
                    lastAlignment = alignment
                    if alignment.cigartuples:
                        (alignment.query_sequence,
                         alignment.query_qualities, _) = _hardClip(
                             alignment.query_sequence,
                             alignment.query_qualities,
                             alignment.cigartuples)

                if ((filterRead is None or
                     filterRead(Read(alignment.query_name,
                                     alignment.query_sequence,
                                     alignment.qual))) and
                    not (
                        (referenceIds and
                         alignment.reference_name not in referenceIds) or
                        (alignment.is_unmapped and dropUnmapped) or
                        (alignment.is_secondary and dropSecondary) or
                        (alignment.is_supplementary and dropSupplementary) or
                        (alignment.is_duplicate and dropDuplicates) or
                        (alignment.is_qcfail and not keepQCFailures))):
                    yield alignment

        self.alignmentCount = count
671134
public Alignment<S> move(int offset) {
        return new Alignment<>(sequence1, mutations, sequence1Range, sequence2Range.move(offset), score);
    }