Merge pull request #90 from J35P312/master

TIDDIT 3.3.0

Author: J35P312

README.md

@@ -10,6 +10,7 @@ INSTALLATION
 ==============
 TIDDIT requires python3 (=> 3.8), cython, pysam, and Numpy.
 
+
 By default, tiddit will require, bwa, fermi2 and ropebwt2 for local assembly; local assembly may be disabled through the "--skip_assembly" parameter.
 
 Installation
@@ -74,6 +75,7 @@ TIDDIT may be fine-tuned by altering these optional parameters:
 	-d	expected reads orientations, possible values "innie" (-> <-) or "outtie" (<- ->). Default: major orientation within the dataset
 	-p	Minimum number of supporting pairs in order to call a variant (default 3)
 	-r	Minimum number of supporting split reads to call a variant (default 3)
+	--threads	Number of threads (default 1)
 	-q	Minimum mapping quality to consider an alignment (default 5)
 	-n	the ploidy of the organism,(default = 2)
 	-e	clustering distance parameter, discordant pairs closer than this distance are considered to belong to the same variant(default = sqrt(insert-size*2)*12)

setup.py

@@ -20,7 +20,8 @@
 
 setup(
     name = 'tiddit',
-    version = '3.2.1',
+    version = '3.3.0',
+
     url = "https://github.com/SciLifeLab/TIDDIT",
     author = "Jesper Eisfeldt",
     author_email= "[email protected]",

tiddit/__main__.py

@@ -17,7 +17,7 @@
 import tiddit.tiddit_contig_analysis as tiddit_contig_analysis
 
 def main():
-	version="3.2.1"
+	version="3.3.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")
@@ -32,6 +32,7 @@ def main():
 		parser.add_argument('-i', type=int, help="paired reads maximum allowed insert size. Pairs aligning on the same chr at a distance higher than this are considered candidates for SV (default= 99.9th percentile of insert size)")
 		parser.add_argument('-d', type=str,help="expected reads orientations, possible values \"innie\" (-> <-) or \"outtie\" (<- ->). Default: major orientation within the dataset")
 		parser.add_argument('-p', type=int,default=3, help="Minimum number of supporting pairs in order to call a variant (default 3)")
+		parser.add_argument('--threads', type=int,default=1, help="Number of threads (default=1)")
 		parser.add_argument('-r', type=int,default=3, help="Minimum number of supporting split reads to call a variant (default 3)")
 		parser.add_argument('-q', type=int,default=5, help="Minimum mapping quality to consider an alignment (default 5)")
 		parser.add_argument('-n', type=int,default=2, help="the ploidy of the organism,(default = 2)")
@@ -90,6 +91,7 @@ def main():
 
 		bam_file_name=args.bam
 		samfile = pysam.AlignmentFile(bam_file_name, "r",reference_filename=args.ref)
+
 		bam_header=samfile.header
 		samfile.close()
 
@@ -110,14 +112,16 @@ def main():
 			contigs.append(contig["SN"])
 			contig_number[contig["SN"]]=i
 			contig_length[ contig["SN"] ]=contig["LN"]
-			i+=0
+			i+=1
 
 		prefix=args.o
 		try:
 			os.mkdir( "{}_tiddit".format(prefix) )
 			os.mkdir("{}_tiddit/clips".format(prefix) )
 		except:
 			print("Folder already exists")
+
+		pysam.index("-c","-m","6","-@",str(args.threads),bam_file_name,"{}_tiddit/{}.csi".format(args.o,sample_id))
 
 		min_mapq=args.q
 		max_ins_len=100000
@@ -131,7 +135,7 @@ def main():
 
 
 		t=time.time()
-		coverage_data=tiddit_signal.main(bam_file_name,args.ref,prefix,min_mapq,max_ins_len,sample_id)
+		coverage_data=tiddit_signal.main(bam_file_name,args.ref,prefix,min_mapq,max_ins_len,sample_id,args.threads,args.min_contig)
 		print("extracted signals in:")
 		print(t-time.time())
 
@@ -163,7 +167,6 @@ def main():
 		f.write(vcf_header+"\n")
 
 		t=time.time()
-		#print(sv_clusters)
 		variants=tiddit_variant.main(bam_file_name,sv_clusters,args,library,min_mapq,samples,coverage_data,contig_number,max_ins_len)
 		print("analyzed clusters in")
 		print(time.time()-t)
@@ -203,7 +206,12 @@ def main():
 			t=time.time()
 			if read.mapq >= args.q:
 				n_reads+=1
-				coverage_data[read.reference_name]=tiddit_coverage.update_coverage(read,args.z,coverage_data[read.reference_name],args.q,end_bin_size[read.reference_name])
+
+				read_position=read.reference_start
+				read_end=read.reference_end
+				read_reference_name=read.reference_name
+
+				coverage_data[read_reference_name]=tiddit_coverage.update_coverage(read_position,read_end,args.z,coverage_data[read_reference_name],end_bin_size[read_reference_name])
 
 		if args.w:
 			tiddit_coverage.print_coverage(coverage_data,bam_header,args.z,"wig",args.o +".wig")

tiddit/tiddit_contig_analysis.pyx

@@ -28,7 +28,14 @@ def read_contigs(aligned_contigs,prefix,sample_id,min_size):
 			continue
 
 		if read.has_tag("SA") and not (read.is_supplementary or read.is_secondary):
-			split_contigs=tiddit_signal.SA_analysis(read,-2,split_contigs,"SA")
+			split=tiddit_signal.SA_analysis(read,-2,"SA",read.reference_name)
+
+			if split:
+				if not split[2] in split_contigs[split[0]][split[1]]:
+					split_contigs[split[0]][split[1]][split[2]]=[]
+				split_contigs[split[0]][split[1]][split[2]]+=split[3:]
+
+			
 		elif read.has_tag("XA") and not (read.is_supplementary or read.is_secondary):
 			XA=read.get_tag("XA")
 			if XA.count(";") == 1:
@@ -44,7 +51,12 @@ def read_contigs(aligned_contigs,prefix,sample_id,min_size):
 					XA=",".join(xa_list)
 
 				read.set_tag("XA",XA)
-				split_contigs=tiddit_signal.SA_analysis(read,-2,split_contigs,"XA")
+				split=tiddit_signal.SA_analysis(read,-2,"XA",read.reference_name)
+
+				if split:
+					if not split[2] in split_contigs[split[0]][split[1]]:
+						split_contigs[split[0]][split[1]][split[2]]=[]
+					split_contigs[split[0]][split[1]][split[2]]+=split[3:]
 
 		elif not (read.is_supplementary or read.is_secondary) and len(read.cigartuples) > 2:
 
@@ -114,9 +126,9 @@ def main(prefix,sample_id,library,contigs,coverage_data,args):
 	f.close()
 	del clips
 
-	os.system("{} -dNCr {}_tiddit/clips.fa | {} assemble -l 81 - > {}_tiddit/clips.fa.assembly.mag".format(args.ropebwt2,prefix,args.fermi2,prefix))
+	os.system("{} -dNCr {}_tiddit/clips.fa | {} assemble -t {} -l 81 - > {}_tiddit/clips.fa.assembly.mag".format(args.ropebwt2,prefix,args.fermi2,args.threads,prefix))
 	os.system("{} simplify -COS -d 0.8 {}_tiddit/clips.fa.assembly.mag 1> {}_tiddit/clips.fa.assembly.clean.mag 2> /dev/null".format(args.fermi2,prefix,prefix))
-	os.system("{} mem -x intractg {} {}_tiddit/clips.fa.assembly.clean.mag  1> {}_tiddit/clips.sam 2> /dev/null".format(args.bwa,args.ref,prefix,prefix))
+	os.system("{} mem -t {} -x intractg {} {}_tiddit/clips.fa.assembly.clean.mag  1> {}_tiddit/clips.sam 2> /dev/null".format(args.bwa,args.threads,args.ref,prefix,prefix))
 
 	read_contigs("{}_tiddit/clips.sam".format(prefix) , prefix, sample_id, args.z)
 

tiddit/tiddit_coverage.pyx

@@ -1,22 +1,24 @@
 import sys
+import time
 cimport numpy
 import numpy
 import math
 cimport cython
 @cython.boundscheck(False)
 @cython.wraparound(False)
 
-def create_coverage(bam_header,bin_size):
+def create_coverage(bam_header,bin_size,c="all"):
 	coverage_data={}
 	end_bin_size={}
 
 	for contig in bam_header["SQ"]:
-		bins= int(math.ceil(contig["LN"]/float(bin_size)))
-		coverage_data[ contig["SN"] ]=numpy.zeros(bins)
-		end_bin_size[contig["SN"]]=contig["LN"]-(bins-1)*bin_size
-
+		if c == "all" or contig["SN"] == c:
+			bins= int(math.ceil(contig["LN"]/float(bin_size)))
+			coverage_data[ contig["SN"] ]=numpy.zeros(bins)
+			end_bin_size[contig["SN"]]=contig["LN"]-(bins-1)*bin_size
+			if c != "all":
+				return(coverage_data[ contig["SN"] ],end_bin_size[contig["SN"]])
 	return(coverage_data,end_bin_size)
-
 def print_coverage(coverage_data,bam_header,bin_size,file_type,outfile):
 	f=open(outfile,"w",buffering=819200)
 
@@ -43,37 +45,34 @@ def print_coverage(coverage_data,bam_header,bin_size,file_type,outfile):
 	f.close()
 
 ctypedef numpy.double_t DTYPE_t
-def update_coverage(read,int bin_size,numpy.ndarray[DTYPE_t, ndim=1] coverage_data,int min_q,int end_bin_size):
-
-	cdef long ref_start=read.reference_start
-	cdef long ref_end=read.reference_end
+def update_coverage(long ref_start,long ref_end,int bin_size,numpy.ndarray[DTYPE_t, ndim=1] coverage_data,int end_bin_size):
 
 	cdef int first_bin=ref_start//bin_size
-	cdef int end_bin=int(ref_end-1)//bin_size
+	cdef int end_bin=(ref_end-1)//bin_size
 
-	cdef int bases_first_bin
+	cdef float bases_first_bin
 
 	if end_bin == first_bin:
 		bases_first_bin=ref_end-ref_start
-		coverage_data[first_bin]=float(bases_first_bin)/bin_size+coverage_data[first_bin]
+		coverage_data[first_bin]=bases_first_bin/bin_size+coverage_data[first_bin]
 
 		return(coverage_data)
 
 	bases_first_bin=((first_bin+1)*bin_size)-ref_start
-	coverage_data[first_bin]=float(bases_first_bin)/bin_size+coverage_data[first_bin]	
-	cdef int bases_last_bin=(ref_end-1)-end_bin*bin_size
+	coverage_data[first_bin]=bases_first_bin/bin_size+coverage_data[first_bin]
+	cdef float bases_last_bin=(ref_end-1)-end_bin*bin_size
+
 
 	if end_bin < len(coverage_data)-1:
-		coverage_data[end_bin]+=float(bases_last_bin)/bin_size
+		coverage_data[end_bin]=bases_last_bin/bin_size+coverage_data[end_bin]
 	else:
-		coverage_data[end_bin]+=float(bases_last_bin)/end_bin_size
+		coverage_data[end_bin]=bases_last_bin/end_bin_size+coverage_data[end_bin]
 
 	for i in range(first_bin+1,end_bin):
-		coverage_data[i]+=1.0
+		coverage_data[i]=1.0+coverage_data[i]
 
 	return(coverage_data)
 
-
 #bam_file_name=sys.argv[1]
 
 #samfile = pysam.AlignmentFile(bam_file_name, "r")

tiddit/tiddit_coverage_analysis.pyx

@@ -75,6 +75,9 @@ def determine_ploidy(dict coverage_data,contigs,dict library,int ploidy,str pref
 		library["avg_coverage"]=c
 
 	for chromosome in contigs:
+		if not chromosome in coverage_data:
+			continue
+
 		avg_coverage_contig=library[ "avg_coverage_{}".format(chromosome) ]
 		library["contig_ploidy_{}".format(chromosome)]=int(round(ploidy*avg_coverage_contig/library["avg_coverage"]))
 		f.write("{}\t{}\t{}\t{}\n".format(chromosome,avg_coverage_contig/library["avg_coverage"]*ploidy,library["contig_ploidy_{}".format(chromosome)],avg_coverage_contig))