annotate SMART/DiffExpAnal/tophat_parallel_unSQL.py @ 33:aa0420172fc6

Deleted selected files

1 
2 #!/usr/bin/env python
3 
4 import optparse, os, shutil, subprocess, sys, tempfile, fileinput, tarfile, glob
5 from commons.core.launcher.Launcher import Launcher
6 from commons.core.sql.TableJobAdaptatorFactory import TableJobAdaptatorFactory
7 from commons.core.utils.FileUtils import FileUtils
8 
9 def stop_err( msg ):
10     sys.stderr.write( "%s\n" % msg )
11     sys.exit()
12 
13 def toTar(tarFileName, accepted_hits_outputNames):
14     tfile = tarfile.open(tarFileName + ".tmp.tar", "w")
15     currentPath = os.getcwd()
16     os.chdir(dir)
17     for file in accepted_hits_outputNames:
18         relativeFileName = os.path.basename(file)
19         tfile.add(relativeFileName)
20     os.system("mv %s %s" % (tarFileName + ".tmp.tar", tarFileName))
21     tfile.close()
22     os.chdir(currentPath)
23     
24 def splitFastQ(fileName, nbOfSeqPerBatch):
25     nbOfLinesPerFile = nbOfSeqPerBatch * 4
26     lOutput = []
27     filePrefix, fileExt = os.path.splitext(os.path.basename(fileName))
28     resDir = os.path.dirname(fileName)
29     with open(fileName) as inF:
30         fileNb = 1
31         line = inF.readline()
32         if not line or nbOfLinesPerFile == 0:
33             outFileName = "%s/%s-%s%s" %(resDir, filePrefix, fileNb, fileExt)
34             lOutput.append(outFileName)
35             f = open(outFileName, "wb")
36             shutil.copyfileobj(open(fileName, "rb"), f)
37             f.close()
38         else:
39             while line:
40                 outFileName = "%s/%s-%s%s" %(resDir, filePrefix, fileNb, fileExt)
41                 lOutput.append(outFileName)
42                 with open(outFileName, "w") as outF:
43                     lineNb = 1
44                     while lineNb <= nbOfLinesPerFile and line:
45                         outF.write(line)
46                         line = inF.readline()
47                         lineNb += 1
48                 fileNb += 1
49     return lOutput
50 
51 def joinBAM(dCutOut2Out):
52     for key in dCutOut2Out.keys():
53         fh = open(key, "w") 
54         fh.close()
55         nbFile = 0
56         cmd = "samtools merge -f %s" % key
57         for fileName in dCutOut2Out[key]:
58             nbFile = nbFile + 1
59             if nbFile < 225:
60                 cmd += " %s" % fileName
61             else:
62                 nbFile = 0
63                 cmd += ";mv %s tmpBAM;" % (key)
64                 cmd += "samtools merge -f %s tmpBAM %s" %  (key, fileName)
65         proc = subprocess.Popen( args=cmd , shell=True)
66         returncode = proc.wait()
67 
68         
69 def _map(iLauncher, cmd, cmdStart, cmdFinish ):
70     lCmds = []
71     lCmds.extend(cmd)
72     lCmdStart = []
73     lCmdStart.extend(cmdStart)
74     lCmdFinish = []
75     lCmdFinish.extend(cmdFinish)
76     return(iLauncher.prepareCommands_withoutIndentation(lCmds, lCmdStart, lCmdFinish))
77 
78 def _createTopHatCommand(iLauncher, options, index_paths, inputFileNames, inputRevFilesNames, space):
79     lArgs = []
80     lArgs.append('-p %s %s' % ( options.num_threads, space ))
81     if options.single_paired == 'paired':
82         lArgs.append('-r %s ' % options.mate_inner_dist)
83     if options.settings == 'preSet':
84         lArgs.append(index_paths)
85         lArgs.append(inputFileNames)
86         if options.input2:
87             lArgs.append(inputRevFilesNames)
88         return iLauncher.getSystemCommand("tophat", lArgs)
89     else:
90         if int( options.min_anchor_length ) >= 3:
91             lArgs.append('-a %s ' % options.min_anchor_length)
92         else:
93             raise Exception, 'Minimum anchor length must be 3 or greater'
94         lArgs.append('-m %s ' % options.splice_mismatches)
95         lArgs.append('-i %s ' % options.min_intron_length)
96         lArgs.append('-I %s ' % options.max_intron_length)
97         if float( options.junction_filter ) != 0.0:
98             lArgs.append('-F %s ' % options.junction_filter)
99         lArgs.append('-g %s ' % options.max_multihits)
00         # Custom junctions options.
01         if options.gene_model_annotations:
02             lArgs.append('-G %s ' % options.gene_model_annotations)
03         if options.raw_juncs:
04             lArgs.append('-j %s ' % options.raw_juncs)
05         if options.no_novel_juncs:
06             lArgs.append('--no-novel-juncs ')
07         if options.library_type:
08             lArgs.append('--library-type %s ' % options.library_type)
09         if options.no_novel_indels:
10             lArgs.append('--no-novel-indels ')
11         else:
12             if options.max_insertion_length:
13                 lArgs.append('--max-insertion-length %i ' % int( options.max_insertion_length ))
14             if options.max_deletion_length:
15                 lArgs.append('--max-deletion-length %i ' % int( options.max_deletion_length ))
16             # Max options do not work for Tophat v1.2.0, despite documentation to the contrary. (Fixed in version 1.3.1)
17             # need to warn user of this fact
18             #sys.stdout.write( "Max insertion length and max deletion length options don't work in Tophat v1.2.0\n" )
19 
20         # Search type options.
21         if options.coverage_search:
22             lArgs.append('--coverage-search --min-coverage-intron %s --max-coverage-intron %s ' % ( options.min_coverage_intron, options.max_coverage_intron ))
23         else:
24             lArgs.append('--no-coverage-search ')
25         if options.closure_search:
26             lArgs.append('--closure-search --min-closure-exon %s --min-closure-intron %s --max-closure-intron %s '  % ( options.min_closure_exon, options.min_closure_intron, options.max_closure_intron )) 
27         else:
28             lArgs.append('--no-closure-search ')
29         if options.microexon_search:
30             lArgs.append('--microexon-search ')
31         if options.single_paired == 'paired':
32             lArgs.append('--mate-std-dev %s ' % options.mate_std_dev)
33         if options.initial_read_mismatches:
34             lArgs.append('--initial-read-mismatches %d ' % int( options.initial_read_mismatches ))
35         if options.seg_mismatches:
36             lArgs.append('--segment-mismatches %d ' % int( options.seg_mismatches ))
37         if options.seg_length:
38             lArgs.append('--segment-length %d ' % int( options.seg_length ))
39         if options.min_segment_intron:
40             lArgs.append('--min-segment-intron %d ' % int( options.min_segment_intron ))
41         if options.max_segment_intron:
42             lArgs.append('--max-segment-intron %d ' % int( options.max_segment_intron ))
43         lArgs.append(index_paths)
44         lArgs.append(inputFileNames)
45         if options.input2:
46             lArgs.append(inputRevFilesNames)
47         return iLauncher.getSystemCommand("tophat", lArgs)
48         
49 
50 
51 def __main__():
52     #Parse Command Line
53     parser = optparse.OptionParser()
54     parser.add_option('-o', '--outputTxtFile', dest='outputTxtFile', help='for Differential expression analysis pipeline, new output option gives a txt output containing the list of mapping results.')
55     parser.add_option('-t', '--tar', dest='outputTar', default=None, help='output all accepted hits results in a tar file.' )
56     parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' )
57     parser.add_option( '-C', '--color-space', dest='color_space', action='store_true', help='This indicates color-space data' )
58     parser.add_option( '-J', '--junctions-output', dest='junctions_output_file', default='junctions_output.bed', help='Junctions output file; formate is BED.' )
59     parser.add_option( '-H', '--hits-output', dest='accepted_hits_output_file', default='hits_output.bam', help='Accepted hits output file; formate is BAM.' )
60     parser.add_option( '', '--own-file', dest='own_file', help='' )
61     parser.add_option( '-D', '--indexes-path', dest='index_path', help='Indexes directory; location of .ebwt and .fa files.' )
62     parser.add_option( '-r', '--mate-inner-dist', dest='mate_inner_dist', help='This is the expected (mean) inner distance between mate pairs. \
63                                                                                 For, example, for paired end runs with fragments selected at 300bp, \
64                                                                                 where each end is 50bp, you should set -r to be 200. There is no default, \
65                                                                                 and this parameter is required for paired end runs.')
66     parser.add_option( '', '--mate-std-dev', dest='mate_std_dev', help='Standard deviation of distribution on inner distances between male pairs.' )
67     parser.add_option( '-a', '--min-anchor-length', dest='min_anchor_length', 
68                         help='The "anchor length". TopHat will report junctions spanned by reads with at least this many bases on each side of the junction.' )
69     parser.add_option( '-m', '--splice-mismatches', dest='splice_mismatches', help='The maximum number of mismatches that can appear in the anchor region of a spliced alignment.' )
70     parser.add_option( '-i', '--min-intron-length', dest='min_intron_length', 
71                         help='The minimum intron length. TopHat will ignore donor/acceptor pairs closer than this many bases apart.' )
72     parser.add_option( '-I', '--max-intron-length', dest='max_intron_length', 
73                         help='The maximum intron length. When searching for junctions ab initio, TopHat will ignore donor/acceptor pairs farther than this many bases apart, except when such a pair is supported by a split segment alignment of a long read.' )
74     parser.add_option( '-F', '--junction_filter', dest='junction_filter', help='Filter out junctions supported by too few alignments (number of reads divided by average depth of coverage)' )
75     parser.add_option( '-g', '--max_multihits', dest='max_multihits', help='Maximum number of alignments to be allowed' )
76     parser.add_option( '', '--initial-read-mismatches', dest='initial_read_mismatches', help='Number of mismatches allowed in the initial read mapping' )
77     parser.add_option( '', '--seg-mismatches', dest='seg_mismatches', help='Number of mismatches allowed in each segment alignment for reads mapped independently' )
78     parser.add_option( '', '--seg-length', dest='seg_length', help='Minimum length of read segments' )
79     parser.add_option( '', '--library-type', dest='library_type', help='TopHat will treat the reads as strand specific. Every read alignment will have an XS attribute tag. Consider supplying library type options below to select the correct RNA-seq protocol.' )
80     parser.add_option( '', '--allow-indels', action="store_true", help='Allow indel search. Indel search is disabled by default.(Not used since version 1.3.0)' )
81     parser.add_option( '', '--max-insertion-length', dest='max_insertion_length', help='The maximum insertion length. The default is 3.' )
82     parser.add_option( '', '--max-deletion-length', dest='max_deletion_length', help='The maximum deletion length. The default is 3.' )
83 
84     # Options for supplying own junctions
85     parser.add_option( '-G', '--GTF', dest='gene_model_annotations', help='Supply TopHat with a list of gene model annotations. \
86                                                                            TopHat will use the exon records in this file to build \
87                                                                            a set of known splice junctions for each gene, and will \
88                                                                            attempt to align reads to these junctions even if they \
89                                                                            would not normally be covered by the initial mapping.')
90     parser.add_option( '-j', '--raw-juncs', dest='raw_juncs', help='Supply TopHat with a list of raw junctions. Junctions are \
91                                                                     specified one per line, in a tab-delimited format. Records \
92                                                                     look like: <chrom> <left> <right> <+/-> left and right are \
93                                                                     zero-based coordinates, and specify the last character of the \
94                                                                     left sequenced to be spliced to the first character of the right \
95                                                                     sequence, inclusive.')
96     parser.add_option( '', '--no-novel-juncs', action="store_true", dest='no_novel_juncs', help="Only look for junctions indicated in the \
97                                                                                             supplied GFF file. (ignored without -G)")
98     parser.add_option( '', '--no-novel-indels', action="store_true", dest='no_novel_indels', help="Skip indel search. Indel search is enabled by default.")
99     # Types of search.
00     parser.add_option( '', '--microexon-search', action="store_true", dest='microexon_search', help='With this option, the pipeline will attempt to find alignments incident to microexons. Works only for reads 50bp or longer.')
01     parser.add_option( '', '--closure-search', action="store_true", dest='closure_search', help='Enables the mate pair closure-based search for junctions. Closure-based search should only be used when the expected inner distance between mates is small (<= 50bp)')
02     parser.add_option( '', '--no-closure-search', action="store_false", dest='closure_search' )
03     parser.add_option( '', '--coverage-search', action="store_true", dest='coverage_search', help='Enables the coverage based search for junctions. Use when coverage search is disabled by default (such as for reads 75bp or longer), for maximum sensitivity.')
04     parser.add_option( '', '--no-coverage-search', action="store_false", dest='coverage_search' )
05     parser.add_option( '', '--min-segment-intron', dest='min_segment_intron', help='Minimum intron length that may be found during split-segment search' )
06     parser.add_option( '', '--max-segment-intron', dest='max_segment_intron', help='Maximum intron length that may be found during split-segment search' )
07     parser.add_option( '', '--min-closure-exon', dest='min_closure_exon', help='Minimum length for exonic hops in potential splice graph' )
08     parser.add_option( '', '--min-closure-intron', dest='min_closure_intron', help='Minimum intron length that may be found during closure search' )
09     parser.add_option( '', '--max-closure-intron', dest='max_closure_intron', help='Maximum intron length that may be found during closure search' )
10     parser.add_option( '', '--min-coverage-intron', dest='min_coverage_intron', help='Minimum intron length that may be found during coverage search' )
11     parser.add_option( '', '--max-coverage-intron', dest='max_coverage_intron', help='Maximum intron length that may be found during coverage search' )
12 
13     # Wrapper options.
14     parser.add_option( '-1', '--input1', dest='input1', help='A list of the (forward or single-end) reads files of Sanger FASTQ format, txt format' )
15     parser.add_option( '-2', '--input2', dest='input2', help='The list of reverse reads file in Sanger FASTQ format' )
16     parser.add_option( '', '--single-paired', dest='single_paired', help='' )
17     parser.add_option( '', '--settings', dest='settings', help='' )
18 
19     (options, args) = parser.parse_args()
20 
21     # output version # of tool
22     try:
23         tmp_files = []
24         tmp = tempfile.NamedTemporaryFile().name
25         tmp_files.append(tmp)
26         tmp_stdout = open( tmp, 'wb' )
27         proc = subprocess.Popen( args='tophat -v', shell=True, stdout=tmp_stdout )
28         tmp_stdout.close()
29         returncode = proc.wait()
30         stdout = open( tmp_stdout.name, 'rb' ).readline().strip()
31         if stdout:
32             sys.stdout.write( '%s\n' % stdout )
33         else:
34             raise Exception
35     except:
36         sys.stdout.write( 'Could not determine Tophat version\n' )
37 
38     # Color or base space
39     space = ''
40     if options.color_space:
41         space = '-C'
42 
43 
44     #reads = options.input1
45     file = open(options.input1,"r")
46     lines = file.readlines()
47     inputFileNames = []
48     accepted_hits_outputNames = []
49     outputName = options.outputTxtFile
50     resDirName = os.path.dirname(outputName) + '/'
51     out = open(outputName, "w")
52     for line in lines:
53         tab = line.split()
54         inputFileNames.append(tab[1])
55         aHitOutName = resDirName + tab[0] + '_' + options.accepted_hits_output_file
56         accepted_hits_outputNames.append(aHitOutName) 
57         out.write(tab[0] + '\t' + aHitOutName + '\n')
58     file.close()
59     out.close()
60     
61     if options.input2:
62         revFile = open(options.input2,"r")
63         lines = revFile.readlines()
64         inputRevFileNames = []
65         for line in lines:
66             revTab = line.split()
67             inputRevFileNames.append(revTab[1])
68         revFile.close()
69 
70     
71     # Creat bowtie index if necessary.
72     tmp_index_dirs = []
73     index_paths = []
74     tmp_index_dir = tempfile.mkdtemp(dir="%s" % os.getcwd())
75     tmp_index_dirs.append(tmp_index_dir)
76     if options.own_file:
77         index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( options.own_file )[1].split( '.' )[:-1] ) )
78         index_paths.append(index_path)
79         try:
80             os.link( options.own_file, index_path + '.fa' )
81         except:
82             # Tophat prefers (but doesn't require) fasta file to be in same directory, with .fa extension
83             pass
84         lCmdsTuples =[]
85         acronym = "tophat_index"
86         jobdb = TableJobAdaptatorFactory.createJobInstance()
87         iLauncher = Launcher(jobdb, os.getcwd(), "", "", os.getcwd(), os.getcwd(), "jobs", "", acronym, acronym, False, True)
88         cmd_index = iLauncher.getSystemCommand("bowtie-build", [space, "-f %s" % options.own_file, index_path])
89         cmd2Launch = []
90         cmdStart = []
91         cmdFinish = []
92         cmd2Launch.append(cmd_index)
93         lCmdsTuples.append(_map(iLauncher, cmd2Launch, cmdStart, cmdFinish)) 
94         iLauncher.runLauncherForMultipleJobs(acronym, lCmdsTuples, True)
95     else:
96         for file in inputFileNames:
97             tmp_index_dir = tempfile.mkdtemp()
98             index_path = tmp_index_dir + '/' + os.path.basename(file).split('.')[0]
99             index_paths.append(index_path)
00             tmp_index_dirs.append(tmp_index_dir)
01 
02     
03     
04     acronym = "tophat"
05     jobdb = TableJobAdaptatorFactory.createJobInstance()
06     iLauncher = Launcher(jobdb, os.getcwd(), "", "", os.getcwd(), os.getcwd(), "jobs", "", acronym, acronym, False, True)
07     lCmdsTuples = []
08     dCutOut2Out = {}
09     lAllFile2remove = []
10     # for inputFileName in inputFileNames:
11     for i in range(len(inputFileNames)):
12         lCutOutput = []
13         lCutInputFile = splitFastQ(inputFileNames[i], 20000)
14         lAllFile2remove.extend(lCutInputFile)
15         if options.input2:
16             lCutPairInputFile = splitFastQ(inputRevFileNames[i], 20000)
17             lAllFile2remove.extend(lCutPairInputFile)
18         for j in range(len(lCutInputFile)):
19             cutOutput = "%s_out" %  lCutInputFile[j]
20             lCutOutput.append(cutOutput)
21             lAllFile2remove.extend(lCutOutput)
22             cmd2Launch = []
23             if options.input2: 
24                 inputRevFile = lCutPairInputFile[j]
25             else:
26                 inputRevFile = ""
27             if options.own_file:
28                 cmd2Launch.append(_createTopHatCommand(iLauncher, options, index_paths[0], lCutInputFile[j], inputRevFile, space))
29             else:
30                 cmd2Launch.append(_createTopHatCommand(iLauncher, options, index_paths[i], lCutInputFile[j], inputRevFile, space))
31             cmdStart = []
32             cmdFinish = ["shutil.copyfile( os.path.join( 'tophat_out', 'accepted_hits.bam' ), '%s')" % cutOutput]
33             lCmdsTuples.append(_map(iLauncher, cmd2Launch, cmdStart, cmdFinish))    
34         dCutOut2Out[accepted_hits_outputNames[i]] = lCutOutput
35     iLauncher.runLauncherForMultipleJobs(acronym, lCmdsTuples, True)
36     
37     joinBAM(dCutOut2Out) 
38     FileUtils.removeFilesFromListIfExist(lAllFile2remove)   
39                  
40     if options.outputTar != None:
41         toTar(options.outputTar, accepted_hits_outputNames)
42 
43     
44     # Clean up temp dirs
45     for tmp_index_dir in tmp_index_dirs:
46         if os.path.exists( tmp_index_dir ):
47             shutil.rmtree( tmp_index_dir )
48 
49     for tmp in tmp_files:
50         os.remove(tmp)
51 
52 
53 if __name__=="__main__": __main__()