From 08e0d5e5329530954669ee1803173413e977dc9a Mon Sep 17 00:00:00 2001 From: Greg Von Kuster Date: Fri, 4 Dec 2009 19:51:46 -0500 Subject: [PATCH] Fixes for Bam data type's set_meta() and sam_to_bam tool, indexes will now be properly created for bam datasets. Fixes also for uploading Bam files, they will no longer be uncompressed on upload. --- lib/galaxy/datatypes/binary.py | 49 ++++----- lib/galaxy/tools/util/maf_utilities.py | 11 +- test-data/1.bam | Bin 1495 -> 3592 bytes test-data/sam_to_bam_out1.bam | Bin 0 -> 3592 bytes tools/data_source/data_source.py | 2 + tools/data_source/upload.py | 57 +++++++---- tools/samtools/sam_to_bam.py | 134 +++++++++++++++---------- tools/samtools/sam_to_bam.xml | 44 ++++---- tools/sr_mapping/bowtie_wrapper.xml | 4 +- tools/sr_mapping/lastz_wrapper.xml | 31 +++--- 10 files changed, 190 insertions(+), 142 deletions(-) create mode 100644 test-data/sam_to_bam_out1.bam diff --git a/lib/galaxy/datatypes/binary.py b/lib/galaxy/datatypes/binary.py index 361719b49f9..074c5f1d735 100644 --- a/lib/galaxy/datatypes/binary.py +++ b/lib/galaxy/datatypes/binary.py @@ -7,7 +7,7 @@ from galaxy.datatypes.metadata import MetadataElement from galaxy.datatypes import metadata from galaxy.datatypes.sniff import * from urllib import urlencode, quote_plus -import zipfile +import zipfile, gzip import os, subprocess, tempfile log = logging.getLogger(__name__) @@ -54,32 +54,35 @@ class Bam( Binary ): def init_meta( self, dataset, copy_from=None ): Binary.init_meta( self, dataset, copy_from=copy_from ) - """ - GVK 12/2/09: just noticed this - not good and doesn't work, so commenting out for now. - def set_meta( self, dataset, overwrite = True, **kwd ): - # Sets index for BAM file. - index_file = dataset.metadata.bam_index - if not index_file: - index_file = dataset.metadata.spec['bam_index'].param.new_file( dataset = dataset ) + def set_meta( self, dataset, overwrite = True, **kwd ): + """ Sets index for BAM file. """ + # These metadata values are not accessible by users, always overwrite + index_file = dataset.metadata.bam_index + if not index_file: + index_file = dataset.metadata.spec['bam_index'].param.new_file( dataset = dataset ) + try: + # Using a symlink from ~/database/files/dataset_XX.dat, create a temporary file + # to store the indexex generated from samtools, something like ~/tmp/dataset_XX.dat.bai tmp_dir = tempfile.gettempdir() - tmpf1 = tempfile.NamedTemporaryFile( dir=tmp_dir ) - tmpf1bai = '%s.bai' % tmpf1.name - try: - os.system( 'cd %s' % tmp_dir ) - os.system( 'cp %s %s' % ( dataset.file_name, tmpf1.name ) ) - os.system( 'samtools index %s' % tmpf1.name ) - os.system( 'cp %s %s' % ( tmpf1bai, index_file.file_name ) ) - except Exception, ex: - sys.stderr.write( 'There was a problem creating the index for the BAM file\n%s\n' + str( ex ) ) - tmpf1.close() - if os.path.exists( tmpf1bai ): - os.remove( tmpf1bai ) - dataset.metadata.bam_index = index_file - """ + tmp_file_path = os.path.join( tmp_dir, os.path.basename( dataset.file_name ) ) + # Here tmp_file_path looks something like /tmp/dataset_XX.dat + os.symlink( dataset.file_name, tmp_file_path ) + command = 'samtools index %s' % tmp_file_path + proc = subprocess.Popen( args=command, shell=True ) + proc.wait() + except: + err_msg = 'Error creating index file (%s) for BAM file (%s)' % ( str( tmp_file_path ), str( dataset.file_name ) ) + log.exception( err_msg ) + sys.stderr.write( err_msg ) + # Move the temporary index file ~/tmp/dataset_XX.dat.bai to be ~/database/files/_metadata_files/dataset_XX.dat + shutil.move( '%s.bai' % ( tmp_file_path ), index_file.file_name ) + os.unlink( tmp_file_path ) + dataset.metadata.bam_index = index_file def sniff( self, filename ): + # BAM is compressed in the BGZF format, and must not be uncompressed in Galaxy. # The first 4 bytes of any bam file is 'BAM\1', and the file is binary. try: - header = open( filename ).read(4) + header = gzip.open( filename ).read(4) if binascii.b2a_hex( header ) == binascii.hexlify( 'BAM\1' ): return True return False diff --git a/lib/galaxy/tools/util/maf_utilities.py b/lib/galaxy/tools/util/maf_utilities.py index 46500eb1046..cdcae1f54ab 100644 --- a/lib/galaxy/tools/util/maf_utilities.py +++ b/lib/galaxy/tools/util/maf_utilities.py @@ -191,13 +191,16 @@ def open_or_build_maf_index( maf_file, index_filename, species = None ): def build_maf_index_species_chromosomes( filename, index_species = None ): species = [] species_chromosomes = {} - indexes = bx.interval_index_file.Indexes() + indexes = bx.interval_index_file.Indexes() + blocks = 0 try: maf_reader = bx.align.maf.Reader( open( filename ) ) while True: pos = maf_reader.file.tell() block = maf_reader.next() - if block is None: break + if block is None: + break + blocks += 1 for c in block.components: spec = c.src chrom = None @@ -225,11 +228,11 @@ def build_maf_index_species_chromosomes( filename, index_species = None ): #most likely a bad MAF log.debug( 'Building MAF index on %s failed: %s' % ( filename, e ) ) return ( None, [], {} ) - return ( indexes, species, species_chromosomes ) + return ( indexes, species, species_chromosomes, blocks ) #builds and returns ( index, index_filename ) for specified maf_file def build_maf_index( maf_file, species = None ): - indexes, found_species, species_chromosomes = build_maf_index_species_chromosomes( maf_file, species ) + indexes, found_species, species_chromosomes, blocks = build_maf_index_species_chromosomes( maf_file, species ) if indexes is not None: fd, index_filename = tempfile.mkstemp() out = os.fdopen( fd, 'w' ) diff --git a/test-data/1.bam b/test-data/1.bam index d49d771903658f48d1990240c00995866893a00e..95c65de857205099e47815d35fa2ae4861d9fba4 100644 GIT binary patch literal 3592 zcmV+j4)^gNiwFb&00000{{{d;LjnNn4V_v|Z)8^y^{`+OkVrQp;UGk4kE1-0&`gi# z(+fmR1y!<`fJ8||*+h;PW`~4&!PeFwc}14|1y=ks7KvEQ4(>hY)a&~u3C!3tw)^AV zck9;2Ij5@q%%0ymXkY(m-`{QDzdHZppP%;^4o=#4q-A=hdAeJgnOpxI(O!&bFSf_+ zNB?S9TI9Uc`_4?;UBLY{z;1-?zu03)+C`_I9sAr)gBKV18_Kb$db4-6k)94ul z#_&vAb~|yPV5%)V%#!Lqd!-m`Le>~E3y#OTO5=rx8MQ!6X(FslW~8M&hQb}|HMI>R z+ieTPiTicSok2h@+oEz&XJ@^bH}wn+ zc(Er@HxM3UIk{R)6g zxPQpuKoAuogaCQvbg`>W8FCvT97xJ~1HSi%fbxz-EY}{~o^LG%-NkzkpEch$SDq_0 z{rOoU-OH(95B^^QOp;%ajm|xzkdrCZPuZ_=*5k-I@CZysN`Wwp@)H!;Elos092U`kFED#*`0t}u`2m^6)vy3b(- zrug1@Hc>_C1Vr1#FT+`J#1!B$MgWGScRlI7Y`J#PQ*VzFyI%k2t)zQ8o}to=oEcvF z0o4&f_8b$DzEIt$EFp-L9RxHKosAG4uE4@2uEj+`#9}9Nko5+X=!IW#0g5tj;JR+2 zIEUhDL6QbhRz3zn(xF6kZW(+Rb-PH3R7vO|6Be$uG6e3Dt4F$GSqe)OykN&a*0=8g}=MVutqeNOYe~Zc??olhNuZW=XMskEOkp(y!jCO^pRQzCqSb;HHOi}a# zU<#6gy`458Ig)Xe$&AiEC2)pdX0wkgn(gE)ex}^yi(_29BoiW;w49ktg%Io`jF>Q_ zKYZi0vIhgZOR6cmQk})I^5D-yiGXcW1dZbMva=%qR_rUko-N^J35CHC;-Bmj}13zAH_ke$GT8@6_;B| z@A=4dQBwd@(4`InXi>wG;DI}B4v8o{U_Yn<2HQzuKOi1t$nJ8162;x<&VuqyL1p)e zl%cy`(f$P$J&fTrW((vEQh$7WT(_1=N|N1at4-kIae+lvpPn22{o!`gsE>#VZoFcs5!j2 zksXjmksLJyPs*~D;}`WBIWv??97&X@Ib*lOR~9G%VN|ojIru@D<&_{PsK|ws28|;y z1yEG*gdA|ztfBfvgPTkmL+z7D71IeqlGJ&SSQc%Ay5gK?o<1Y5#x}vxmA#^1yu`0m z7x0i|N?y&Hj3Sv9lOL;D#N05TR_6MyKFcK#1`l^2%aYkOhf8gvX&=MGj2w3jA#cqP z#I+&#FPb#f)uQOiJX5OWzMUnHsc9!M4Z&gHT)dVX&q5P1MR*RchHv240;V~1^iC8c z*32Wd#}Z66vZ>al)T<84TZh+KMuXX8l-#q@;ie==NS4EH-f2KHEiSzpXY5^sBFw`sPn;KW|44D{kV%{ops9y1{aT& zBh=J1wHSJ#isml%l_oT6=2NXQI|VOaN&bf>42wi`j6NLsoAVr}x@ApEL+#*_bJW5( zb(0XMyGpQv?FoqfwFvZH^>NF1R+n(9 z9fD3q)nsTUAwU>+b$E6q=DP(`UG`8U!c=i{N-@I7K|BYEi1NGLV`wXaQ@X0fJXc zraOsrC8SdYJmxUxOVBwr3S&a2WQk@4PR<;$2{{|3C9Nh?G5eg%gQugLYE(G)(`Bzy z=MnNk;I<%1uQ^<3zQ)hgJZi|G23`ZcC0%*xSM^U>pG@v2#w=qFX%^zgeO@3Am}`jH zXd)5h&Vc@81hW#Y!Q!tWn|m>qW_=5i+?S4feSlMw4b!K0HV+@hF}iESUR)z351Td1 z5kVu*UlH)d?1dpdQ@HXoDN+Hn`&|POhYnN8FgfHkR8?dZkNeVK6VKHrF)2@JKw&j= zuw}=g&y}Rytr`TO<2><|gBL^I{Wuw;qoMyy0($lGUf~-dJk;6pz54Vg_{lVtvC5%OCdK+u&D0M7rTZCIg1#*KD=kMPpzc^ zs4YY-Wj-@O-9k_!a(tuubBWF|hl%pg3#kkd6Vu-`o}| zFk+KNjpCP@l@VoJnyDoIVvNp?FG{Hh$ysByo2#36B_&W=Q4_e5lRzrAx}sT8=VC&j zO(dsXIwFN|Rv&JKROn(<1=o6J?9YSNH&_^v1*EL35r^kaM#G-mi#IY*!5kw$;4p#_ z1+#2bIhACugGiVB8o@N?k7FoKtuAE1*acJjX>__@>JobVY`la!GaWfxjzr}}j4Hfx zE!04|&m*9$*^1|^T*r@sLY~o|dL;k2s3Md6blx(m0~nU-^bA9Bbol-VfSO3XyWjPE zcd@(ZA1%K2?|JE)E$hF(VOMZ;@g>qOrG({4&EJOLxD6QY^#@T!In|rGrW!ep(XBkM zze6gtXOg8xL;G<6mVWPmGgB_aOJB6_|Fo}vzIAYL=bew%w+`MuIM!?Zt6%fsA3pl! zZ|wI^x1Vi4JGuAx$rH0rJ~`d?e~+F%IeoZ&bnxr%92{IFZy&Zt20t&;8aXTfg=EgHIoSeDCA?k7tv9_k$ambka) z4_AjrM_;_0uGZM|CqG*nvw2;U{{Evf>EKOFdc1vf`e1tklk95q?pN2BSF6p@<@IFO z)77gU`KGULH|ej(H!|tb$%E7T+bvDHz4`6Grx%Co7wh%)@s%B4{vVTCSns^|?IwNt z-i-qHaC>t4c>540O*XAJoAu>tb9MRR?w708dV1^o)AZBn&dqXZ_5RIrX?ySV=~H@X zI=XyubZNKkaCJGo`{i_aHNCT%j;D8~GU|tKWYS-Lej}3}Z?~uSA3TLg>&t|E@ zzWwa(%k}lq|45F9dBS>OHYcsgn}ygyApzM+%;_2@=Uy1#wW zJjGO6uU1oo$JOd^wVCdIdDyeNZ~FEdn)F{EsySUOEC2u>iwFb&00000{{{d;LjnLB O00RI3000000001Z^!VBU literal 1495 zcmbW%J&)5c7zc0%iKRl~Y$B1O6NlYY_ke^NH}9En?vU<=79<9$TTX(74e2IUTS3Z1 z$NLt1GlrfRx`CZEy#~~(%@f7R(+>XWZ^yZqSV)%{9_st^mzxQ)1wwQ4tFE6tga-Xf zVxc4C>i*BF%zqajFJHeEA5wnJuYE605+P=@jGHS;lPu&Bdh-CGQc|Geyre_{jFy-Z zaKMyZjqSm*R{eBxX3hASH95bEr`~iNuSxG8cOvtMK@bTw>92koVh`&rDxiUa}(r9@(?WHDM` zOrZA&;5iKbp_q2f+)l!$ok%|SvNWwz0wxP4g)JN`;bZ}hf<6Wp`Z|mn_2k~R2kG0( zPNbOgY#gLjN`zu6NDlh8^U;!&>|XHS+a9E^uR4)}lxLA2SENEG#27AfS8w$BGMD*) zZhmdUP1g?64|DM7`s5+RI-EhY)a&~u3C!3tw)^AV zck9;2Ij5@q%%0ymXkY(m-`{QDzdHZppP%;^4o=#4q-A=hdAeJgnOpxI(O!&bFSf_+ zNB?S9TI9Uc`_4?;UBLY{z;1-?zu03)+C`_I9sAr)gBKV18_Kb$db4-6k)94ul z#_&vAb~|yPV5%)V%#!Lqd!-m`Le>~E3y#OTO5=rx8MQ!6X(FslW~8M&hQb}|HMI>R z+ieTPiTicSok2h@+oEz&XJ@^bH}wn+ zc(Er@HxM3UIk{R)6g zxPQpuKoAuogaCQvbg`>W8FCvT97xJ~1HSi%fbxz-EY}{~o^LG%-NkzkpEch$SDq_0 z{rOoU-OH(95B^^QOp;%ajm|xzkdrCZPuZ_=*5k-I@CZysN`Wwp@)H!;Elos092U`kFED#*`0t}u`2m^6)vy3b(- zrug1@Hc>_C1Vr1#FT+`J#1!B$MgWGScRlI7Y`J#PQ*VzFyI%k2t)zQ8o}to=oEcvF z0o4&f_8b$DzEIt$EFp-L9RxHKosAG4uE4@2uEj+`#9}9Nko5+X=!IW#0g5tj;JR+2 zIEUhDL6QbhRz3zn(xF6kZW(+Rb-PH3R7vO|6Be$uG6e3Dt4F$GSqe)OykN&a*0=8g}=MVutqeNOYe~Zc??olhNuZW=XMskEOkp(y!jCO^pRQzCqSb;HHOi}a# zU<#6gy`458Ig)Xe$&AiEC2)pdX0wkgn(gE)ex}^yi(_29BoiW;w49ktg%Io`jF>Q_ zKYZi0vIhgZOR6cmQk})I^5D-yiGXcW1dZbMva=%qR_rUko-N^J35CHC;-Bmj}13zAH_ke$GT8@6_;B| z@A=4dQBwd@(4`InXi>wG;DI}B4v8o{U_Yn<2HQzuKOi1t$nJ8162;x<&VuqyL1p)e zl%cy`(f$P$J&fTrW((vEQh$7WT(_1=N|N1at4-kIae+lvpPn22{o!`gsE>#VZoFcs5!j2 zksXjmksLJyPs*~D;}`WBIWv??97&X@Ib*lOR~9G%VN|ojIru@D<&_{PsK|ws28|;y z1yEG*gdA|ztfBfvgPTkmL+z7D71IeqlGJ&SSQc%Ay5gK?o<1Y5#x}vxmA#^1yu`0m z7x0i|N?y&Hj3Sv9lOL;D#N05TR_6MyKFcK#1`l^2%aYkOhf8gvX&=MGj2w3jA#cqP z#I+&#FPb#f)uQOiJX5OWzMUnHsc9!M4Z&gHT)dVX&q5P1MR*RchHv240;V~1^iC8c z*32Wd#}Z66vZ>al)T<84TZh+KMuXX8l-#q@;ie==NS4EH-f2KHEiSzpXY5^sBFw`sPn;KW|44D{kV%{ops9y1{aT& zBh=J1wHSJ#isml%l_oT6=2NXQI|VOaN&bf>42wi`j6NLsoAVr}x@ApEL+#*_bJW5( zb(0XMyGpQv?FoqfwFvZH^>NF1R+n(9 z9fD3q)nsTUAwU>+b$E6q=DP(`UG`8U!c=i{N-@I7K|BYEi1NGLV`wXaQ@X0fJXc zraOsrC8SdYJmxUxOVBwr3S&a2WQk@4PR<;$2{{|3C9Nh?G5eg%gQugLYE(G)(`Bzy z=MnNk;I<%1uQ^<3zQ)hgJZi|G23`ZcC0%*xSM^U>pG@v2#w=qFX%^zgeO@3Am}`jH zXd)5h&Vc@81hW#Y!Q!tWn|m>qW_=5i+?S4feSlMw4b!K0HV+@hF}iESUR)z351Td1 z5kVu*UlH)d?1dpdQ@HXoDN+Hn`&|POhYnN8FgfHkR8?dZkNeVK6VKHrF)2@JKw&j= zuw}=g&y}Rytr`TO<2><|gBL^I{Wuw;qoMyy0($lGUf~-dJk;6pz54Vg_{lVtvC5%OCdK+u&D0M7rTZCIg1#*KD=kMPpzc^ zs4YY-Wj-@O-9k_!a(tuubBWF|hl%pg3#kkd6Vu-`o}| zFk+KNjpCP@l@VoJnyDoIVvNp?FG{Hh$ysByo2#36B_&W=Q4_e5lRzrAx}sT8=VC&j zO(dsXIwFN|Rv&JKROn(<1=o6J?9YSNH&_^v1*EL35r^kaM#G-mi#IY*!5kw$;4p#_ z1+#2bIhACugGiVB8o@N?k7FoKtuAE1*acJjX>__@>JobVY`la!GaWfxjzr}}j4Hfx zE!04|&m*9$*^1|^T*r@sLY~o|dL;k2s3Md6blx(m0~nU-^bA9Bbol-VfSO3XyWjPE zcd@(ZA1%K2?|JE)E$hF(VOMZ;@g>qOrG({4&EJOLxD6QY^#@T!In|rGrW!ep(XBkM zze6gtXOg8xL;G<6mVWPmGgB_aOJB6_|Fo}vzIAYL=bew%w+`MuIM!?Zt6%fsA3pl! zZ|wI^x1Vi4JGuAx$rH0rJ~`d?e~+F%IeoZ&bnxr%92{IFZy&Zt20t&;8aXTfg=EgHIoSeDCA?k7tv9_k$ambka) z4_AjrM_;_0uGZM|CqG*nvw2;U{{Evf>EKOFdc1vf`e1tklk95q?pN2BSF6p@<@IFO z)77gU`KGULH|ej(H!|tb$%E7T+bvDHz4`6Grx%Co7wh%)@s%B4{vVTCSns^|?IwNt z-i-qHaC>t4c>540O*XAJoAu>tb9MRR?w708dV1^o)AZBn&dqXZ_5RIrX?ySV=~H@X zI=XyubZNKkaCJGo`{i_aHNCT%j;D8~GU|tKWYS-Lej}3}Z?~uSA3TLg>&t|E@ zzWwa(%k}lq|45F9dBS>OHYcsgn}ygyApzM+%;_2@=Uy1#wW zJjGO6uU1oo$JOd^wVCdIdDyeNZ~FEdn)F{EsySUOEC2u>iwFb&00000{{{d;LjnLB O00RI3000000001Z^!VBU literal 0 HcmV?d00001 diff --git a/tools/data_source/data_source.py b/tools/data_source/data_source.py index f48754d2cc5..71f1e5a5bbb 100644 --- a/tools/data_source/data_source.py +++ b/tools/data_source/data_source.py @@ -12,6 +12,7 @@ def stop_err( msg ): sys.exit() def check_gzip( filename ): + # TODO: This needs to check for BAM files since they are compressed and must remain so ( see upload.py ) temp = open( filename, "U" ) magic_check = temp.read( 2 ) temp.close() @@ -66,6 +67,7 @@ def __main__(): out.write( chunk ) out.close() if check_gzip( filename ): + # TODO: This needs to check for BAM files since they are compressed and must remain so ( see upload.py ) fd, uncompressed = tempfile.mkstemp() gzipped_file = gzip.GzipFile( filename ) while 1: diff --git a/tools/data_source/upload.py b/tools/data_source/upload.py index cc99a82f833..db2d7f91fe5 100644 --- a/tools/data_source/upload.py +++ b/tools/data_source/upload.py @@ -4,7 +4,7 @@ # WARNING: Changes in this tool (particularly as related to parsing) may need # to be reflected in galaxy.web.controllers.tool_runner and galaxy.tools -import urllib, sys, os, gzip, tempfile, shutil, re, gzip, zipfile, codecs +import urllib, sys, os, gzip, tempfile, shutil, re, gzip, zipfile, codecs, binascii from galaxy import eggs # need to import model before sniff to resolve a circular import dependency import galaxy.model @@ -18,7 +18,6 @@ assert sys.version_info[:2] >= ( 2, 4 ) def stop_err( msg, ret=1 ): sys.stderr.write( msg ) sys.exit( ret ) - def file_err( msg, dataset, json_file ): json_file.write( to_json_string( dict( type = 'dataset', ext = 'data', @@ -28,7 +27,6 @@ def file_err( msg, dataset, json_file ): os.remove( dataset.path ) except: pass - def safe_dict(d): """ Recursively clone json structure with UTF-8 dictionary keys @@ -40,7 +38,6 @@ def safe_dict(d): return [safe_dict(x) for x in d] else: return d - def check_html( temp_name, chunk=None ): if chunk is None: temp = open(temp_name, "U") @@ -64,7 +61,6 @@ def check_html( temp_name, chunk=None ): if chunk is None: temp.close() return False - def check_binary( temp_name, chunk=None ): if chunk is None: temp = open( temp_name, "U" ) @@ -85,21 +81,42 @@ def check_binary( temp_name, chunk=None ): if chunk is None: temp.close() return False - def check_gzip( temp_name ): + # This is sort of hacky. BAM is compressed in the BGZF format, and must + # not be uncompressed in upon upload ( it will be detected as gzipped ). + # The tuple we're returning from here contains boolean values for + # ( is_compressed, is_valid, is_bam ). temp = open( temp_name, "U" ) magic_check = temp.read( 2 ) temp.close() if magic_check != util.gzip_magic: - return ( False, False ) + return ( False, False, False ) CHUNK_SIZE = 2**15 # 32Kb gzipped_file = gzip.GzipFile( temp_name ) chunk = gzipped_file.read( CHUNK_SIZE ) gzipped_file.close() - if check_html( temp_name, chunk=chunk ) or check_binary( temp_name, chunk=chunk ): - return( True, False ) - return ( True, True ) - + if check_html( temp_name, chunk=chunk ): + return ( True, False, False ) + if check_binary( temp_name, chunk=chunk ): + # We do support some binary data types, so check if the compressed binary file is valid + # We currently only check for [ 'sff', 'bam' ] + # TODO: this should be fixed to more easily support future-supported binary data types. + # This is currently just copied from the sniff methods. + # The first 4 bytes of any bam file is 'BAM\1', and the file is binary. + try: + header = gzip.open( temp_name ).read(4) + if binascii.b2a_hex( header ) == binascii.hexlify( 'BAM\1' ): + return ( True, True, True ) + except: + pass + try: + header = gzip.open( temp_name ).read(4) + if binascii.b2a_hex( header ) == binascii.hexlify( '.sff' ): + return ( True, True, False ) + except: + pass + return ( True, False, False ) + return ( True, True, False ) def check_zip( temp_name ): if not zipfile.is_zipfile( temp_name ): return ( False, False, None ) @@ -116,14 +133,12 @@ def check_zip( temp_name ): if ext != test_ext: return ( True, False, test_ext ) return ( True, True, test_ext ) - def parse_outputs( args ): rval = {} for arg in args: id, files_path, path = arg.split( ':', 2 ) rval[int( id )] = ( path, files_path ) return rval - def add_file( dataset, json_file, output_path ): data_type = None line_count = None @@ -153,15 +168,19 @@ def add_file( dataset, json_file, output_path ): ext = sniff.guess_ext( dataset.path, is_multi_byte=True ) else: # See if we have a gzipped file, which, if it passes our restrictions, we'll uncompress - is_gzipped, is_valid = check_gzip( dataset.path ) + is_gzipped, is_valid, is_bam = check_gzip( dataset.path ) if is_gzipped and not is_valid: file_err( 'The uploaded file contains inappropriate content', dataset, json_file ) return - elif is_gzipped and is_valid: - # We need to uncompress the temp_name file + elif is_gzipped and is_valid and is_bam: + ext = 'bam' + data_type = 'bam' + elif is_gzipped and is_valid and not is_bam: + # We need to uncompress the temp_name file, but BAM files must remain compressed + # in order for samtools to function on them CHUNK_SIZE = 2**20 # 1Mb - fd, uncompressed = tempfile.mkstemp( prefix='data_id_%s_upload_gunzip_' % dataset.dataset_id, dir=os.path.dirname( dataset.path ) ) - gzipped_file = gzip.GzipFile( dataset.path ) + fd, uncompressed = tempfile.mkstemp( prefix='data_id_%s_upload_gunzip_' % dataset.dataset_id, dir=os.path.dirname( dataset.path ), text=False ) + gzipped_file = gzip.GzipFile( dataset.path, 'rb' ) while 1: try: chunk = gzipped_file.read( CHUNK_SIZE ) @@ -229,7 +248,7 @@ def add_file( dataset, json_file, output_path ): if check_html( dataset.path ): file_err( 'The uploaded file contains inappropriate content', dataset, json_file ) return - if data_type != 'binary' and data_type != 'zip': + if data_type != 'bam' and data_type != 'binary' and data_type != 'zip': if dataset.space_to_tab: line_count = sniff.convert_newlines_sep2tabs( dataset.path ) else: diff --git a/tools/samtools/sam_to_bam.py b/tools/samtools/sam_to_bam.py index 7f7458282ce..114d1b76b07 100644 --- a/tools/samtools/sam_to_bam.py +++ b/tools/samtools/sam_to_bam.py @@ -1,31 +1,27 @@ #! /usr/bin/python - """ -Converts SAM data to BAM format. - -usage: %prog [options] - -i, --input1=i: SAM file to be converted - -d, --dbkey=d: dbkey value - -r, --ref_file=r: Reference file if choosing from history - -o, --output1=o: BAM output - -x, --index_dir=x: Index directory - -usage: %prog input_file dbkey ref_list output_file +Converts SAM data to sorted BAM data. +usage: sam_to_bam.py [options] + --input1: SAM file to be converted + --dbkey: dbkey value + --ref_file: Reference file if choosing from history + --output1: output dataset in bam format + --index_dir: GALAXY_DATA_INDEX_DIR """ -import os, sys, tempfile +import optparse, os, sys, subprocess, tempfile, shutil, gzip from galaxy import eggs import pkg_resources; pkg_resources.require( "bx-python" ) from bx.cookbook import doc_optparse +from galaxy import util def stop_err( msg ): sys.stderr.write( "%s\n" % msg ) sys.exit() -def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ): - seq_file = "%s/sam_fa_indices.loc" % GALAXY_DATA_INDEX_DIR +def check_seq_file( dbkey, cached_seqs_pointer_file ): seq_path = '' - for line in open( seq_file ): + for line in open( cached_seqs_pointer_file ): line = line.rstrip( '\r\n' ) if line and not line.startswith( "#" ) and line.startswith( 'index' ): fields = line.split( '\t' ) @@ -38,48 +34,80 @@ def check_seq_file( dbkey, GALAXY_DATA_INDEX_DIR ): def __main__(): #Parse Command Line - options, args = doc_optparse.parse( __doc__ ) - seq_path = check_seq_file( options.dbkey, options.index_dir ) + parser = optparse.OptionParser() + parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' ) + parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' ) + parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) + parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' ) + parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' ) + ( options, args ) = parser.parse_args() + + cached_seqs_pointer_file = "%s/sam_fa_indices.loc" % options.index_dir + if not os.path.exists( cached_seqs_pointer_file ): + stop_err( "The required file (%s) does not exist." % cached_seqs_pointer_file ) + # If found for the dbkey, seq_path will look something like /depot/data2/galaxy/equCab2/sam_index/equCab2.fa, + # and the equCab2.fa file will contain fasta sequences. + seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file ) tmp_dir = tempfile.gettempdir() - os.chdir(tmp_dir) - tmpf1 = tempfile.NamedTemporaryFile(dir=tmp_dir) - tmpf1fai = '%s.fai' % tmpf1.name - tmpf2 = tempfile.NamedTemporaryFile(dir=tmp_dir) - tmpf3 = tempfile.NamedTemporaryFile(dir=tmp_dir) - tmpf3bam = '%s.bam' % tmpf3.name if options.ref_file == "None": - full_path = "%s.fai" % seq_path - if not os.path.exists( full_path ): - stop_err( "No sequences are available for '%s', request them by reporting this error." % options.dbkey ) - cmd1 = "cp %s %s; cp %s %s" % (seq_path, tmpf1.name, full_path, tmpf1fai) + # We're using locally cached reference sequences( e.g., /depot/data2/galaxy/equCab2/sam_index/equCab2.fa ). + # The indexes for /depot/data2/galaxy/equCab2/sam_index/equCab2.fa will be contained in + # a file named /depot/data2/galaxy/equCab2/sam_index/equCab2.fa.fai + fai_index_file_path = "%s.fai" % seq_path + if not os.path.exists( fai_index_file_path ): + stop_err( "No sequences are available for build (%s), request them by reporting this error." % options.dbkey ) else: - cmd1 = "cp %s %s; samtools faidx %s 2>/dev/null" % (options.ref_file, tmpf1.name, tmpf1.name) - cmd2 = "samtools view -bt %s -o %s %s 2>/dev/null" % (tmpf1fai, tmpf2.name, options.input1) - cmd3 = "samtools sort %s %s 2>/dev/null" % (tmpf2.name, tmpf3.name) - cmd4 = "cp %s %s" % (tmpf3bam, options.output1) - # either create index based on fa file or copy provided index to temp directory + try: + # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will: + # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence + # - if no region is specified, faidx will index the file and create .fai on the disk + # - if regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format + # - the input file can be compressed in the RAZF format. + # IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset + # every time we run this tool. It would be nice if we could somehow keep track of user's specific + # index files so they could be re-used. + fai_index_file_path = os.path.join( tmp_dir, os.path.basename( options.ref_file ) ) + # At this point, fai_index_file_path will look something like /tmp/dataset_13.dat + os.symlink( options.ref_file, fai_index_file_path ) + command = "samtools faidx %s 2>/dev/null" % fai_index_file_path + proc = subprocess.Popen( args=command, shell=True ) + proc.wait() + except Exception, e: + stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) ) try: - os.system(cmd1) - except Exception, eq: - stop_err("Error creating the reference list index.\n" + str(eq)) - # create original bam file + # Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ). + tmp_aligns_file = tempfile.NamedTemporaryFile() + tmp_aligns_file_name = tmp_aligns_file.name + tmp_aligns_file.close() + # IMPORTANT NOTE: for some reason the samtools view command gzips the resulting bam file without warning, + # and the docs do not currently state that this occurs ( very bad ). + command = "samtools view -bt %s -o %s %s 2>/dev/null" % ( fai_index_file_path, tmp_aligns_file_name, options.input1 ) + proc = subprocess.Popen( args=command, shell=True ) + proc.wait() + except Exception, e: + stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) ) try: - os.system(cmd2) - except Exception, eq: - stop_err("Error running view command.\n" + str(eq)) - # sort original bam file to produce sorted output bam file - try: - os.system(cmd3) - os.system(cmd4) - except Exception, eq: - stop_err("Error sorting data and creating output file.\n" + str(eq)) - # cleanup temp files - tmpf1.close() - tmpf2.close() - tmpf3.close() - if os.path.exists(tmpf1fai): - os.remove(tmpf1fai) - if os.path.exists(tmpf3bam): - os.remove(tmpf3bam) + # Sort alignments by leftmost coordinates. File .bam will be created. This command + # may also create temporary files .%d.bam when the whole alignment cannot be fitted + # into memory ( controlled by option -m ). + tmp_sorted_aligns_file = tempfile.NamedTemporaryFile() + tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name + tmp_sorted_aligns_file.close() + command = "samtools sort %s %s 2>/dev/null" % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name ) + proc = subprocess.Popen( args=command, shell=True ) + proc.wait() + except Exception, e: + stop_err( 'Error sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) ) + # Move tmp_aligns_file_name to our output dataset location + sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name + shutil.move( sorted_bam_file, options.output1 ) + if options.ref_file != "None": + # Remove the symlink from /tmp/dataset_13.dat to ~/database/files/000/dataset_13.dat + os.unlink( fai_index_file_path ) + # Remove the index file + index_file_name = '%s.fai' % fai_index_file_path + os.unlink( index_file_name ) + # Remove the tmp_aligns_file_name + os.unlink( tmp_aligns_file_name ) if __name__=="__main__": __main__() diff --git a/tools/samtools/sam_to_bam.xml b/tools/samtools/sam_to_bam.xml index 0c17ed261e2..d8302ed4c8b 100644 --- a/tools/samtools/sam_to_bam.xml +++ b/tools/samtools/sam_to_bam.xml @@ -1,32 +1,29 @@ converts SAM format to BAM format - sam_to_bam.py - --input1=$source.input1 - --dbkey=${input1.metadata.dbkey} - #if $source.indexSource == "history": - --ref_file=$ref_file - #else - --ref_file="None" - #end if - --output1=$output1 - --index_dir=${GALAXY_DATA_INDEX_DIR} +sam_to_bam.py --input1=$source.input1 --dbkey=${input1.metadata.dbkey} +#if $source.index_source == "history": +--ref_file=$source.ref_file +#else +--ref_file="None" +#end if +--output1=$output1 --index_dir=${GALAXY_DATA_INDEX_DIR} - - + + - + - - + + @@ -34,19 +31,16 @@ + - - - + + + - diff --git a/tools/sr_mapping/bowtie_wrapper.xml b/tools/sr_mapping/bowtie_wrapper.xml index 501f7326262..c74346e1ad6 100644 --- a/tools/sr_mapping/bowtie_wrapper.xml +++ b/tools/sr_mapping/bowtie_wrapper.xml @@ -338,7 +338,7 @@ - + @@ -349,7 +349,7 @@ - + diff --git a/tools/sr_mapping/lastz_wrapper.xml b/tools/sr_mapping/lastz_wrapper.xml index 2335bb458b5..4222b899aa8 100644 --- a/tools/sr_mapping/lastz_wrapper.xml +++ b/tools/sr_mapping/lastz_wrapper.xml @@ -125,10 +125,13 @@ lastz - + + @@ -156,10 +159,13 @@ - + + @@ -173,14 +179,7 @@ - +