update the pip version

.gitignore

@@ -33,3 +33,7 @@
 
 # Custom rules (everything added below won't be overriden by 'Generate .gitignore File' if you use 'Update' option)
 ./tests
+setup.sh
+./build
+./docs
+./dist

build/lib/searcHPV/alignment.py

@@ -13,7 +13,8 @@
 #outputDir: output directory
 #multi: if fastq file is in gz format: default = True
 #index: if True, index the reference files; if False, not index the references files
-def alignment(fq1, fq2, humRef, virRef, outputDir, index, gz):
+#memory: memory size allocated
+def alignment(fq1, fq2, humRef, virRef, outputDir, index, gz, memory,thread):
     #make output dir
     outputDir = os.path.abspath(outputDir)
     mkdir(outputDir)
@@ -26,8 +27,8 @@ def alignment(fq1, fq2, humRef, virRef, outputDir, index, gz):
 
     alignmentFile = scriptDir + "/orignal.alignment.sh"
     indelFile = scriptDir + "/indel.alignment.sh"
-    generate_alignment_bash(alignmentFile,ref,fq1,fq2,scriptDir,gz)
-    generate_indel_alignment_bash(indelFile,ref,scriptDir)
+    generate_alignment_bash(alignmentFile,ref,fq1,fq2,scriptDir,memory,thread,gz)
+    generate_indel_alignment_bash(indelFile,ref,scriptDir, memory,thread)
     check_file(alignmentFile)
     check_file(indelFile)
     bashFile = scriptDir + f"/alignment.sh"
@@ -51,7 +52,7 @@ def alignment(fq1, fq2, humRef, virRef, outputDir, index, gz):
     if index:
     ##generate mkdup alignment bash file
         mkdupFile = scriptDir + "/mkdup.alignment.sh"
-        generate_mkdup_bash(mkdupFile,scriptDir)
+        generate_mkdup_bash(mkdupFile,scriptDir,thread)
 
         with open(bashFile,'w') as output:
             output.write(f'''#!/bin/bash
@@ -65,7 +66,7 @@ def alignment(fq1, fq2, humRef, virRef, outputDir, index, gz):
     else:
         ##generate mkdup alignment bash file
         mkdupFile = scriptDir + "/mkdup.alignment.sh"
-        generate_mkdup_bash(mkdupFile,scriptDir)
+        generate_mkdup_bash(mkdupFile,scriptDir,thread)
 
         with open(bashFile,'w') as output:
             output.write(f'''#!/bin/bash

build/lib/searcHPV/assemble.py

@@ -9,7 +9,7 @@
 #out_dir: output directory for searcHPV
 #virRef: virus reference genome
 #window: the length of region searching for informative reads, default=300
-def assemble(fq1, fq2, out_dir,virRef,gz,window):
+def assemble(fq1, fq2, out_dir,virRef,gz,window,memory,thread):
     bam = f'{out_dir}/alignment/alignment.RG.indelre.mkdup.sort.bam'
     check_file(bam)
     assemble_out_dir = f'{out_dir}/assemble/'
@@ -20,7 +20,7 @@ def assemble(fq1, fq2, out_dir,virRef,gz,window):
     subprocess.call(script_read_seq)
 
     #preprocessForPear(assemble_out_dir)
-    script_pear = pear(assemble_out_dir)
+    script_pear = pear(assemble_out_dir,memory,thread)
     os.system(f'chmod +x {script_pear}')
     subprocess.call(script_pear)
 

build/lib/searcHPV/generate_alignment.py

@@ -42,21 +42,21 @@ def indexRef(bash_file,ref):
 #fq2: fastq2 file
 #out_dir: outputPath
 #gz: if fastq file is in gz format: default = True
-def generate_alignment_bash(bash_file,ref,fq1,fq2,out_dir,gz = True):
+def generate_alignment_bash(bash_file,ref,fq1,fq2,out_dir,memory,thread,gz = True):
     with open(bash_file,'w') as output:
         if gz:
             output.write(f'''
-    bwa mem {ref} '<zcat {fq1}' '<zcat {fq2}' > {out_dir}/alignment.sam
-    samtools view -bhS {out_dir}/alignment.sam > {out_dir}/alignment.bam
-    samtools sort {out_dir}/alignment.bam -o {out_dir}/alignment.sort.bam
+    bwa mem -t {thread} {ref} '<zcat {fq1}' '<zcat {fq2}' > {out_dir}/alignment.sam
+    samtools view -@ {thread} -bhS {out_dir}/alignment.sam > {out_dir}/alignment.bam
+    samtools sort -@ {thread} {out_dir}/alignment.bam -o {out_dir}/alignment.sort.bam
     rm {out_dir}/alignment.sam
     echo \'alignment done\'
     ''')
         else:
             output.write(f'''
-    bwa mem {ref} '<cat {fq1}' '<cat {fq2}' > {out_dir}/alignment.sam
-    samtools view -bhS {out_dir}/alignment.sam > {out_dir}/alignment.bam
-    samtools sort {out_dir}/alignment.bam -o {out_dir}/alignment.sort.bam
+    bwa mem -t {thread} {ref} '<cat {fq1}' '<cat {fq2}' > {out_dir}/alignment.sam
+    samtools view -@ {thread} -bhS {out_dir}/alignment.sam > {out_dir}/alignment.bam
+    samtools sort -@ {thread} {out_dir}/alignment.bam -o {out_dir}/alignment.sort.bam
     rm {out_dir}/alignment.sam
     echo \'alignment done\'
     ''')
@@ -78,7 +78,7 @@ def generate_alignment_bash(bash_file,ref,fq1,fq2,out_dir,gz = True):
 #picard: full path of picard
 #gatk: full path of GATK
 #java: full path of java
-def generate_indel_alignment_bash(bash_File,ref,out_dir):
+def generate_indel_alignment_bash(bash_File,ref,out_dir,memory,thread):
     with open(bash_File,'w') as output:
         output.write(f'''picard \
 AddOrReplaceReadGroups \
@@ -89,35 +89,37 @@ def generate_indel_alignment_bash(bash_File,ref,out_dir):
 RGPU=NA \
 RGSM=sample \
 RGLB=sample
-samtools sort {out_dir}/alignment.RG.bam -o {out_dir}/alignment.RG.sort.bam
-samtools index {out_dir}/alignment.RG.sort.bam
+samtools sort -@ {thread} {out_dir}/alignment.RG.bam -o {out_dir}/alignment.RG.sort.bam
+samtools index -@ {thread} {out_dir}/alignment.RG.sort.bam
 GenomeAnalysisTK \
+-Xmx{memory} \
 -T RealignerTargetCreator \
 -R {ref} \
 -I {out_dir}/alignment.RG.sort.bam \
 -o {out_dir}/alignment.RG.intervals
 GenomeAnalysisTK \
+-Xmx{memory} \
 -T IndelRealigner \
 -R {ref} \
 -I {out_dir}/alignment.RG.sort.bam \
 -targetIntervals {out_dir}/alignment.RG.intervals \
 -o {out_dir}/alignment.RG.indelre.bam
-samtools index {out_dir}/alignment.RG.indelre.bam
+samtools index -@ {thread} {out_dir}/alignment.RG.indelre.bam
 echo \'indel alignment done\'''')
 
 
 
-def generate_mkdup_bash(bash_File,out_dir):
+def generate_mkdup_bash(bash_File,out_dir,thread):
     with open(bash_File,'w') as output:
         output.write(f'''
-samtools sort -n {out_dir}/alignment.RG.indelre.bam -o {out_dir}/alignment.RG.indelre.sortbyQ.bam
+samtools sort -@ {thread} -n {out_dir}/alignment.RG.indelre.bam -o {out_dir}/alignment.RG.indelre.sortbyQ.bam
 picard MarkDuplicates \
 I={out_dir}/alignment.RG.indelre.sortbyQ.bam \
 O={out_dir}/alignment.RG.indelre.mkdup.bam \
 M={out_dir}/alignment.RG.indelre.mkdup.txt \
 TAGGING_POLICY=All ASSUME_SORT_ORDER=queryname
-samtools sort {out_dir}/alignment.RG.indelre.mkdup.bam -o {out_dir}/alignment.RG.indelre.mkdup.sort.bam
-samtools index {out_dir}/alignment.RG.indelre.mkdup.sort.bam
+samtools sort -@ {thread} {out_dir}/alignment.RG.indelre.mkdup.bam -o {out_dir}/alignment.RG.indelre.mkdup.sort.bam
+samtools index -@ {thread} {out_dir}/alignment.RG.indelre.mkdup.sort.bam
 echo \'indel alignment done\'''')
 
 

build/lib/searcHPV/generate_assemble.py

@@ -3,6 +3,7 @@
 import re
 import pysam
 from searcHPV.general import *
+import string
 
 ########################
 #get informative reads (SP + PE) from bam file
@@ -67,7 +68,9 @@ def extract_read_name(bam,out_dir,virRef,window):
             if '>' in each:
                 virus_chrm_list.append(each.split()[0].replace('>',''))
     for virus_chrm in virus_chrm_list:
-        fusionRes = f'{out_dir}/call_fusion/{virus_chrm}.all.filtered.clustered.result'
+        invalidCharacter = re.escape(string.punctuation)
+        virus_chrm_file_name = re.sub(r'['+invalidCharacter+']',"_",virus_chrm)
+        fusionRes = f'{out_dir}/call_fusion/{virus_chrm_file_name}.all.filtered.clustered.result'
         check_file(fusionRes)
         with open(fusionRes) as candidate_in:
             candidate_in = candidate_in.read().rstrip()
@@ -172,7 +175,7 @@ def preprocessForPear(out_dir):
 #Function:generate bash script for PEAR
 #out_dir:output directory for assemble
 #return:path of bash script
-def pear(out_dir):
+def pear(out_dir,memory,thread):
     listSites = os.listdir(out_dir)
     with open(f'{out_dir}/pear.sh','w') as output:
         output.write('#!/bin/bash\n')
@@ -184,6 +187,8 @@ def pear(out_dir):
                     os.mkdir(outputPath)
                 output.write(f'''
     pear \
+    -j {thread} \
+    -y {memory} \
     -f {fqPath}/{site}.informativeReads.1.fq \
     -r {fqPath}/{site}.informativeReads.2.fq \
     -o {outputPath}/{site}''')

build/lib/searcHPV/generate_call_fusion.py

@@ -2,6 +2,8 @@
 import sys
 from searcHPV.general import *
 import re
+import string
+import re
 
 ##################
 #find the end of each cigar string slot
@@ -153,12 +155,13 @@ def define_fusion(bam,virus_chrm,out_dir):
                 if read.next_reference_name == virus_chrm:
                     paired_evidence_count += 1
         b[i] += (str(paired_evidence_count),)
-
-    with open(f'{out_dir}/{virus_chrm}.genome_fusion.txt','w') as output:
+    invalidCharacter = re.escape(string.punctuation)
+    virus_chrm_file_name = re.sub(r'['+invalidCharacter+']',"_",virus_chrm)
+    with open(f'{out_dir}/{virus_chrm_file_name}.genome_fusion.txt','w') as output:
         output.write('chrm\tpos\tsingle_evidence\tpaired_evidence\n')
         for i in b:
             output.write(f'{i[0]}\t{i[1]}\t{i[2]}\t{i[3]}\n')
-    return f'{out_dir}/{virus_chrm}.genome_fusion.txt'
+    return f'{out_dir}/{virus_chrm_file_name}.genome_fusion.txt'
 
 ##################
 #cluster fusion points within certain base pair

build/lib/searcHPV/generate_call_fusion_virus.py

@@ -8,7 +8,7 @@
 #out_dir:out_dir for searcHPV
 #humRef: human reference genome
 #return: path of bash script 
-def mapToRef(out_dir):
+def mapToRef(out_dir,thread):
     newRef = f'{out_dir}/hg_hpv.fa'
     sitesPath = f'{out_dir}/assemble/'
     listSites = os.listdir(sitesPath)
@@ -23,11 +23,11 @@ def mapToRef(out_dir):
                     contigPath = f'{sitesPath}/{site}/pearOutput/'            
                     sitePath = f'{out_dir}/call_fusion_virus/{site}/'
                     mkdir(sitePath)
-                    bashFile.write( f'''bwa mem -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
+                    bashFile.write( f'''bwa mem -t {thread} -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
 {newRef} {contigPath}/{site}.all.fa.cap.contigs > {sitePath}/{site}.contig.sam;
-samtools view -bhS {sitePath}/{site}.contig.sam > {sitePath}/{site}.contig.bam;
-samtools sort {sitePath}/{site}.contig.bam -o {sitePath}/{site}.contig.sort.bam;
-samtools index {sitePath}/{site}.contig.sort.bam;
+samtools view -@ {thread} -bhS {sitePath}/{site}.contig.sam > {sitePath}/{site}.contig.bam;
+samtools sort -@ {thread} {sitePath}/{site}.contig.bam -o {sitePath}/{site}.contig.sort.bam;
+samtools index -@ {thread} {sitePath}/{site}.contig.sort.bam;
 rm {sitePath}/{site}.contig.bam;\n''')
     return f'{outputPath}/alignContigsToRef.sh'
 
@@ -39,7 +39,7 @@ def mapToRef(out_dir):
 #out_dir:out_dir for searcHPV
 #humRef: human reference genome
 #return: path of bash script 
-def mapToHgRef(out_dir,humRef):
+def mapToHgRef(out_dir,humRef,thread):
     sitesPath = f'{out_dir}/assemble/'
     listSites = os.listdir(sitesPath)
     outputPath = f'{out_dir}/call_fusion_virus/'
@@ -53,11 +53,11 @@ def mapToHgRef(out_dir,humRef):
                     contigPath = f'{sitesPath}/{site}/pearOutput/'            
                     sitePath = f'{out_dir}/call_fusion_virus/{site}/'
                     mkdir(sitePath)
-                    bashFile.write( f'''bwa mem -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
+                    bashFile.write( f'''bwa mem -t {thread} -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
 {humRef} {contigPath}/{site}.all.fa.cap.contigs > {sitePath}/{site}.contigToGenome.sam;
-samtools view -bhS {sitePath}/{site}.contigToGenome.sam > {sitePath}/{site}.contigToGenome.bam;
-samtools sort {sitePath}/{site}.contigToGenome.bam -o {sitePath}/{site}.contigToGenome.sort.bam;
-samtools index {sitePath}/{site}.contigToGenome.sort.bam;
+samtools view -@ {thread} -bhS {sitePath}/{site}.contigToGenome.sam > {sitePath}/{site}.contigToGenome.bam;
+samtools sort -@ {thread} {sitePath}/{site}.contigToGenome.bam -o {sitePath}/{site}.contigToGenome.sort.bam;
+samtools index -@ {thread} {sitePath}/{site}.contigToGenome.sort.bam;
 rm {sitePath}/{site}.contigToGenome.bam;\n''')
     return f'{outputPath}/alignContigsToGenome.sh'
 
@@ -66,7 +66,7 @@ def mapToHgRef(out_dir,humRef):
 #out_dir:out_dir for searcHPV
 #virRef: virus reference genome
 #return: path of bash script 
-def mapToVirRef(out_dir,virRef):
+def mapToVirRef(out_dir,virRef,thread):
     sitesPath = f'{out_dir}/assemble/'
     listSites = os.listdir(sitesPath)
     outputPath = f'{out_dir}/call_fusion_virus/'
@@ -79,12 +79,12 @@ def mapToVirRef(out_dir,virRef):
                 if ".sh" not in site:
                     contigPath = f'{sitesPath}/{site}/pearOutput/'
                     sitePath = f'{out_dir}/call_fusion_virus/{site}/'
-                    bashFile.write( f'''bwa mem -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
+                    bashFile.write( f'''bwa mem -t {thread} -R \'@RG\\tID:hpv\\tSM:hpv\\tLB:hpv\\tPL:ILLUMINA\' -M -t 8 \
 {virRef} {contigPath}/{site}.all.fa.cap.contigs > {contigPath}/{site}.contigToVirus.sam;
-samtools view -bhS {contigPath}/{site}.contigToVirus.sam >  {contigPath}/{site}.contigToVirus.bam;
-samtools sort {contigPath}/{site}.contigToVirus.bam -o {sitePath}/{site}.contigToVirus.sort.bam;
-samtools index {sitePath}/{site}.contigToVirus.sort.bam;
-samtools faidx {contigPath}/{site}.all.fa.cap.contigs;
+samtools view -@ {thread} -bhS {contigPath}/{site}.contigToVirus.sam >  {contigPath}/{site}.contigToVirus.bam;
+samtools sort -@ {thread} {contigPath}/{site}.contigToVirus.bam -o {sitePath}/{site}.contigToVirus.sort.bam;
+samtools index -@ {thread} {sitePath}/{site}.contigToVirus.sort.bam;
+samtools faidx -@ {thread} {contigPath}/{site}.all.fa.cap.contigs;
 rm {contigPath}/{site}.contigToVirus.bam;\n''')
     return f'{outputPath}/alignContigsToVirus.sh'
 

build/lib/searcHPV/genome_fusion.py

@@ -2,6 +2,7 @@
 import sys
 from searcHPV.general import *
 from searcHPV.generate_call_fusion import *
+import string
 
 def genomeFusion(window,out_dir,virRef):
 
@@ -24,17 +25,20 @@ def genomeFusion(window,out_dir,virRef):
 
         #filter and cluster fusion points
         ##sort result
-        os.system(f'(head -n 1 {res} && tail -n +2 {res} | sort -k3,3rn) > {out_dir}/{virus_chrm}.genome_fusion.sort.txt')
+        #if there are invalid characters in the chromosome name:
+        invalidCharacter = re.escape(string.punctuation)
+        virus_chrm_file_name = re.sub(r'['+invalidCharacter+']',"_",virus_chrm)
+        os.system(f'(head -n 1 {res} && tail -n +2 {res} | sort -k3,3rn) > {out_dir}/{virus_chrm_file_name}.genome_fusion.sort.txt')
 
 
         #change format for cluster
         chrm_li = range(1,22)
         chrm_li = list(map(str,chrm_li))
         chrm_li+=['X','Y']
 
-        with open(f'{out_dir}/{virus_chrm}.all.result','w') as output:   
+        with open(f'{out_dir}/{virus_chrm_file_name}.all.result','w') as output:   
                 fusion_li = []
-                with open(f'{out_dir}/{virus_chrm}.genome_fusion.sort.txt') as res:
+                with open(f'{out_dir}/{virus_chrm_file_name}.genome_fusion.sort.txt') as res:
                     res.readline()
                     for line in res.read().rstrip().split('\n'):
                         elements = line.rstrip().split('\t')
@@ -51,12 +55,12 @@ def genomeFusion(window,out_dir,virRef):
                 output.write(f'{to_print}')
 
         ##cluster the events within 100bp from each other, maybe becasue of SVs or CNVs
-        cluster_result(f'{out_dir}/{virus_chrm}.all.result',f'{out_dir}/{virus_chrm}.all.clustered.result',window)
+        cluster_result(f'{out_dir}/{virus_chrm_file_name}.all.result',f'{out_dir}/{virus_chrm_file_name}.all.clustered.result',window)
 
 
         ##filter for sites with at least 2 split read  and 2 pairs of read support(high cutoff) and their summation greater than 5
-        with open(f'{out_dir}/{virus_chrm}.all.clustered.result') as inf:
-            with open(f'{out_dir}/{virus_chrm}.all.filtered.clustered.result','w') as outf:
+        with open(f'{out_dir}/{virus_chrm_file_name}.all.clustered.result') as inf:
+            with open(f'{out_dir}/{virus_chrm_file_name}.all.filtered.clustered.result','w') as outf:
                 for line in inf.read().rstrip().split('\n'):
                     elements = line.rstrip()
                     if elements == "":

build/lib/searcHPV/main.py

@@ -61,6 +61,13 @@ def main():
                     default = False, dest ='index', 
                     help ="index the original human and virus reference files, default=False")
 
+    parser.add_argument('-memory', type = str, default="10G",dest = 'memory',
+                        help ='specify the maximum size, default "10G"', 
+                        )
+    parser.add_argument('-thread', type = int, default=1,dest = 'thread',
+                        help ='number of threads, default 1', 
+                        )
+
 
 
 
@@ -72,19 +79,19 @@ def main():
 
 
     if args.alignment:
-        alignment(fq1 = args.fq1, fq2 = args.fq2, humRef = args.humRef, virRef = args.virRef, index = args.index, outputDir = args.outputDir, gz = args.gz)
+        alignment(fq1 = args.fq1, fq2 = args.fq2, humRef = args.humRef, virRef = args.virRef, index = args.index, outputDir = args.outputDir, gz = args.gz, memory = args.memory,thread = args.thread)
     elif args.genomeFusion: 
         genomeFusion(args.clusterWindow,args.outputDir,args.virRef)
     elif args.assemble: 
         #check result from genomeFusion
-        assemble(args.fq1, args.fq2, args.outputDir,args.virRef,args.gz,args.window)
+        assemble(args.fq1, args.fq2, args.outputDir,args.virRef,args.gz,args.window, memory = args.memory,thread = args.thread)
     elif args.hpvFusion:
         #check result from assemble
-        virus_fusion(args.humRef,args.virRef,args.outputDir,args.n)
+        virus_fusion(args.humRef,args.virRef,args.outputDir,args.n,thread = args.thread)
     else:
-        alignment(fq1 = args.fq1, fq2 = args.fq2, humRef = args.humRef, virRef = args.virRef, index = args.index, outputDir = args.outputDir, gz = args.gz)
+        alignment(fq1 = args.fq1, fq2 = args.fq2, humRef = args.humRef, virRef = args.virRef, index = args.index, outputDir = args.outputDir, gz = args.gz, memory = args.memory,thread = args.thread)
         genomeFusion(args.clusterWindow,args.outputDir,args.virRef)
-        assemble(args.fq1, args.fq2, args.outputDir,args.virRef,args.gz,args.window)
-        virus_fusion(args.humRef,args.virRef,args.outputDir,args.n)
+        assemble(args.fq1, args.fq2, args.outputDir,args.virRef,args.gz,args.window, memory = args.memory,thread = args.thread)
+        virus_fusion(args.humRef,args.virRef,args.outputDir,args.n,thread = args.thread)
 
 

build/lib/searcHPV/virus_fusion.py

@@ -9,10 +9,10 @@
 #virRef: virus reference genome
 #out_dir: output directory for seacHPV
 #n: poly(n), n*d(A/T/C/G), will report low confidence if contig contains poly(n)
-def virus_fusion(humRef,virRef,out_dir,n):
-    script_map = mapToRef(out_dir)
-    script_mapHg = mapToHgRef(out_dir,humRef)
-    script_mapVir = mapToVirRef(out_dir,virRef)
+def virus_fusion(humRef,virRef,out_dir,n,thread):
+    script_map = mapToRef(out_dir,thread)
+    script_mapHg = mapToHgRef(out_dir,humRef,thread)
+    script_mapVir = mapToVirRef(out_dir,virRef,thread)
     os.system(f'chmod +x {script_map}')
     subprocess.call(script_map)
     os.system(f'chmod +x {script_mapHg}')