Python pysam.view函数代码示例

    def execute(self, inBam, exclude, readList, outBam, picardOptions=None, JVMmemory=None):    # pylint: disable=W0221
        picardOptions = picardOptions or []

        if tools.samtools.SamtoolsTool().isEmpty(inBam):
            # Picard FilterSamReads cannot deal with an empty input BAM file
            shutil.copyfile(inBam, outBam)
        elif os.path.getsize(readList) == 0:
            # Picard FilterSamReads cannot deal with an empty READ_LIST_FILE
            if exclude:
                shutil.copyfile(inBam, outBam)
            else:
                tmpf = util.file.mkstempfname('.sam')
                if inBam.endswith('.sam'):
                    # output format (sam/bam) is inferred by samtools based on file extension
                    header = pysam.view('-o', tmpf, '-H', '-S', inBam, catch_stdout=False)
                else:
                    header = pysam.view('-o', tmpf, '-H', inBam, catch_stdout=False)
                # pysam.AlignmentFile cannot write an empty file
                # samtools cannot convert SAM -> BAM on an empty file
                # but Picard SamFormatConverter can deal with empty files
                opts = ['INPUT=' + tmpf, 'OUTPUT=' + outBam, 'VERBOSITY=ERROR']
                PicardTools.execute(self, 'SamFormatConverter', opts, JVMmemory='50m')
        else:
            opts = [
                'INPUT=' + inBam, 'OUTPUT=' + outBam, 'READ_LIST_FILE=' + readList,
                'FILTER=' + (exclude and 'excludeReadList' or 'includeReadList'), 'WRITE_READS_FILES=false'
            ]
            PicardTools.execute(self, self.subtoolName, opts + picardOptions, JVMmemory)

    def execute(self, inBam, exclude, readList, outBam, picardOptions=None, JVMmemory=None):
        picardOptions = picardOptions or []

        if os.path.getsize(readList) == 0:
            # Picard FilterSamReads cannot deal with an empty READ_LIST_FILE
            if exclude:
                shutil.copyfile(inBam, outBam)
            else:
                tmpf = util.file.mkstempfname('.sam')
                with open(tmpf, 'wt') as outf:
                    if inBam.endswith('.sam'):
                        header = pysam.view('-H', '-S', inBam)
                    else:
                        header = pysam.view('-H', inBam)
                    for line in header:
                        outf.write(line)
                # pysam.AlignmentFile cannot write an empty file
                # samtools cannot convert SAM -> BAM on an empty file
                # but Picard SamFormatConverter can deal with empty files
                opts = ['INPUT=' + tmpf, 'OUTPUT=' + outBam, 'VERBOSITY=ERROR']
                PicardTools.execute(self, 'SamFormatConverter', opts, JVMmemory='50m')
        else:
            opts = ['INPUT=' + inBam, 'OUTPUT=' + outBam, 'READ_LIST_FILE=' + readList, 'FILTER=' +
                    (exclude and 'excludeReadList' or 'includeReadList'), 'WRITE_READS_FILES=false']
            PicardTools.execute(self, self.subtoolName, opts + picardOptions, JVMmemory)

def checkSamtoolsViewEqual(filename1, filename2,
                           without_header=False):
    '''return true if the two files are equal in their
    content through samtools view.
    '''

    # strip MD and NM tags, as not preserved in CRAM files
    args = ["-x", "MD", "-x", "NM"]
    if not without_header:
        args.append("-h")

    lines1 = pysam.view(*(args + [filename1]))
    lines2 = pysam.view(*(args + [filename2]))

    if len(lines1) != len(lines2):
        return False

    if lines1 != lines2:
        # line by line comparison
        # sort each line, as tags get rearranged between
        # BAM/CRAM
        for n, pair in enumerate(zip(lines1, lines2)):
            l1, l2 = pair
            l1 = sorted(l1[:-1].split("\t"))
            l2 = sorted(l2[:-1].split("\t"))
            if l1 != l2:
                print "mismatch in line %i" % n
                print l1
                print l2
                return False
        else:
            return False

    return True

def filter_bam(in_fpath, out_fpath, min_mapq=0, required_flag_tags=None,
               filtering_flag_tags=None, regions=None):
    cmd = ['-bh']

    # The following line:
    cmd.append('-o' + out_fpath)
    # should be
    # cmd.extend(['-o', out_fpath])
    # but it is a workaround, take a look at:
    # https://groups.google.com/forum/#!msg/pysam-user-group/ooHgIiNVe4c/CcY06d45rzQJ

    if min_mapq:
        cmd.extend(['-q', str(min_mapq)])

    if required_flag_tags:
        flag = create_flag(required_flag_tags)
        cmd.extend(['-f', str(flag)])

    if filtering_flag_tags:
        flag = create_flag(filtering_flag_tags)
        cmd.extend(['-F', str(flag)])

    cmd.extend([in_fpath])

    if regions:
        regions = ['{0}:{1}-{2}'.format(*s) for s in regions.segments]
        cmd.extend(regions)

    pysam.view(*cmd)

 def _generate_bam_file(self, sam_content, file_prefix):
     sam_file = "%s.sam" % file_prefix
     bam_file = "%s.bam" % file_prefix
     sam_fh = open(sam_file, "w")
     sam_fh.write(sam_content)
     sam_fh.close()
     pysam.view("-Sb", "-o%s" % bam_file, sam_file)
     pysam.index(bam_file)

def convert_bam_to_sam(in_file):
    if not is_bam(in_file):
        raise ValueError("Non BAM file passed to convert_sam_to_bam: "
                         "%s" % (in_file))
    out_file = replace_suffix(in_file, ".sam")
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tmp_out_file:
        pysam.view("-h", "-o" + tmp_out_file, in_file)
    return out_file

def bam2sam(in_file):
    """
    converts a bam file to a sam file
    bam2sam("file.bam") -> "file.sam"
    """
    out_file = replace_suffix(in_file, ".sam")
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tmp_out_file:
        pysam.view("-h", "-o" + tmp_out_file, in_file)
    return out_file

 def good_header(self):
   try:#Test file integrity
     header = pysam.view("-H",self.inputFilePath)
     content = pysam.view(self.inputFilePath)
     outFile = open(self.outputFileRoot+".header","w+")
     outFile.write(''.join(header))
     outFile.write(''.join(content))
     outFile.close()
     return True
   except Exception as e:
     print str(e)
     print >> sys.stderr, "Cannot read binary header, please check BAM file.)"
     return False

def _bam_to_sam(local_name, temp_name):
    temp_local = tempfile.NamedTemporaryFile(suffix='.sam', prefix='local_bam_converted_to_sam_')
    fd, temp_temp = tempfile.mkstemp(suffix='.sam', prefix='history_bam_converted_to_sam_')
    os.close(fd)
    try:
        pysam.view('-h', '-o%s' % temp_local.name, local_name)
    except Exception as e:
        raise Exception("Converting local (test-data) BAM to SAM failed: %s" % e)
    try:
        pysam.view('-h', '-o%s' % temp_temp, temp_name)
    except Exception as e:
        raise Exception("Converting history BAM to SAM failed: %s" % e)
    os.remove(temp_name)
    return temp_local, temp_temp

def combine_samfiles(multi=False, clipped=False):
	#Seperate out clipped and unclipped!
	#Look at naming!
	if multi:
		sam1 = "unclipped_multimap.sam"
		sam2 = "clipped_multimap.sam"
		bam1 = "unclipped_multimap.bam"
		bam2 = "clipped_multimap.bam"
		out = open("multi_mapped.sam", "w")
	else:
		sam1 = "unclipped_unique.sam"
		sam2 = "clipped_unique.sam"
		bam1 = "unclipped_unique.bam"
		bam2 = "clipped_unique.bam"
		out = open("unique_mapped.sam", "w")
	#Convert unclipped sam to bam

	#Converts sam to bam
	bam1_o = open(bam1, "w")
	a = pysam.view("-bS", sam1)
	for r in a:                                     
		bam1_o.write(r)
	bam1_o.close()
	#Converts clipped sam to bam
	if clipped == True:
		if os.stat(sam2).st_size > 0: #Checking file is not empty
			try:
				bam2_o = open(bam2, "w")
				b = pysam.view("-bS", sam2)
				for r in b:                                     
					bam2_o.write(r)
				bam2_o.close()
			except:
				print "Samtools raised error, will assume Sam file is empty!"
			#Merge clipped and unclipped
			input_filenames = ["-f", bam1, bam2]
			output_filename = "tmp1.bam"
			merge_parameters = [output_filename] + input_filenames
			pysam.merge(*merge_parameters)
			pysam.sort("-n", "tmp1.bam", "tmp2" )
			subprocess.call(["rm", sam2, bam2])
	else:
		#If no clipped bam, just sort 
		pysam.sort("-n", bam1, "tmp2" )
	#Converts file to sam
	d = pysam.view("-h", "tmp2.bam")
	for r in d:                                     
		out.write(r)
	subprocess.call(["rm", "tmp2.bam", "tmp1.bam", sam1, bam1])

def bed_from_sam(samIN, name):
    file = open(name, 'wa')
    pysam.view("-bS", "-o"+name, samIN)
    pysam.sort(name, name)
    Bamname = name+".bam"
    pysam.index(Bamname)
    #delete Sam file
    os.remove(samIN)
    bamfile = pysam.AlignmentFile(Bamname, "rb")
    for read in bamfile.fetch():
        if read.mapq > mapQ:
            #this may eliminate reads that aligned more than once, to cout how may use the XS flag
            line = "%s\t%i\t%i\n" % (str(read).split("\t")[0], read.reference_start, read.reference_end)
            file.write(line)
    file.close()

def check_inputs(file):
	"""check that input files exist and identify appropriate naming convention for filtering reads by chromosome"""
	try:
		samfile = pysam.Samfile(file, "rb")
	except IOError:
		print 'file not found'
		pass
		
	try:
		pysam.view("-X", file, "chr11:1000-1100")
		head = "chr"
	except NameError:
		head = ""
	
	return samfile, head

def create_bam(filename):
    """
    Function that create a BAM file from a SAM file.

    Args :
        filename [STR] = SAM filename

    Returns:
        bamfile [STR] = BAM filename
    """
    # name of the bam file to create
    bamfile = os.path.dirname(filename)[:-3] + "bam/" + os.path.basename(filename)[:-3] + "bam" 
    # convert sam to bam using pysam
    pysam.view('-Sb',filename, '-o', bamfile, catch_stdout=False)

    return bamfile

    def test_bam_extract_01(self):
        TEST_DIR, T_TEST_DIR = self.__get_temp_dirs()

        input_file = TEST_DIR + "test_terg_02.bam"
        output_file = T_TEST_DIR + "test_terg_02.filtered.bam"
        output_file_s = T_TEST_DIR + "test_terg_02.filtered.sam"
        test_file = TEST_DIR + "test_terg_02.filtered.sam"

        # c = BAMExtract(input_file)
        # c.extract("chr21:39000000-40000000", "chr5:1-2", output_file)
        command = ["bin/dr-disco",
                   "bam-extract",
                   "chr21:39000000-40000000",
                   "chr5:1-2",
                   output_file,
                   input_file]

        self.assertEqual(subprocess.call(command), 0)

        # Bam2Sam
        fhq = open(output_file_s, "w")
        fhq.write(pysam.view(output_file))
        fhq.close()

        if not filecmp.cmp(output_file_s, test_file):
            print 'diff \'' + output_file_s + '\' \'' + test_file + '\''

        self.assertTrue(filecmp.cmp(output_file_s, test_file))

def _get_sort_order(in_bam, config):
    for line in pysam.view("-H", in_bam).split("\r\n"):
        if line.startswith("@HD"):
            for keyval in line.split()[1:]:
                key, val = keyval.split(":")
                if key == "SO":
                    return val