[2] | 1 | """ |
---|
| 2 | released under the terms of the LGPL |
---|
| 3 | copyright ross lazarus August 2007 |
---|
| 4 | for the rgenetics project |
---|
| 5 | |
---|
| 6 | Special galaxy tool for the camp2007 data |
---|
| 7 | Allows grabbing genotypes from an arbitrary region and estimating |
---|
| 8 | ld using haploview |
---|
| 9 | |
---|
| 10 | stoopid haploview won't allow control of dest directory for plots - always end |
---|
| 11 | up where the data came from - need to futz to get it where it belongs |
---|
| 12 | |
---|
| 13 | Needs a mongo results file in the location hardwired below or could be passed in as |
---|
| 14 | a library parameter - but this file must have a very specific structure |
---|
| 15 | rs chrom offset float1...floatn |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | """ |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | import sys, array, os, string, tempfile, shutil, subprocess, glob |
---|
| 22 | from rgutils import galhtmlprefix |
---|
| 23 | |
---|
| 24 | progname = os.path.split(sys.argv[0])[1] |
---|
| 25 | |
---|
| 26 | javabin = 'java' |
---|
| 27 | #hvbin = '/usr/local/bin/Haploview.jar' |
---|
| 28 | #hvbin = '/home/universe/linux-i686/haploview/Haploview.jar' |
---|
| 29 | # get this from tool as a parameter - can use |
---|
| 30 | |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | atrandic = {'A':'1','C':'2','G':'3','T':'4','N':'0','-':'0','1':'1','2':'2','3':'3','4':'4','0':'0'} |
---|
| 34 | |
---|
| 35 | class NullDevice: |
---|
| 36 | """ """ |
---|
| 37 | def write(self, s): |
---|
| 38 | pass |
---|
| 39 | |
---|
| 40 | def ld(): |
---|
| 41 | """ ### Sanity check the arguments |
---|
| 42 | |
---|
| 43 | <command interpreter="python"> |
---|
| 44 | rgHaploView.py "$ucsc_region" "$rslist" "$title" "$out_file1" |
---|
| 45 | "$lhistIn.extra_files_path" "$lhistIn.metadata.base_name" |
---|
| 46 | "$minmaf" "$maxdist" "$ldtype" "$hires" "$memsize" "$out_file1.files_path" |
---|
| 47 | "$infoTrack" "$tagr2" "$hmpanel" ${GALAXY_DATA_INDEX_DIR}/rg/bin/haploview.jar |
---|
| 48 | </command> |
---|
| 49 | |
---|
| 50 | remember to figure out chromosome and complain if > 1? |
---|
| 51 | and use the -chromosome <1-22,X,Y> parameter to haploview |
---|
| 52 | skipcheck? |
---|
| 53 | """ |
---|
| 54 | progname = os.path.split(sys.argv[0])[-1] |
---|
| 55 | ts = '%s%s' % (string.punctuation,string.whitespace) |
---|
| 56 | ptran = string.maketrans(ts,'_'*len(ts)) |
---|
| 57 | if len(sys.argv) < 16: |
---|
| 58 | s = '##!%s: Expected 16 params in sys.argv, got %d (%s)' % (progname,len(sys.argv), sys.argv) |
---|
| 59 | print s |
---|
| 60 | sys.exit(1) |
---|
| 61 | |
---|
| 62 | ### Figure out what genomic region we are interested in |
---|
| 63 | region = sys.argv[1] |
---|
| 64 | orslist = sys.argv[2].replace('X',' ').lower() # galaxy replaces newlines with XX - go figure |
---|
| 65 | title = sys.argv[3].translate(ptran) |
---|
| 66 | # for outputs |
---|
| 67 | outfile = sys.argv[4] |
---|
| 68 | logfn = 'Log_%s.txt' % (title) |
---|
| 69 | histextra = sys.argv[5] |
---|
| 70 | base_name = sys.argv[6] |
---|
| 71 | pedFileBase = os.path.join(histextra,base_name) |
---|
| 72 | print 'pedfilebase=%s' % pedFileBase |
---|
| 73 | minMaf=sys.argv[7] |
---|
| 74 | if minMaf: |
---|
| 75 | try: |
---|
| 76 | minMaf = float(minMaf) |
---|
| 77 | except: |
---|
| 78 | minMaf = 0.0 |
---|
| 79 | maxDist=sys.argv[8] or None |
---|
| 80 | ldType=sys.argv[9] or 'RSQ' |
---|
| 81 | hiRes = (sys.argv[10].lower() == 'hi') |
---|
| 82 | memSize= sys.argv[11] or '1000' |
---|
| 83 | memSize = int(memSize) |
---|
| 84 | outfpath = sys.argv[12] |
---|
| 85 | infotrack = False # note that otherwise this breaks haploview in headless mode |
---|
| 86 | #infotrack = sys.argv[13] == 'info' |
---|
| 87 | # this fails in headless mode as at april 2010 with haploview 4.2 |
---|
| 88 | tagr2 = sys.argv[14] or '0.8' |
---|
| 89 | hmpanels = sys.argv[15] # eg "['CEU','YRI']" |
---|
| 90 | if hmpanels: |
---|
| 91 | hmpanels = hmpanels.replace('[','') |
---|
| 92 | hmpanels = hmpanels.replace(']','') |
---|
| 93 | hmpanels = hmpanels.replace("'",'') |
---|
| 94 | hmpanels = hmpanels.split(',') |
---|
| 95 | hvbin = sys.argv[16] # added rml june 2008 |
---|
| 96 | bindir = os.path.split(hvbin)[0] |
---|
| 97 | # jan 2010 - always assume utes are on path to avoid platform problems |
---|
| 98 | pdfjoin = 'pdfjoin' # os.path.join(bindir,'pdfjoin') |
---|
| 99 | pdfnup = 'pdfnup' # os.path.join(bindir,'pdfnup') |
---|
| 100 | mogrify = 'mogrify' # os.path.join(bindir,'mogrify') |
---|
| 101 | convert = 'convert' # os.path.join(bindir,'convert') |
---|
| 102 | log_file = os.path.join(outfpath,logfn) |
---|
| 103 | MAP_FILE = '%s.map' % pedFileBase |
---|
| 104 | DATA_FILE = '%s.ped' % pedFileBase |
---|
| 105 | try: |
---|
| 106 | os.makedirs(outfpath) |
---|
| 107 | s = '## made new path %s\n' % outfpath |
---|
| 108 | except: |
---|
| 109 | pass |
---|
| 110 | lf = file(log_file,'w') |
---|
| 111 | s = 'PATH=%s\n' % os.environ.get('PATH','?') |
---|
| 112 | lf.write(s) |
---|
| 113 | hlogf = os.path.join(outfpath,'%s.log' % pedFileBase) |
---|
| 114 | chromosome = '' |
---|
| 115 | spos = epos = -9 |
---|
| 116 | rslist = [] |
---|
| 117 | rsdict = {} |
---|
| 118 | hlog = [] |
---|
| 119 | |
---|
| 120 | if region > '': |
---|
| 121 | useRs = [] |
---|
| 122 | useRsdict={} |
---|
| 123 | try: # TODO make a regexp? |
---|
| 124 | c,rest = region.split(':') |
---|
| 125 | chromosome = c.replace('chr','') |
---|
| 126 | rest = rest.replace(',','') # remove commas |
---|
| 127 | spos,epos = rest.split('-') |
---|
| 128 | spos = int(spos) |
---|
| 129 | epos = int(epos) |
---|
| 130 | s = '## %s parsing chrom %s from %d to %d\n' % (progname,chromosome,spos,epos) |
---|
| 131 | lf.write(s) |
---|
| 132 | lf.write('\n') |
---|
| 133 | print >> sys.stdout, s |
---|
| 134 | except: |
---|
| 135 | s = '##! %s unable to parse region %s - MUST look like "chr8:10,000-100,000\n' % (progname,region) |
---|
| 136 | print >> sys.stdout, s |
---|
| 137 | lf.write(s) |
---|
| 138 | lf.write('\n') |
---|
| 139 | lf.close() |
---|
| 140 | sys.exit(1) |
---|
| 141 | else: |
---|
| 142 | useRs = orslist.split() # galaxy replaces newlines with XX - go figure |
---|
| 143 | useRsdict = dict(zip(useRs,useRs)) |
---|
| 144 | useTemp = False |
---|
| 145 | try: |
---|
| 146 | dfile = open(DATA_FILE, 'r') |
---|
| 147 | except: # bad input file name? |
---|
| 148 | s = '##! RGeno unable to open file %s\n' % (DATA_FILE) |
---|
| 149 | lf.write(s) |
---|
| 150 | lf.write('\n') |
---|
| 151 | lf.close() |
---|
| 152 | print >> sys.stdout, s |
---|
| 153 | raise |
---|
| 154 | sys.exit(1) |
---|
| 155 | try: |
---|
| 156 | mfile = open(MAP_FILE, 'r') |
---|
| 157 | except: # bad input file name? |
---|
| 158 | s = '##! RGeno unable to open file %s' % (MAP_FILE) |
---|
| 159 | lf.write(s) |
---|
| 160 | lf.write('\n') |
---|
| 161 | lf.close() |
---|
| 162 | print >> sys.stdout, s |
---|
| 163 | raise |
---|
| 164 | sys.exit(1) |
---|
| 165 | if len(useRs) > 0 or spos <> -9 : # subset region |
---|
| 166 | useTemp = True |
---|
| 167 | ### Figure out which markers are in this region |
---|
| 168 | markers = [] |
---|
| 169 | snpcols = {} |
---|
| 170 | chroms = {} |
---|
| 171 | minpos = 2**32 |
---|
| 172 | maxpos = 0 |
---|
| 173 | for lnum,row in enumerate(mfile): |
---|
| 174 | line = row.strip() |
---|
| 175 | if not line: continue |
---|
| 176 | chrom, snp, genpos, abspos = line.split() |
---|
| 177 | try: |
---|
| 178 | ic = int(chrom) |
---|
| 179 | except: |
---|
| 180 | ic = None |
---|
| 181 | if ic and ic <= 23: |
---|
| 182 | try: |
---|
| 183 | abspos = int(abspos) |
---|
| 184 | if abspos > maxpos: |
---|
| 185 | maxpos = abspos |
---|
| 186 | if abspos < minpos: |
---|
| 187 | minpos = abspos |
---|
| 188 | except: |
---|
| 189 | abspos = epos + 999999999 # so next test fails |
---|
| 190 | if useRsdict.get(snp,None) or (spos <> -9 and chrom == chromosome and (spos <= abspos <= epos)): |
---|
| 191 | if chromosome == '': |
---|
| 192 | chromosome = chrom |
---|
| 193 | chroms.setdefault(chrom,chrom) |
---|
| 194 | markers.append((chrom,abspos,snp)) # decorate for sort into genomic |
---|
| 195 | snpcols[snp] = lnum # so we know which col to find genos for this marker |
---|
| 196 | markers.sort() |
---|
| 197 | rslist = [x[2] for x in markers] # drop decoration |
---|
| 198 | rsdict = dict(zip(rslist,rslist)) |
---|
| 199 | if len(rslist) == 0: |
---|
| 200 | s = '##! %s found no rs numbers in %s' % (progname,sys.argv[1:3]) |
---|
| 201 | lf.write(s) |
---|
| 202 | lf.write('\n') |
---|
| 203 | lf.close() |
---|
| 204 | print >> sys.stdout, s |
---|
| 205 | sys.exit(1) |
---|
| 206 | if spos == -9: |
---|
| 207 | spos = minpos |
---|
| 208 | epos = maxpos |
---|
| 209 | s = '## %s looking for %d rs (%s)' % (progname,len(rslist),rslist[:5]) |
---|
| 210 | lf.write(s) |
---|
| 211 | print >> sys.stdout, s |
---|
| 212 | wewant = [(6+(2*snpcols[x])) for x in rslist] # |
---|
| 213 | # column indices of first geno of each marker pair to get the markers into genomic |
---|
| 214 | ### ... and then parse the rest of the ped file to pull out |
---|
| 215 | ### the genotypes for all subjects for those markers |
---|
| 216 | # /usr/local/galaxy/data/rg/1/lped/ |
---|
| 217 | tempMapName = os.path.join(outfpath,'%s.info' % title) |
---|
| 218 | tempMap = file(tempMapName,'w') |
---|
| 219 | tempPedName = os.path.join(outfpath,'%s.ped' % title) |
---|
| 220 | tempPed = file(tempPedName,'w') |
---|
| 221 | pngpath = '%s.LD.PNG' % tempPedName |
---|
| 222 | map = ['%s\t%s' % (x[2],x[1]) for x in markers] # snp,abspos in genomic order for haploview |
---|
| 223 | tempMap.write('%s\n' % '\n'.join(map)) |
---|
| 224 | tempMap.close() |
---|
| 225 | nrows = 0 |
---|
| 226 | for line in dfile: |
---|
| 227 | line = line.strip() |
---|
| 228 | if not line: |
---|
| 229 | continue |
---|
| 230 | fields = line.split() |
---|
| 231 | preamble = fields[:6] |
---|
| 232 | g = ['%s %s' % (fields[snpcol], fields[snpcol+1]) for snpcol in wewant] |
---|
| 233 | g = ' '.join(g) |
---|
| 234 | g = g.split() # we'll get there |
---|
| 235 | g = [atrandic.get(x,'0') for x in g] # numeric alleles... |
---|
| 236 | tempPed.write('%s %s\n' % (' '.join(preamble), ' '.join(g))) |
---|
| 237 | nrows += 1 |
---|
| 238 | tempPed.close() |
---|
| 239 | s = '## %s: wrote %d markers, %d subjects for region %s\n' % (progname,len(rslist),nrows,region) |
---|
| 240 | lf.write(s) |
---|
| 241 | lf.write('\n') |
---|
| 242 | print >> sys.stdout,s |
---|
| 243 | else: # even if using all, must set up haploview info file instead of map |
---|
| 244 | markers = [] |
---|
| 245 | chroms = {} |
---|
| 246 | spos = sys.maxint |
---|
| 247 | epos = -spos |
---|
| 248 | for lnum,row in enumerate(mfile): |
---|
| 249 | line = row.strip() |
---|
| 250 | if not line: continue |
---|
| 251 | chrom, snp, genpos, abspos = line.split() |
---|
| 252 | try: |
---|
| 253 | ic = int(chrom) |
---|
| 254 | except: |
---|
| 255 | ic = None |
---|
| 256 | if ic and ic <= 23: |
---|
| 257 | if chromosome == '': |
---|
| 258 | chromosome = chrom |
---|
| 259 | chroms.setdefault(chrom,chrom) |
---|
| 260 | try: |
---|
| 261 | p = int(abspos) |
---|
| 262 | if p < spos and p <> 0: |
---|
| 263 | spos = p |
---|
| 264 | if p > epos and p <> 0: |
---|
| 265 | epos = p |
---|
| 266 | except: |
---|
| 267 | pass |
---|
| 268 | markers.append('%s %s' % (snp,abspos)) # no sort - pass |
---|
| 269 | # now have spos and epos for hapmap if hmpanels |
---|
| 270 | tempMapName = os.path.join(outfpath,'%s.info' % title) |
---|
| 271 | tempMap = file(tempMapName,'w') |
---|
| 272 | tempMap.write('\n'.join(markers)) |
---|
| 273 | tempMap.close() |
---|
| 274 | tempPedName = os.path.join(outfpath,'%s.ped' % title) |
---|
| 275 | try: # will fail on winblows! |
---|
| 276 | os.symlink(DATA_FILE,tempPedName) |
---|
| 277 | except: |
---|
| 278 | shutil.copy(DATA_FILE,tempPedName) # wasteful but.. |
---|
| 279 | nchroms = len(chroms) # if > 1 can't really do this safely |
---|
| 280 | if nchroms > 1: |
---|
| 281 | s = '## warning - multiple chromosomes found in your map file - %s\n' % ','.join(chroms.keys()) |
---|
| 282 | lf.write(s) |
---|
| 283 | print >> sys.stdout,s |
---|
| 284 | sys.exit(1) |
---|
| 285 | DATA_FILE = tempPedName # for haploview |
---|
| 286 | INFO_FILE = tempMapName |
---|
| 287 | dfile.close() |
---|
| 288 | mfile.close() |
---|
| 289 | fblog,blog = tempfile.mkstemp() |
---|
| 290 | ste = open(blog,'w') # to catch the blather |
---|
| 291 | # if no need to rewrite - set up names for haploview call |
---|
| 292 | vcl = [javabin,'-jar',hvbin,'-n','-memory','%d' % memSize,'-pairwiseTagging', |
---|
| 293 | '-pedfile',DATA_FILE,'-info',INFO_FILE,'-tagrsqcounts', |
---|
| 294 | '-tagrsqcutoff',tagr2, |
---|
| 295 | '-ldcolorscheme',ldType] |
---|
| 296 | if minMaf: |
---|
| 297 | vcl += ['-minMaf','%f' % minMaf] |
---|
| 298 | if maxDist: |
---|
| 299 | vcl += ['-maxDistance',maxDist] |
---|
| 300 | if hiRes: |
---|
| 301 | vcl.append('-png') |
---|
| 302 | else: |
---|
| 303 | vcl.append('-compressedpng') |
---|
| 304 | if nchroms == 1: |
---|
| 305 | vcl += ['-chromosome',chromosome] |
---|
| 306 | if infotrack: |
---|
| 307 | vcl.append('-infoTrack') |
---|
| 308 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=ste,stdout=lf) |
---|
| 309 | retval = p.wait() |
---|
| 310 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 311 | lf.write(s) |
---|
| 312 | vcl = [mogrify, '-resize 800x400!', '*.PNG'] |
---|
| 313 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 314 | retval = p.wait() |
---|
| 315 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 316 | lf.write(s) |
---|
| 317 | inpng = '%s.LD.PNG' % DATA_FILE # stupid but necessary - can't control haploview name mangle |
---|
| 318 | inpng = inpng.replace(' ','') |
---|
| 319 | inpng = os.path.split(inpng)[-1] |
---|
| 320 | tmppng = '%s.tmp.png' % title |
---|
| 321 | tmppng = tmppng.replace(' ','') |
---|
| 322 | outpng = '1_%s.png' % title |
---|
| 323 | outpng = outpng.replace(' ','') |
---|
| 324 | outpng = os.path.split(outpng)[-1] |
---|
| 325 | vcl = [convert, '-resize 800x400!', inpng, tmppng] |
---|
| 326 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 327 | retval = p.wait() |
---|
| 328 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 329 | lf.write(s) |
---|
| 330 | s = "text 10,300 '%s'" % title[:40] |
---|
| 331 | vcl = ['convert', '-pointsize 25','-fill maroon', |
---|
| 332 | '-draw "%s"' % s, tmppng, outpng] |
---|
| 333 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 334 | retval = p.wait() |
---|
| 335 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 336 | lf.write(s) |
---|
| 337 | try: |
---|
| 338 | os.remove(os.path.join(outfpath,tmppng)) |
---|
| 339 | except: |
---|
| 340 | pass # label all the plots then delete all the .PNG files before munging |
---|
| 341 | fnum=1 |
---|
| 342 | if hmpanels: |
---|
| 343 | sp = '%d' % (spos/1000.) # hapmap wants kb |
---|
| 344 | ep = '%d' % (epos/1000.) |
---|
| 345 | for panel in hmpanels: |
---|
| 346 | if panel > '' and panel.lower() <> 'none': # in case someone checks that option too :) |
---|
| 347 | ptran = panel.strip() |
---|
| 348 | ptran = ptran.replace('+','_') |
---|
| 349 | fnum += 1 # preserve an order or else we get sorted |
---|
| 350 | vcl = [javabin,'-jar',hvbin,'-n','-memory','%d' % memSize, |
---|
| 351 | '-chromosome',chromosome, '-panel',panel.strip(), |
---|
| 352 | '-hapmapDownload','-startpos',sp,'-endpos',ep, |
---|
| 353 | '-ldcolorscheme',ldType] |
---|
| 354 | if minMaf: |
---|
| 355 | vcl += ['-minMaf','%f' % minMaf] |
---|
| 356 | if maxDist: |
---|
| 357 | vcl += ['-maxDistance',maxDist] |
---|
| 358 | if hiRes: |
---|
| 359 | vcl.append('-png') |
---|
| 360 | else: |
---|
| 361 | vcl.append('-compressedpng') |
---|
| 362 | if infotrack: |
---|
| 363 | vcl.append('-infoTrack') |
---|
| 364 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=ste,stdout=lf) |
---|
| 365 | retval = p.wait() |
---|
| 366 | inpng = 'Chromosome%s%s.LD.PNG' % (chromosome,panel) |
---|
| 367 | inpng = inpng.replace(' ','') # mysterious spaces! |
---|
| 368 | outpng = '%d_HapMap_%s_%s.png' % (fnum,ptran,chromosome,) |
---|
| 369 | # hack for stupid chb+jpt |
---|
| 370 | outpng = outpng.replace(' ','') |
---|
| 371 | tmppng = '%s.tmp.png' % title |
---|
| 372 | tmppng = tmppng.replace(' ','') |
---|
| 373 | outpng = os.path.split(outpng)[-1] |
---|
| 374 | vcl = [convert, '-resize 800x400!', inpng, tmppng] |
---|
| 375 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 376 | retval = p.wait() |
---|
| 377 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 378 | lf.write(s) |
---|
| 379 | s = "text 10,300 'HapMap %s'" % ptran.strip() |
---|
| 380 | vcl = [convert, '-pointsize 25','-fill maroon', |
---|
| 381 | '-draw "%s"' % s, tmppng, outpng] |
---|
| 382 | p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 383 | retval = p.wait() |
---|
| 384 | s = '## executing %s returned %d\n' % (' '.join(vcl),retval) |
---|
| 385 | lf.write(s) |
---|
| 386 | try: |
---|
| 387 | os.remove(os.path.join(outfpath,tmppng)) |
---|
| 388 | except: |
---|
| 389 | pass |
---|
| 390 | nimages = len(glob.glob(os.path.join(outfpath,'*.png'))) # rely on HaploView shouting - PNG @! |
---|
| 391 | lf.write('### nimages=%d\n' % nimages) |
---|
| 392 | if nimages > 0: # haploview may fail? |
---|
| 393 | vcl = '%s -format pdf -resize 800x400! *.png' % mogrify |
---|
| 394 | p=subprocess.Popen(vcl,shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 395 | retval = p.wait() |
---|
| 396 | lf.write('## executing %s returned %d\n' % (vcl,retval)) |
---|
| 397 | vcl = '%s "*.pdf" --fitpaper true --outfile alljoin.pdf' % pdfjoin |
---|
| 398 | p=subprocess.Popen(vcl,shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 399 | retval = p.wait() |
---|
| 400 | lf.write('## executing %s returned %d\n' % (vcl,retval)) |
---|
| 401 | vcl = '%s alljoin.pdf --nup 1x%d --outfile allnup.pdf' % (pdfnup,nimages) |
---|
| 402 | p=subprocess.Popen(vcl,shell=True,cwd=outfpath,stderr=lf,stdout=lf) |
---|
| 403 | retval = p.wait() |
---|
| 404 | lf.write('## executing %s returned %d\n' % (vcl,retval)) |
---|
| 405 | #vcl = ['convert', 'allnup.pdf', 'allnup.png'] # this fails - bad pdf? |
---|
| 406 | #p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath) |
---|
| 407 | #retval = p.wait() |
---|
| 408 | #lf.write('## executing %s returned %d\n' % (vcl,retval)) |
---|
| 409 | lf.write('\n'.join(hlog)) |
---|
| 410 | ste.close() # temp file used to catch haploview blather |
---|
| 411 | hblather = open(blog,'r').readlines() # to catch the blather |
---|
| 412 | os.unlink(blog) |
---|
| 413 | if len(hblather) > 0: |
---|
| 414 | lf.write('## In addition, Haploview complained:') |
---|
| 415 | lf.write(''.join(hblather)) |
---|
| 416 | lf.write('\n') |
---|
| 417 | lf.close() |
---|
| 418 | flist = glob.glob(os.path.join(outfpath, '*')) |
---|
| 419 | flist.sort() |
---|
| 420 | ts = '!"#$%&\'()*+,-/:;<=>?@[\\]^_`{|}~' + string.whitespace |
---|
| 421 | ftran = string.maketrans(ts,'_'*len(ts)) |
---|
| 422 | outf = file(outfile,'w') |
---|
| 423 | outf.write(galhtmlprefix % progname) |
---|
| 424 | s = '<h4>rgenetics for Galaxy %s, wrapping HaploView</h4>' % (progname) |
---|
| 425 | outf.write(s) |
---|
| 426 | """ |
---|
| 427 | as at ashg 2009, convert fails on allnup.pdf - probably too complex... |
---|
| 428 | mainthumb = 'allnup.png' |
---|
| 429 | mainpdf = 'allnup.pdf' |
---|
| 430 | if os.path.exists(mainpdf): |
---|
| 431 | if not os.path.exists(mainthumb): |
---|
| 432 | outf.write('<table><tr><td colspan="3"><a href="%s">Main combined LD plot</a></td></tr></table>\n' % (mainpdf)) |
---|
| 433 | else: |
---|
| 434 | outf.write('<table><tr><td><a href="%s"><img src="%s" alt="Main combined LD image" hspace="10" align="middle">') |
---|
| 435 | outf.write('</td><td>Click this thumbnail to display the main combined LD image</td></tr></table>\n' % (mainpdf,mainthumb)) |
---|
| 436 | else: |
---|
| 437 | outf.write('(No main image was generated - this usually means a Haploview error connecting to Hapmap site - please try later)<br/>\n') |
---|
| 438 | outf.write('## Called as %s' % sys.argv) |
---|
| 439 | """ |
---|
| 440 | outf.write('<br><div><hr><ul>\n') |
---|
| 441 | for i, data in enumerate( flist ): |
---|
| 442 | dn = os.path.split(data)[-1] |
---|
| 443 | if dn[:3] <> 'all': |
---|
| 444 | continue |
---|
| 445 | newdn = dn.translate(ftran) |
---|
| 446 | if dn <> newdn: |
---|
| 447 | os.rename(os.path.join(outfpath,dn),os.path.join(outfpath,newdn)) |
---|
| 448 | dn = newdn |
---|
| 449 | dnlabel = dn |
---|
| 450 | ext = dn.split('.')[-1] |
---|
| 451 | if dn == 'allnup.pdf': |
---|
| 452 | dnlabel = 'All pdf plots on a single page' |
---|
| 453 | elif dn == 'alljoin.pdf': |
---|
| 454 | dnlabel = 'All pdf plots, each on a separate page' |
---|
| 455 | outf.write('<li><a href="%s">%s - %s</a></li>\n' % (dn,dn,dnlabel)) |
---|
| 456 | for i, data in enumerate( flist ): |
---|
| 457 | dn = os.path.split(data)[-1] |
---|
| 458 | if dn[:3] == 'all': |
---|
| 459 | continue |
---|
| 460 | newdn = dn.translate(ftran) |
---|
| 461 | if dn <> newdn: |
---|
| 462 | os.rename(os.path.join(outfpath,dn),os.path.join(outfpath,newdn)) |
---|
| 463 | dn = newdn |
---|
| 464 | dnlabel = dn |
---|
| 465 | ext = dn.split('.')[-1] |
---|
| 466 | if dn == 'allnup.pdf': |
---|
| 467 | dnlabel = 'All pdf plots on a single page' |
---|
| 468 | elif dn == 'alljoin.pdf': |
---|
| 469 | dnlabel = 'All pdf plots, each on a separate page' |
---|
| 470 | elif ext == 'info': |
---|
| 471 | dnlabel = '%s map data for Haploview input' % title |
---|
| 472 | elif ext == 'ped': |
---|
| 473 | dnlabel = '%s genotype data for Haploview input' % title |
---|
| 474 | elif dn.find('CEU') <> -1 or dn.find('YRI') <> -1 or dn.find('CHB_JPT') <> -1: # is hapmap |
---|
| 475 | dnlabel = 'Hapmap data' |
---|
| 476 | if ext == 'TAGS' or ext == 'TESTS' or ext == 'CHAPS': |
---|
| 477 | dnlabel = dnlabel + ' Tagger output' |
---|
| 478 | outf.write('<li><a href="%s">%s - %s</a></li>\n' % (dn,dn,dnlabel)) |
---|
| 479 | outf.write('</ol><br>') |
---|
| 480 | outf.write("</div><div><hr>Job Log follows below (see %s)<pre>" % logfn) |
---|
| 481 | s = file(log_file,'r').readlines() |
---|
| 482 | s = '\n'.join(s) |
---|
| 483 | outf.write('%s</pre><hr></div>' % s) |
---|
| 484 | outf.write('</body></html>') |
---|
| 485 | outf.close() |
---|
| 486 | if useTemp: |
---|
| 487 | try: |
---|
| 488 | os.unlink(tempMapName) |
---|
| 489 | os.unlink(tempPedName) |
---|
| 490 | except: |
---|
| 491 | pass |
---|
| 492 | |
---|
| 493 | if __name__ == "__main__": |
---|
| 494 | ld() |
---|
| 495 | |
---|