Mercurial > repos > tyty > structurefold
comparison Iterative_mapping/iterative_map.py @ 56:9d26c2e4953e draft
Uploaded
| author | tyty |
|---|---|
| date | Tue, 18 Nov 2014 01:00:33 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 55:20b8f6ab0e05 | 56:9d26c2e4953e |
|---|---|
| 1 #!/usr/bin/env python | |
| 2 # -*- coding: utf-8 -*- | |
| 3 | |
| 4 import sys | |
| 5 import os | |
| 6 from read_file import * | |
| 7 from read_s_file import * | |
| 8 import random | |
| 9 import string | |
| 10 | |
| 11 type_input = sys.argv[1] | |
| 12 seq_file = sys.argv[2] | |
| 13 ref_file = sys.argv[3] | |
| 14 shift = sys.argv[4] | |
| 15 length = sys.argv[5] | |
| 16 t_end = sys.argv[6] | |
| 17 map_type = sys.argv[7] | |
| 18 output_file = sys.argv[8] | |
| 19 | |
| 20 | |
| 21 if map_type!="default": | |
| 22 s = "" | |
| 23 s = s+"-v "+sys.argv[9] | |
| 24 s = s+" -5 "+sys.argv[10] | |
| 25 s = s+" -3 "+sys.argv[11] | |
| 26 s = s+" -k "+sys.argv[12] | |
| 27 if sys.argv[13]: | |
| 28 s = s+" -a" | |
| 29 if int(sys.argv[14])>=1: | |
| 30 s = s+" -m "+sys.argv[14] | |
| 31 if sys.argv[15]: | |
| 32 s = s+" --best --strata " | |
| 33 | |
| 34 else: | |
| 35 s = "-v 3 -a --best --strata " | |
| 36 | |
| 37 ospath = os.path.realpath(sys.argv[0]) | |
| 38 ost = ospath.split('/') | |
| 39 syspath = "" | |
| 40 for i in range(len(ost)-1): | |
| 41 syspath = syspath+ost[i].strip() | |
| 42 syspath = syspath+'/' | |
| 43 | |
| 44 rs = ''.join(random.sample(string.ascii_letters + string.digits, 8)) | |
| 45 os.system("mkdir "+syspath+"output_"+rs) | |
| 46 | |
| 47 syspathrs = syspath+"output_"+rs+'/' | |
| 48 | |
| 49 os.system("bowtie-build -f "+ref_file+" "+syspathrs+"ref > "+syspathrs+"log.txt") | |
| 50 | |
| 51 os.system("cp "+seq_file+" "+syspathrs+"seq0.fa") | |
| 52 | |
| 53 if type_input == "fasta": | |
| 54 tp = 'fasta' | |
| 55 if type_input == "fastq": | |
| 56 tp = 'fastq' | |
| 57 | |
| 58 k = 0 | |
| 59 print(type_input) | |
| 60 while(True): | |
| 61 if type_input == "fasta": | |
| 62 os.system("bowtie "+s+"-f "+syspathrs+"ref"+" "+syspathrs+"seq"+str(k)+".fa --quiet -S > "+syspathrs+"map"+str(k)+".sam") | |
| 63 if type_input == "fastq": | |
| 64 os.system("bowtie "+s+"-q "+syspathrs+"ref"+" "+syspathrs+"seq"+str(k)+".fa --quiet -S > "+syspathrs+"map"+str(k)+".sam") | |
| 65 os.system("samtools view -Sb -F 0xfff "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"mapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get mapped reads | |
| 66 os.system("samtools view -Sb -f 0x4 "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"umapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get unmapped reads | |
| 67 os.system("samtools view -Sb -f 0x10 "+syspathrs+"map"+str(k)+".sam > "+syspathrs+"rmapped"+str(k)+".bam 2>"+syspathrs+"log.txt") #get reversed mapped reads | |
| 68 os.system("samtools merge -f "+syspathrs+"unmapped"+str(k)+".bam "+syspathrs+"umapped"+str(k)+".bam "+syspathrs+"rmapped"+str(k)+".bam") #get reversed mapped reads | |
| 69 os.system("samtools view -h -o "+syspathrs+"unmapped"+str(k)+".sam "+syspathrs+"unmapped"+str(k)+".bam") #get reversed mapped reads | |
| 70 if k>0: | |
| 71 os.system("samtools view -h -o "+syspathrs+"mapped"+str(k)+".sam "+syspathrs+"mapped"+str(k)+".bam") #get reversed mapped reads | |
| 72 os.system("cut -f 1 "+syspathrs+"unmapped"+str(k)+".sam > "+syspathrs+"unmapped"+str(k)+".txt") | |
| 73 os.system("cut -f 1 "+syspathrs+"mapped"+str(k)+".sam > "+syspathrs+"mapped"+str(k)+".txt") | |
| 74 os.system("python "+syspath+"remove_map.py "+syspathrs+"unmapped"+str(k)+".txt "+syspathrs+"mapped"+str(k)+".txt "+syspathrs+"runmapped"+str(k)+".txt") | |
| 75 os.system("rm "+syspathrs+"mapped"+str(k)+".sam") | |
| 76 os.system("rm "+syspathrs+"mapped"+str(k)+".txt") | |
| 77 os.system("rm "+syspathrs+"unmapped"+str(k)+".txt") | |
| 78 else: | |
| 79 os.system("cut -f 1 "+syspathrs+"unmapped"+str(k)+".sam > "+syspathrs+"runmapped"+str(k)+".txt") | |
| 80 | |
| 81 os.system("rm "+syspathrs+"unmapped"+str(k)+".bam") | |
| 82 os.system("rm "+syspathrs+"umapped"+str(k)+".bam") | |
| 83 os.system("rm "+syspathrs+"rmapped"+str(k)+".bam") | |
| 84 os.system("python "+syspath+"seq_track.py "+syspathrs+"runmapped"+str(k)+".txt "+syspathrs+"seq"+str(k)+".fa "+syspathrs+"unmap_seq"+str(k)+".fa "+tp) #get unmapped sequence | |
| 85 os.system("python "+syspath+"truncate.py "+syspathrs+"unmap_seq"+str(k)+".fa "+shift+" "+syspathrs+"seq"+str(k+1)+".fa "+length+" "+t_end) #truncate unmapped sequence | |
| 86 os.system("rm "+syspathrs+"seq"+str(k)+".fa") #Remove sequences being mapped | |
| 87 os.system("rm "+syspathrs+"map"+str(k)+".sam") #Remove mapping file | |
| 88 os.system("rm "+syspathrs+"unmap_seq"+str(k)+".fa") #Remove unmapped sequnce | |
| 89 os.system("rm "+syspathrs+"runmapped"+str(k)+".txt") | |
| 90 os.system("rm "+syspathrs+"unmapped"+str(k)+".sam") | |
| 91 | |
| 92 os.system("wc -l "+syspathrs+"seq"+str(k+1)+".fa > "+syspathrs+"count"+str(k+1)+".txt") | |
| 93 c = read_sp_file(syspathrs+"count"+str(k+1)+".txt") | |
| 94 if c[0][0] == '0': #If no reads is in the sequence file, stop | |
| 95 os.system("rm "+syspathrs+"count"+str(k+1)+".txt") | |
| 96 os.system("rm "+syspathrs+"seq"+str(k+1)+".fa") | |
| 97 break | |
| 98 os.system("rm "+syspathrs+"count"+str(k+1)+".txt") | |
| 99 k = k+1 | |
| 100 | |
| 101 ss = "" | |
| 102 for i in range(0,k+1): | |
| 103 ss = ss+" "+syspathrs+"mapped"+str(i)+".bam" | |
| 104 | |
| 105 | |
| 106 os.system("samtools merge -f "+output_file+" "+ss) | |
| 107 #print("samtools merge mapped_all.bam"+ss) | |
| 108 os.system("rm "+syspathrs+"mapped*.bam") | |
| 109 os.system("rm "+syspathrs+"ref*") | |
| 110 os.system("rm -r "+syspathrs) | |
| 111 | |
| 112 |