Merge pull request #105 from J35P312/master

3.6.0

Author: J35P312

setup.py

@@ -20,7 +20,8 @@
 
 setup(
     name = 'tiddit',
-    version = '3.5.2',
+    version = '3.6.0',
+
 
     url = "https://github.com/SciLifeLab/TIDDIT",
     author = "Jesper Eisfeldt",

tiddit/__main__.py

@@ -17,7 +17,7 @@
 import tiddit.tiddit_contig_analysis as tiddit_contig_analysis
 
 def main():
-	version="3.5.2"
+	version="3.6.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")
@@ -155,7 +155,7 @@ def main():
 			args.e=int(library["avg_insert_size"]/2.0)
 
 		t=time.time()
-		sv_clusters=tiddit_cluster.main(prefix,contigs,contig_length,samples,library["mp"],args.e,args.l,max_ins_len,args.min_contig,args.skip_assembly)
+		sv_clusters=tiddit_cluster.main(prefix,contigs,contig_length,samples,library["mp"],args.e,args.l,max_ins_len,args.min_contig,args.skip_assembly,args.r)
 
 		print("generated clusters in")
 		print(time.time()-t)

tiddit/tiddit_cluster.pyx

@@ -37,7 +37,7 @@ def find_discordant_pos(fragment,is_mp):
 
 	return(posA,posB)
 
-def main(prefix,chromosomes,contig_length,samples,is_mp,epsilon,m,max_ins_len,min_contig,skip_assembly):
+def main(prefix,chromosomes,contig_length,samples,is_mp,epsilon,m,max_ins_len,min_contig,skip_assembly,min_reads):
 
 	discordants={}
 	contigs=set([])
@@ -265,13 +265,18 @@ def main(prefix,chromosomes,contig_length,samples,is_mp,epsilon,m,max_ins_len,mi
 				candidates[chrA][chrB][candidate]["N_splits"]=len(candidates[chrA][chrB][candidate]["splits"])
 				candidates[chrA][chrB][candidate]["N_contigs"]=len(candidates[chrA][chrB][candidate]["contigs"])
 
-				if candidates[chrA][chrB][candidate]["N_contigs"]:
+
+				if candidates[chrA][chrB][candidate]["N_splits"] and min_reads <= candidates[chrA][chrB][candidate]["N_splits"]:
+					candidates[chrA][chrB][candidate]["posA"]=mode(candidates[chrA][chrB][candidate]["positions_A"]["splits"])
+					candidates[chrA][chrB][candidate]["posB"]=mode(candidates[chrA][chrB][candidate]["positions_B"]["splits"])
+
+				elif candidates[chrA][chrB][candidate]["N_contigs"]:
 					candidates[chrA][chrB][candidate]["posA"]=mode(candidates[chrA][chrB][candidate]["positions_A"]["contigs"])
 					candidates[chrA][chrB][candidate]["posB"]=mode(candidates[chrA][chrB][candidate]["positions_B"]["contigs"])
+
 				elif candidates[chrA][chrB][candidate]["N_splits"]:
 					candidates[chrA][chrB][candidate]["posA"]=mode(candidates[chrA][chrB][candidate]["positions_A"]["splits"])
-					candidates[chrA][chrB][candidate]["posB"]=mode(candidates[chrA][chrB][candidate]["positions_B"]["splits"])
-
+					candidates[chrA][chrB][candidate]["posB"]=mode(candidates[chrA][chrB][candidate]["positions_B"]["splits"])				
 
 				else:
 					reverse_A = candidates[chrA][chrB][candidate]["positions_A"]["orientation_discordants"].count("True")
@@ -328,10 +333,6 @@ def main(prefix,chromosomes,contig_length,samples,is_mp,epsilon,m,max_ins_len,mi
 						candidates[chrA][chrB][candidate]["posB"]=mode(candidates[chrA][chrB][candidate]["positions_B"]["discordants"])
 
 
-
-					candidates[chrA][chrB][candidate]["posA"]
-					candidates[chrA][chrB][candidate]["posB"]
-
 				candidates[chrA][chrB][candidate]["startB"]=min(candidates[chrA][chrB][candidate]["positions_B"]["start"])
 				candidates[chrA][chrB][candidate]["endB"]=max(candidates[chrA][chrB][candidate]["positions_B"]["end"])
 

tiddit/tiddit_signal.pyx

@@ -65,7 +65,13 @@ def SA_analysis(read,min_q,tag,reference_name):
 
 	cdef long read_query_alignment_end=read.query_alignment_end
 
-	if (read.query_alignment_start ) < (read.query_length - read_query_alignment_end):
+	clip_before=False
+
+	supplementry_alignment=find_SA_query_range(SA_data)
+	if supplementry_alignment.query_alignment_start < read.query_alignment_start:
+		clip_before=True
+
+	if not clip_before:
 		if read.is_reverse:
 
 
@@ -78,7 +84,6 @@ def SA_analysis(read,min_q,tag,reference_name):
 		else:
 			split_pos=read.reference_start+1
 
-	supplementry_alignment=find_SA_query_range(SA_data)
 	SA_chr=SA_data[0]
 
 	startA=read.reference_start+1
@@ -87,8 +92,7 @@ def SA_analysis(read,min_q,tag,reference_name):
 	startB=supplementry_alignment.reference_start
 	endB=supplementry_alignment.reference_end
 
-
-	if (supplementry_alignment.query_alignment_start ) < (supplementry_alignment.query_length - read_query_alignment_end):
+	if clip_before:
 		if SA_data[2] == "-":
 
 			SA_split_pos=supplementry_alignment.reference_start
@@ -192,6 +196,9 @@ def worker(str chromosome, str bam_file_name,str ref,str prefix,int min_q,int ma
 		if read.mate_is_unmapped:
 			continue
 
+		if not read.is_paired:
+			continue
+
 
 		if ( abs(read.isize) > max_ins or mate_chromosome != read_chromosome ):
 			read_query_name=read.query_name