Mercurial > repos > devteam > cuffmerge

comparison cuffmerge_wrapper.py @ 4:9cfc0b970015

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge heads.
author Dave Bouvier <dave@bx.psu.edu>
date Wed, 08 Jan 2014 09:16:08 -0500
parents 5b285b6e4ee3
children b6e3849293b1
comparison
equal deleted inserted replaced
3:aa209d311ece 4:9cfc0b970015
6 6
7 def stop_err( msg ): 7 def stop_err( msg ):
8 sys.stderr.write( '%s\n' % msg ) 8 sys.stderr.write( '%s\n' % msg )
9 sys.exit() 9 sys.exit()
10 10
11 # Copied from sam_to_bam.py:
12 def check_seq_file( dbkey, cached_seqs_pointer_file ):
13 seq_path = ''
14 for line in open( cached_seqs_pointer_file ):
15 line = line.rstrip( '\r\n' )
16 if line and not line.startswith( '#' ) and line.startswith( 'index' ):
17 fields = line.split( '\t' )
18 if len( fields ) < 3:
19 continue
20 if fields[1] == dbkey:
21 seq_path = fields[2].strip()
22 break
23 return seq_path
24
25 def __main__(): 11 def __main__():
26 #Parse Command Line 12 #Parse Command Line
27 parser = optparse.OptionParser() 13 parser = optparse.OptionParser()
28 parser.add_option( '-g', dest='ref_annotation', help='An optional "reference" annotation GTF. Each sample is matched against this file, and sample isoforms are tagged as overlapping, matching, or novel where appropriate. See the refmap and tmap output file descriptions below.' ) 14 parser.add_option( '-g', dest='ref_annotation', help='An optional "reference" annotation GTF. Each sample is matched against this file, and sample isoforms are tagged as overlapping, matching, or novel where appropriate. See the refmap and tmap output file descriptions below.' )
29 parser.add_option( '-s', dest='use_seq_data', action="store_true", help='Causes cuffmerge to look into for fasta files with the underlying genomic sequences (one file per contig) against which your reads were aligned for some optional classification functions. For example, Cufflinks transcripts consisting mostly of lower-case bases are classified as repeats. Note that <seq_dir> must contain one fasta file per reference chromosome, and each file must be named after the chromosome, and have a .fa or .fasta extension.') 15 parser.add_option( '-s', dest='use_seq_data', action="store_true", help='Causes cuffmerge to look into for fasta files with the underlying genomic sequences (one file per contig) against which your reads were aligned for some optional classification functions. For example, Cufflinks transcripts consisting mostly of lower-case bases are classified as repeats. Note that <seq_dir> must contain one fasta file per reference chromosome, and each file must be named after the chromosome, and have a .fa or .fasta extension.')
30 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' ) 16 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' )
31 17
32 18
33 # Wrapper / Galaxy options. 19 # Wrapper / Galaxy options.
34 parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' ) 20 parser.add_option( '', '--index', dest='index', help='The path of the reference genome' )
35 parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' )
36 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) 21 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' )
37 22
38 # Outputs. 23 # Outputs.
39 parser.add_option( '', '--merged-transcripts', dest='merged_transcripts' ) 24 parser.add_option( '', '--merged-transcripts', dest='merged_transcripts' )
40 25
59 except: 44 except:
60 sys.stdout.write( 'Could not determine Cuffmerge version\n' ) 45 sys.stdout.write( 'Could not determine Cuffmerge version\n' )
61 46
62 # Set/link to sequence file. 47 # Set/link to sequence file.
63 if options.use_seq_data: 48 if options.use_seq_data:
64 if options.ref_file != 'None': 49 if options.ref_file:
65 # Sequence data from history. 50 # Sequence data from history.
66 # Create symbolic link to ref_file so that index will be created in working directory. 51 # Create symbolic link to ref_file so that index will be created in working directory.
67 seq_path = "ref.fa" 52 seq_path = "ref.fa"
68 os.symlink( options.ref_file, seq_path ) 53 os.symlink( options.ref_file, seq_path )
69 else: 54 else:
70 # Sequence data from loc file. 55 if not os.path.exists( options.index ):
71 cached_seqs_pointer_file = os.path.join( options.index_dir, 'sam_fa_indices.loc' ) 56 stop_err( 'Reference genome %s not present, request it by reporting this error.' % options.index )
72 if not os.path.exists( cached_seqs_pointer_file ): 57 seq_path = options.index
73 stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
74 # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa,
75 # and the equCab2.fa file will contain fasta sequences.
76 seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
77 if seq_path == '':
78 stop_err( 'No sequence data found for dbkey %s, so sequence data cannot be used.' % options.dbkey )
79 58
80 # Build command. 59 # Build command.
81 60
82 # Base. 61 # Base.
83 cmd = "cuffmerge -o cm_output " 62 cmd = "cuffmerge -o cm_output "