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

def main(BAM):
#   retreive the region from filename
    geo=BAM.split('.')[4]
#   create two pipe files.
    bampipe =  BAM.rsplit('.',1)[0] + '.pipe'
    bampipetohadoop =  BAM.rsplit('.',1)[0] + '.hadooppipe'
    os.mkfifo(bampipe)
    os.mkfifo(bampipetohadoop)
#   Start 2 subprocesses, one to download file fron swift and one to pipe filterd data to hdfs
    command = 'swift download GenomeData ' + BAM + ' -o -  > '  +  bampipe
    p = subprocess.Popen(command, shell=True)
    command2 = '/usr/local/hadoop/bin/hadoop fs -put -f  - < ' + bampipetohadoop +  '  /genome/' + BAM
    p2 = subprocess.Popen(command2, shell=True)
#   open the hadoop pipe
    f=open(bampipetohadoop,'w')
#   read the swift pipe
    rows = pysam.view('-B', bampipe)
#   call the filter function
    for r in rows:
        fline(r,f,geo)
    f.close()
#   remove pipe files.
    os.remove(bampipe)
    os.remove(bampipetohadoop)
#   create and empty file so that the job can be restarted without processing all files again.
    open( BAM,  'w').close()

    def test_01(self):
        basename = 'multilength_fragments_per_position_001'

        if not os.path.exists('tmp/' + basename + '.bam'):
            fhq = open('tmp/' + basename + '.bam', "wb")
            fhq.write(pysam.view('-bS', 'tests/data/' + basename + ".sam"))
            fhq.close()

        args = CLI(['tmp/' + basename + '.bam', '--verbose'])

        args.parameters.left_padding = 0
        args.parameters.right_padding = 0

        flaimapper = FlaiMapper(args)
        i = 0
        for region in flaimapper.regions():
            self.assertEqual(region.region[0], 'SNORD78')
            for result in region:
                if i == 0:
                    self.assertEqual(region.region[1] + result.start, 11)
                    self.assertEqual(region.region[1] + result.stop, 11 + 61)
                elif i == 1:
                    self.assertEqual(region.region[1] + result.start, 44)
                    self.assertEqual(region.region[1] + result.stop, 44 + 28)
                i += 1
        self.assertEqual(i, 2)

def convert_sam_to_bam():
    """
    This method should take a newly create .sam file from alignment and
        - convert it to .bam
        - sort .bam
        - index .bam
    """
    ids = generate_ids()
    for id in ids:
        start_time = time()
        print 'converting: %s'%id
        base_path = os.path.join(SAMPLE_DIR, id)
        sam_path = os.path.join(base_path, id+'-bwape.sam')
        bam_path = os.path.join(base_path, id+'-bwape.bam')

        bam_content = pysam.view('-bS', sam_path)
        bam_file = open(bam_path, 'w+')
        bam_file.writelines(bam_content)
        bam_file.close()

        pysam.sort(bam_path, bam_path+'_sorted')
        pysam.index(bam_path+'_sorted.bam')

        # indexing creates file.bam.bam. Move it to file.bam
        bam_call = "mv {0} {1}".format(bam_path+'_sorted.bam', bam_path)
        index_call = "mv {0} {1}".format(bam_path+'_sorted.bam.bai',
                                         bam_path+'.bam.bai')
        subprocess.call(bam_call, shell=True)
        subprocess.call(index_call, shell=True)
        end_time = time()
        print 'completed: %.3fs'%(end_time-start_time)

def convert_bam_bed(bam, name, paired, outdir):
	count = 0
	print "==> Converting bam to bed...\n"
#	if aligner=="T":
	outbam = open("{}/{}.unique.bam".format(outdir, name), "wb")
	filtered_bam = pysam.view( "-bq 50", bam) ##Filters for uniquely aligned reads!
	for read in filtered_bam:
		count += 1 
		outbam.write(read)

	inbam = pybedtools.BedTool("{}/{}.unique.bam".format(outdir, name))
	bed = inbam.bam_to_bed(split=True)
	bed.saveas("{}/{}.BED".format(outdir, name))
	#STAR conversion
	#elif aligner=="S":
#		samfile = pysam.Samfile(name+".bam", "rb")
#		for alignedread in samfile.fetch():
#			count += 1
#		samfile.close()
#		inbam = pybedtools.BedTool(name+".bam")
#		bed = inbam.bam_to_bed(split=True)
#		bed.saveas(name+".BED")
	if paired:
		count /= 2
	return count

def only_mapped(filename):
    """
    Function that keep only mapped reads in the BAM file

    Args:
        filename [STR] = BAM file, containing all alignments

    Returns:
        mappedfile [STR] = BAM file with only the mapped reads
    """
    # new bamfile name
    mappedfile = '{0}/filtered_{1}'.format(os.path.dirname(filename),
                                       os.path.basename(filename)[7:])
    # get only mapped reads
    pysam.view('-b', '-F', '4', filename, '-o', mappedfile, catch_stdout=False)
    
    return mappedfile