Mercurial > repos > jjohnson > gmap
changeset 6:3be0e0a858fe
refactor and update README
| author | Jim Johnson <jj@umn.edu> |
|---|---|
| date | Tue, 08 Nov 2011 13:22:34 -0600 |
| parents | f4b4c1712e39 |
| children | 561503a442f0 |
| files | README gmap/README gmap/gmap_indices.loc.sample gmap/lib/galaxy/datatypes/gmap.py lib/galaxy/datatypes/gmap.py tool-data/gmap_indices.loc.sample |
| diffstat | 6 files changed, 553 insertions(+), 542 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README Tue Nov 08 13:22:34 2011 -0600 @@ -0,0 +1,71 @@ +GMAP applications and citation info are available from: http://research-pub.gene.com/gmap/ + + + Installation instructions are in the README file in the download, + and online: http://research-pub.gene.com/gmap/src/README + + These tools were consistent with gmap version: 2011-10-16 + + +GMAP and GSNAP use added datatypes: + + add datatype definition file: lib/galaxy/datatypes/gmap.py + + add the following import line to: lib/galaxy/datatypes/registry.py + import gmap # added for gmap tools + + add to datatypes_conf.xml + <!-- Start GMAP Datatypes --> + <datatype extension="gmapdb" type="galaxy.datatypes.gmap:GmapDB" display_in_upload="False"/> + <datatype extension="gmapsnpindex" type="galaxy.datatypes.gmap:GmapSnpIndex" display_in_upload="False"/> + <datatype extension="iit" type="galaxy.datatypes.gmap:IntervalIndexTree" display_in_upload="True"/> + <datatype extension="splicesites.iit" type="galaxy.datatypes.gmap:SpliceSitesIntervalIndexTree" display_in_upload="True"/> + <datatype extension="introns.iit" type="galaxy.datatypes.gmap:IntronsIntervalIndexTree" display_in_upload="True"/> + <datatype extension="snps.iit" type="galaxy.datatypes.gmap:SNPsIntervalIndexTree" display_in_upload="True"/> + <datatype extension="tally.iit" type="galaxy.datatypes.gmap:TallyIntervalIndexTree" display_in_upload="True"/> + <datatype extension="gmap_annotation" type="galaxy.datatypes.gmap:IntervalAnnotation" display_in_upload="False"/> + <datatype extension="gmap_splicesites" type="galaxy.datatypes.gmap:SpliceSiteAnnotation" display_in_upload="True"/> + <datatype extension="gmap_introns" type="galaxy.datatypes.gmap:IntronAnnotation" display_in_upload="True"/> + <datatype extension="gmap_snps" type="galaxy.datatypes.gmap:SNPAnnotation" display_in_upload="True"/> + <datatype extension="gsnap_tally" type="galaxy.datatypes.gmap:TallyAnnotation" display_in_upload="True"/> + <datatype extension="gsnap" type="galaxy.datatypes.gmap:GsnapResult" display_in_upload="True"/> + <!-- End GMAP Datatypes --> + +Tools: + GMAP_Build - create a GmapDB set of index files for a reference sequence and optional set of annotations + GMAP - map sequences to a reference sequence GmapDB index + GSNAP - align sequences to a reference and detect splicing + + Add to tool_conf.xml ( probably in the "NGS: Mapping" section ) + <tool file="gmap/gmap.xml" /> + <tool file="gmap/gsnap.xml" /> + <tool file="gmap/gmap_build.xml" /> + <tool file="gmap/snpindex.xml" /> + <tool file="gmap/iit_store.xml" /> + +Admin built cached gmapdb indexes defined in tool-data/gmap_indices.loc + + +TODO: + + + Add classes to gmap.py + CmetIndex - an index created by cmetindex + AtoiIndex - an index created by atoiindex + + Add tally creation + gsnap default output -> gsnap_tally -> iit_store + + Add goby support + Should add separate tools and datatypes for goby + GSNAP goby output relies on goby input, might be better to have a separate gsnap tool for goby + + Possibly add Tools: + get_genome - retrieves from a gmapdb + cmetindex - create methylcytosine index + atoiindex - create A-to-I RNA editing index + + + + +
--- a/gmap/README Tue Nov 08 13:07:25 2011 -0600 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -GMAP applications and citation info are available from: http://research-pub.gene.com/gmap/ - - - Installation instructions are in the README file in the download, - and online: http://research-pub.gene.com/gmap/src/README - - These galaxy tools assume gmap is configured and built with goby support. - These tools were consistent with gmap version: 2011-10-07 - - -GMAP and GSNAP use added datatypes: - - add datatype definition file: lib/galaxy/datatypes/gmap.py - - add the following import line to: lib/galaxy/datatypes/registry.py - import gmap # added for gmap tools - - add to datatypes_conf.xml - <!-- Start GMAP Datatypes --> - <datatype extension="gmapdb" type="galaxy.datatypes.gmap:GmapDB" display_in_upload="False"/> - <datatype extension="gmapsnpindex" type="galaxy.datatypes.gmap:GmapSnpIndex" display_in_upload="False"/> - <datatype extension="iit" type="galaxy.datatypes.gmap:IntervalIndexTree" display_in_upload="True"/> - <datatype extension="splicesites.iit" type="galaxy.datatypes.gmap:SpliceSitesIntervalIndexTree" display_in_upload="True"/> - <datatype extension="introns.iit" type="galaxy.datatypes.gmap:IntronsIntervalIndexTree" display_in_upload="True"/> - <datatype extension="snps.iit" type="galaxy.datatypes.gmap:SNPsIntervalIndexTree" display_in_upload="True"/> - <datatype extension="gmap_annotation" type="galaxy.datatypes.gmap:IntervalAnnotation" display_in_upload="False"/> - <datatype extension="gmap_splicesites" type="galaxy.datatypes.gmap:SpliceSiteAnnotation" display_in_upload="True"/> - <datatype extension="gmap_introns" type="galaxy.datatypes.gmap:IntronAnnotation" display_in_upload="True"/> - <datatype extension="gmap_snps" type="galaxy.datatypes.gmap:SNPAnnotation" display_in_upload="True"/> - <!-- End GMAP tools --> - -Tools: - GMAP_Build - create a GmapDB set of index files for a reference sequence and optional set of annotations - GMAP - map sequences to a reference sequence GmapDB index - GSNAP - align sequences to a reference and detect splicing - - Add to tool_conf.xml ( probably in the "NGS: Mapping" section ) - <tool file="gmap/gmap.xml" /> - <tool file="gmap/gsnap.xml" /> - <tool file="gmap/gmap_build.xml" /> - <tool file="gmap/snpindex.xml" /> - <tool file="gmap/iit_store.xml" /> - -Admin built cached gmapdb indexes defined in tool-data/gmap_indices.loc - - -TODO: - Add classes to gmap.py - CmetIndex - an index created by cmetindex - AtoiIndex - an index created by atoiindex - - Possibly add Tools: - get_genome - retrieves from a gmapdb - cmetindex - create methylcytosine index - atoiindex - create A-to-I RNA editing index - - - - -
--- a/gmap/gmap_indices.loc.sample Tue Nov 08 13:07:25 2011 -0600 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -#This is a sample file distributed with Galaxy that enables tools -#to use a directory of GMAPDB indexed sequences data files. You will need -#to create these data files using gmap_build and then create a gmap_indices.loc file -#similar to this one (store it in this directory) that points to -#the directories in which those files are stored. The gmap_indices.loc -#file has this format (white space characters are TAB characters): -# -#<unique_build_id> <dbkey> <display_name> <kmers> <map,map> <snp,snp> <file_base_path> -#hg18 hg18 hg18 (cmet atoi) 12,13,14,15 splicesites,introns snps /depot/data2/galaxy/gmap/hg18 -#hg19 hg19 hg19 (cmet atoi) 12,13,14,15 splicesites,introns,snps snps,dbsnp /depot/data2/galaxy/gmap/hg19
--- a/gmap/lib/galaxy/datatypes/gmap.py Tue Nov 08 13:07:25 2011 -0600 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,472 +0,0 @@ -""" -GMAP indexes -""" -import logging -import os,os.path,re -import data -from data import Text -from galaxy import util -from metadata import MetadataElement - -log = logging.getLogger(__name__) - -class GmapDB( Text ): - """ - A GMAP DB for indexes - """ - MetadataElement( name="db_name", desc="The db name for this index set", default='unknown', set_in_upload=True, readonly=True ) - MetadataElement( name="basesize", default="12", desc="The basesize for offsetscomp", visible=True, readonly=True ) - MetadataElement( name="kmers", default=[''], desc="The kmer sizes for indexes", visible=True, no_value=[''], readonly=True ) - MetadataElement( name="map_dir", desc="The maps directory", default='unknown', set_in_upload=True, readonly=True ) - MetadataElement( name="maps", default=[''], desc="The names of maps stored for this gmap gmapdb", visible=True, no_value=[''], readonly=True ) - MetadataElement( name="snps", default=[''], desc="The names of SNP indexes stored for this gmapdb", visible=True, no_value=[''], readonly=True ) - MetadataElement( name="cmet", default=False, desc="Has a cmet index", visible=True, readonly=True ) - MetadataElement( name="atoi", default=False, desc="Has a atoi index", visible=True, readonly=True ) - - file_ext = 'gmapdb' - is_binary = True - composite_type = 'auto_primary_file' - allow_datatype_change = False - - def generate_primary_file( self, dataset = None ): - """ - This is called only at upload to write the html file - cannot rename the datasets here - they come with the default unfortunately - """ - return '<html><head></head><body>AutoGenerated Primary File for Composite Dataset</body></html>' - - def regenerate_primary_file(self,dataset): - """ - cannot do this until we are setting metadata - """ - bn = dataset.metadata.db_name - log.info( "GmapDB regenerate_primary_file %s" % (bn)) - rval = ['<html><head><title>GMAPDB %s</title></head><p/><H3>GMAPDB %s</H3><p/>cmet %s<br>atoi %s<H4>Maps:</H4><ul>' % (bn,bn,dataset.metadata.cmet,dataset.metadata.atoi)] - for i,name in enumerate(dataset.metadata.maps): - rval.append( '<li>%s' % name) - rval.append( '</ul></html>' ) - f = file(dataset.file_name,'w') - f.write("\n".join( rval )) - f.write('\n') - f.close() - - def set_peek( self, dataset, is_multi_byte=False ): - log.info( "GmapDB set_peek %s" % (dataset)) - if not dataset.dataset.purged: - dataset.peek = "GMAPDB index %s\n cmet %s\n atoi %s\n maps %s" % ( dataset.metadata.db_name,dataset.metadata.cmet,dataset.metadata.atoi,dataset.metadata.maps ) - dataset.blurb = "GMAPDB %s" % ( dataset.metadata.db_name ) - else: - dataset.peek = 'file does not exist' - dataset.blurb = 'file purged from disk' - def display_peek( self, dataset ): - try: - return dataset.peek - except: - return "GMAP index file" - - def sniff( self, filename ): - return False - def set_meta( self, dataset, overwrite = True, **kwd ): - """ - Expecting: - extra_files_path/<db_name>/db_name>.ref<basesize><kmer>3<index> - extra_files_path/db_name/db_name.ref1[2345]1[2345]3offsetscomp - extra_files_path/db_name/db_name.ref1[2345]1[2345]3positions - extra_files_path/db_name/db_name.ref1[2345]1[2345]3gammaptrs - index maps: - extra_files_path/db_name/db_name.maps/*.iit - """ - log.info( "GmapDB set_meta %s %s" % (dataset,dataset.extra_files_path)) - pat = '(.*)\.((ref)|(met)[atgc][atgc]|(a2i)[atgc][atgc])((\d\d)(\d\d))?3positions(\.(.+))?' - efp = dataset.extra_files_path - flist = os.listdir(efp) - for i,fname in enumerate(flist): - log.info( "GmapDB set_meta %s %s" % (i,fname)) - fpath = os.path.join(efp,fname) - if os.path.isdir(fpath): - ilist = os.listdir(fpath) - kmers = {'':'default'} # HACK '' empty key added so user has default choice when selecting kmer from metadata - for j,iname in enumerate(ilist): - log.info( "GmapDB set_meta file %s %s" % (j,iname)) - ipath = os.path.join(fpath,iname) - if os.path.isdir(ipath): # find maps - dataset.metadata.map_dir = iname - for mapfile in os.listdir(ipath): - mapname = mapfile.replace('.iit','') - log.info( "GmapDB set_meta map %s %s" % (mapname,mapfile)) - dataset.metadata.maps.append(mapname) - else: - m = re.match(pat,iname) - if m: - log.info( "GmapDB set_meta m %s %s " % (iname, m)) - assert len(m.groups()) == 10 - dataset.metadata.db_name = fname - if m.groups()[2] == 'ref': - if m.groups()[-1] != None: - dataset.metadata.snps.append(m.groups()[-1]) - else: - if m.groups()[-3] != None: - k = int(m.groups()[-3]) - kmers[k] = k - if m.groups()[-4] != None: - dataset.metadata.basesize = int( m.groups()[-4]) - elif m.groups()[3] == 'met': - dataset.metadata.cmet = True - elif m.groups()[4] == 'a2i': - dataset.metadata.atoi = True - dataset.metadata.kmers = kmers.keys() - -class GmapSnpIndex( Text ): - """ - A GMAP SNP index created by snpindex - """ - MetadataElement( name="db_name", desc="The db name for this index set", default='unknown', set_in_upload=True, readonly=True ) - MetadataElement( name="snps_name", default='snps', desc="The name of SNP index", visible=True, no_value='', readonly=True ) - - file_ext = 'gmapsnpindex' - is_binary = True - composite_type = 'auto_primary_file' - allow_datatype_change = False - - def generate_primary_file( self, dataset = None ): - """ - This is called only at upload to write the html file - cannot rename the datasets here - they come with the default unfortunately - """ - return '<html><head></head><body>AutoGenerated Primary File for Composite Dataset</body></html>' - - def regenerate_primary_file(self,dataset): - """ - cannot do this until we are setting metadata - """ - bn = dataset.metadata.db_name - log.info( "GmapDB regenerate_primary_file %s" % (bn)) - rval = ['<html><head><title>GMAPDB %s</title></head><p/><H3>GMAPDB %s</H3><p/>cmet %s<br>atoi %s<H4>Maps:</H4><ul>' % (bn,bn,dataset.metadata.cmet,dataset.metadata.atoi)] - for i,name in enumerate(dataset.metadata.maps): - rval.append( '<li>%s' % name) - rval.append( '</ul></html>' ) - f = file(dataset.file_name,'w') - f.write("\n".join( rval )) - f.write('\n') - f.close() - def set_peek( self, dataset, is_multi_byte=False ): - log.info( "GmapSnpIndex set_peek %s" % (dataset)) - if not dataset.dataset.purged: - dataset.peek = "GMAP SNPindex %s on %s\n" % ( dataset.metadata.snps_name,dataset.metadata.db_name) - dataset.blurb = "GMAP SNPindex %s on %s\n" % ( dataset.metadata.snps_name,dataset.metadata.db_name) - else: - dataset.peek = 'file does not exist' - dataset.blurb = 'file purged from disk' - def display_peek( self, dataset ): - try: - return dataset.peek - except: - return "GMAP SNP index" - - def sniff( self, filename ): - return False - def set_meta( self, dataset, overwrite = True, **kwd ): - """ - Expecting: - extra_files_path/snp_name.iit - extra_files_path/db_name/db_name.ref1[2345]1[2345]3offsetscomp.snp_name - extra_files_path/db_name/db_name.ref1[2345]1[2345]3positions.snp_name - extra_files_path/db_name/db_name.ref1[2345]1[2345]3gammaptrs.snp_name - """ - log.info( "GmapSnpIndex set_meta %s %s" % (dataset,dataset.extra_files_path)) - pat = '(.*)\.(ref((\d\d)(\d\d))?3positions)\.(.+)?' - efp = dataset.extra_files_path - flist = os.listdir(efp) - for i,fname in enumerate(flist): - m = re.match(pat,fname) - if m: - assert len(m.groups()) == 6 - dataset.metadata.db_name = m.groups()[0] - dataset.metadata.snps_name = m.groups()[-1] - - - - -class IntervalIndexTree( Text ): - """ - A GMAP Interval Index Tree Map - created by iit_store - (/path/to/map)/(mapname).iit - """ - file_ext = 'iit' - is_binary = True - -class SpliceSitesIntervalIndexTree( IntervalIndexTree ): - """ - A GMAP Interval Index Tree Map - created by iit_store - """ - file_ext = 'splicesites.iit' - -class IntronsIntervalIndexTree( IntervalIndexTree ): - """ - A GMAP Interval Index Tree Map - created by iit_store - """ - file_ext = 'introns.iit' - -class SNPsIntervalIndexTree( IntervalIndexTree ): - """ - A GMAP Interval Index Tree Map - created by iit_store - """ - file_ext = 'snps.iit' - -class TallyIntervalIndexTree( IntervalIndexTree ): - """ - A GMAP Interval Index Tree Map - created by iit_store - """ - file_ext = 'tally.iit' - -class IntervalAnnotation( Text ): - """ - Class describing a GMAP Interval format: - >label coords optional_tag - optional_annotation (which may be zero, one, or multiple lines) - The coords should be of the form: - chr:position - chr:startposition..endposition - """ - file_ext = 'gmap_annotation' - """Add metadata elements""" - MetadataElement( name="annotations", default=0, desc="Number of interval annotations", readonly=True, optional=True, visible=False, no_value=0 ) - - def set_meta( self, dataset, **kwd ): - """ - Set the number of annotations and the number of data lines in dataset. - """ - data_lines = 0 - annotations = 0 - for line in file( dataset.file_name ): - line = line.strip() - if line and line.startswith( '>' ): - annotations += 1 - data_lines +=1 - else: - data_lines += 1 - dataset.metadata.data_lines = data_lines - dataset.metadata.annotations = annotations - def set_peek( self, dataset, is_multi_byte=False ): - if not dataset.dataset.purged: - dataset.peek = data.get_file_peek( dataset.file_name, is_multi_byte=is_multi_byte ) - if dataset.metadata.annotations: - dataset.blurb = "%s annotations" % util.commaify( str( dataset.metadata.annotations ) ) - else: - dataset.blurb = data.nice_size( dataset.get_size() ) - else: - dataset.peek = 'file does not exist' - dataset.blurb = 'file purged from disk' - - def sniff( self, filename ): - """ - Determines whether the file is a gmap annotation file - Format: - >label coords optional_tag - optional_annotation (which may be zero, one, or multiple lines) - For example, the label may be an EST accession, with the coords - representing its genomic position. Labels may be duplicated if - necessary. - The coords should be of the form - chr:position - chr:startposition..endposition - The term "chr:position" is equivalent to "chr:position..position". If - you want to indicate that the interval is on the minus strand or - reverse direction, then <endposition> may be less than <startposition>. - """ - try: - pat = '>(\S+)\s((\S+):(\d+)(\.\.(\d+))?(\s.(.+))?$' #>label chr:position[..endposition][ optional_tag] - fh = open( filename ) - count = 0 - while True and count < 10: - line = fh.readline() - if not line: - break #EOF - line = line.strip() - if line: #first non-empty line - if line.startswith( '>' ): - count += 1 - if re.match(pat,line) == None: # Failed to match - return False - finally: - fh.close() - return False - -class SpliceSiteAnnotation(IntervalAnnotation): - file_ext = 'gmap_splicesites' - """ - Example: - >NM_004448.ERBB2.exon1 17:35110090..35110091 donor 6678 - >NM_004448.ERBB2.exon2 17:35116768..35116769 acceptor 6678 - >NM_004448.ERBB2.exon2 17:35116920..35116921 donor 1179 - >NM_004448.ERBB2.exon3 17:35118099..35118100 acceptor 1179 - >NM_004449.ERG.exon1 21:38955452..38955451 donor 783 - >NM_004449.ERG.exon2 21:38878740..38878739 acceptor 783 - >NM_004449.ERG.exon2 21:38878638..38878637 donor 360 - >NM_004449.ERG.exon3 21:38869542..38869541 acceptor 360 - Each line must start with a ">" character, then be followed by an - identifier, which may have duplicates and can have any format, with - the gene name or exon number shown here only as a suggestion. Then - there should be the chromosomal coordinates which straddle the - exon-intron boundary, so one coordinate is on the exon and one is on - the intron. (Coordinates are all 1-based, so the first character of a - chromosome is number 1.) Finally, there should be the splice type: - "donor" or "acceptor". You may optionally store the intron distance - at the end. GSNAP can use this intron distance, if it is longer than - its value for --localsplicedist, to look for long introns at that - splice site. The same splice site may have different intron distances - in the database; GSNAP will use the longest intron distance reported - in searching for long introns. - """ - def sniff( self, filename ): # TODO - """ - Determines whether the file is a gmap splice site annotation file - """ - try: - pat = '>(\S+\.intron\d+)\s((\S+):(\d+)\.\.(\d+))\s(donor|acceptor)(\s(\d+))?$' #>label chr:position..position donor|acceptor[ intron_dist] - fh = open( filename ) - count = 0 - while True and count < 10: - line = fh.readline() - if not line: - break #EOF - line = line.strip() - if line: #first non-empty line - count += 1 - if re.match(pat,line) == None: # Failed to match - return False - finally: - fh.close() - return False - -class IntronAnnotation(IntervalAnnotation): - file_ext = 'gmap_introns' - """ - Example: - >NM_004448.ERBB2.intron1 17:35110090..35116769 - >NM_004448.ERBB2.intron2 17:35116920..35118100 - >NM_004449.ERG.intron1 21:38955452..38878739 - >NM_004449.ERG.intron2 21:38878638..38869541 - The coordinates are 1-based, and specify the exon coordinates - surrounding the intron, with the first coordinate being from the donor - exon and the second one being from the acceptor exon. - """ - def sniff( self, filename ): # TODO - """ - Determines whether the file is a gmap Intron annotation file - """ - try: - pat = '>(\S+\.intron\d+)\s((\S+):(\d+)\.\.(\d+)(\s(.)+)?$' #>label chr:position - fh = open( filename ) - count = 0 - while True and count < 10: - line = fh.readline() - if not line: - break #EOF - line = line.strip() - if line: #first non-empty line - count += 1 - if re.match(pat,line) == None: # Failed to match - return False - finally: - fh.close() - return False - -class SNPAnnotation(IntervalAnnotation): - file_ext = 'gmap_snps' - """ - Example: - >rs62211261 21:14379270 CG - >rs62211262 21:14379281 AT - >rs62211263 21:14379298 WN - Each line must start with a ">" character, then be followed by an - identifier (which may have duplicates). Then there should be the - chromosomal coordinate of the SNP. (Coordinates are all 1-based, so - the first character of a chromosome is number 1.) Finally, there - should be the two possible alleles. (Previous versions required that - these be in alphabetical order: "AC", "AG", "AT", "CG", "CT", or "GT", - but that is no longer a requirement.) These alleles must correspond - to the possible nucleotides on the plus strand of the genome. If the - one of these two letters does not match the allele in the reference - sequence, that SNP will be ignored in subsequent processing as a - probable error. - - GSNAP also supports the idea of a wildcard SNP. A wildcard SNP allows - all nucleotides to match at that position, not just a given reference - and alternate allele. It is essentially as if an "N" were recorded at - that genomic location, although the index files still keep track of - the reference allele. To indicate that a position has a wildcard SNP, - you can indicate the genotype as "WN", where "W" is the reference - allele. Another indication of a wildcard SNP is to provide two - separate lines at that position with the genotypes "WX" and "WY", - where "W" is the reference allele and "X" and "Y" are two different - alternate alleles. - """ - def sniff( self, filename ): - """ - Determines whether the file is a gmap SNP annotation file - """ - try: - pat = '>(\S+)\s((\S+):(\d+)\s([TACGW][TACGN])$' #>label chr:position ATCG - fh = open( filename ) - count = 0 - while True and count < 10: - line = fh.readline() - if not line: - break #EOF - line = line.strip() - if line: #first non-empty line - count += 1 - if re.match(pat,line) == None: # Failed to match - return False - finally: - fh.close() - return False - - -class TallyAnnotation(IntervalAnnotation): - file_ext = 'gsnap_tally' - """ - Output produced by gsnap_tally - Example: - >144 chr20:57268791..57268935 - G0 - A1(1@7|1Q-3) - A2(1@36,1@1|1Q2,1Q-8) - C2 0.889,0.912,0.889,0.889,0.933,0.912,0.912,0.889,0.889,0.889 -2.66,-2.89,-2.66,-2.66,-3.16,-2.89,-2.89,-2.66,-2.66,-2.66 - C1 T1 0.888,0.9,0.888,0.9,0.913,0.9,0.911,0.888,0.9,0.913 -2.66,-2.78,-2.66,-2.78,-2.91,-2.78,-2.89,-2.66,-2.78,-2.91 - """ - def sniff( self, filename ): # TODO - """ - Determines whether the file is a gmap splice site annotation file - """ - try: - pat = '^>(\d+)\s((\S+):(\d+)\.\.(\d+))$' #>total chr:position..position - pat2 = '^[GATCN]\d.*$' #BaseCountDeatails - fh = open( filename ) - count = 0 - while True and count < 10: - line = fh.readline() - if not line: - break #EOF - line = line.strip() - if line: #first non-empty line - count += 1 - if re.match(pat,line) == None and re.match(pat2,line) == None: # Failed to match - return False - finally: - fh.close() - return False - -class GsnapResult( Text ): - """ - The default output format for gsnap. Can be used as input for gsnap_tally. - """ - file_ext = 'gsnap' - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib/galaxy/datatypes/gmap.py Tue Nov 08 13:22:34 2011 -0600 @@ -0,0 +1,472 @@ +""" +GMAP indexes +""" +import logging +import os,os.path,re +import data +from data import Text +from galaxy import util +from metadata import MetadataElement + +log = logging.getLogger(__name__) + +class GmapDB( Text ): + """ + A GMAP DB for indexes + """ + MetadataElement( name="db_name", desc="The db name for this index set", default='unknown', set_in_upload=True, readonly=True ) + MetadataElement( name="basesize", default="12", desc="The basesize for offsetscomp", visible=True, readonly=True ) + MetadataElement( name="kmers", default=[''], desc="The kmer sizes for indexes", visible=True, no_value=[''], readonly=True ) + MetadataElement( name="map_dir", desc="The maps directory", default='unknown', set_in_upload=True, readonly=True ) + MetadataElement( name="maps", default=[''], desc="The names of maps stored for this gmap gmapdb", visible=True, no_value=[''], readonly=True ) + MetadataElement( name="snps", default=[''], desc="The names of SNP indexes stored for this gmapdb", visible=True, no_value=[''], readonly=True ) + MetadataElement( name="cmet", default=False, desc="Has a cmet index", visible=True, readonly=True ) + MetadataElement( name="atoi", default=False, desc="Has a atoi index", visible=True, readonly=True ) + + file_ext = 'gmapdb' + is_binary = True + composite_type = 'auto_primary_file' + allow_datatype_change = False + + def generate_primary_file( self, dataset = None ): + """ + This is called only at upload to write the html file + cannot rename the datasets here - they come with the default unfortunately + """ + return '<html><head></head><body>AutoGenerated Primary File for Composite Dataset</body></html>' + + def regenerate_primary_file(self,dataset): + """ + cannot do this until we are setting metadata + """ + bn = dataset.metadata.db_name + log.info( "GmapDB regenerate_primary_file %s" % (bn)) + rval = ['<html><head><title>GMAPDB %s</title></head><p/><H3>GMAPDB %s</H3><p/>cmet %s<br>atoi %s<H4>Maps:</H4><ul>' % (bn,bn,dataset.metadata.cmet,dataset.metadata.atoi)] + for i,name in enumerate(dataset.metadata.maps): + rval.append( '<li>%s' % name) + rval.append( '</ul></html>' ) + f = file(dataset.file_name,'w') + f.write("\n".join( rval )) + f.write('\n') + f.close() + + def set_peek( self, dataset, is_multi_byte=False ): + log.info( "GmapDB set_peek %s" % (dataset)) + if not dataset.dataset.purged: + dataset.peek = "GMAPDB index %s\n cmet %s\n atoi %s\n maps %s" % ( dataset.metadata.db_name,dataset.metadata.cmet,dataset.metadata.atoi,dataset.metadata.maps ) + dataset.blurb = "GMAPDB %s" % ( dataset.metadata.db_name ) + else: + dataset.peek = 'file does not exist' + dataset.blurb = 'file purged from disk' + def display_peek( self, dataset ): + try: + return dataset.peek + except: + return "GMAP index file" + + def sniff( self, filename ): + return False + def set_meta( self, dataset, overwrite = True, **kwd ): + """ + Expecting: + extra_files_path/<db_name>/db_name>.ref<basesize><kmer>3<index> + extra_files_path/db_name/db_name.ref1[2345]1[2345]3offsetscomp + extra_files_path/db_name/db_name.ref1[2345]1[2345]3positions + extra_files_path/db_name/db_name.ref1[2345]1[2345]3gammaptrs + index maps: + extra_files_path/db_name/db_name.maps/*.iit + """ + log.info( "GmapDB set_meta %s %s" % (dataset,dataset.extra_files_path)) + pat = '(.*)\.((ref)|(met)[atgc][atgc]|(a2i)[atgc][atgc])((\d\d)(\d\d))?3positions(\.(.+))?' + efp = dataset.extra_files_path + flist = os.listdir(efp) + for i,fname in enumerate(flist): + log.info( "GmapDB set_meta %s %s" % (i,fname)) + fpath = os.path.join(efp,fname) + if os.path.isdir(fpath): + ilist = os.listdir(fpath) + kmers = {'':'default'} # HACK '' empty key added so user has default choice when selecting kmer from metadata + for j,iname in enumerate(ilist): + log.info( "GmapDB set_meta file %s %s" % (j,iname)) + ipath = os.path.join(fpath,iname) + if os.path.isdir(ipath): # find maps + dataset.metadata.map_dir = iname + for mapfile in os.listdir(ipath): + mapname = mapfile.replace('.iit','') + log.info( "GmapDB set_meta map %s %s" % (mapname,mapfile)) + dataset.metadata.maps.append(mapname) + else: + m = re.match(pat,iname) + if m: + log.info( "GmapDB set_meta m %s %s " % (iname, m)) + assert len(m.groups()) == 10 + dataset.metadata.db_name = fname + if m.groups()[2] == 'ref': + if m.groups()[-1] != None: + dataset.metadata.snps.append(m.groups()[-1]) + else: + if m.groups()[-3] != None: + k = int(m.groups()[-3]) + kmers[k] = k + if m.groups()[-4] != None: + dataset.metadata.basesize = int( m.groups()[-4]) + elif m.groups()[3] == 'met': + dataset.metadata.cmet = True + elif m.groups()[4] == 'a2i': + dataset.metadata.atoi = True + dataset.metadata.kmers = kmers.keys() + +class GmapSnpIndex( Text ): + """ + A GMAP SNP index created by snpindex + """ + MetadataElement( name="db_name", desc="The db name for this index set", default='unknown', set_in_upload=True, readonly=True ) + MetadataElement( name="snps_name", default='snps', desc="The name of SNP index", visible=True, no_value='', readonly=True ) + + file_ext = 'gmapsnpindex' + is_binary = True + composite_type = 'auto_primary_file' + allow_datatype_change = False + + def generate_primary_file( self, dataset = None ): + """ + This is called only at upload to write the html file + cannot rename the datasets here - they come with the default unfortunately + """ + return '<html><head></head><body>AutoGenerated Primary File for Composite Dataset</body></html>' + + def regenerate_primary_file(self,dataset): + """ + cannot do this until we are setting metadata + """ + bn = dataset.metadata.db_name + log.info( "GmapDB regenerate_primary_file %s" % (bn)) + rval = ['<html><head><title>GMAPDB %s</title></head><p/><H3>GMAPDB %s</H3><p/>cmet %s<br>atoi %s<H4>Maps:</H4><ul>' % (bn,bn,dataset.metadata.cmet,dataset.metadata.atoi)] + for i,name in enumerate(dataset.metadata.maps): + rval.append( '<li>%s' % name) + rval.append( '</ul></html>' ) + f = file(dataset.file_name,'w') + f.write("\n".join( rval )) + f.write('\n') + f.close() + def set_peek( self, dataset, is_multi_byte=False ): + log.info( "GmapSnpIndex set_peek %s" % (dataset)) + if not dataset.dataset.purged: + dataset.peek = "GMAP SNPindex %s on %s\n" % ( dataset.metadata.snps_name,dataset.metadata.db_name) + dataset.blurb = "GMAP SNPindex %s on %s\n" % ( dataset.metadata.snps_name,dataset.metadata.db_name) + else: + dataset.peek = 'file does not exist' + dataset.blurb = 'file purged from disk' + def display_peek( self, dataset ): + try: + return dataset.peek + except: + return "GMAP SNP index" + + def sniff( self, filename ): + return False + def set_meta( self, dataset, overwrite = True, **kwd ): + """ + Expecting: + extra_files_path/snp_name.iit + extra_files_path/db_name/db_name.ref1[2345]1[2345]3offsetscomp.snp_name + extra_files_path/db_name/db_name.ref1[2345]1[2345]3positions.snp_name + extra_files_path/db_name/db_name.ref1[2345]1[2345]3gammaptrs.snp_name + """ + log.info( "GmapSnpIndex set_meta %s %s" % (dataset,dataset.extra_files_path)) + pat = '(.*)\.(ref((\d\d)(\d\d))?3positions)\.(.+)?' + efp = dataset.extra_files_path + flist = os.listdir(efp) + for i,fname in enumerate(flist): + m = re.match(pat,fname) + if m: + assert len(m.groups()) == 6 + dataset.metadata.db_name = m.groups()[0] + dataset.metadata.snps_name = m.groups()[-1] + + + + +class IntervalIndexTree( Text ): + """ + A GMAP Interval Index Tree Map + created by iit_store + (/path/to/map)/(mapname).iit + """ + file_ext = 'iit' + is_binary = True + +class SpliceSitesIntervalIndexTree( IntervalIndexTree ): + """ + A GMAP Interval Index Tree Map + created by iit_store + """ + file_ext = 'splicesites.iit' + +class IntronsIntervalIndexTree( IntervalIndexTree ): + """ + A GMAP Interval Index Tree Map + created by iit_store + """ + file_ext = 'introns.iit' + +class SNPsIntervalIndexTree( IntervalIndexTree ): + """ + A GMAP Interval Index Tree Map + created by iit_store + """ + file_ext = 'snps.iit' + +class TallyIntervalIndexTree( IntervalIndexTree ): + """ + A GMAP Interval Index Tree Map + created by iit_store + """ + file_ext = 'tally.iit' + +class IntervalAnnotation( Text ): + """ + Class describing a GMAP Interval format: + >label coords optional_tag + optional_annotation (which may be zero, one, or multiple lines) + The coords should be of the form: + chr:position + chr:startposition..endposition + """ + file_ext = 'gmap_annotation' + """Add metadata elements""" + MetadataElement( name="annotations", default=0, desc="Number of interval annotations", readonly=True, optional=True, visible=False, no_value=0 ) + + def set_meta( self, dataset, **kwd ): + """ + Set the number of annotations and the number of data lines in dataset. + """ + data_lines = 0 + annotations = 0 + for line in file( dataset.file_name ): + line = line.strip() + if line and line.startswith( '>' ): + annotations += 1 + data_lines +=1 + else: + data_lines += 1 + dataset.metadata.data_lines = data_lines + dataset.metadata.annotations = annotations + def set_peek( self, dataset, is_multi_byte=False ): + if not dataset.dataset.purged: + dataset.peek = data.get_file_peek( dataset.file_name, is_multi_byte=is_multi_byte ) + if dataset.metadata.annotations: + dataset.blurb = "%s annotations" % util.commaify( str( dataset.metadata.annotations ) ) + else: + dataset.blurb = data.nice_size( dataset.get_size() ) + else: + dataset.peek = 'file does not exist' + dataset.blurb = 'file purged from disk' + + def sniff( self, filename ): + """ + Determines whether the file is a gmap annotation file + Format: + >label coords optional_tag + optional_annotation (which may be zero, one, or multiple lines) + For example, the label may be an EST accession, with the coords + representing its genomic position. Labels may be duplicated if + necessary. + The coords should be of the form + chr:position + chr:startposition..endposition + The term "chr:position" is equivalent to "chr:position..position". If + you want to indicate that the interval is on the minus strand or + reverse direction, then <endposition> may be less than <startposition>. + """ + try: + pat = '>(\S+)\s((\S+):(\d+)(\.\.(\d+))?(\s.(.+))?$' #>label chr:position[..endposition][ optional_tag] + fh = open( filename ) + count = 0 + while True and count < 10: + line = fh.readline() + if not line: + break #EOF + line = line.strip() + if line: #first non-empty line + if line.startswith( '>' ): + count += 1 + if re.match(pat,line) == None: # Failed to match + return False + finally: + fh.close() + return False + +class SpliceSiteAnnotation(IntervalAnnotation): + file_ext = 'gmap_splicesites' + """ + Example: + >NM_004448.ERBB2.exon1 17:35110090..35110091 donor 6678 + >NM_004448.ERBB2.exon2 17:35116768..35116769 acceptor 6678 + >NM_004448.ERBB2.exon2 17:35116920..35116921 donor 1179 + >NM_004448.ERBB2.exon3 17:35118099..35118100 acceptor 1179 + >NM_004449.ERG.exon1 21:38955452..38955451 donor 783 + >NM_004449.ERG.exon2 21:38878740..38878739 acceptor 783 + >NM_004449.ERG.exon2 21:38878638..38878637 donor 360 + >NM_004449.ERG.exon3 21:38869542..38869541 acceptor 360 + Each line must start with a ">" character, then be followed by an + identifier, which may have duplicates and can have any format, with + the gene name or exon number shown here only as a suggestion. Then + there should be the chromosomal coordinates which straddle the + exon-intron boundary, so one coordinate is on the exon and one is on + the intron. (Coordinates are all 1-based, so the first character of a + chromosome is number 1.) Finally, there should be the splice type: + "donor" or "acceptor". You may optionally store the intron distance + at the end. GSNAP can use this intron distance, if it is longer than + its value for --localsplicedist, to look for long introns at that + splice site. The same splice site may have different intron distances + in the database; GSNAP will use the longest intron distance reported + in searching for long introns. + """ + def sniff( self, filename ): # TODO + """ + Determines whether the file is a gmap splice site annotation file + """ + try: + pat = '>(\S+\.intron\d+)\s((\S+):(\d+)\.\.(\d+))\s(donor|acceptor)(\s(\d+))?$' #>label chr:position..position donor|acceptor[ intron_dist] + fh = open( filename ) + count = 0 + while True and count < 10: + line = fh.readline() + if not line: + break #EOF + line = line.strip() + if line: #first non-empty line + count += 1 + if re.match(pat,line) == None: # Failed to match + return False + finally: + fh.close() + return False + +class IntronAnnotation(IntervalAnnotation): + file_ext = 'gmap_introns' + """ + Example: + >NM_004448.ERBB2.intron1 17:35110090..35116769 + >NM_004448.ERBB2.intron2 17:35116920..35118100 + >NM_004449.ERG.intron1 21:38955452..38878739 + >NM_004449.ERG.intron2 21:38878638..38869541 + The coordinates are 1-based, and specify the exon coordinates + surrounding the intron, with the first coordinate being from the donor + exon and the second one being from the acceptor exon. + """ + def sniff( self, filename ): # TODO + """ + Determines whether the file is a gmap Intron annotation file + """ + try: + pat = '>(\S+\.intron\d+)\s((\S+):(\d+)\.\.(\d+)(\s(.)+)?$' #>label chr:position + fh = open( filename ) + count = 0 + while True and count < 10: + line = fh.readline() + if not line: + break #EOF + line = line.strip() + if line: #first non-empty line + count += 1 + if re.match(pat,line) == None: # Failed to match + return False + finally: + fh.close() + return False + +class SNPAnnotation(IntervalAnnotation): + file_ext = 'gmap_snps' + """ + Example: + >rs62211261 21:14379270 CG + >rs62211262 21:14379281 AT + >rs62211263 21:14379298 WN + Each line must start with a ">" character, then be followed by an + identifier (which may have duplicates). Then there should be the + chromosomal coordinate of the SNP. (Coordinates are all 1-based, so + the first character of a chromosome is number 1.) Finally, there + should be the two possible alleles. (Previous versions required that + these be in alphabetical order: "AC", "AG", "AT", "CG", "CT", or "GT", + but that is no longer a requirement.) These alleles must correspond + to the possible nucleotides on the plus strand of the genome. If the + one of these two letters does not match the allele in the reference + sequence, that SNP will be ignored in subsequent processing as a + probable error. + + GSNAP also supports the idea of a wildcard SNP. A wildcard SNP allows + all nucleotides to match at that position, not just a given reference + and alternate allele. It is essentially as if an "N" were recorded at + that genomic location, although the index files still keep track of + the reference allele. To indicate that a position has a wildcard SNP, + you can indicate the genotype as "WN", where "W" is the reference + allele. Another indication of a wildcard SNP is to provide two + separate lines at that position with the genotypes "WX" and "WY", + where "W" is the reference allele and "X" and "Y" are two different + alternate alleles. + """ + def sniff( self, filename ): + """ + Determines whether the file is a gmap SNP annotation file + """ + try: + pat = '>(\S+)\s((\S+):(\d+)\s([TACGW][TACGN])$' #>label chr:position ATCG + fh = open( filename ) + count = 0 + while True and count < 10: + line = fh.readline() + if not line: + break #EOF + line = line.strip() + if line: #first non-empty line + count += 1 + if re.match(pat,line) == None: # Failed to match + return False + finally: + fh.close() + return False + + +class TallyAnnotation(IntervalAnnotation): + file_ext = 'gsnap_tally' + """ + Output produced by gsnap_tally + Example: + >144 chr20:57268791..57268935 + G0 + A1(1@7|1Q-3) + A2(1@36,1@1|1Q2,1Q-8) + C2 0.889,0.912,0.889,0.889,0.933,0.912,0.912,0.889,0.889,0.889 -2.66,-2.89,-2.66,-2.66,-3.16,-2.89,-2.89,-2.66,-2.66,-2.66 + C1 T1 0.888,0.9,0.888,0.9,0.913,0.9,0.911,0.888,0.9,0.913 -2.66,-2.78,-2.66,-2.78,-2.91,-2.78,-2.89,-2.66,-2.78,-2.91 + """ + def sniff( self, filename ): # TODO + """ + Determines whether the file is a gmap splice site annotation file + """ + try: + pat = '^>(\d+)\s((\S+):(\d+)\.\.(\d+))$' #>total chr:position..position + pat2 = '^[GATCN]\d.*$' #BaseCountDeatails + fh = open( filename ) + count = 0 + while True and count < 10: + line = fh.readline() + if not line: + break #EOF + line = line.strip() + if line: #first non-empty line + count += 1 + if re.match(pat,line) == None and re.match(pat2,line) == None: # Failed to match + return False + finally: + fh.close() + return False + +class GsnapResult( Text ): + """ + The default output format for gsnap. Can be used as input for gsnap_tally. + """ + file_ext = 'gsnap' + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool-data/gmap_indices.loc.sample Tue Nov 08 13:22:34 2011 -0600 @@ -0,0 +1,10 @@ +#This is a sample file distributed with Galaxy that enables tools +#to use a directory of GMAPDB indexed sequences data files. You will need +#to create these data files using gmap_build and then create a gmap_indices.loc file +#similar to this one (store it in this directory) that points to +#the directories in which those files are stored. The gmap_indices.loc +#file has this format (white space characters are TAB characters): +# +#<unique_build_id> <dbkey> <display_name> <kmers> <map,map> <snp,snp> <file_base_path> +#hg18 hg18 hg18 (cmet atoi) 12,13,14,15 splicesites,introns snps /depot/data2/galaxy/gmap/hg18 +#hg19 hg19 hg19 (cmet atoi) 12,13,14,15 splicesites,introns,snps snps,dbsnp /depot/data2/galaxy/gmap/hg19