| 
2
 | 
     1 #!/usr/bin/python
 | 
| 
 | 
     2 
 | 
| 
 | 
     3 import argparse, optparse, os, shutil, subprocess, sys, tempfile, shlex, time
 | 
| 
 | 
     4 
 | 
| 
 | 
     5 
 | 
| 
 | 
     6 parser = argparse.ArgumentParser(description='')
 | 
| 
 | 
     7 
 | 
| 
 | 
     8 parser.add_argument ( '-i1', dest='input_bam', help='the bam input file' )
 | 
| 
 | 
     9 parser.add_argument ( '-o1', dest='output_raw', help='the output file' )
 | 
| 
 | 
    10 parser.add_argument ( '-o2', dest='output_vcf', help='the output file' )
 | 
| 
 | 
    11 
 | 
| 
 | 
    12 
 | 
| 
 | 
    13 tmp_dir = tempfile.mkdtemp()
 | 
| 
 | 
    14 errorFile = tmp_dir+"/errorLog"
 | 
| 
 | 
    15 
 | 
| 
 | 
    16 bam2cfg_path = "/home/jeremie/Tools/breakdancer-1.1.2/perl/bam2cfg.pl"
 | 
| 
 | 
    17 breakdancer_path = "/usr/bin/breakdancer-max"
 | 
| 
 | 
    18 breakdancer2vcf_path = "/home/jeremie/Tools/breakdancer-1.1.2/breakdancer2vcf.py"
 | 
| 
 | 
    19 
 | 
| 
 | 
    20 
 | 
| 
 | 
    21 def bam2cfg(args):
 | 
| 
 | 
    22 	config = tmp_dir+"/breakdancer_config"
 | 
| 
 | 
    23 	cmd = 'perl %s %s' % (bam2cfg_path, args.input_bam)
 | 
| 
 | 
    24 	execute(cmd, output=config)
 | 
| 
 | 
    25 	print ("\ncmd = %s  \n" ) %(cmd)
 | 
| 
 | 
    26 	return config
 | 
| 
 | 
    27 
 | 
| 
 | 
    28 
 | 
| 
 | 
    29 def breakdancer(args, config):
 | 
| 
 | 
    30 	cmd = '%s %s' % (breakdancer_path, config)
 | 
| 
 | 
    31 	execute(cmd, output=args.output_raw)
 | 
| 
 | 
    32 
 | 
| 
 | 
    33 
 | 
| 
 | 
    34 def breakdancer2vcf(args):
 | 
| 
 | 
    35 	cmd = "python %s -i %s -o %s" % ( breakdancer2vcf_path, args.output_raw, args.output_vcf )
 | 
| 
 | 
    36 	execute(cmd)
 | 
| 
 | 
    37 
 | 
| 
 | 
    38 
 | 
| 
 | 
    39 def execute( cmd, output="" ):
 | 
| 
 | 
    40 	try: 
 | 
| 
 | 
    41 		err = open(tmp_dir+"/errorLog", 'a')
 | 
| 
 | 
    42 		if output != "":
 | 
| 
 | 
    43 			out = open(output, 'w')
 | 
| 
 | 
    44 		else:
 | 
| 
 | 
    45 			out = subprocess.PIPE
 | 
| 
 | 
    46 		process = subprocess.Popen( args=shlex.split(cmd), stdout=out, stderr=err )
 | 
| 
 | 
    47 		process.wait()
 | 
| 
 | 
    48 		err.close()
 | 
| 
 | 
    49 		if out != subprocess.PIPE:
 | 
| 
 | 
    50 			out.close()
 | 
| 
 | 
    51 	except Exception, e:
 | 
| 
 | 
    52 		sys.stderr.write("problem doing : %s\n" %(cmd))
 | 
| 
 | 
    53 		sys.stderr.write( '%s\n\n' % str(e) )
 | 
| 
 | 
    54 
 | 
| 
 | 
    55 
 | 
| 
 | 
    56 def check( output ):
 | 
| 
 | 
    57 	if (os.path.getsize(output)>0):
 | 
| 
 | 
    58 		return True
 | 
| 
 | 
    59 	else:
 | 
| 
 | 
    60 		sys.stderr.write('The output file is empty : %s\n' % (output))
 | 
| 
 | 
    61 
 | 
| 
 | 
    62 
 | 
| 
 | 
    63 def getLine(file):
 | 
| 
 | 
    64 	try:
 | 
| 
 | 
    65 		f = open(file, 'r')
 | 
| 
 | 
    66 		lignes = f.readlines()
 | 
| 
 | 
    67 		n=0
 | 
| 
 | 
    68 		for ligne in lignes:
 | 
| 
 | 
    69 			if ligne.strip()[0]!="#":
 | 
| 
 | 
    70 				n+=1
 | 
| 
 | 
    71 		# n = len(lignes)
 | 
| 
 | 
    72 		f.close()
 | 
| 
 | 
    73 	except Exception, e:
 | 
| 
 | 
    74 		sys.stderr.write( '%s\n' % str(e) )
 | 
| 
 | 
    75 	return n
 | 
| 
 | 
    76 
 | 
| 
 | 
    77 
 | 
| 
 | 
    78 def compare( output1, output2 ):
 | 
| 
 | 
    79 	# compare le nombre de ligne entre 2 fichiers
 | 
| 
 | 
    80 	# pour verifier qu'aucune ligne n'a ete saute
 | 
| 
 | 
    81 	num_raw = getLine(output1)
 | 
| 
 | 
    82 	num_vcf = getLine(output2)
 | 
| 
 | 
    83 	if (num_raw==num_vcf):
 | 
| 
 | 
    84 		return True
 | 
| 
 | 
    85 	else:
 | 
| 
 | 
    86 		sys.stderr.write('Not the same number of variant between the raw file and the vcf file : %d vs %d\n' % (num_raw, num_vcf))
 | 
| 
 | 
    87 
 | 
| 
 | 
    88 
 | 
| 
 | 
    89 def __main__():
 | 
| 
 | 
    90 
 | 
| 
 | 
    91 	time0 = time.time()
 | 
| 
 | 
    92 	args = parser.parse_args()
 | 
| 
 | 
    93 
 | 
| 
 | 
    94 	try:
 | 
| 
 | 
    95 		config = bam2cfg(args)
 | 
| 
 | 
    96 		breakdancer(args, config)
 | 
| 
 | 
    97 		breakdancer2vcf(args)
 | 
| 
 | 
    98 
 | 
| 
 | 
    99 		#	quelques tests
 | 
| 
 | 
   100 		# verifier que les fichiers de sorties ne sont pas vides
 | 
| 
 | 
   101 		check(args.output_raw)
 | 
| 
 | 
   102 		check(args.output_vcf)
 | 
| 
 | 
   103 		# comparer le nombre de ligne entre les 2 fichiers
 | 
| 
 | 
   104 		# pour etre sur que toute les variations ont ete prises en compte
 | 
| 
 | 
   105 		compare(args.output_raw, args.output_vcf)
 | 
| 
 | 
   106 			
 | 
| 
 | 
   107 		sys.stdout.write( '\nDone in %d seconds\n' %(int(time.time()-time0)))
 | 
| 
 | 
   108 
 | 
| 
 | 
   109 		if ( os.path.getsize(errorFile)>0 ):
 | 
| 
 | 
   110 			sys.stdout.write( 'At least one non fatal error was send, check the file : %s\n' %(errorFile) )
 | 
| 
 | 
   111 			try:
 | 
| 
 | 
   112 				err = open(errorFile, 'r')
 | 
| 
 | 
   113 				errors = err.read()
 | 
| 
 | 
   114 				err.close()
 | 
| 
 | 
   115 				sys.stdout.write( "Errors :\n%s" %(errors))
 | 
| 
 | 
   116 			except Exception, e:
 | 
| 
 | 
   117 				sys.stderr.write( '%s\n' % str(e) )
 | 
| 
 | 
   118 		else:
 | 
| 
 | 
   119 			sys.stdout.write( 'BreakDancer successful' )
 | 
| 
 | 
   120 
 | 
| 
 | 
   121 	except Exception, e:
 | 
| 
 | 
   122 		sys.stdout.write( 'BreakDancer fail\n' )
 | 
| 
 | 
   123 		sys.stderr.write( '%s\n' % str(e) )
 | 
| 
 | 
   124 		sys.exit()
 | 
| 
 | 
   125 
 | 
| 
 | 
   126 	finally:
 | 
| 
 | 
   127 		if os.path.exists( tmp_dir ):
 | 
| 
 | 
   128 			shutil.rmtree( tmp_dir )
 | 
| 
 | 
   129 
 | 
| 
 | 
   130 if __name__=="__main__":
 | 
| 
 | 
   131 	__main__()
 |