Added tool version directories to regVariation tools directory.

This commit is contained in:
Greg Von Kuster
2008-02-22 20:09:25 +00:00
parent 218684f610
commit b75509fa17
26 changed files with 14 additions and 882 deletions
+6 -6
View File
@@ -117,13 +117,13 @@
<tool file="visualization/build_ucsc_custom_track.xml" />
</section>
<section name="Regional Variation" id="regVar">
<tool file="regVariation/windowSplitter.xml" />
<tool file="regVariation/featureCounter.xml" />
<tool file="regVariation/quality_filter.xml" />
<tool file="regVariation/maf_cpg_filter.xml" />
<tool file="regVariation/winSplitter/1.0.0/windowSplitter.xml" />
<tool file="regVariation/featureCoverage1/1.0.0/featureCounter.xml" />
<tool file="regVariation/qualityFilter/1.0.0/quality_filter.xml" />
<tool file="regVariation/cpgFilter/1.0.0/maf_cpg_filter.xml" />
<!--
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/getIndels_2way/1.0.0/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way/1.0.0/getIndels_3way.xml" />
-->
</section>
<section name="Evolution: HyPhy" id="hyphy">
+6 -6
View File
@@ -122,12 +122,12 @@
<tool file="visualization/build_ucsc_custom_track.xml" />
</section>
<section name="Regional Variation" id="regVar">
<tool file="regVariation/windowSplitter.xml" />
<tool file="regVariation/featureCounter.xml" />
<tool file="regVariation/quality_filter.xml" />
<tool file="regVariation/maf_cpg_filter.xml" />
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/winSplitter/1.0.0/windowSplitter.xml" />
<tool file="regVariation/featureCoverage1/1.0.0/featureCounter.xml" />
<tool file="regVariation/qualityFilter/1.0.0/quality_filter.xml" />
<tool file="regVariation/cpgFilter/1.0.0/maf_cpg_filter.xml" />
<tool file="regVariation/getIndels_2way/1.0.0/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way/1.0.0/getIndels_3way.xml" />
</section>
<section name="Evolution: HyPhy" id="hyphy">
<tool file="hyphy/hyphy_branch_lengths_wrapper1/1.0.0/hyphy_branch_lengths_wrapper.xml" />
-53
View File
@@ -1,53 +0,0 @@
#Adapted from bx/intervals/operations/coverage.py
"""
Determine amount of each interval in one set covered by the intervals of
another set. Adds two columns to the first input, giving number of bases
covered and percent coverage on the second input.
"""
import pkg_resources
pkg_resources.require( "bx-python" )
import psyco_full
import traceback
import fileinput
from warnings import warn
from bx.intervals.io import *
from bx.intervals.operations import *
def coverage(readers, comments=True):
# Read all but first into bitsets and union to one
primary = readers[0]
intersect = readers[1:]
bitsets = intersect[0].binned_bitsets()
intersect = intersect[1:]
for andset in intersect:
bitset2 = andset.binned_bitsets()
for chrom in bitsets:
if chrom not in bitset2: continue
bitsets[chrom].ior(bitset2[chrom])
intersect = intersect[1:]
# Read remaining intervals and give coverage
for interval in primary:
if type( interval ) is Header:
yield interval
if type( interval ) is Comment and comments:
yield interval
elif type( interval ) == GenomicInterval:
chrom = interval.chrom
start = int(interval.start)
end = int(interval.end)
if start > end: warn( "Interval start after end!" )
if chrom not in bitsets:
bases_covered = 0
percent = 0.0
else:
bases_covered = bitsets[ chrom ].count_range( start, end-start )
if (end - start) == 0: percent = 0
else: percent = float(bases_covered) / float(end - start)
interval.fields.append(str(bases_covered))
interval.fields.append(str(percent))
yield interval
-136
View File
@@ -1,136 +0,0 @@
#!/usr/bin/env python
"""
Calculate coverage of one query on another, and append the coverage to
the last two columns as bases covered and percent coverage.
usage: %prog bed_file_1 bed_file_2 out_file
-1, --cols1=N,N,N,N: Columns for start, end, strand in first file
-2, --cols2=N,N,N,N: Columns for start, end, strand in second file
"""
import pkg_resources
pkg_resources.require( "bx-python" )
import sys
import traceback
import fileinput
from warnings import warn
#from bx.intervals import *
from bx.intervals.io import *
#from bx.intervals.operations.coverage import *
from bx.cookbook import doc_optparse
from featureCounter_code import *
from galaxyops import *
def main():
upstream_pad = 0
downstream_pad = 0
options, args = doc_optparse.parse( __doc__ )
try:
chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
in_fname, in2_fname, out_fname = args
except:
doc_optparse.exception()
g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
chrom_col=chr_col_1,
start_col=start_col_1,
end_col=end_col_1,
strand_col=strand_col_1,
fix_strand=True)
g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
chrom_col=chr_col_2,
start_col=start_col_2,
end_col=end_col_2,
strand_col=strand_col_2,
fix_strand=True)
out_file = open( out_fname, "w" )
chr_list = []
start_list = []
end_list = []
for line in open(in2_fname, 'r'):
line = line.rstrip("\r\n")
if not(line) or line == "" or line[0:1] == '#':
continue
elems = line.split('\t')
chr_list.append(elems[chr_col_2])
start_list.append(elems[start_col_2])
end_list.append(elems[end_col_2])
for line1 in open(in_fname, 'r'):
fp = 0 #number of features fully present
fp_p = [] #% covered by fp features
pp_l = 0 #partly present on left side of the window
pp_l_p = [] #% covered by pp_l features
pp_r = 0 #partly present on right side of the window
pp_r_p = [] #% covered by pp_r features
np = 0 #number of features present
chr1 = line1.split('\t')[chr_col_1]
start1 = int(line1.split('\t')[start_col_1])
end1 = int(line1.split('\t')[end_col_1])
if not(chr1 in chr_list):
continue
for i, chr2 in enumerate(chr_list):
#chr2 = line2.split('\t')[chr_col_2]
#start2 = int(line2.split('\t')[start_col_2])
#end2 = int(line2.split('\t')[end_col_2])
start2 = int(start_list[i])
end2 = int(end_list[i])
if chr2 == chr1:
"""
if end2 < start2:
tmp = start2
start2 = end2
end2 = tmp
"""
if (start2 < start1 and end2 < start1) or (start2 >end1 and end2 >end1):
continue
if start2 in range(start1,end1):
if end2 in range(start1,end1):
fp+=1
fp_p.append("%2.2f"%(100.0*abs(end2 - start2)/abs(end1-start1)))
else:
pp_r = pp_r + (1.0*abs(end1 - start2)/abs(end2-start2))
pp_r_p.append("%2.2f"%(100.0*abs(end1 - start2)/abs(end1-start1)))
#pp_r_p = pp_r_p + abs(end1 - start2)
elif end2 in range(start1,end1):
pp_l = pp_l + (1.0*abs(end2 - start1)/abs(end2 - start2))
pp_l_p.append("%2.2f"%(100.0*abs(end2 - start1)/abs(end1-start1)))
#pp_l_p = pp_l_p + abs(end2 - start1)
else:
if (start1 in range(start2,end2)) and (end1 in range(start2,end2)):
fp = fp + (1.0*abs(end1 - start1)/abs(end2 - start2))
fp_p.append('100.00')
print >>out_file, "%s\t%2.2f\t%s\t%2.2f\t%s\t%2.2f\t%s" %(line1.strip(),fp,fp_p,pp_l,pp_l_p,pp_r,pp_r_p)
#print >>out_file, "%s\t%d\t%2.2f\t%d\t%2.2f\t%d\t%2.2f" %(line1.strip(),fp,float(fp_p)/abs(end1-start1),pp_l,float(pp_l_p)/abs(end1-start1),pp_r,float(pp_r_p)/abs(end1-start1))
"""
fo=open("pfct","w")
print >>fo, g1
print >>fo, g2
try:
for line in mycoverage([g1,g2]):
if type( line ) is GenomicInterval:
print >> out_file, "\t".join( line.fields )
else:
print >> out_file, line
except ParseError, exc:
print >> sys.stderr, "Invalid file format: ", str( exc )
if g1.skipped > 0:
print skipped( g1, filedesc=" of 1st dataset" )
if g2.skipped > 0:
print skipped( g2, filedesc=" of 2nd dataset" )
"""
if __name__ == "__main__":
main()
-127
View File
@@ -1,127 +0,0 @@
#!/usr/bin/env python
"""
Calculate coverage of one query on another, and append the coverage to
the last two columns as bases covered and percent coverage.
usage: %prog bed_file_1 bed_file_2 out_file
-1, --cols1=N,N,N,N: Columns for start, end, strand in first file
-2, --cols2=N,N,N,N: Columns for start, end, strand in second file
"""
import pkg_resources
pkg_resources.require( "bx-python" )
import sys
import traceback
import fileinput
from warnings import warn
#from bx.intervals import *
from bx.intervals.io import *
#from bx.intervals.operations.coverage import *
from bx.cookbook import doc_optparse
from featureCounter_code import *
from galaxyops import *
def main():
upstream_pad = 0
downstream_pad = 0
options, args = doc_optparse.parse( __doc__ )
try:
chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
in_fname, in2_fname, out_fname = args
except:
doc_optparse.exception()
g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
chrom_col=chr_col_1,
start_col=start_col_1,
end_col=end_col_1,
strand_col=strand_col_1,
fix_strand=True)
g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
chrom_col=chr_col_2,
start_col=start_col_2,
end_col=end_col_2,
strand_col=strand_col_2,
fix_strand=True)
out_file = open( out_fname, "w" )
chr_list = []
start_list = []
end_list = []
for line in open(in2_fname, 'r'):
line = line.rstrip("\r\n")
if not(line) or line == "" or line[0:1] == '#':
continue
elems = line.split('\t')
chr_list.append(elems[chr_col_2])
start_list.append(elems[start_col_2])
end_list.append(elems[end_col_2])
for line1 in open(in_fname, 'r'):
fp = 0 #number of features fully present
fp_p = [] #% covered by fp features
pp_l = 0 #partly present on left side of the window
pp_l_p = [] #% covered by pp_l features
pp_r = 0 #partly present on right side of the window
pp_r_p = [] #% covered by pp_r features
np = 0 #number of features present
chr1 = line1.split('\t')[chr_col_1]
start1 = int(line1.split('\t')[start_col_1])
end1 = int(line1.split('\t')[end_col_1])
for line2 in open(in2_fname, 'r'):
chr2 = line2.split('\t')[chr_col_2]
start2 = int(line2.split('\t')[start_col_2])
end2 = int(line2.split('\t')[end_col_2])
if chr2 == chr1:
if end2 < start2:
tmp = start2
start2 = end2
end2 = tmp
if start2 in range(start1,end1):
if end2 in range(start1,end1):
fp+=1
fp_p.append(round(100.0*abs(end2 - start2)/abs(end1-start1)))
else:
pp_r+=1
pp_r_p.append(round(100.0*abs(end1 - start2)/abs(end1-start1)))
#pp_r_p = pp_r_p + abs(end1 - start2)
elif end2 in range(start1,end1):
pp_l+=1
pp_l_p.append(round(100.0*abs(end2 - start1)/abs(end1-start1)))
#pp_l_p = pp_l_p + abs(end2 - start1)
else:
if (start1 in range(start2,end2)) and (end1 in range(start2,end2)):
fp+=1
fp_p.append(100.0)
print >>out_file, "%s\t%d\t%s\t%d\t%s\t%d\t%s" %(line1.strip(),fp,fp_p,pp_l,pp_l_p,pp_r,pp_r_p)
#print >>out_file, "%s\t%d\t%2.2f\t%d\t%2.2f\t%d\t%2.2f" %(line1.strip(),fp,float(fp_p)/abs(end1-start1),pp_l,float(pp_l_p)/abs(end1-start1),pp_r,float(pp_r_p)/abs(end1-start1))
"""
fo=open("pfct","w")
print >>fo, g1
print >>fo, g2
try:
for line in mycoverage([g1,g2]):
if type( line ) is GenomicInterval:
print >> out_file, "\t".join( line.fields )
else:
print >> out_file, line
except ParseError, exc:
print >> sys.stderr, "Invalid file format: ", str( exc )
if g1.skipped > 0:
print skipped( g1, filedesc=" of 1st dataset" )
if g2.skipped > 0:
print skipped( g2, filedesc=" of 2nd dataset" )
"""
if __name__ == "__main__":
main()
-134
View File
@@ -1,134 +0,0 @@
#!/usr/bin/env python
"""
Calculate coverage of one query on another, and append the coverage to
the last two columns as bases covered and percent coverage.
usage: %prog bed_file_1 bed_file_2 out_file
-1, --cols1=N,N,N,N: Columns for start, end, strand in first file
-2, --cols2=N,N,N,N: Columns for start, end, strand in second file
"""
import pkg_resources
pkg_resources.require( "bx-python" )
import sys
import traceback
import fileinput
from warnings import warn
#from bx.intervals import *
from bx.intervals.io import *
#from bx.intervals.operations.coverage import *
from bx.cookbook import doc_optparse
from featureCounter_code import *
from galaxyops import *
def main():
upstream_pad = 0
downstream_pad = 0
options, args = doc_optparse.parse( __doc__ )
try:
chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
in_fname, in2_fname, out_fname = args
except:
doc_optparse.exception()
g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
chrom_col=chr_col_1,
start_col=start_col_1,
end_col=end_col_1,
strand_col=strand_col_1,
fix_strand=True)
g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
chrom_col=chr_col_2,
start_col=start_col_2,
end_col=end_col_2,
strand_col=strand_col_2,
fix_strand=True)
out_file = open( out_fname, "w" )
chr_list = []
start_list = []
end_list = []
for line in open(in2_fname, 'r'):
line = line.rstrip("\r\n")
if not(line) or line == "" or line[0:1] == '#':
continue
elems = line.split('\t')
chr_list.append(elems[chr_col_2])
start_list.append(elems[start_col_2])
end_list.append(elems[end_col_2])
for line1 in open(in_fname, 'r'):
fp = 0 #number of features fully present
fp_p = [] #% covered by fp features
pp_l = 0 #partly present on left side of the window
pp_l_p = [] #% covered by pp_l features
pp_r = 0 #partly present on right side of the window
pp_r_p = [] #% covered by pp_r features
np = 0 #number of features present
chr1 = line1.split('\t')[chr_col_1]
start1 = int(line1.split('\t')[start_col_1])
end1 = int(line1.split('\t')[end_col_1])
if not(chr1 in chr_list):
continue
for i, chr2 in enumerate(chr_list):
#chr2 = line2.split('\t')[chr_col_2]
#start2 = int(line2.split('\t')[start_col_2])
#end2 = int(line2.split('\t')[end_col_2])
start2 = int(start_list[i])
end2 = int(end_list[i])
if chr2 == chr1:
"""
if end2 < start2:
tmp = start2
start2 = end2
end2 = tmp
"""
if start2 in range(start1,end1):
if end2 in range(start1,end1):
fp+=1
fp_p.append("%2.2f"%(100.0*abs(end2 - start2)/abs(end1-start1)))
else:
pp_r+=1
pp_r_p.append("%2.2f"%(100.0*abs(end1 - start2)/abs(end1-start1)))
#pp_r_p = pp_r_p + abs(end1 - start2)
elif end2 in range(start1,end1):
pp_l+=1
pp_l_p.append("%2.2f"%(100.0*abs(end2 - start1)/abs(end1-start1)))
#pp_l_p = pp_l_p + abs(end2 - start1)
else:
if (start1 in range(start2,end2)) and (end1 in range(start2,end2)):
fp+=1
fp_p.append(100.00)
print >>out_file, "%s\t%d\t%s\t%d\t%s\t%d\t%s" %(line1.strip(),fp,fp_p,pp_l,pp_l_p,pp_r,pp_r_p)
#print >>out_file, "%s\t%d\t%2.2f\t%d\t%2.2f\t%d\t%2.2f" %(line1.strip(),fp,float(fp_p)/abs(end1-start1),pp_l,float(pp_l_p)/abs(end1-start1),pp_r,float(pp_r_p)/abs(end1-start1))
"""
fo=open("pfct","w")
print >>fo, g1
print >>fo, g2
try:
for line in mycoverage([g1,g2]):
if type( line ) is GenomicInterval:
print >> out_file, "\t".join( line.fields )
else:
print >> out_file, line
except ParseError, exc:
print >> sys.stderr, "Invalid file format: ", str( exc )
if g1.skipped > 0:
print skipped( g1, filedesc=" of 1st dataset" )
if g2.skipped > 0:
print skipped( g2, filedesc=" of 2nd dataset" )
"""
if __name__ == "__main__":
main()
-73
View File
@@ -1,73 +0,0 @@
"""
Determine amount of each interval in one set covered by the intervals of
another set. Adds two columns to the first input, giving number of bases
covered and percent coverage on the second input.
"""
import pkg_resources
pkg_resources.require( "bx-python" )
import psyco_full
import traceback
import fileinput
from warnings import warn
from bx.intervals.io import *
from bx.intervals.operations import *
fo = open("pfc","w")
def mycoverage(readers, comments=True):
# Read all but first into bitsets and union to one
primary = readers[0]
intersect = readers[1:]
bitsets = intersect[0].binned_bitsets()
#print >>fo, primary
#print >>fo, intersect
#print >>fo, bitsets
#print >>fo, primary.keys()
#print >>fo, primary.values()
#print >>fo, intersect.keys()
#print >>fo, intersect.values()
print >>fo, bitsets.keys()
print >>fo, bitsets.values()
intersect = intersect[1:]
i=j=0
for andset in intersect:
print >>fo, "inside j for"
print >>fo, andset
j = j+1
bitset2 = andset.binned_bitsets()
for chrom in bitsets:
i+=1
if chrom not in bitset2: continue
bitsets[chrom].ior(bitset2[chrom])
print >>fo, "inside i for"
print >>fo, i
intersect = intersect[1:]
print >>fo, j
total_features = 0
# Read remaining intervals and give coverage
for interval in primary:
if type( interval ) is Header:
yield interval
if type( interval ) is Comment and comments:
yield interval
elif type( interval ) == GenomicInterval:
chrom = interval.chrom
start = int(interval.start)
end = int(interval.end)
if start > end: warn( "Interval start after end!" )
if chrom not in bitsets:
bases_covered = 0
percent = 0.0
else:
total_features += 1
bases_covered = bitsets[ chrom ].count_range( start, end-start )
if (end - start) == 0: percent = 0
else: percent = float(bases_covered) / float(end - start)
interval.fields.append(str(total_features))
interval.fields.append(str(percent))
yield interval
@@ -22,7 +22,7 @@ from bx.intervals.io import *
from bx.intervals.operations.merge import *
from bx.cookbook import doc_optparse
from galaxyops import *
from galaxy.tools.util.galaxyops import *
def stop_err(msg):
sys.stderr.write(msg)
-34
View File
@@ -1,34 +0,0 @@
"""
Utility functions for galaxyops
"""
from bx.bitset import *
from bx.intervals.io import *
import sys
def warn( msg ):
print >> sys.stderr, msg
def fail( msg ):
print >> sys.stderr, msg
sys.exit( 1 )
# Default chrom, start, end, stran cols for a bed file
BED_DEFAULT_COLS = 0, 1, 2, 5
def parse_cols_arg( cols ):
"""Parse a columns command line argument into a four-tuple"""
if cols:
return map( lambda x: int( x ) - 1, cols.split(",") )
else:
return BED_DEFAULT_COLS
def default_printer( stream, exc, obj ):
print >> stream, "%d: %s" % ( obj.linenum, obj.current_line )
print >> stream, "\tError: %s" % ( str(exc) )
def skipped( reader, filedesc="" ):
first_line, line_contents = reader.skipped_lines[0]
return 'Data issue: skipped %d invalid lines%s starting at line #%d which is "%s"' \
% ( reader.skipped, filedesc, first_line, line_contents )
-126
View File
@@ -1,126 +0,0 @@
#!/usr/bin/python2.4
"""
Estimate INDELs for pait-wise alignments.
usage: %prog maf_input out_file1 out_file2
"""
from __future__ import division
import pkg_resources
pkg_resources.require( "bx-python" )
pkg_resources.require( "lrucache" )
try:
pkg_resources.require("numpy")
pkg_resources.require( "python-lzo" )
except:
pass
import psyco_full
import sys
import os, os.path
from UserDict import DictMixin
import bx.wiggle
from bx.binned_array import BinnedArray, FileBinnedArray
from bx.bitset import *
from bx.bitset_builders import *
from fpconst import isNaN
from bx.cookbook import doc_optparse
from galaxy.tools.exception_handling import *
import bx.align.maf
def main():
# Parsing Command Line here
options, args = doc_optparse.parse( __doc__ )
try:
inp_file, out_file1 = args
except:
print >> sys.stderr, "Tool initialization error."
sys.exit()
try:
fin = open(inp_file,'r')
except:
print >> sys.stderr, "Unable to open input file"
sys.exit()
try:
fout1 = open(out_file1,'w')
#fout2 = open(out_file2,'w')
except:
print >> sys.stderr, "Unable to open output file"
sys.exit()
try:
maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
except:
print >>sys.stderr, "Your MAF file appears to be malformed."
sys.exit()
maf_count = 0
print >>fout1, "#Block\tSource\tStart\tEnd\tIndel_length"
for block_ind, block in enumerate(maf_reader):
if len(block.components) < 2:
continue
seq1 = block.components[0].text
src1 = block.components[0].src
start1 = block.components[0].start
if len(block.components) == 2:
seq2 = block.components[1].text
src2 = block.components[1].src
start2 = block.components[1].start
#for pos in range(len(seq1)):
nt_pos1 = start1-1 #position of the nucleotide (without counting gaps)
nt_pos2 = start2-1
pos = 0 #character column position
gaplen1 = 0
gaplen2 = 0
prev_pos_gap1 = 0
prev_pos_gap2 = 0
while pos < len(seq1):
if prev_pos_gap1 == 0:
gaplen1 = 0
if prev_pos_gap2 == 0:
gaplen2 = 0
if seq1[pos] == '-':
if seq2[pos] != '-':
nt_pos2 += 1
gaplen1 += 1
prev_pos_gap1 = 1
#write 2
if prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
if pos == len(seq1)-1:
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1,nt_pos1+1,gaplen1)
else:
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
elif prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
else:
nt_pos1 += 1
if seq2[pos] != '-':
nt_pos2 += 1
#write both
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
elif prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
else:
gaplen2 += 1
prev_pos_gap2 = 1
#write 1
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
if pos == len(seq1)-1:
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2,nt_pos2+1,gaplen2)
pos += 1
if __name__ == "__main__":
main()
-131
View File
@@ -1,131 +0,0 @@
#!/usr/bin/python2.4
"""
Estimate INDEL rates.
usage: %prog maf_input out_file1 out_file2
"""
from __future__ import division
import pkg_resources
pkg_resources.require( "bx-python" )
pkg_resources.require( "lrucache" )
try:
pkg_resources.require("numpy")
pkg_resources.require( "python-lzo" )
except:
pass
import psyco_full
import sys
import os, os.path
from UserDict import DictMixin
import bx.wiggle
from bx.binned_array import BinnedArray, FileBinnedArray
from bx.bitset import *
from bx.bitset_builders import *
from fpconst import isNaN
from bx.cookbook import doc_optparse
from galaxy.tools.exception_handling import *
import bx.align.maf
def main():
# Parsing Command Line here
options, args = doc_optparse.parse( __doc__ )
try:
#chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols )
inp_file, out_file1 = args
except:
doc_optparse.exception()
try:
fin = open(inp_file,'r')
except:
print >> sys.stderr, "Unable to open input file"
sys.exit()
try:
fout1 = open(out_file1,'w')
#fout2 = open(out_file2,'w')
except:
print >> sys.stderr, "Unable to open output file"
sys.exit()
try:
maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
except:
print >>sys.stderr, "Your MAF file appears to be malformed."
sys.exit()
maf_count = 0
print >>fout1, "#Block\tSource\tStart\tEnd\tEvent"
for block_ind, block in enumerate(maf_reader):
if len(block.components) != 2:
continue
seqs = []
srcs = []
starts = []
gaplens = []
gapstatus = []
starts.append(block.components[0].start)
for seq_num in range(len(block.components)):
seqs.append(block.components[seq_num].text)
srcs.append(block.components[seq_num].src)
starts.append(block.components[seq_num].start)
gaplens.append(0)
gapstatus.append(0)
pos = 0 #character column position
nt_pos1 = 0 #nt positions
nt_pos2 = 0
while pos < len(seqs[0]):
for j,elem in enumerate(seqs):
if gapstatus[j] == 0:
gaplens[j] = 0
if seqs[0][pos] == '-':
next = pos+1
leng = 1
while next < len(seqs[0]):
if seqs[0][next] == '-':
leng += 1
if seq2[pos] != '-':
nt_pos2 += 1
gaplen1 += 1
prev_pos_gap1 = 1
#write 2
if prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
if pos == len(seq1)-1:
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1,nt_pos1+1,gaplen1)
else:
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
elif prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
else:
nt_pos1 += 1
if seq2[pos] != '-':
nt_pos2 += 1
#write both
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
elif prev_pos_gap2 == 1:
prev_pos_gap2 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2-1,nt_pos2,gaplen2)
else:
gaplen2 += 1
prev_pos_gap2 = 1
#write 1
if prev_pos_gap1 == 1:
prev_pos_gap1 = 0
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src1,nt_pos1-1,nt_pos1,gaplen1)
if pos == len(seq1)-1:
print >>fout1,"%d\t%s\t%s\t%s\t%s" %(block_ind+1,src2,nt_pos2,nt_pos2+1,gaplen2)
pos += 1
if __name__ == "__main__":
main()
-24
View File
@@ -1,24 +0,0 @@
#!/usr/bin/env python
import sys, re, os, tempfile
qual_dir = "/depot/data2/galaxy/mm8/align/multiz17way"
fout = open("ploc","w")
os.chdir(qual_dir)
tmpfile = tempfile.NamedTemporaryFile()
cmdline = "ls " + "*.lzo | cat >> " + tmpfile.name
os.system (cmdline)
fstr = "17-way multiZ (mm8)\t17_WAY_MULTIZ_mm8\tmm8\t"
for j,qual_file in enumerate(tmpfile.readlines()):
if j!=0:
fstr = fstr + ',' + qual_dir + '/' + qual_file.strip()
else:
fstr = fstr + qual_dir + '/' + qual_file.strip()
print >>fout, fstr
print fstr
os.system("echo '%s' | cat >> /depot/data2/galaxy/maf_index.loc" %(fstr))
-30
View File
@@ -1,30 +0,0 @@
#!/usr/bin/env python
import sys, re, os, tempfile
qual_dir = "/home/gua110/Desktop/rhesus_quality_scores/chr"
qual_file= open("/home/gua110/Desktop/rhesus_quality_scores/rheMac2.qual.qa", "r")
qual_file_contents = qual_file.read(1000000000)
while qual_file_contents != "":
contents_list = qual_file_contents.split(">")
os.chdir(qual_dir)
print "len_contents_list", len(contents_list)
if len(contents_list) == 1:
os.system("echo %s | cat >> %s.qa" %(contents_list, prev_file))
else:
#print len(contents_list[0])
#print len(contents_list[1])
for elem in contents_list:
if elem == '':
continue
if elem.startswith("chr"):
elems = elem.replace("\r","\n").split("\n")
cmdline = "echo %s | cat >> %s.qa" %(elems[1:], elems[0])
os.system("echo >%s\n%s | cat >> %s.qa" %(elems[0], elems[1:], elems[0]))
else:
os.system("echo %s | cat >> %s.qa" %(elem, prev_file))
prev_file = elems[0]
qual_file_contents = qual_file.read(1000000000)
@@ -11,7 +11,7 @@ import sys, sets, re, os
import pkg_resources; pkg_resources.require( "bx-python" )
from bx.cookbook import doc_optparse
from galaxyops import *
from galaxy.tools.util.galaxyops import *
def main():
# Parsing Command Line here