Mercurial > repos > tyty > structurefold
view Iterative_mapping/iterative_map.py @ 53:52538ebd3c4b draft
Uploaded
author | tyty |
---|---|
date | Mon, 17 Nov 2014 15:14:21 -0500 |
parents | d56631911cc1 |
children |
line wrap: on
line source
#!/usr/bin/env python # -*- coding: utf-8 -*- import sys import os from read_file import * from read_s_file import * type_input = sys.argv[1] seq_file = sys.argv[2] ref_file = sys.argv[3] shift = sys.argv[4] length = sys.argv[5] map_type = sys.argv[6] output_file = sys.argv[7] if map_type!="default": s = "" s = s+"-v "+sys.argv[8] s = s+" -5 "+sys.argv[9] s = s+" -3 "+sys.argv[10] s = s+" -k "+sys.argv[11] if sys.argv[12]: s = s+" -a" if int(sys.argv[13])>=1: s = s+" -m "+sys.argv[13] if sys.argv[14]: s = s+" --best --strata " else: s = "-v 3 -a --best --strata " ospath = os.path.realpath(sys.argv[0]) ost = ospath.split('/') syspath = "" for i in range(len(ost)-1): syspath = syspath+ost[i].strip() syspath = syspath+'/' os.system("bowtie-build -f "+ref_file+" "+syspath+"ref > "+syspath+"log.txt") os.system("cp "+seq_file+" "+syspath+"seq0.fa") if type_input == "fasta": tp = 'fasta' if type_input == "fastq": tp = 'fastq' k = 0 print(type_input) while(True): if type_input == "fasta": os.system("bowtie "+s+"-f "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam") if type_input == "fastq": os.system("bowtie "+s+"-q "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam") os.system("samtools view -Sb -F 0xfff "+syspath+"map"+str(k)+".sam > "+syspath+"mapped"+str(k)+".bam 2>"+syspath+"log.txt") #get mapped reads os.system("samtools view -Sb -f 0x4 "+syspath+"map"+str(k)+".sam > "+syspath+"umapped"+str(k)+".bam 2>"+syspath+"log.txt") #get unmapped reads os.system("samtools view -Sb -f 0x10 "+syspath+"map"+str(k)+".sam > "+syspath+"rmapped"+str(k)+".bam 2>"+syspath+"log.txt") #get reversed mapped reads os.system("samtools merge -f "+syspath+"unmapped"+str(k)+".bam "+syspath+"umapped"+str(k)+".bam "+syspath+"rmapped"+str(k)+".bam") #get reversed mapped reads os.system("samtools view -h -o "+syspath+"unmapped"+str(k)+".sam "+syspath+"unmapped"+str(k)+".bam") #get reversed mapped reads if k>0: os.system("samtools view -h -o "+syspath+"mapped"+str(k)+".sam "+syspath+"mapped"+str(k)+".bam") #get reversed mapped reads os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"unmapped"+str(k)+".txt") os.system("cut -f 1 "+syspath+"mapped"+str(k)+".sam > "+syspath+"mapped"+str(k)+".txt") os.system("python "+syspath+"remove_map.py "+syspath+"unmapped"+str(k)+".txt "+syspath+"mapped"+str(k)+".txt "+syspath+"runmapped"+str(k)+".txt") os.system("rm "+syspath+"mapped"+str(k)+".sam") os.system("rm "+syspath+"mapped"+str(k)+".txt") os.system("rm "+syspath+"unmapped"+str(k)+".txt") else: os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"runmapped"+str(k)+".txt") os.system("rm "+syspath+"unmapped"+str(k)+".bam") os.system("rm "+syspath+"umapped"+str(k)+".bam") os.system("rm "+syspath+"rmapped"+str(k)+".bam") os.system("python "+syspath+"seq_track.py "+syspath+"runmapped"+str(k)+".txt "+syspath+"seq"+str(k)+".fa "+syspath+"unmap_seq"+str(k)+".fa "+tp) #get unmapped sequence os.system("python "+syspath+"truncate.py "+syspath+"unmap_seq"+str(k)+".fa "+shift+" "+syspath+"seq"+str(k+1)+".fa "+length) #truncate unmapped sequence os.system("rm "+syspath+"seq"+str(k)+".fa") #Remove sequences being mapped os.system("rm "+syspath+"map"+str(k)+".sam") #Remove mapping file os.system("rm "+syspath+"unmap_seq"+str(k)+".fa") #Remove unmapped sequnce os.system("rm "+syspath+"runmapped"+str(k)+".txt") os.system("rm "+syspath+"unmapped"+str(k)+".sam") os.system("wc -l "+syspath+"seq"+str(k+1)+".fa > "+syspath+"count"+str(k+1)+".txt") c = read_sp_file(syspath+"count"+str(k+1)+".txt") if c[0][0] == '0': #If no reads is in the sequence file, stop os.system("rm "+syspath+"count"+str(k+1)+".txt") os.system("rm "+syspath+"seq"+str(k+1)+".fa") break os.system("rm "+syspath+"count"+str(k+1)+".txt") k = k+1 ss = "" for i in range(0,k+1): ss = ss+" "+syspath+"mapped"+str(i)+".bam" os.system("samtools merge -f "+output_file+" "+ss) #print("samtools merge mapped_all.bam"+ss) os.system("rm "+syspath+"mapped*.bam") os.system("rm "+syspath+"ref*")