annotate SMART/DiffExpAnal/gsnap_parallel_unSQL.py @ 32:3441fe98a2ba

Deleted selected files

1 #!/usr/bin/env python
2 
3 import optparse, os, shutil, subprocess, sys, tempfile, fileinput, tarfile, glob 
4 import time
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 from optparse import OptionParser
9 
10 def stop_err( msg ):
11     sys.stderr.write( "%s\n" % msg )
12     sys.exit()
13 
14 def toTar(tarFileName, accepted_hits_outputNames):
15     tfile = tarfile.open(tarFileName + ".tmp.tar", "w")
16     currentPath = os.getcwd()
17     os.chdir(dir)
18     for file in accepted_hits_outputNames:
19         relativeFileName = os.path.basename(file)
20         tfile.add(relativeFileName)
21     os.system("mv %s %s" % (tarFileName + ".tmp.tar", tarFileName))
22     tfile.close()
23     os.chdir(currentPath)
24     
25 def joinSAM(dCutOut2Out):
26     for key in dCutOut2Out.keys():
27         FileUtils.catFilesFromList(dCutOut2Out[key],key, False)
28         
29 def _map(iLauncher, cmd, cmdStart, cmdFinish ):
30     lCmds = []
31     lCmds.extend(cmd)
32     lCmdStart = []
33     lCmdStart.extend(cmdStart)
34     lCmdFinish = []
35     lCmdFinish.extend(cmdFinish)
36     return(iLauncher.prepareCommands_withoutIndentation(lCmds, lCmdStart, lCmdFinish))
37 
38 def _createGsnapSplicingOptions(options):
39     lArgs = []
40     lArgs.append("-N %s" % options.novelsplicing)
41     if options.useSplicing:
42         lArgs.append("-s %s" % options.useSplicing)
43     lArgs.append("-w %s" % options.localsplicedist)
44     lArgs.append("-e %s" % options.localSplicePenality)
45     lArgs.append("-E %s" % options.distantSplicePenality)
46     lArgs.append("-K %s" % options.distantSpliceEndlength)
47     lArgs.append("-l %s" % options.shortendSpliceEndlength)
48     
49     
50     return lArgs
51 
52 def _createGsnapPairedEndOptions(options):
53     lArgs = []
54     if not(options.useSplicing or options.pairedEndFile):
55         lArgs.append("--pairmax-dna %s" % options.pairmaxRna)
56     if options.useSplicing or options.pairedEndFile:
57         lArgs.append("--pairmax-rna %s" % options.pairmaxRna)
58     lArgs.append("--pairexpect=%s" % options.pairexpect)
59     lArgs.append("--pairdev=%s" % options.pairedev)
60     
61     
62 
63 def _createGsnapCommand(iLauncher, options, workingDir, inputFileNames, inputRevFilesNames, outputFileName, batchNumber, numberOfBatch):
64     lArgs = []
65     lArgs.append("-d %s" % options.genomeName)
66     lArgs.append("-k %s" % options.kmer)
67     lArgs.append("-D %s" % workingDir)
68     lArgs.append("-A %s" % options.outputFormat)
69     lArgs.append("-q %s/%s" % (batchNumber, numberOfBatch))
70     lArgs.append("--no-sam-headers")
71     lArgs.append(inputFileNames)
72     print 'N option: %s, pairedEndFile option: %s' %(options.novelsplicing, options.pairedEndFile)
73     if options.pairedEndFile:
74         lArgs.append(inputRevFilesNames)
75     if options.novelsplicing == '1':
76         lArgs.extend(_createGsnapSplicingOptions(options))
77     elif options.pairedEndFile:
78         lArgs.extend(_createGsnapPairedEndOptions(options))
79     
80     lArgs.append("> %s" % outputFileName)  
81     return iLauncher.getSystemCommand("gsnap", lArgs)   
82 
83 def __main__():
84     #Parse Command Line
85     description = "GMAP/GSNAP version:2012-12-20."
86     parser = OptionParser(description = description)
87     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.')
88     parser.add_option("-q", "--inputTxt",         dest="inputTxt",         action="store",                 type="string",     help="input, a txt file for a list of input reads files [compulsory]")
89     parser.add_option('-t', '--tar', dest='outputTar', default=None, help='output all accepted hits results in a tar file.' )
90     parser.add_option("-d", "--genomeName", dest="genomeName", help="Define the reference genome name.[compulsory]")
91 #    parser.add_option("-o", "--outputFile", dest="outputfile", help="output[compulsory]")
92     parser.add_option("-k", "--kmer", dest="kmer", default=12, help="Choose kmer value (<=16), a big kmer value can take more RAM(4Go).[compulsory]")
93     parser.add_option("-i", "--inputFasta", dest="inputFastaFile", help="Reference genome file, fasta format.[compulsory]")
94     parser.add_option("-p", "--pairedEnd", dest="pairedEndFile", default=None, help="Input paired-end fastq file.")
95     parser.add_option("-A", "--outputFormat", dest="outputFormat", default="sam", help="Choose an output format [sam, goby (need to re-compile with appropriate options)].")
96     #Splicing options for RNA-Seq
97     parser.add_option("-N","--novelsplicing", dest="novelsplicing", default=0, help="Look for novel splicing (0=no (default), 1=yes)")
98     parser.add_option("-s","--use-splicing", dest="useSplicing",action="store", type="string", help="Look for splicing involving known sites or known introns (in <STRING>.iit), at short or long distances")
99     parser.add_option("-w","--localsplicedist", dest="localsplicedist", default=200000, help="Definition of local novel splicing event (default 200000)")
00     parser.add_option("-e","--local-splice-penality", dest="localSplicePenality", default=0, help="Penalty for a local splice (default 0).  Counts against mismatches allowed")
01     parser.add_option("-E","--distant-splice-penality", dest="distantSplicePenality", default=1, help="Penalty for a distant splice (default 1).  A distant splice is one where the intron length exceeds the value of -w, or --localsplicedist, or is an inversion, scramble, or translocation between two different chromosomes Counts against mismatches allowed")
02     parser.add_option("-K","--distant-splice-endlength", dest="distantSpliceEndlength", default=16, help="Minimum length at end required for distant spliced alignments (default 16, min allowed is the value of -k, or kmer size)")
03     parser.add_option("-l","--shortend-splice-endlength", dest="shortendSpliceEndlength", default=2, help="Minimum length at end required for short-end spliced alignments (default 2, but unless known splice sites are provided with the -s flag, GSNAP may still need the end length to be the value of -k, or kmer size to find a given splice")
04 
05     #Specific paired-ends reads
06     parser.add_option("--pairmax-dna", dest="pairmaxDna", default=1000, help="Max total genomic length for DNA-Seq paired reads, or other reads without splicing (default 1000).")
07     parser.add_option("--pairmax-rna", dest="pairmaxRna", default=2000, help="Max total genomic length for RNA-Seq paired reads, or other reads that could have a splice (default 200000).")
08     parser.add_option("--pairexpect", dest="pairexpect", default=200, help="Expected paired-end length, used for calling splices in medial part of paired-end reads (default 200)")
09     parser.add_option("--pairdev", dest="pairdev", default=25, help="Allowable deviation from expected paired-end length, used for calling splices in medial part of paired-end reads (default 25)")
10     
11     (options, args) = parser.parse_args()    
12 
13     workingDir = os.path.dirname(options.inputFastaFile)
14     
15     file = open(options.inputTxt,"r")
16     lines = file.readlines()
17     inputFileNames = []
18     gsnapOutputNames = []
19     outputName = options.outputTxtFile
20     resDirName = os.path.dirname(outputName) + '/'
21     out = open(outputName, "w")
22     for line in lines:
23         timeId = time.strftime("%Y%m%d%H%M%S")
24         tab = line.split()
25         inputFileNames.append(tab[1])
26         OutputName = resDirName + tab[0] + '_samOutput_%s.sam' % timeId
27         gsnapOutputNames.append(OutputName) 
28         out.write(tab[0] + '\t' + OutputName + '\n')
29     file.close()
30     out.close()
31     
32     if options.pairedEndFile:
33         revFile = open(options.pairedEndFile,"r")
34         lines = revFile.readlines()
35         inputRevFileNames = []
36         for line in lines:
37             revTab = line.split()
38             inputRevFileNames.append(revTab[1])
39         revFile.close()
40 
41     #Create gsnap make 
42     lCmdsTuples =[]
43     acronym = "gsnap_make"
44     jobdb = TableJobAdaptatorFactory.createJobInstance()
45     iLauncher = Launcher(jobdb, os.getcwd(), "", "", os.getcwd(), os.getcwd(), "jobs", "", acronym, acronym, False, True)
46     cmds = []
47     cmd_setup = "gmap_setup -d %s -D %s -k %s %s;" % (options.genomeName, workingDir, options.kmer, options.inputFastaFile)
48     cmds.append(cmd_setup)
49     cmd_make_coords = "make -f Makefile.%s coords;" % options.genomeName 
50     cmds.append(cmd_make_coords)
51     cmd_make_gmapdb = "make -f Makefile.%s gmapdb;" % options.genomeName
52     cmds.append(cmd_make_gmapdb)
53     cmd_make_install = "make -f Makefile.%s install;" % options.genomeName
54     cmds.append(cmd_make_install)
55     cmd_index = iLauncher.getSystemCommand("", cmds)
56     cmd2Launch = []
57     cmdStart = []
58     cmdFinish = []
59     cmd2Launch.append(cmd_index)
60     lCmdsTuples.append(_map(iLauncher, cmd2Launch, cmdStart, cmdFinish)) 
61     iLauncher.runLauncherForMultipleJobs(acronym, lCmdsTuples, True)    
62     
63     acronym = "gsnap"
64     jobdb = TableJobAdaptatorFactory.createJobInstance()
65     iLauncher = Launcher(jobdb, os.getcwd(), "", "", os.getcwd(), os.getcwd(), "jobs", "", acronym, acronym, False, True)
66     lCmdsTuples = []
67     dCutOut2Out = {}
68     lAllFile2remove = []
69     numberOfBatch = 20  #usually for testing, working on to find a value for default launch on galaxy
70     for i in range(len(inputFileNames)):
71         lCutOutput = []
72         for j in range(numberOfBatch):
73             cutOutput = "%s_out_%s" % (inputFileNames[i], j)
74             lCutOutput.append(cutOutput)
75             lAllFile2remove.extend(lCutOutput)
76             cmd2Launch = []
77             if options.pairedEndFile: 
78                 inputRevFile = inputRevFileNames[i]
79             else:
80                 inputRevFile = ""
81             cmd2Launch.append(_createGsnapCommand(iLauncher, options, workingDir, inputFileNames[i], inputRevFile, cutOutput, j, numberOfBatch))
82             cmdStart = []
83             cmdFinish = []
84             lCmdsTuples.append(_map(iLauncher, cmd2Launch, cmdStart, cmdFinish))    
85             dCutOut2Out[gsnapOutputNames[i]] = lCutOutput
86     iLauncher.runLauncherForMultipleJobs(acronym, lCmdsTuples, True)
87     
88     joinSAM(dCutOut2Out) 
89     FileUtils.removeFilesFromListIfExist(lAllFile2remove)   
90                  
91     if options.outputTar != None:
92         toTar(options.outputTar, gsnapOutputNames)
93 
94 
95 if __name__=="__main__": __main__()