# HG changeset patch # User Jim Johnson # Date 1319215135 18000 # Node ID 52da588232b0fac11b295514b4e6befaafd8ce4b # Parent 30d42bb409b8808190b044c11d9957c633f430e4 Add datatypes for maps and snpindex, add iit_store and snpindex tools, update GMAP and GSNAP to use these. diff -r 30d42bb409b8 -r 52da588232b0 gmap/README --- a/gmap/README Tue Oct 18 11:51:15 2011 -0500 +++ b/gmap/README Fri Oct 21 11:38:55 2011 -0500 @@ -9,9 +9,16 @@ add to datatypes_conf.xml - Tools: GMAP_Build - create a GmapDB set of index files for a reference sequence and optional set of annotations @@ -19,28 +26,22 @@ GSNAP - align sequences to a reference and detect splicing Add to tool_conf.xml ( probably in the "NGS: Mapping" section ) - + + + + +Admin built cached gmapdb indexes defined in tool-data/gmap_indices.loc TODO: Add classes to gmap.py - IntervalIndexTree - datatype for the iit_store result - IntervalAnnotation - generic class for gmap annotation formats - SpliceSiteAnnotation - class for gmap splicesite annotation format - IntronAnnotation - class for gmap introns annotation format - SNPAnnotation - class for gmap snp annotation format - GmapIndex - generic class for a gmap index - SnpIndex - an index created by snpindex CmetIndex - an index created by cmetindex AtoiIndex - an index created by atoiindex Possibly add Tools: - iit_store - creates a mapping index .iit from: gtf,gff3,gmapannotation, UCSC refGene table - iit_get - retrieves from .iit store get_genome - retrieves from a gmapdb - snpindex - create a SNPIndex cmetindex - create methylcytosine index atoiindex - create A-to-I RNA editing index diff -r 30d42bb409b8 -r 52da588232b0 gmap/gmap.xml --- a/gmap/gmap.xml Tue Oct 18 11:51:15 2011 -0500 +++ b/gmap/gmap.xml Fri Oct 21 11:38:55 2011 -0500 @@ -1,7 +1,13 @@ Genomic Mapping and Alignment Program for mRNA and EST sequences - gmap + gmap + + gmapdb + gmap_annotation + gmap_splicesites + gmap_introns + gmap_snps gmap --version @@ -113,8 +119,20 @@ #end if + + + + + + + + + + + + - + @@ -155,8 +173,6 @@ - @@ -174,18 +190,89 @@ - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + @@ -246,92 +333,20 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + - + (split_output == False) @@ -339,10 +354,10 @@ - + + @@ -352,6 +367,10 @@ + + + + @@ -361,6 +380,10 @@ + + + + @@ -370,6 +393,10 @@ + + + + @@ -379,6 +406,10 @@ + + + + diff -r 30d42bb409b8 -r 52da588232b0 gmap/gmap_build.xml --- a/gmap/gmap_build.xml Tue Oct 18 11:51:15 2011 -0500 +++ b/gmap/gmap_build.xml Fri Oct 21 11:38:55 2011 -0500 @@ -1,16 +1,23 @@ - a GMAP DB Index + a database genome index for GMAP and GSNAP gmap_build + + gmapdb + gmap_snps gmap --version /bin/bash $shscript 2>1 1> $output - - - + + + + + + + @@ -79,7 +86,6 @@ - @@ -89,6 +95,11 @@ #set $gt = chr(62) #set $lt = chr(60) #set $ad = chr(38) +## #set $ref_files = '' +## #for $i in $inputs: + ## #set $ref_files = $ref_files $i.input +## #end for +## echo $ref_files #import os.path #set $gmapdb = $output.extra_files_path #set $mapsdir = $os.path.join($os.path.join($gmapdb,str($refname)), str($refname) + '.maps') @@ -96,7 +107,7 @@ ## export GMAPDB required for cmetindex and atoiindex export GMAPDB=$gmapdb #for $k in $kmer.__str__.split(','): -gmap_build -D $gmapdb -d $refname -s numeric-alpha -k $k $input +gmap_build -D $gmapdb -d $refname -s numeric-alpha -k $k #for i in $inputs# ${i.input}#end for# #end for get-genome -D $gmapdb -d '?' | sed 's/^Available .*/gmap db: /' echo "kmers: " $kmer @@ -116,23 +127,26 @@ cat $splicesite.gff3Genes | gff3_introns | iit_store -o $os.path.join($mapsdir,'introns') #end if #end if +#if $dbsnp.snp_source != 'none' and $dbsnp.snps.__str__ != 'None': #if $dbsnp.snp_source == 'snpTable': -#if $dbsnp.snps.__str__ != 'None': #if $dbsnp.snpsex.__str__ != 'None': cat $dbsnp.snps | dbsnp_iit -w $dbsnp.weight -e $dbsnp.snpsex | iit_store -o $os.path.join($mapsdir,'snps') #else: cat $dbsnp.snps | dbsnp_iit -w $dbsnp.weight | iit_store -o $os.path.join($mapsdir,'snps') #end if +#else: +cat $dbsnp.snps | iit_store -o $os.path.join($mapsdir,'snps') +#end if snpindex -d $refname -v snps -#end if +echo "snpindex" -d $refname -v snps #end if #if $cmetindex.__str__ == 'yes': cmetindex -d $refname -echo "cmetindex" +echo "cmetindex" -d $refname #end if #if $atoiindex.__str__ == 'yes': atoiindex -d $refname -echo "atoiindex" +echo "atoiindex" -d $refname #end if get-genome -D $gmapdb -d $refname -m '?' | sed 's/^Available maps .*/maps: /' @@ -146,7 +160,7 @@ **GMAP Build** -GMAP Build creates an index of a genomic sequence for alignments using GMAP_ (Genomic Mapping and Alignment Program for mRNA and EST sequences) and GSNAP_ (Genomic Short-read Nucleotide Alignment Program). +GMAP Build creates an index of a genomic sequence for mapping and alignment using GMAP_ (Genomic Mapping and Alignment Program for mRNA and EST sequences) and GSNAP_ (Genomic Short-read Nucleotide Alignment Program). (GMAP Build uses GMSP commands: gmap_build, iit_store, psl_splicesites, psl_introns, gtf_splicesites, gtf_introns, gff3_splicesites, gff3_introns, dbsnp_iit, snpindex, cmetindex, and atoiindex.) You will want to read the README_ diff -r 30d42bb409b8 -r 52da588232b0 gmap/gsnap.xml --- a/gmap/gsnap.xml Tue Oct 18 11:51:15 2011 -0500 +++ b/gmap/gsnap.xml Fri Oct 21 11:38:55 2011 -0500 @@ -2,31 +2,43 @@ Genomic Short-read Nucleotide Alignment Program gsnap + + gmapdb + gmapsnpindex + splicesites.iit + introns.iit gsnap --version #import os.path, re gsnap --nthreads="4" --ordered - #if $refGenomeSource.genomeSource == "history": - --gseg=$refGenomeSource.ownFile - #elif $refGenomeSource.genomeSource == "gmapdb": + #if $refGenomeSource.genomeSource == "gmapdb": #set $gmapdb = $os.listdir($refGenomeSource.gmapdb.extra_files_path)[0] - --dir=$refGenomeSource.gmapdb.extra_files_path --db=$gmapdb - #if $refGenomeSource.kmer != None and len($refGenomeSource.kmer.__str__) == 2: - --kmer=$refGenomeSource.kmer - #end if - #if $refGenomeSource.splicemap != None and len($refGenomeSource.splicemap.__str__) == 2: - --use-splices=$refGenomeSource.splicemap - #end if - #if $refGenomeSource.snpindex != None and len($refGenomeSource.snpindex.__str__) == 2: - --use-snps=$refGenomeSource.snpindex - #end if + --dir=$refGenomeSource.gmapdb.extra_files_path --db=$refGenomeSource.gmapdb.metadata.db_name #else: --dir=$os.path.dirname($refGenomeSource.gmapindex.value) --db=$os.path.basename($refGenomeSource.gmapindex.value) - #if $refGenomeSource.kmer != None and len($refGenomeSource.kmer.__str__) == 2: - --kmer=$refGenomeSource.kmer + #end if + #if $refGenomeSource.kmer != None and len($refGenomeSource.kmer.__str__) == 2: + --kmer=$refGenomeSource.kmer + #end if + #if $refGenomeSource.use_splicing.src == 'gmapdb': + #if $refGenomeSource.use_splicing.splicemap != None and len($refGenomeSource.use_splicing.splicemap.__str__) > 0: + -s $refGenomeSource.use_splicing.splicemap.value #end if + #elif $refGenomeSource.use_splicing.src == 'history': + #if $refGenomeSource.use_splicing.splicemap != None and len($refGenomeSource.use_splicing.splicemap.__str__) > 0: + -S $os.path.dirname($refGenomeSource.use_splicing.splicemap) -s $os.path.basename($refGenomeSource.use_splicing.splicemap) + #end if + #end if + #if $refGenomeSource.use_snps.src == 'gmapdb': + #if $refGenomeSource.use_snps.snpindex != None and len($refGenomeSource.use_snps.snpindex.__str__) > 0: + -v $refGenomeSource.use_snps.snpindex.value + #end if + #elif $refGenomeSource.use_snps.src == 'history': + #if $refGenomeSource.use_snps.snpindex != None and len($refGenomeSource.use_snps.snpindex.__str__) > 0: + -V $refGenomeSource.use_snps.snpindex.extra_files_path -v $refGenomeSource.use_snps.snpindex.metadata.snps_name + #end if #end if #if $mode.__str__ != '': --mode=$mode @@ -181,86 +193,9 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + @@ -308,6 +243,7 @@ + @@ -315,9 +251,147 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + @@ -352,8 +426,7 @@ - - + @@ -376,7 +449,7 @@ - + @@ -384,7 +457,7 @@ - + @@ -398,7 +471,7 @@ - + @@ -439,7 +512,7 @@ help="Separate outputs for: nomapping, halfmapping_uniq, halfmapping_mult, unpaired_uniq, unpaired_mult, paired_uniq, paired_mult, concordant_uniq, and concordant_mult results"/> - + (split_output == False) diff -r 30d42bb409b8 -r 52da588232b0 gmap/iit_store.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gmap/iit_store.xml Fri Oct 21 11:38:55 2011 -0500 @@ -0,0 +1,188 @@ + + Create a map store for known genes or SNPs + + iit_store + + gmap_annotation + gmap_snps + iit + splicesites.iit + introns.iit + snps.iit + + iit_store --version + /bin/bash $shscript 2> $log + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + (map['type'] == 'genes' and 'splicesites' in map['maps']) + + + (map['type'] == 'genes' and 'introns' in map['maps']) + + + (map['type'] == 'snps') + + + (map['type'] == 'gmap') + + + + +#!/bin/bash +#set $catcmd = 'gzcat -f' +#set $catcmd = 'cat' +#set $ds = chr(36) +#set $gt = chr(62) +#set $lt = chr(60) +#set $ad = chr(38) +#set $ep = chr(33) +#set $toerr = ''.join([$gt,$ad,'2']) +#import os.path +#if $map.type == 'genes': +if [ $ep -e $map.src.genes ]; then echo "$map.src.genes does not exist" $toerr; exit 1; fi +if [ $ep -s $map.src.genes ]; then echo "$map.src.genes is empty" $toerr; exit 2; fi + #if $map.src.src_format == 'refGeneTable': + #if 'splicesites' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | psl_splicesites -s $map.src.col_skip | iit_store -o $splicesites_iit + #end if + #if 'introns' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | psl_introns -s $map.src.col_skip | iit_store -o $introns_iit + #end if + #elif $map.src.src_format == 'gtf': + #if 'splicesites' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | gtf_splicesites | iit_store -o $splicesites_iit + #end if + #if 'introns' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | gtf_introns | iit_store -o $introns_iit + #end if + #elif $map.src.src_format == 'gff3': + #if 'splicesites' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | gff3_splicesites | iit_store -o $splicesites_iit + #end if + #if 'introns' in [ $map.maps.__str__ ]: + $catcmd $map.src.genes | gff3_introns | iit_store -o $introns_iit + #end if + #end if +#elif $map.type == 'snps': +if [ $ep -s $map.src.snps ]; then echo "$map.src.snps is empty" $toerr; exit 2; fi + #if $map.src.snpsex.__str__ != 'None': + $catcmd $map.src.snps | dbsnp_iit -w $map.src.weight -e $map.src.snpsex | iit_store -o $snps_iit + #else: + $catcmd $map.src.snps | dbsnp_iit -w $map.src.weight | iit_store -o $snps_iit + #end if +#else: + $catcmd $map.src.snps | iit_store -o $map_iit +#end if + + + + + + + + + +**iit_store** + +GMAP IIT creates an Interval Index Tree map of known splice sites, introns, or SNPs (it uses iit_store described in the GMAP documentation). The maps can be used in GMAP_ (Genomic Mapping and Alignment Program for mRNA and EST sequences) and GSNAP_ (Genomic Short-read Nucleotide Alignment Program). Maps are typically used for known splice sites, introns, or SNPs. + +You will want to read the README_ + +Publication_ citation: Thomas D. Wu, Colin K. Watanabe Bioinformatics 2005 21(9):1859-1875; doi:10.1093/bioinformatics/bti310 + +.. _GMAP: http://research-pub.gene.com/gmap/ +.. _GSNAP: http://research-pub.gene.com/gmap/ +.. _README: http://research-pub.gene.com/gmap/src/README +.. _Publication: http://bioinformatics.oxfordjournals.org/cgi/content/full/21/9/1859 + + +**inputs** + + + + diff -r 30d42bb409b8 -r 52da588232b0 gmap/lib/galaxy/datatypes/gmap.py --- a/gmap/lib/galaxy/datatypes/gmap.py Tue Oct 18 11:51:15 2011 -0500 +++ b/gmap/lib/galaxy/datatypes/gmap.py Fri Oct 21 11:38:55 2011 -0500 @@ -3,7 +3,9 @@ """ 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__) @@ -61,6 +63,7 @@ return dataset.peek except: return "GMAP index file" + def sniff( self, filename ): return False def set_meta( self, dataset, overwrite = True, **kwd ): @@ -113,57 +116,309 @@ dataset.metadata.atoi = True dataset.metadata.kmers = kmers.keys() -## class IntervalIndexTree( Text ): -## """ -## A GMAP Interval Index Tree Map -## created by iit_store -## (/path/to/map)/(mapname).iit -## """ -## MetadataElement( name="map_name", desc="The map name for this index set", default='unknown', set_in_upload=True, readonly=False ) -## file_ext = 'iit' -## is_binary = True -## composite_type = 'auto_primary_file' -## allow_datatype_change = False -## -## class IntervalAnnotation(data.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 = 'gmapannotation' -## -## class SpliceSiteAnnotation(IntervalAnnotation): -## file_ext = 'gmapsplicesites' -## """ -## 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 -## """ -## -## class IntronAnnotation(IntervalAnnotation): -## file_ext = 'gmapintrons' -## """ -## 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 -## """ -## -## class SNPAnnotation(IntervalAnnotation): -## file_ext = 'gmapsnps' -## """ -## Example: -## >rs62211261 21:14379270 CG -## >rs62211262 21:14379281 CG -## """ +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 'AutoGenerated Primary File for Composite Dataset' + + 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 = ['GMAPDB %s

GMAPDB %s

cmet %s
atoi %s

Maps:

    ' % (bn,bn,dataset.metadata.cmet,dataset.metadata.atoi)] + for i,name in enumerate(dataset.metadata.maps): + rval.append( '
  • %s' % name) + rval.append( '
' ) + 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 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 may be less than . + """ + 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 + + diff -r 30d42bb409b8 -r 52da588232b0 gmap/snpindex.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gmap/snpindex.xml Fri Oct 21 11:38:55 2011 -0500 @@ -0,0 +1,141 @@ + + build index files for known SNPs + + snpindex + + gmapsnpindex + gmapdb + gmap_snps + snps.iit + + snpindex --version + /bin/bash $shscript 2>1 1> $output + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +#!/bin/bash +#set $ds = chr(36) +#set $gt = chr(62) +#set $lt = chr(60) +#set $ad = chr(38) +#import os.path +#if $refGenomeSource.genomeSource == "gmapdb": +#set $gmapdb = $refGenomeSource.gmapdb.extra_files_path +#set $refname = $refGenomeSource.gmapdb.metadata.db_name +#else: +#set $gmapdb = $os.path.dirname($refGenomeSource.gmapindex.value) +$refname = $os.path.basename($refGenomeSource.gmapindex.value) +#end if +#set $gmapsnpdir = $output.extra_files_path +mkdir -p $gmapsnpdir +#set $snpsname = $snps_name.__str__ +#set $snpsiit = '.'.join([$snpsname,'iit']) +#set $pathsnps = $os.path.join($gmapsnpdir,$snpsname) +#set $pathsnpsiit = $os.path.join($gmapsnpdir,$snpsiit) +#if $dbsnp.snp_source != 'none' and $dbsnp.snps.__str__ != 'None': +#if $dbsnp.snp_source == 'snpTable': +#if $dbsnp.snpsex.__str__ != 'None': +cat $dbsnp.snps | dbsnp_iit -w $dbsnp.weight -e $dbsnp.snpsex | iit_store -o $pathsnps +#else: +cat $dbsnp.snps | dbsnp_iit -w $dbsnp.weight | iit_store -o $pathsnps +#end if +#elif $dbsnp.snp_source == 'snpFile': +cat $dbsnp.snps | iit_store -o $pathsnps +#elif $dbsnp.snp_source == 'snpIIT': +cat $dbsnp.snps > $pathsnpsiit +#end if +snpindex -D $gmapdb -d $refname -V $output.extra_files_path -v $snpsname $pathsnpsiit +echo snpindex -D $gmapdb -d $refname -V $output.extra_files_path -v $snpsname $pathsnpsiit +#end if + + + + + + + + + +**GMAP SNP Index** + +GMAP SNP Index (snpindex in the GMAP documentaion) creates an index for known SNPs allowing for SNP tolerant mapping and alignment when using GMAP_ (Genomic Mapping and Alignment Program for mRNA and EST sequences) and GSNAP_ (Genomic Short-read Nucleotide Alignment Program). + +You will want to read the README_ + +Publication_ citation: Thomas D. Wu, Colin K. Watanabe Bioinformatics 2005 21(9):1859-1875; doi:10.1093/bioinformatics/bti310 + +.. _GMAP: http://research-pub.gene.com/gmap/ +.. _GSNAP: http://research-pub.gene.com/gmap/ +.. _README: http://research-pub.gene.com/gmap/src/README +.. _Publication: http://bioinformatics.oxfordjournals.org/cgi/content/full/21/9/1859 + + + + +