opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Adding tools to compute Substitution rates.

Browse Source
This commit is contained in:
Guruprasad Anada
2008-09-21 17:36:28 -04:00
parent cb74e79f21
commit e2f4f46a75
5 changed files with 306 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tool_conf.xml.sample
+2
View File
@@ -128,6 +128,8 @@
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/getIndelRates_3way.xml" />
<tool file="regVariation/substitutions.xml" />
<tool file="regVariation/substitution_rates.xml" />
</section>
<section name="Multiple regression" id="multReg">
<tool file="regVariation/linear_regression.xml" />
tools/regVariation/substitution_rates.py
+118
View File
@@ -0,0 +1,118 @@
#! /usr/bin/python
#guruprasad Ananda
"""
Estimates substitution rates from pairwise alignments using JC69 model.
"""
from galaxy import eggs
from galaxy.tools.util.galaxyops import *
from galaxy.tools.util import maf_utilities
import bx.align.maf
import sys, fileinput
def stop_err(msg):
sys.stderr.write(msg)
sys.exit()
if len(sys.argv) < 3:
stop_err("Incorrect number of arguments.")
inp_file = sys.argv[1]
out_file = sys.argv[2]
fout = open(out_file, 'w')
int_file = sys.argv[3]
if int_file != "None": #The user has specified an interval file
dbkey_i = sys.argv[4]
chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] )
def rateEstimator(block):
global alignlen, mismatches
src1 = block.components[0].src
sequence1 = block.components[0].text
start1 = block.components[0].start
end1 = block.components[0].end
len1 = int(end1)-int(start1)
len1_withgap = len(sequence1)
mismatch = 0.0
for seq in range (1,len(block.components)):
src2 = block.components[seq].src
sequence2 = block.components[seq].text
start2 = block.components[seq].start
end2 = block.components[seq].end
len2 = int(end2)-int(start2)
for nt in range(len1_withgap):
if sequence1[nt] not in '-#$^*?' and sequence2[nt] not in '-#$^*?': #Not a gap or masked character
if sequence1[nt].upper() != sequence2[nt].upper():
mismatch += 1
if int_file == "None":
p = mismatch/min(len1,len2)
print >>fout, "%s\t%s\t%s\t%s\t%s\t%s\t%d\t%d\t%.4f" %(src1,start1,end1,src2,start2,end2,min(len1,len2),mismatch,p)
else:
mismatches += mismatch
alignlen += min(len1,len2)
def main():
skipped = 0
not_pairwise = 0
if int_file == "None":
try:
maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
except:
stop_err("Your MAF file appears to be malformed.")
print >>fout, "#Seq1\tStart1\tEnd1\tSeq2\tStart2\tEnd2\tL\tN\tp"
for block in maf_reader:
if len(block.components) != 2:
not_pairwise += 1
continue
try:
rateEstimator(block)
except:
skipped += 1
else:
index, index_filename = maf_utilities.build_maf_index( inp_file, species = [dbkey_i] )
if index is None:
print >> sys.stderr, "Your MAF file appears to be malformed."
sys.exit()
win = NiceReaderWrapper( fileinput.FileInput( int_file ),
chrom_col=chr_col_i,
start_col=start_col_i,
end_col=end_col_i,
strand_col=strand_col_i,
fix_strand=True)
species=None
mincols = 0
global alignlen, mismatches
for interval in win:
alignlen = 0
mismatches = 0.0
src = "%s.%s" % ( dbkey_i, interval.chrom )
for block in maf_utilities.get_chopped_blocks_for_region( index, src, interval, species, mincols ):
if len(block.components) != 2:
not_pairwise += 1
continue
try:
rateEstimator(block)
except:
skipped += 1
if alignlen:
p = mismatches/alignlen
else:
p = 'NA'
interval.fields.append(str(alignlen))
interval.fields.append(str(mismatches))
interval.fields.append(str(p))
print >>fout, "\t".join(interval.fields)
#num_blocks += 1
if not_pairwise:
print "Skipped %d non-pairwise blocks" %(not_pairwise)
if skipped:
print "Skipped %d blocks as invalid" %(skipped)
if __name__ == "__main__":
main()
tools/regVariation/substitution_rates.xml
+61
View File
@@ -0,0 +1,61 @@
<tool id="subRate1" name="Estimate substitution rates " version="1.0.0">
<description> for non-coding regions</description>
<command interpreter="python">
substitution_rates.py
$input
$out_file1
#if $region.type == "win":
${region.input2} ${region.input2.dbkey} ${region.input2.metadata.chromCol},$region.input2.metadata.startCol,$region.input2.metadata.endCol,$region.input2.metadata.strandCol
#else:
"None"
#end if
</command>
<inputs>
<param format="maf" name="input" type="data" label="Select pair-wise alignment data"/>
<conditional name="region">
<param name="type" type="select" label="Estimate rates corresponding to" multiple="false">
<option value="align">Alignment block</option>
<option value="win">Intervals in your history</option>
</param>
<when value="win">
<param format="interval" name="input2" type="data" label="Choose intervals">
<validator type="unspecified_build" />
</param>
</when>
<when value="align" />
</conditional>
</inputs>
<outputs>
<data format="tabular" name="out_file1" metadata_source="input"/>
</outputs>
<tests>
<test>
<param name="input" value="Interval2Maf_pairwise_out.maf"/>
<param name="type" value="align"/>
<output name="out_file1" file="subRates1.out"/>
</test>
</tests>
<help>
.. class:: infomark
**What it does**
This tool takes a pairwise MAF file as input and estimates substitution rate according to Jukes-Cantor JC69 model. The 3 new columns appended to the output are explanied below:
- L: number of nucleotides compared
- N: number of different nucleotides
- p = N/L
-----
.. class:: warningmark
**Note**
Any block/s not containing exactly two sequences, will be omitted.
</help>
</tool>
tools/regVariation/substitutions.py
+87
View File
@@ -0,0 +1,87 @@
#! /usr/bin/python
#Guruprasad ANanda
"""
Fetches substitutions from pairwise alignments.
"""
from galaxy import eggs
from galaxy.tools.util import maf_utilities
import bx.align.maf
import sys
import os, fileinput
def stop_err(msg):
sys.stderr.write(msg)
sys.exit()
if len(sys.argv) < 3:
stop_err("Incorrect number of arguments.")
inp_file = sys.argv[1]
out_file = sys.argv[2]
fout = open(out_file, 'w')
def fetchSubs(block):
src1 = block.components[0].src
sequence1 = block.components[0].text
start1 = block.components[0].start
end1 = block.components[0].end
len1 = int(end1)-int(start1)
len1_withgap = len(sequence1)
for seq in range (1,len(block.components)):
src2 = block.components[seq].src
sequence2 = block.components[seq].text
start2 = block.components[seq].start
end2 = block.components[seq].end
len2 = int(end2)-int(start2)
sub_begin = None
sub_end = None
begin = False
for nt in range(len1_withgap):
if sequence1[nt] not in '-#$^*?' and sequence2[nt] not in '-#$^*?': #Not a gap or masked character
if sequence1[nt].upper() != sequence2[nt].upper():
if not(begin):
sub_begin = nt
begin = True
sub_end = nt
else:
if begin:
print >>fout, "%s\t%s\t%s" %(src1,start1+sub_begin-sequence1[0:sub_begin].count('-'),start1+sub_end-sequence1[0:sub_end].count('-'))
print >>fout, "%s\t%s\t%s" %(src2,start2+sub_begin-sequence2[0:sub_begin].count('-'),start2+sub_end-sequence2[0:sub_end].count('-'))
begin = False
else:
if begin:
print >>fout, "%s\t%s\t%s" %(src1,start1+sub_begin-sequence1[0:sub_begin].count('-'),end1+sub_end-sequence1[0:sub_end].count('-'))
print >>fout, "%s\t%s\t%s" %(src2,start2+sub_begin-sequence2[0:sub_begin].count('-'),end2+sub_end-sequence2[0:sub_end].count('-'))
begin = False
ended = False
def main():
skipped = 0
not_pairwise = 0
try:
maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
except:
stop_err("Your MAF file appears to be malformed.")
print >>fout, "#Chr\tStart\tEnd"
for block in maf_reader:
if len(block.components) != 2:
not_pairwise += 1
continue
try:
fetchSubs(block)
except:
skipped += 1
if not_pairwise:
print "Skipped %d non-pairwise blocks" %(not_pairwise)
if skipped:
print "Skipped %d blocks" %(skipped)
if __name__ == "__main__":
main()
tools/regVariation/substitutions.xml
+38
View File
@@ -0,0 +1,38 @@
<tool id="substitutions1" name="Fetch substitutions " version="1.0.0">
<description> from pairwise alignments</description>
<command interpreter="python">
substitutions.py
$input
$out_file1
</command>
<inputs>
<param format="maf" name="input" type="data" label="Select pair-wise alignment data"/>
</inputs>
<outputs>
<data format="tabular" name="out_file1" metadata_source="input"/>
</outputs>
<tests>
<test>
<param name="input" value="Interval2Maf_pairwise_out.maf"/>
<output name="out_file1" file="subs.out"/>
</test>
</tests>
<help>
.. class:: infomark
**What it does**
This tool takes a pairwise MAF file as input and fetches substitutions per alignment block.
-----
.. class:: warningmark
**Note**
Any block/s not containing exactly two sequences, will be omitted.
</help>
</tool>