Extracting Genomic DNA can now result in an interval file with the

sequence appended as a new column.
This commit is contained in:
Daniel Blankenberg
2007-09-25 16:12:37 +00:00
parent 1c3c0d2659
commit a3b137d94b
3 changed files with 22 additions and 5 deletions
+10 -4
View File
@@ -18,7 +18,7 @@ my %seen = ();
my @locFields = ();
my %seqLocation = ();
die "Cannot fetch sequences for unspecified genome\n" unless @ARGV == 11;
die "Cannot fetch sequences for unspecified genome\n" unless @ARGV == 13;
die "Please specify genome build by clicking on pencil icon in the original dataset\n" if $ARGV[10] =~ m/\?/;
# Read /depot/data2/galaxy/alignseq.loc
@@ -47,7 +47,7 @@ if ($ARGV[8] > 0) {
$ARGV[8] = 1000000;
}
$ARGV[12] = int($ARGV[12]);
open (BED, "<$ARGV[1]") or die "Cannot open $ARGV[1] for reading :$!\n";
open (FASTA,">$ARGV[3]") or die "Cannot open $ARGV[3] for writing :$!\n";
@@ -64,13 +64,19 @@ while (<BED>) {
if (-e "$seqLocation{$ARGV[10]}$columns[$ARGV[5]].nib") {
$call = "fasta-subseq $seqLocation{$ARGV[10]}$columns[$ARGV[5]].nib ".($columns[$ARGV[6]]+1)." $columns[$ARGV[7]] $columns[$ARGV[8]]"; #+1 to start position to fix coordinate system
open (GET_SUBSEQ, "$call |") or die "Cannot start fasta-subseq:$!\n";
my $seq = "";
while (<GET_SUBSEQ>) {
if (!m/\>/) { # unless header print seq
print FASTA;
if ($ARGV[12]) {
chomp $_;
$seq = $seq . $_; }
else {
print FASTA; }
} else { # if header replace with seqdata
print FASTA ">$ARGV[10]_$columns[$ARGV[5]]_$columns[$ARGV[6]]_$columns[$ARGV[7]]_$columns[$ARGV[8]]\n";
if (!($ARGV[12])) { print FASTA ">$ARGV[10]_$columns[$ARGV[5]]_$columns[$ARGV[6]]_$columns[$ARGV[7]]_$columns[$ARGV[8]]\n"; }
}
}
if ($ARGV[12]) { print FASTA join("\t",@columns) . "\t$seq\n"; }
close (GET_SUBSEQ);
} else {
push (@errors, "Sequence $columns[$ARGV[5]] was not found for genome build $ARGV[10]\nMost likely your data lists wrong chromosome number for this organism\nCheck your genome build selection");
+7 -1
View File
@@ -1,15 +1,21 @@
<tool id="Extract genomic DNA 1" name="Extract genomic DNA">
<description>using coordinates from assembled genomes</description>
<command interpreter="perl">fasta-subseq-wrapper.pl -i $input -o $out_file1 -p $input_chromCol $input_startCol $input_endCol $input_strandCol -g $dbkey</command>
<command interpreter="perl">fasta-subseq-wrapper.pl -i $input -o $out_file1 -p $input_chromCol $input_startCol $input_endCol $input_strandCol -g $dbkey -f $out_format</command>
<inputs>
<param format="interval" name="input" type="data" label="Fetch sequences corresponding to Query"/>
<param name="out_format" type="select" label="Output Type">
<option value="0">FASTA</option>
<option value="1">Interval</option>
</param>
</inputs>
<outputs>
<data format="fasta" name="out_file1" />
</outputs>
<code file="fasta-subseq-wrapper_code.py" />
<tests>
<test>
<param name="input" value="1.bed" dbkey="hg17" ftype="bed" />
<param name="out_format" value="0"/>
<output name="out_file1" file="fsa_extract_genomic_dna.dat" />
</test>
</tests>
@@ -0,0 +1,5 @@
#Change format from FASTA to Interval if needed; use metadata from input file
def exec_before_job(app, inp_data, out_data, param_dict, tool):
if param_dict['out_format'] == "1":
out_data['out_file1'].change_datatype('interval')
out_data['out_file1'].init_meta( copy_from=inp_data['input'] )