annotate pyPRADA_1.2/prada-preprocess-bi @ 3:f17965495ec9 draft default tip

Uploaded
author siyuan
date Tue, 11 Mar 2014 12:14:01 -0400
parents acc2ca1a3ba4
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
1 #!/usr/bin/env python
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
2
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
3 #This is a program to implement the PRADA pipeline "process subsection", based on the work of Rahul Vegesna and Wandaliz Torres-Garcia
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
4 #It merely generates a PBS file, depending on the user input entry point (step)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
5 #Author: Siyuan Zheng, szheng2@mdanderson.org
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
6 #Copy right belongs to Roel Verhaak's lab from MD Anderson Cancer Center, Department of Bioinformatics and Computational Biology.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
7 #Last revision: 02/17/2014
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
8
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
9 import subprocess
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
10 import os,os.path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
11 import sys
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
12 import time
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
13 import ioprada
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
14 import re
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
15
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
16 ########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
17 args=sys.argv
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
18
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
19 help_menu='''\nPipeline for RNAseq Data Analaysis - preprocessing pipeline (PRADA).
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
20 \t**Usage**:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
21 \tprada-preprocess-bi -conf xx.txt -inputdir .. -sample XX -tag TCGA-XX -platform illumina -step 1_1 -intermediate no -pbs xxx -outdir ... -submit no
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
22 \t**Parameters**:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
23 \t-h print help message
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
24 \t-step_info print complete steps curated in the module.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
25 \t-inputdir the dir where the input bam or fastq can be found.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
26 \t-sample input sample name. PRADA searches for sample.bam or sample.end1.fastq/sample.end2.fastq, etc, depending on the
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
27 \t initiating step number. See step_info for more information.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
28 \t-conf config file for references and parameters. Default is conf.txt in py-PRADA installation folder.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
29 \t-tag a tag to describe the sample, likely sample ID, such as TCGA-LGG-01; no default.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
30 \t-platform only illumina at present (default).
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
31 \t-step values: 1_1/2,2_e1/2_1/2/3/4,3_e1/2_1/2,4_1/2,5,6_1/2,7,8; example 2_e1_1; no default.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
32 \t-outdir output dir. Default is the directory where the input bam is.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
33 \t-pbs name for output pbs file and log file. Default (time-stamp) is used if no input.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
34 \t-intermediate values:yes/no; if intermediate files should be kept. Default is not.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
35 \t-submit if submit the job to HPC, default is no. If yes, ppn is set to 12.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
36 \t-v print version information.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
37 '''
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
38
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
39 steps_info='''
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
40 Command orders (sample XX)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
41 step 1_1 --> XX.sorted.bam [sort input bam by name]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
42 step 1_2 --> XX.end1/2.fastq [extract reads from bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
43 step 2_e1_1 --> XX.end1.sai [realign end1 reads to composite reference]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
44 step 2_e1_2 --> XX.end1.sam [generate sam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
45 step 2_e1_3 --> XX.end1.bam [generate bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
46 step 2_e1_4 --> XX.end1.sorted.bam [sort bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
47 step 2_e2_1 --> XX.end2.sai [realign end2 reads to composite reference]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
48 step 2_e2_2 --> XX.end2.sam [generate sam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
49 step 2_e2_3 --> XX.end2.bam [generate bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
50 step 2_e2_4 --> XX.end2.sorted.bam [sort bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
51 step 3_e1_1 --> XX.end1.remapped.bam [remap end1 to genome]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
52 step 3_e1_2 --> XX.end1.remapped.sorted.bam [sort remapped bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
53 step 3_e2_1 --> XX.end2.remapped.bam [remap end2 to genome]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
54 step 3_e2_2 --> XX.end2.remapped.sorted.bam [sort remapped bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
55 step 4_1 --> XX.paired.bam [pair end1 and end1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
56 step 4_2 --> XX.paired.sorted.bam [sort paired bam]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
57 step 5 --> XX.withRG.paired.sorted.bam [add read group]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
58 step 6_1 --> XX.orig.csv [prepair recalibration table]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
59 step 6_2 --> XX.withRG.GATKRecalibrated.bam [recalibration]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
60 step 7 --> XX.withRG.GATKRecalibrated.flagged.bam [flag duplication reads]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
61 step 8 --> folder XX for gene expression, QC metrics etc. [generate QC and expression]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
62 '''
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
63
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
64 if '-h' in args or '-help' in args or len(args)==1:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
65 print help_menu
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
66 sys.exit(0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
67 if '-step_info' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
68 print steps_info
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
69 sys.exit(0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
70 if '-v' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
71 import version
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
72 print version.version
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
73 sys.exit(0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
74
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
75 if '-sample' not in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
76 sys.exit('ERROR: Sample name is needed')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
77 if '-step' not in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
78 sys.exit('ERROR: Step number is needed')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
79 if '-tag' not in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
80 sys.exit('ERROR: A tag is needed')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
81
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
82 i=args.index('-sample')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
83 sample=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
84 if '-inputdir' not in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
85 inputpath=os.path.abspath('./')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
86 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
87 i=args.index('-inputdir')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
88 inputpath=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
89 bampath=os.path.abspath(inputpath)+'/%s.bam'%sample
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
90 fq1path=os.path.abspath(inputpath)+'/%s.end1.fastq'%sample
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
91 fq2path=os.path.abspath(inputpath)+'/%s.end2.fastq'%sample
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
92
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
93 if '-outdir' not in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
94 outpath=os.path.dirname(bampath)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
95 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
96 i=args.index('-outdir')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
97 outpath=os.path.abspath(args[i+1])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
98 if not os.path.exists(outpath):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
99 os.mkdir(outpath)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
100
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
101 #bam=os.path.basename(bampath)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
102 #sample=bam[:-4]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
103 i=args.index('-step')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
104 step=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
105 i=args.index('-tag')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
106 tag=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
107
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
108 prada_path=os.path.dirname(os.path.abspath(__file__)) ####
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
109 ref_search_path=[prada_path,os.getcwd()] #search path for ref file if not specified in command line
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
110
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
111 if '-conf' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
112 i=args.index('-conf')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
113 reffile=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
114 if os.path.exists(reffile):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
115 pass
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
116 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
117 for pth in ref_search_path:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
118 new_reffile='%s/%s'%(pth, os.path.basename(reffile))
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
119 if os.path.exists(new_reffile):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
120 reffile=new_reffile
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
121 break
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
122 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
123 sys.exit('ERROR: conf file %s not found'%reffile)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
124 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
125 reffile='%s/conf.txt'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
126 if not os.path.exists(reffile):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
127 sys.exit('ERROR: No default conf.txt found and none specified')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
128
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
129 if '-platform' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
130 i=args.index('-platform')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
131 plat=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
132 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
133 plat='illumina'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
134 if '-submit' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
135 i=args.index('-submit')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
136 submit=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
137 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
138 submit='no'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
139 if '-intermediate' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
140 i=args.index('-intermediate')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
141 keepmed=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
142 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
143 keepmed='False'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
144
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
145 if keepmed in ['False','FALSE','false','F','NO','No','no','n']:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
146 keepmed_flag='no'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
147 elif keepmed in ['True','TRUE','true','T','YES','Yes','yes','y']:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
148 keepmed_flag='yes'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
149 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
150 sys.exit('ERROR: -intermediate value not recognized')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
151
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
152 if '-pbs' in args:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
153 i=args.index('-pbs')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
154 docstr=args[i+1]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
155 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
156 a=time.ctime().split()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
157 b=time.time()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
158 timestamp='_'.join([a[-1],a[1],a[2]])+'.'+str(b)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
159 docstr='prada_prep_'+timestamp
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
160 logfilename=docstr+'.log'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
161 pbsfilename=docstr+'.pbs'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
162 pbspath=outpath+'/'+pbsfilename
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
163 logpath=outpath+'/'+logfilename
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
164 ########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
165
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
166 ########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
167 #underlying utilities, automatically detected
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
168 samtools='%s/tools/samtools-0.1.16/samtools'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
169 bwa='%s/tools/bwa-0.5.7-mh/bwa'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
170 gatk='%s/tools/GATK/'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
171 picard='%s/tools/Picard/'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
172 seqc='%s/tools/RNA-SeQC_v1.1.7.jar'%prada_path
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
173 #Default uses 12 nodes in HPC
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
174 #########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
175
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
176 #########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
177 #reference files
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
178 refdict=ioprada.read_conf(reffile)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
179 genome_gtf=refdict['--REF--']['genome_gtf']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
180 compdb_fasta=refdict['--REF--']['compdb_fasta']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
181 compdb_map=refdict['--REF--']['compdb_map']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
182 genome_fasta=refdict['--REF--']['genome_fasta']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
183 dbsnp_vcf=refdict['--REF--']['dbsnp_vcf']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
184 select_tx=refdict['--REF--']['select_tx']
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
185 pat=re.compile('ppn=(\d*)')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
186 parallel_n=pat.search(refdict['--PBS--']['-l']).groups()[0]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
187 #########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
188
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
189 #########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
190 #pipeline command lines.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
191 ##Cleaning up steps: if -intermediate is yes, none will be executed.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
192 step_1_1_cmd=['%s sort -n -m 1000000000 %s %s.sorted'%(samtools,bampath,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
193 post_1_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
194 step_1_2_cmd=['java -Xmx8g -jar %s/SamToFastq.jar INPUT=%s.sorted.bam FASTQ=%s.end1.fastq SECOND_END_FASTQ=%s.end2.fastq INCLUDE_NON_PF_READS=true VALIDATION_STRINGENCY=SILENT TMP_DIR=tmp/'%(picard,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
195 post_1_2_clean=['rm -f %s.sorted.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
196 alnparams=' '.join([' '.join(x) for x in refdict['--BWA aln--'].items()])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
197 samseparams=' '.join([' '.join(x) for x in refdict['--BWA samse--'].items()])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
198 if step == '2_e1_1' or step == '2_e2_1':
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
199 step_2_e1_1_cmd=['%s aln %s %s %s > %s.end1.sai'%(bwa,alnparams,compdb_fasta,fq1path,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
200 post_2_e1_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
201 step_2_e1_2_cmd=['%s samse -s %s %s %s.end1.sai %s > %s.end1.sam'%(bwa,samseparams,compdb_fasta,sample,fq1path,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
202 post_2_e1_2_clean=['rm -f %s.end1.sai'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
203 step_2_e2_1_cmd=['%s aln %s %s %s > %s.end2.sai'%(bwa,alnparams,compdb_fasta,fq2path,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
204 post_2_e2_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
205 step_2_e2_2_cmd=['%s samse -s %s %s %s.end2.sai %s > %s.end2.sam'%(bwa,samseparams,compdb_fasta,sample,fq2path,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
206 post_2_e2_2_clean=['rm -f %s.end2.sai'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
207 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
208 step_2_e1_1_cmd=['%s aln %s %s %s.end1.fastq > %s.end1.sai'%(bwa,alnparams,compdb_fasta,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
209 post_2_e1_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
210 step_2_e1_2_cmd=['%s samse -s %s %s %s.end1.sai %s.end1.fastq > %s.end1.sam'%(bwa,samseparams,compdb_fasta,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
211 post_2_e1_2_clean=['rm -f %s.end1.sai'%sample,'rm -f %s.end1.fastq'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
212 step_2_e2_1_cmd=['%s aln %s %s %s.end2.fastq > %s.end2.sai'%(bwa,alnparams,compdb_fasta,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
213 post_2_e2_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
214 step_2_e2_2_cmd=['%s samse -s %s %s %s.end2.sai %s.end2.fastq > %s.end2.sam'%(bwa,samseparams,compdb_fasta,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
215 post_2_e2_2_clean=['rm -f %s.end2.sai'%sample,'rm -f %s.end2.fastq'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
216 step_2_e1_3_cmd=['%s view -bS -o %s.end1.bam %s.end1.sam'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
217 post_2_e1_3_clean=['rm -f %s.end1.sam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
218 step_2_e1_4_cmd=['%s sort -n -m 1000000000 %s.end1.bam %s.end1.sorted'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
219 post_2_e1_4_clean=['rm -f %s.end1.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
220 step_2_e2_3_cmd=['%s view -bS -o %s.end2.bam %s.end2.sam'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
221 post_2_e2_3_clean=['rm -f %s.end2.sam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
222 step_2_e2_4_cmd=['%s sort -n -m 1000000000 %s.end2.bam %s.end2.sorted'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
223 post_2_e2_4_clean=['rm -f %s.end2.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
224 step_3_e1_1_cmd=['java -Djava.io.tmpdir=tmp/ -cp %s/RemapAlignments.jar -Xmx8g org.broadinstitute.cga.tools.gatk.rna.RemapAlignments M=%s IN=%s.end1.sorted.bam OUT=%s.end1.remapped.bam R=%s REDUCE=TRUE'%(gatk,compdb_map,sample,sample,genome_fasta)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
225 post_3_e1_1_clean=['rm -f %s.end1.sorted.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
226 step_3_e1_2_cmd=['%s sort -n -m 1000000000 %s.end1.remapped.bam %s.end1.remapped.sorted'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
227 post_3_e1_2_clean=['rm -f %s.end1.remapped.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
228 step_3_e2_1_cmd=['java -Djava.io.tmpdir=tmp/ -cp %s/RemapAlignments.jar -Xmx8g org.broadinstitute.cga.tools.gatk.rna.RemapAlignments M=%s IN=%s.end2.sorted.bam OUT=%s.end2.remapped.bam R=%s REDUCE=TRUE'%(gatk,compdb_map,sample,sample,genome_fasta)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
229 post_3_e2_1_clean=['rm -f %s.end2.sorted.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
230 step_3_e2_2_cmd=['%s sort -n -m 1000000000 %s.end2.remapped.bam %s.end2.remapped.sorted'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
231 post_3_e2_2_clean=['rm -f %s.end2.remapped.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
232 step_4_1_cmd=['java -Djava.io.tmpdir=tmp/ -Xmx8g -jar %s/PairMaker.jar IN1=%s.end1.remapped.sorted.bam IN2=%s.end2.remapped.sorted.bam OUTPUT=%s.paired.bam TMP_DIR=tmp/'%(gatk,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
233 post_4_1_clean=['rm -f %s.end1.remapped.sorted.bam'%sample,'rm -f %s.end2.remapped.sorted.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
234 step_4_2_cmd=['%s sort -m 1000000000 %s.paired.bam %s.paired.sorted'%(samtools,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
235 post_4_2_clean=['rm -f %s.paired.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
236 step_5_cmd=['java -Xmx8g -jar %s/AddOrReplaceReadGroups.jar I=%s.paired.sorted.bam O=%s.withRG.paired.sorted.bam RGLB=%s RGPL=%s RGPU=%s RGSM=%s'%(picard,sample,sample,tag,plat,tag,tag)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
237 post_5_clean=['rm -f %s.paired.sorted.bam'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
238 step_6_1_cmd=['%s index %s.withRG.paired.sorted.bam'%(samtools,sample),'java -Xmx8g -jar %s/GenomeAnalysisTK.jar -l INFO -R %s --default_platform %s --knownSites %s -I %s.withRG.paired.sorted.bam --downsample_to_coverage 10000 -T CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate -nt %s -recalFile %s.orig.csv'%(gatk,genome_fasta,plat,dbsnp_vcf,sample,parallel_n,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
239 post_6_1_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
240 step_6_2_cmd=['java -Xmx8g -jar %s/GenomeAnalysisTK.jar -l INFO -R %s --default_platform %s -I %s.withRG.paired.sorted.bam -T TableRecalibration --out %s.withRG.GATKRecalibrated.bam -recalFile %s.orig.csv'%(gatk,genome_fasta,plat,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
241 post_6_2_clean=['rm -f %s.withRG.paired.sorted.bam'%sample,'rm -f %s.withRG.paired.sorted.bam.bai'%sample,'rm -f %s.orig.csv'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
242 step_7_cmd=['java -Xmx8g -jar %s/MarkDuplicates.jar I=%s.withRG.GATKRecalibrated.bam O=%s.withRG.GATKRecalibrated.flagged.bam METRICS_FILE=%s.Duplicates_metrics.txt VALIDATION_STRINGENCY=SILENT TMP_DIR=tmp/'%(picard,sample,sample,sample),'%s index %s.withRG.GATKRecalibrated.flagged.bam'%(samtools,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
243 post_7_clean=['rm -f %s.withRG.GATKRecalibrated.bam'%sample,'rm -f %s.withRG.GATKRecalibrated.bai'%sample,'rm -f %s.Duplicates_metrics.txt'%sample]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
244 step_8_cmd=["java -Xmx8g -jar %s -ttype 2 -t %s -r %s -s '%s|%s.withRG.GATKRecalibrated.flagged.bam|Disc' -o %s/"%(seqc,genome_gtf,genome_fasta,sample,sample,sample)]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
245 post_8_clean=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
246
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
247 cmdset=[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
248 for item in globals().keys():
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
249 if item.endswith('_cmd'):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
250 cmdset.append(item)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
251 cmdset.sort()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
252
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
253 #########################################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
254 #write PBS file
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
255
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
256 ###############Entry Point
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
257 try:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
258 cmd_entry=cmdset.index('step_'+step+'_cmd') ##entry point
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
259 except ValueError:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
260 sys.exit('ERROR: STEP not recognized')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
261
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
262 def _parsecmd(cmd):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
263 '''pase cmd str for step information'''
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
264 info=cmd.split('_')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
265 cstep='_'.join(info[1:-1])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
266 return cstep
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
267
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
268 ###########################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
269 ##headers
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
270 outfile=open(pbspath,'w')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
271 outfile.write('#! /bin/sh\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
272 outfile.write('#PBS -V\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
273 outfile.write('#PBS -N %s\n'%tag)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
274 outfile.write('#PBS -j oe\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
275 outfile.write('#PBS -o %s\n'%logpath)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
276 for item in refdict['--PBS--']:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
277 outfile.write('#PBS %s %s\n'%(item,refdict['--PBS--'][item]))
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
278 outfile.write('#PBS -d %s\n'%outpath)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
279
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
280 #commands
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
281 outfile.write('echo "Job start: `date`"\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
282
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
283 for i in range(cmd_entry,len(cmdset)):
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
284 cmd=cmdset[i]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
285 cstep=_parsecmd(cmd)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
286 clean_flag=1
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
287 outfile.write('echo "step %s start: `date`"\n'%cstep)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
288 outfile.write('if\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
289 outfile.write('\t%s\n'%('\n'.join(eval(cmd))))
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
290 outfile.write('then\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
291 outfile.write('\techo "step %s done: `date`"\n'%cstep)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
292 outfile.write('else\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
293 outfile.write('\techo "step %s ERROR"\n'%cstep)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
294 outfile.write('\texit\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
295 outfile.write('fi\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
296
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
297 if keepmed_flag=='yes': #if so, keep files by force
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
298 clean_flag=0
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
299
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
300 if clean_flag==1:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
301 outfile.write('\n'.join(eval('post_%s_clean'%cstep)))
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
302 outfile.write('\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
303
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
304 outfile.write('echo "PIPELINE FINISHED"\n')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
305 outfile.close()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
306 ######################################################################
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
307
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
308 if submit in ['False','FALSE','false','F','NO','No','no','n']:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
309 jid='Not_Submitted'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
310 logpath='None'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
311 elif submit in ['True','TRUE','true','T','YES','Yes','yes','y']:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
312 cmdstr='qsub %s'%pbspath
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
313 cmd=cmdstr.split()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
314 cmdout=subprocess.Popen(cmd,stdout=subprocess.PIPE,stderr=subprocess.STDOUT)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
315 jid=cmdout.stdout.read().strip() ##JOB ID
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
316 else:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
317 sys.exit('ERROR: submit parameter not recognized')
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
318
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
319 print '#!#%s'%tag
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
320 print 'BAM\t%s'%bampath
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
321 print 'Entry\t%s'%step
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
322 print 'PBS\t%s'%pbspath
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
323 print 'JOB\t%s'%jid
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
324 print 'LOG\t%s'%logpath
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
325