require that both pairs are mapped in the sam file in order to remove the reads
493763
def filter_paired(list):
	
	pairs = {}
	filtered = []
	for id in list:
		read = id.rsplit('/')[0]
		if read not in pairs:
			pairs[read] = []
		pairs[read].append(id)
	for read in pairs:
		ids = pairs[read]
		if len(ids) == 2:
			filtered.extend(ids)
	return set(filtered)
461773
def _check_paired(files, force_single, separators):
    """"""
    full_name = _check_stems(files)
    if files[0].endswith(".bam"):
        return files
    elif is_gsm(files[0]):
        return files
    return combine_pairs(files, force_single, full_name, separators)
192268
def _run_qsnp_paired(align_bams, items, ref_file, assoc_files,
                     region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-paired-variants.vcf" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        out_file = out_file.replace(".gz", "")
        with file_transaction(config, out_file) as tx_out_file:
            with tx_tmpdir(config) as tmpdir:
                with utils.chdir(tmpdir):
                    paired = get_paired_bams(align_bams, items)
                    qsnp = config_utils.get_program("qsnp", config)
                    resources = config_utils.get_resources("qsnp", config)
                    mem = " ".join(resources.get("jvm_opts", ["-Xms750m -Xmx4g"]))
                    qsnp_log = os.path.join(tmpdir, "qsnp.log")
                    qsnp_init = os.path.join(tmpdir, "qsnp.ini")
                    if region:
                        paired = _create_bam_region(paired, region, tmpdir)
                    _create_input(paired, tx_out_file, ref_file, assoc_files['dbsnp'], qsnp_init)
                    cl = ("{qsnp} {mem} -i {qsnp_init} -log {qsnp_log}")
                    do.run(cl.format(**locals()), "Genotyping paired variants with Qsnp", {})
        out_file = _filter_vcf(out_file)
        out_file = bgzip_and_index(out_file, config)
    return out_file
26289
protected UnmappedReads<Read> getUnmappedMatesOfMappedReads(String readsetId) 
      throws GeneralSecurityException, IOException {
    LOG.info("Collecting unmapped mates of mapped reads for injection");
    final Iterable<Read> unmappedReadsIterable = getUnmappedReadsIterator(readsetId); 
    final UnmappedReads<Read> unmappedReads = createUnmappedReads();
    for (Read read : unmappedReadsIterable) {
      unmappedReads.maybeAddRead(read);
    }
    LOG.info("Finished collecting unmapped mates of mapped reads: " + 
        unmappedReads.getReadCount() + " found.");
    return unmappedReads;
  }
236192
def samse(args, opts):
    
    dbfile, readfile = args
    dbfile = check_index(dbfile)
    saifile = check_aln(dbfile, readfile, cpus=opts.cpus)

    samfile, _, unmapped = get_samfile(readfile, dbfile,
                                       bam=opts.bam, unmapped=opts.unmapped)
    if not need_update((dbfile, saifile), samfile):
        logging.error("`{0}` exists. `bwa samse` already run.".format(samfile))
        return "", samfile

    cmd = "bwa samse {0} {1} {2}".format(dbfile, saifile, readfile)
    cmd += " " + opts.extra
    if opts.uniq:
        cmd += " -n 1"

    return cmd, samfile
257665
def remove_readonly(func, path, _):
    ""
    os.chmod(path, stat.S_IWRITE)
    func(path)
959938
def filter_multimappers(align_file, data):
    
    config = dd.get_config(data)
    type_flag = "" if bam.is_bam(align_file) else "S"
    base, ext = os.path.splitext(align_file)
    out_file = base + ".unique" + ext
    bed_file = dd.get_variant_regions(data) or dd.get_sample_callable(data)
    bed_cmd = '-L {0}'.format(bed_file) if bed_file else " "
    if utils.file_exists(out_file):
        return out_file
    base_filter = '-F "not unmapped {paired_filter} and not duplicate and [XA] == null and [SA] == null and not supplementary " '
    if bam.is_paired(align_file):
        paired_filter = "and paired and proper_pair"
    else:
        paired_filter = ""
    filter_string = base_filter.format(paired_filter=paired_filter)
    sambamba = config_utils.get_program("sambamba", config)
    num_cores = dd.get_num_cores(data)
    with file_transaction(out_file) as tx_out_file:
        cmd = ('{sambamba} view -h{type_flag} '
               '--nthreads {num_cores} '
               '-f bam {bed_cmd} '
               '{filter_string} '
               '{align_file} '
               '> {tx_out_file}')
        message = "Removing multimapped reads from %s." % align_file
        do.run(cmd.format(**locals()), message)
    bam.index(out_file, config)
    return out_file
1034078
def remove_nopairs(in_bam, out_dir, config):
    
    runner = broad.runner_from_config(config)
    out_bam = os.path.join(out_dir, "{}-safepair{}".format(*os.path.splitext(os.path.basename(in_bam))))
    if not utils.file_exists(out_bam):
        read_counts = collections.defaultdict(int)
        with pysam.Samfile(in_bam, "rb") as in_pysam:
            for read in in_pysam:
                if read.is_paired:
                    read_counts[read.qname] += 1
        with pysam.Samfile(in_bam, "rb") as in_pysam:
            with file_transaction(out_bam) as tx_out_bam:
                with pysam.Samfile(tx_out_bam, "wb", template=in_pysam) as out_pysam:
                    for read in in_pysam:
                        if read_counts[read.qname] == 2:
                            out_pysam.write(read)
    return runner.run_fn("picard_sort", out_bam, "queryname")
779385
def unmapped(sam, mates):
    
    for read in sam:
        if read.startswith('@') is True:
            continue
        read = read.strip().split()
        if read[2] == '*' and read[6] == '*':
                yield read
        elif mates is True:
            if read[2] == '*' or read[6] == '*':
                yield read
            for i in read:
                if i == 'YT:Z:UP':
                    yield read
623274
def clean_bam_file(bam_in, mask=None):
    
    seq_obj = defaultdict(int)
    if mask:
        mask_file = op.splitext(bam_in)[0] + "_mask.bam"
        if not file_exists(mask_file):
            pybedtools.BedTool(bam_file).intersect(b=mask, v=True).saveas(mask_file)
        bam_in = mask_file
    out_file = op.splitext(bam_in)[0] + "_rmlw.bam"
    # bam.index(bam_in, {'algorithm':{}})
    run("samtools index %s" % bam_in)
    if not file_exists(bam_in + ".bai"):
        raise IOError("Failed to created bam index of %s. Try to do it manually" % bam_in)
    bam_handle = pysam.AlignmentFile(bam_in, "rb")
    with pysam.AlignmentFile(out_file, "wb", template=bam_handle) as out_handle:
        for read in bam_handle.fetch():
            seq_name = int(read.query_name.replace('seq_', ''))
            match_size = [nts for oper, nts in read.cigartuples if oper == 0]
            subs_size = [nts for oper, nts in read.cigartuples if oper == 4]
            if match_size[0] < 17:
                continue
            if subs_size:
                if subs_size[0] > 3:
                    continue
            try:
                nh = read.get_tag('NH')
            except KeyError:
                nh = 1
            seq_obj[seq_name] = sequence(seq_name)
            seq_obj[seq_name].align = nh
            out_handle.write(read)
    return out_file, seq_obj
1062248
def bam_to_fastq(bam_file, is_paired):
    
    out_files, out_handles = _get_fastq_handles(bam_file,
            is_paired)
    if len(out_handles) > 0:
        in_bam = pysam.Samfile(bam_file, mode='rb')
        for read in in_bam:
            num = 1 if (not read.is_paired or read.is_read1) else 2
            # reverse the sequence and quality if mapped to opposite strand
            if read.is_reverse:
                seq = str(Seq.reverse_complement(Seq.Seq(read.seq)))
                qual = "".join(reversed(read.qual))
            else:
                seq = read.seq
                qual = read.qual
            out_handles[num].write("@%s\n%s\n+\n%s\n" % (read.qname,
                seq, qual))
    [h.close() for h in out_handles.values()]
    return out_files
1223411
def _break_reads(self, contig, position, fout, min_read_length=250):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(contig):
            seqs = []
            if read.pos < position < read.reference_end - 1:
                split_point = position - read.pos
                if split_point - 1 >= min_read_length:
                    sequence = mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out).subseq(0, split_point)
                    sequence.id += '.left'
                    seqs.append(sequence)
                if read.query_length - split_point >= min_read_length:
                    sequence = mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out).subseq(split_point, read.query_length)
                    sequence.id += '.right'
                    seqs.append(sequence)
            else:
                seqs.append(mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out))

            for seq in seqs:
                if read.is_reverse:
                    seq.revcomp()
                print(seq, file=fout)
216287
def _all_reads_from_contig(self, contig, fout):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(contig):
            print(mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out), file=fout)
242501
def filter_primary(bam_file, data):
    
    stem, ext = os.path.splitext(bam_file)
    out_file = stem + ".primary" + ext
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            cores = dd.get_num_cores(data)
            cmd = ("samtools view -@ {cores} -F 2304 -b {bam_file} > {tx_out_file}")
            do.run(cmd.format(**locals()), ("Filtering primary alignments in %s." %
                                            os.path.basename(bam_file)))
    return out_file
912820
def ignore(mapping):
    
    if isinstance(mapping, Mapping):
        return AsDict(mapping)
    elif isinstance(mapping, list):
        return [ignore(item) for item in mapping]
    return mapping
823053
def _discarded_reads2_out_file_name(self):
        """"""
        if self.Parameters['-4'].isOn():
            discarded_reads2 = self._absolute(str(self.Parameters['-4'].Value))
        else:
            raise ValueError(
                "No discarded-reads2 (flag -4) output path specified")
        return discarded_reads2
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()))