Mercurial > repos > jeremie > pindel_test

changeset 14:ee3a2e0cc270 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jeremie
date Wed, 09 Jul 2014 10:17:54 -0400
parents 077a2dbd906a
children
files pindel.py
diffstat 1 files changed, 31 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/pindel.py	Wed Jul 09 10:17:45 2014 -0400
+++ b/pindel.py	Wed Jul 09 10:17:54 2014 -0400
@@ -1,14 +1,13 @@
 
-import argparse, optparse, os, shutil, subprocess, sys, tempfile, time, shlex
+import argparse, os, shutil, subprocess, sys, tempfile, time, shlex
 
 parser = argparse.ArgumentParser(description='')
 parser.add_argument('-r', dest='inputFastaFile', required=True, help='the reference file')
 parser.add_argument('-b', dest='inputBamFile', required=True, help='the bam file')
-parser.add_argument('-s', dest='insert_size', type=int, default='500', required=True, help='the insert size')
-parser.add_argument('-n', dest='sampleTag', default='sample', required=False, help='the sample tag')
+parser.add_argument('-s', dest='insert_size', type=int, default=None, required=False, help='the insert size')
+parser.add_argument('-t', dest='sampleTag', default='default', required=False, help='the sample tag')
 # parser.add_argument('-o1', dest='outputRaw', help='the output raw', required=True)
 parser.add_argument('-o2', dest='outputVcfFile', help='the output vcf', required=True)
-parser.add_argument('--name', dest='name', default='name')
 parser.add_argument('--number_of_threads', dest='number_of_threads', type=int, default='1')
 parser.add_argument('--window_size', dest='window_size', type=int, default='5')
 parser.add_argument('--sequencing_error_rate', dest='sequencing_error_rate', type=float, default='0.01')
@@ -31,6 +30,7 @@
 parser.add_argument('-B', '--balance_cutoff', dest='balance_cutoff', type=int, default='0')
 parser.add_argument('-A', '--anchor_quality', dest='anchor_quality', type=int, default='0')
 parser.add_argument('-M', '--minimum_support_for_event', dest='minimum_support_for_event', type=int, default='3')
+parser.add_argument('-n', '--NM', dest='NM', type=int, default='2')
 parser.add_argument('--detect_DD', dest='detect_DD', action='store_true', default=False)
 parser.add_argument('--MAX_DD_BREAKPOINT_DISTANCE', dest='MAX_DD_BREAKPOINT_DISTANCE', type=int, default='350')
 parser.add_argument('--MAX_DISTANCE_CLUSTER_READS', dest='MAX_DISTANCE_CLUSTER_READS', type=int, default='100')
@@ -39,22 +39,8 @@
 parser.add_argument('--MIN_DD_MAP_DISTANCE', dest='MIN_DD_MAP_DISTANCE', type=int, default='8000')
 parser.add_argument('--DD_REPORT_DUPLICATION_READS', dest='DD_REPORT_DUPLICATION_READS', action='store_true', default=False)
 
+parser.add_argument('-z', '--input_SV_Calls_for_assembly', dest='input_SV_Calls_for_assembly', action='store_true', default=False)
 
-# def execute(cmd, output=""):
-# 	try: 
-# 		err = open(tempDir+"/errorLog", 'a')
-# 		if output != "":
-# 			out = open(output, 'w')
-# 		else:
-# 			out = subprocess.PIPE
-# 		process = subprocess.Popen(args=shlex.split(cmd), stdout=out, stderr=err)
-# 		process.wait()
-# 		err.close()
-# 		if out != subprocess.PIPE:
-# 			out.close()
-# 	except Exception, e:
-# 		sys.stderr.write("problem doing : %s\n" %(cmd))
-# 		sys.stderr.write('%s\n\n' % str(e))
 
 def execute(cmd, output=None):
 	import subprocess, sys, shlex
@@ -82,29 +68,14 @@
 	execute(cmd)
 
 
-def config(inputBamFile, meanInsertSize, name, tempDir):
+def config(inputBamFile, meanInsertSize, tag, tempDir):
 	configFile = tempDir+"/pindel_configFile"
-	# cmd = 'echo %s\t%s\t%s > %s' %(bam, size, name, configFile)
-	# print("commande = "+cmd)
-	# os.system(cmd)
 	fil = open(configFile, 'w')
-	fil.write("%s\t%s\t%s" %(inputBamFile, meanInsertSize, name))
+	fil.write("%s\t%s\t%s" %(inputBamFile, meanInsertSize, tag))
 	fil.close()
 	return configFile
 
 
-# def pindel(reference, configFile, output):
-# 	print("pindel")
-# 	cmd="pindel -f %s -i %s -o %s" %(reference, configFile, output)
-# 	print("commande = "+cmd)
-# 	os.system(cmd)
-
-# 	os.system("rm %s" %(output))
-# 	if not os.path.exists(output):
-# 		os.makedirs(output)
-# 	os.system("mv %s_* %s" %(output, output))
-
-
 def pindel(reference, configFile, args, tempDir):
 	pindelTempDir=tempDir+"/pindel" 
 	cmd = "pindel -f %s -i %s -o %s " %(reference, configFile, pindelTempDir)
@@ -161,23 +132,23 @@
 		execute("cp %s %s" %(pindelTempDir+".vcf", outputVcf))
 
 
-def getMeanInsertSize(bamFile, tempDir):
-	samFile = tempDir+"/sam"
-	histogramFile = tempDir+"/histogram"
-	outputFile = tempDir+"/output"
-	cmd = "samtools view -h -o %s %s" % (bamFile, samFile)
-	execute(cmd, )
-	cmd = "picard-tools CollectInsertSizeMetrics H=%s O=%s I=%s" % (histogramFile, outputFile, samFile)
-	execute(cmd)
-	try:
-		f = open(outputFile, 'r')
-		lignes = f.readlines()
-		f.close()
-		meanInsertSize = lignes[7].split('\t')[0]
-	except:
-		sys.stderr.write("problem while opening file %s " %(outputFile))
-		return 0
-	return meanInsertSize
+# def getMeanInsertSize(bamFile, tempDir):
+# 	samFile = tempDir+"/sam"
+# 	histogramFile = tempDir+"/histogram"
+# 	outputFile = tempDir+"/output"
+# 	cmd = "samtools view -h %s -o %s" % (bamFile, samFile)
+# 	execute(cmd, )
+# 	cmd = "picard-tools CollectInsertSizeMetrics H=%s O=%s I=%s" % (histogramFile, outputFile, samFile)
+# 	execute(cmd)
+# 	try:
+# 		f = open(outputFile, 'r')
+# 		lignes = f.readlines()
+# 		f.close()
+# 		meanInsertSize = lignes[7].split('\t')[0]
+# 	except:
+# 		sys.stderr.write("problem while opening file %s " %(outputFile))
+# 		return 0
+# 	return meanInsertSize
 
 
 def __main__():
@@ -185,24 +156,16 @@
 	args = parser.parse_args()
 	tempDir = tempfile.mkdtemp()
 
-	# try:
-	# 	cmd = 'pindel'
-	# 	cmd += ' -i1 %s -i2 %s -o %s ' %(args.input1, args.input2, args.output1)
-	# except:
-	# 	sys.stdout.write('problem args : %s' %(cmd))
-	# try:
-	# 	os.system(cmd)
-	# except:
-	# 	sys.stdout.write('problem cmd')
-	# pipeline(args.reference, args.bam, args.output_raw, args.outputVcfFile, size, name, args)
-	
 	try:
 		indexBam(args.inputFastaFile, args.inputBamFile)
-		meanInsertSize = int(getMeanInsertSize(args.inputBamFile, tempDir))
-		configFile = config(args.inputBamFile, meanInsertSize, args.name, tempDir)
+		# if args.insert_size==None:
+		# 	meanInsertSize = int(getMeanInsertSize(args.inputBamFile, tempDir))
+		# else:
+		meanInsertSize=args.insert_size
+		configFile = config(args.inputBamFile, meanInsertSize, args.sampleTag, tempDir)
 		pindelTempDir = pindel(args.inputFastaFile, configFile, args, tempDir)
-		pindel2vcf(args.outputVcfFile, args.inputFastaFile, args.name, pindelTempDir)
-		
+		pindel2vcf(args.outputVcfFile, args.inputFastaFile, args.sampleTag, pindelTempDir)
+
 	finally:
 		if os.path.exists(tempDir):
 			shutil.rmtree(tempDir)