Python pysam.Samfile方法代码示例

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def validate_bam_fasta_pairs(bam_path, fasta_sequences, genome):
    """
    Make sure that this BAM is actually aligned to this fasta. Every sequence should be the same length. Sequences
    can exist in the reference that do not exist in the BAM, but not the other way around.
    """
    handle = pysam.Samfile(bam_path, 'rb')
    bam_sequences = {(n, s) for n, s in zip(*[handle.references, handle.lengths])}
    difference = bam_sequences - fasta_sequences
    if len(difference) > 0:
        base_err = 'Error: BAM {} has the following sequence/length pairs not found in the {} fasta: {}.'
        err = base_err.format(bam_path, genome, ','.join(['-'.join(map(str, x)) for x in difference]))
        raise UserException(err)
    missing_seqs = fasta_sequences - bam_sequences
    if len(missing_seqs) > 0:
        base_msg = 'BAM {} does not have the following sequence/length pairs in its header: {}.'
        msg = base_msg.format(bam_path, ','.join(['-'.join(map(str, x)) for x in missing_seqs]))
        logger.warning(msg) 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def get_chrom_sizes(bam_handle):
    """
    return the list of chromosome names and their
    size from the bam file
    The return value is a list of the form
    [('chr1', 2343434), ('chr2', 43432432)]

    >>> test = Tester()
    >>> get_chrom_sizes(pysam.Samfile(test.bam_file_1, 'rb'))
    [('contig-1', 7125), ('contig-2', 3345)]
    """

    # in some cases there are repeated entries in
    # the bam file. Thus, I first convert to dict,
    # then to list.
    list_chrom_sizes = OrderedDict(zip(bam_handle.references,
                                       bam_handle.lengths))
    return list(list_chrom_sizes.items()) 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def main(input_bam, output_bam):
    if input_bam.endswith(".sam") or input_bam.endswith("sam.gz"):
        infile = pysam.Samfile(input_bam, "r")
    else:
        # assume binary BAM file
        infile = pysam.Samfile(input_bam, "rb")

    if output_bam.endswith(".sam"):
        # output in text SAM format
        outfile = pysam.Samfile(output_bam, "w", template=infile)
    elif output_bam.endswith(".bam"):
        # output in binary compressed BAM format
        outfile = pysam.Samfile(output_bam, "wb", template=infile)
    else:
        raise ValueError("name of output file must end with .bam or .sam")

    filter_reads(infile, outfile)

    infile.close()
    outfile.close() 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_get_overlapping_snps_softclip(self):
        """Test that soft-clipped part of read is not used"""
        data = Data()
        data.setup()

        # write a single read with softclipping on left end
        sam_file = open(data.sam_filename, "w")
        data.write_sam_header(sam_file)
        data.write_sam_read(sam_file, cigar="10S20M")
        sam_file.close()
        
        sam_file = pysam.Samfile(data.sam_filename)
        read = next(sam_file)

        snp_tab = snptable.SNPTable()
        snp_tab.read_file(data.snp_filename)
        snp_idx, snp_read_pos, \
            indel_idx, indel_read_pos = snp_tab.get_overlapping_snps(read)

        # check that overlapping SNPs are found and in correct locations
        assert len(snp_idx) == 2
        assert snp_idx[0] == 0
        assert snp_idx[1] == 1
        assert snp_read_pos[0] == 20
        assert snp_read_pos[1] == 30 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_get_overlapping_indel(self):
        """Test that indels can be correctly obtained"""
        data = Data()
        data.snp_list = [(10, "A", "-")]
        data.setup()

        # write a single read with match
        sam_file = open(data.sam_filename, "w")
        data.write_sam_header(sam_file)
        data.write_sam_read(sam_file, cigar="30M")
        sam_file.close()
        
        sam_file = pysam.Samfile(data.sam_filename)
        read = next(sam_file)

        snp_tab = snptable.SNPTable()
        snp_tab.read_file(data.snp_filename)
        snp_idx, snp_read_pos, \
            indel_idx, indel_read_pos = snp_tab.get_overlapping_snps(read)

        # check that overlapping indel found in correct location
        assert len(snp_idx) == 0
        assert len(indel_idx) == 1
        assert indel_idx[0] == 0
        assert indel_read_pos[0] == 10 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def check_bampe_sorted(filename, input_dir):
        infile = pysam.Samfile(input_dir + filename, 'rb')
        count_list = []
        count = 1
        pre_name = ''
        for idx,line in enumerate(infile.fetch(until_eof = True)):
            if line.query_name == pre_name:
                count += 1
            else:
                count_list.append(count)
                count = 1
                pre_name = line.query_name
            if idx == 999:
                break

        count1 = len([i for i in count_list if i==1])
        ratio = float(count1)/sum(count_list)
        #print filename, ratio
        if ratio > 0.8:
            warning("%s may not be sorted by read name. Please check.",filename) 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def parse_bam_for_f_r(filename, input_dir):
    num = 0
    forward = {}
    reverse = {}

    infile =pysam.Samfile(input_dir+filename, 'rb')
    for line in infile.fetch(until_eof = True):
        num += 1
        if num % 10000000 == 0 :
            print ("{0:,} lines processed in {1}".format(num, filename))
        if line.is_unmapped is False:
            chr = infile.getrname(line.tid)
            if line.is_reverse is False:
                try: forward[chr].append(line.reference_start)
                except KeyError:
                    forward[chr] = array.array('i',[line.reference_start])
            else:
                try: reverse[chr].append(line.reference_end)
                except KeyError:
                    reverse[chr] = array.array('i',[line.reference_end])
    return forward,reverse 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def get_chr_info_bam(parameter, filename):
    chr_info = {}
    num = 0
    infile = pysam.Samfile(parameter.input_directory+filename, 'rb')
    for line in infile.fetch(until_eof=True):
        num += 1
        if num %10000000 == 0:
            print("{0:,} lines processed in {1}".format(num, filename))
        if line.is_unmapped is False:
            chr = infile.getrname(line.tid)
            try: 
                chr_info[chr] = max(chr_info[chr], line.pos)
            except KeyError: 
                chr_info[chr] = line.pos
    parameter.chr_info = chr_info
    return 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def get_read_length_from_bam(parameter):
    for filename in parameter.get_filenames():
        with pysam.Samfile(parameter.input_directory+filename, 'rb') as infile:
            length_list = []
            num = 0
            #for idx in range(1000):
            #    line = infile.fetch(until_eof=True).__next__()
            for line in infile.fetch(until_eof=True):
                num += 1
                if num == 1000:
                    break
                try:
                    length_list.append(line.query_length)
                # rlen is deprecated starting from pysam v0.9
                except AttributeError: 
                    length_list.append(line.rlen)
            length = max(length_list)
        parameter.read_length_dict[filename] = length
    return 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def parse_bam_for_f_r(filename):
    num = 0
    forward = {}
    reverse = {}

    infile =pysam.Samfile(filename, 'rb')
    line = infile.fetch(until_eof=True).__next__()
    length = line.alen
    for line in infile.fetch(until_eof = True):
        num += 1
        if num % 1000000 == 0 :
            print ("{0:,} lines processed in {1}".format(num, filename))
        if line.is_unmapped is False:
            chr = infile.getrname(line.tid)
            if line.is_reverse is False:
                try: forward[chr].append(line.pos)
                except KeyError:
                    forward[chr] = array.array('i',[line.pos])
            else:
                try: reverse[chr].append(line.pos)
                except KeyError:
                    reverse[chr] = array.array('i',[line.pos])
    return forward,reverse, length 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_map_reads(self):
        '''test _map_reads'''
        a = assembly.Assembly(contigs_file=os.path.join(data_dir, 'assembly_test.fa'))
        reads_prefix = os.path.join(data_dir, 'assembly_test.to_map')
        out_prefix = 'tmp.assembly_test.out'
        a._map_reads(reads_prefix + '_1.fastq', reads_prefix + '_2.fastq', out_prefix)

        # different smalt version output slightly different BAMs. Some columns
        # should never change, so check just those ones
        def get_sam_columns(bamfile):
            sams = []
            sam_reader = pysam.Samfile(bamfile, "rb")
            for sam in sam_reader.fetch(until_eof=True):
                if sam.is_unmapped:
                    refname = None
                else:
                    refname = sam_reader.getrname(sam.tid)
                sams.append((sam.qname, sam.flag, refname, sam.pos, sam.cigar, sam.seq))
            return sams

        expected = get_sam_columns(os.path.join(data_dir, 'assembly_test.mapped.bam'))
        got = get_sam_columns(out_prefix + '.bam')
        self.assertListEqual(expected, got)
        os.unlink(out_prefix + '.bam') 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_soft_clipped(self):
        '''Test soft_clipped'''
        expected = [
            (5, 0),
            (0, 0),
            (0, 0),
            (0, 5),
            (0, 0),
            None,
            (0, 0),
            (0, 0),
            (2, 0),
            (0, 1),
            None,
            None,
            (1, 1),
            (0, 1)
        ]

        sam_reader = pysam.Samfile(os.path.join(data_dir, 'mapping_test.smalt.out.bam'), "rb")
        i = 0
        for sam in sam_reader.fetch(until_eof=True):
            self.assertEqual(mapping.soft_clipped(sam), expected[i])
            i += 1 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_can_extend(self):
        '''Test can_extend'''
        expected = [
            (True, False),
            (False, False),
            (False, False),
            (False, True),
            (False, False),
            (False, False),
            (False, False),
            (False, False),
            (True, False),
            (False, False),
            (False, False),
            (False, False),
            (False, False),
            (False, False),

        ]

        sam_reader = pysam.Samfile(os.path.join(data_dir, 'mapping_test.smalt.out.bam'), "rb")
        i = 0
        for sam in sam_reader.fetch(until_eof=True):
            self.assertEqual(mapping._can_extend(sam, 190, min_clip=2), expected[i])
            i += 1 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def test_get_pair_type(self):
        '''Test get_pair_type'''
        expected = [
            (mapping.CAN_EXTEND_LEFT, mapping.KEEP),
            (mapping.KEEP, mapping.CAN_EXTEND_RIGHT),
            (mapping.KEEP, mapping.KEEP),
            (mapping.NOT_USEFUL, mapping.NOT_USEFUL),
            (mapping.CAN_EXTEND_LEFT, mapping.KEEP),
            (mapping.BOTH_UNMAPPED, mapping.BOTH_UNMAPPED),
            (mapping.NOT_USEFUL, mapping.NOT_USEFUL)
        ]

        sam_reader = pysam.Samfile(os.path.join(data_dir, 'mapping_test.smalt.out.bam'), "rb")
        previous_sam = None
        i = 0
        for sam in sam_reader.fetch(until_eof=True):
            if previous_sam is None:
                previous_sam = sam
                continue

            types = mapping.get_pair_type(previous_sam, sam, 190, 1000, min_clip=2)
            self.assertEqual(types, expected[i])
            i += 1
            previous_sam = None 

# 需要导入模块: import pysam [as 别名]
# 或者: from pysam import Samfile [as 别名]
def prepare_unmapped_sequences(options):
  start_time = time.time() 
  counter = 0;
  bam = pysam.Samfile(options.bamfile, 'rb')
  fas = open('%s/temp/unmapped.fas' % options.directory, 'wb')
  map = open('%s/temp/unmapped.map' % options.directory, 'wb')
  
  for read in bam:
    if read.is_unmapped and not read.mate_is_unmapped and bam.references[read.rnext].find('chr') == 0:
      fas.write('>read_%d\n%s\n' % (counter, read.seq))
      map.write('%s\tread_%d\n' % (read.qname, counter))
      counter+=1
  fas.close()
  map.close()
  bam.close()
  end_time = time.time() - start_time;
  print "Prepared sequences for searching against HMMs: %fs" % end_time