# HG changeset patch # User tyty # Date 1413783191 14400 # Node ID 74ae6f6a3383c4fcc52f962c687d3a0a91478bea # Parent b10c6790f01e9f891319a2c757b4fc32e1ed89c3 Uploaded diff -r b10c6790f01e -r 74ae6f6a3383 Iterative_mapping/iterative_map.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Iterative_mapping/iterative_map.py Mon Oct 20 01:33:11 2014 -0400 @@ -0,0 +1,103 @@ +#!/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*") + +