modified: TIDDIT.py

	modified:   src/TIDDIT.cpp
	modified:   src/TIDDIT_calling.py
	modified:   src/TIDDIT_coverage.py
	modified:   src/TIDDIT_signals.py

Author: J35P312

TIDDIT.py

@@ -7,7 +7,7 @@
 sys.path.insert(0, '{}/src/'.format(wd))
 import TIDDIT_calling
 
-version = "2.6.0"
+version = "2.7.0"
 parser = argparse.ArgumentParser("""TIDDIT-{}""".format(version),add_help=False)
 parser.add_argument('--sv'       , help="call structural variation", required=False, action="store_true")
 parser.add_argument('--cov'        , help="generate a coverage bed file", required=False, action="store_true")
@@ -37,14 +37,14 @@
 	args= parser.parse_args()
 	args.wd=os.path.dirname(os.path.realpath(__file__))
 	if args.l < 2:
-		print "error, too low --l value!"
+		print ("error, too low --l value!")
 		quit()
 	if args.ref:
 		if not os.path.isfile(args.ref):
-			print "error,  could not find the reference file"
+			print ("error,  could not find the reference file")
 			quit()
 	if not os.path.isfile(args.bam):
-		print "error,  could not find the bam file"
+		print ("error,  could not find the bam file")
 		quit()
 
 	command_str="{}/bin/TIDDIT --sv -b {} -o {} -p {} -r {} -q {} -n {} -s {}".format(args.wd,args.bam,args.o,args.p,args.r,args.q,args.n,args.s)

src/TIDDIT.cpp

@@ -45,7 +45,7 @@ int main(int argc, char **argv) {
 	int min_variant_size= 100;
 	int sample = 100000000;
 	string outputFileHeader ="output";
-	string version = "2.6.0";
+	string version = "2.7.0";
 
 	//collect all options as a vector
 	vector<string> arguments(argv, argv + argc);

src/TIDDIT_calling.py

@@ -212,11 +212,11 @@ def cluster(args):
 	if not args.e:
 		args.e=int(math.sqrt(library_stats["MeanInsertSize"]*2)*12)
 	n=1
-	print "clustering signals on chromosome:"
+	print ("clustering signals on chromosome:")
 	calls={}
 	for chrA in chromosomes:
 		calls[chrA] =[]
-		print "{}".format(chrA)
+		print ("{}".format(chrA))
 		for chrB in chromosomes:
 			signal_data=numpy.array([ [hit[0],hit[1],hit[2],hit[3],hit[4],hit[5],hit[6],hit[7]] for hit in args.c.execute('SELECT posA,posB,forwardA,qualA,forwardB,qualB,resolution,name FROM TIDDITcall WHERE chrA == \'{}\' AND chrB == \'{}\''.format(chrA,chrB)).fetchall()])
 
@@ -292,8 +292,8 @@ def cluster(args):
 			output="\t".join(output)+"\n"
 			outfile.write(output)
 
-	print "variant clustering completed in {}".format(time.time()-start_time)
-	print "Work complete!"
+	print ("variant clustering completed in {}".format(time.time()-start_time))
+	print ("Work complete!")
 	os.remove("{}.db".format(args.o))
 	return()
 

src/TIDDIT_coverage.py

@@ -16,7 +16,7 @@ def coverage(args):
 	coverage_data={}
 	chromosome_index={}
 	for chromosome in chromosome_size:
-		coverage_data[chromosome]=numpy.zeros( (chromosome_size[chromosome],3) )
+		coverage_data[chromosome]=numpy.zeros( (chromosome_size[chromosome],3), dtype=numpy.float32 )
 		chromosome_index[chromosome]=0
 
 	for line in open(args.o+".tab"):
@@ -35,7 +35,8 @@ def determine_ploidy(args,chromosomes,coverage_data,Ncontent,library_stats):
 	library_stats["chr_cov"]={}
 	ploidies={}
 	avg_coverage=[]
-	cov=[]
+
+	cov=numpy.array([],dtype=numpy.float32)
 	for chromosome in chromosomes:		
 		try:
 			if args.ref:
@@ -46,17 +47,16 @@ def determine_ploidy(args,chromosomes,coverage_data,Ncontent,library_stats):
 
 			if len(chr_cov):
 				chromosomal_average=numpy.median(chr_cov)
-				cov+= list(chr_cov)
+				cov = numpy.append(cov,chr_cov)
 			else:
 				chromosomal_average=0
 			library_stats["chr_cov"][chromosome]=chromosomal_average
 
 		except:
-			print "error: reference mismatch!"
-			print "make sure that the contigs of the bam file and the reference match"
+			print ("error: reference mismatch!")
+			print ("make sure that the contigs of the bam file and the reference match")
 			quit()
 
-	cov=numpy.array(cov)
 	if len(cov):
 		coverage_norm=numpy.median(cov)
 	else:
@@ -93,7 +93,7 @@ def determine_ploidy(args,chromosomes,coverage_data,Ncontent,library_stats):
 
 #normalise the coverage based on GC content
 def gc_norm(args,median_coverage,normalising_chromosomes,coverage_data,Ncontent):
-	gc_vals=[0., 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 , 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1]
+	gc_vals=range(0,101)
 	gc_dict={}
 
 	for gc in gc_vals:
@@ -130,7 +130,7 @@ def retrieve_N_content(args):
 	Ncontent={}
 	chromosome_index={}
 	for chromosome in chromosome_size:
-		Ncontent[chromosome]=numpy.zeros( chromosome_size[chromosome] )
+		Ncontent[chromosome]=numpy.zeros( chromosome_size[chromosome],dtype=numpy.int8 )
 		chromosome_index[chromosome]=0
 
 	for line in open(args.o+".gc.tab"):
@@ -139,7 +139,8 @@ def retrieve_N_content(args):
 
 		content=line.strip().split()
 		contig=content[0]
-		gc=round(float(content[-2]),2)
+
+		gc=int(round(float(content[-2])*100))
 		n=float(content[-1])
 		if n > args.n_mask:
 			Ncontent[contig][chromosome_index[contig]]=-1

src/TIDDIT_signals.py

@@ -36,7 +36,7 @@ def sample(args,coverage_data,library_stats,samfile):
 		index=numpy.random.randint(low=1000, high=len(coverage_data[chromosome][:,0])-1000,size=n)
 		for i in range(0,n):
 
-			median_coverage=numpy.median(coverage_data[chromosome][index[i]:index[i]+region/bin_size,0])
+			median_coverage=numpy.median(coverage_data[chromosome][index[i]:index[i]+int( math.floor(region/bin_size) ),0])
 			roi=[index[i]*bin_size,index[i]*bin_size+region]
 			if median_coverage == 0 or median_coverage > library_stats["Coverage"]*5:
 				continue
@@ -239,15 +239,15 @@ def read_cigar(cigar):
 
 		if SC[i*2+1] == "M":
 			length += int( SC[i*2] )
-			bases=range(0,int( SC[i*2] ))
+			bases=list(range(0,int( SC[i*2] )))
 			for j in range(0,len(bases)):
 				bases[j] += current_pos
 
 			aligned_range += bases
 			current_pos += int( SC[i*2] )
 		elif SC[i*2+1] == "I":
 			insertions+=1
-			bases=range(0,int( SC[i*2] ))
+			bases=list(range(0,int( SC[i*2] )))
 			for j in range(0,len(bases)):
 				bases[j] += current_pos
 			aligned_range += bases