diff --git a/tools/rgenetics/rgManQQ.xml b/tools/rgenetics/rgManQQ.xml index 9de5f162318..f5a9c90fc68 100644 --- a/tools/rgenetics/rgManQQ.xml +++ b/tools/rgenetics/rgManQQ.xml @@ -24,7 +24,7 @@ help='Select "None" if offset not available or no Manhattan plot required' dynamic_options="get_phecols(i,True,'offs')" /> + label="Grey scale for Manhattan plot (default is colour"/> @@ -76,7 +76,7 @@ how busy the server is - be patient please. **Summary** This tool will create a qq plot and a Manhattan plot for one or more GWA P value columns from a tabular -dataset. For Manhatten plots, the data must include the chromosome (eg use 23,24,25 for x,y,mt...) and +dataset. For Manhattan plots, the data must include the chromosome (eg use 23,24,25 for x,y,mt...) and offset. Many analysis files contain the required fields but even without chromosome and offset, a qq plot can be created. diff --git a/tools/rgenetics/rgQC.py b/tools/rgenetics/rgQC.py index 664247c37e9..a5ba96a7626 100644 --- a/tools/rgenetics/rgQC.py +++ b/tools/rgenetics/rgQC.py @@ -984,13 +984,18 @@ def subjectRep(froot='cleantest',outfname="srep",newfpath='.',logf = None): f.close() return res,Tops -def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None ): +def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None,maplist=None ): """by marker (hwe = .hwe, missingness=.lmiss, freq = .frq) keep a list of marker order but keep all stats in dicts write out a fake xls file for R or SAS etc kinda clunky, but.. TODO: ensure stable if any file not found? """ + mapdict = {} + if maplist <> None: + rslist = [x[1] for x in maplist] + offset = [x[3] for x in maplist] + mapdict = dict(zip(rslist,offset)) hwefile = '%s.hwe' % froot lmissfile = '%s.lmiss' % froot freqfile = '%s.frq' % froot @@ -1139,7 +1144,7 @@ def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None ): else: logf.write('No %s file - assuming not family data\n' % lmendfile) # now assemble result list - rhead = ['snp','chrom','maf','a1','a2','missfrac','p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel'] + rhead = ['snp','chromosome','offset','maf','a1','a2','missfrac','p_hwe_all','logp_hwe_all','p_hwe_unaff','logp_hwe_unaff','N_Mendel'] res = [] fres = [] for i in xrange(len(markerlist)): # for each snp in found order @@ -1150,7 +1155,8 @@ def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None ): hwe_all = hwedict[rs].get('ALL',('NA','NA')) # hope this doesn't change... hwe_unaff = hwedict[rs].get('UNAFF',('NA','NA')) nmend = lmenddict.get(rs,'NA') - res.append([rs,chrom,maf,a1,a2,f_missing,hwe_all[0],hwe_all[1],hwe_unaff[0],hwe_unaff[1],nmend]) + offset=mapdict.get(rs,'0') + res.append([rs,chrom,offset,maf,a1,a2,f_missing,hwe_all[0],hwe_all[1],hwe_unaff[0],hwe_unaff[1],nmend]) try: fmaf = '%f' % float(maf) except: @@ -1172,7 +1178,7 @@ def markerRep(froot='cleantest',outfname="mrep",newfpath='.',logf=None ): except: ff_missing = 'NA' #fres.append([rs,chrom,fmaf,a1,a2,ff_missing,hwe_all[0],hwe_all[1],hwe_unaff[0],fhwe,inmend]) - arow = [rs,chrom,fmaf,a1,a2,ff_missing,hwe_all[0],fhweall,hwe_unaff[0],fhweunaff,inmend] + arow = [rs,chrom,offset,fmaf,a1,a2,ff_missing,hwe_all[0],fhweall,hwe_unaff[0],fhweunaff,inmend] fres.append(arow) ntokeep = max(10,len(res)/keepfrac) for i,col in enumerate(Tsorts): @@ -1256,6 +1262,16 @@ if __name__ == "__main__": pass ofn = basename bfn = options.infile + try: + mapf = '%s.bim' % bfn + maplist = file(mapf,'r').readlines() + maplist = [x.split() for x in maplist] + except: + maplist = None + alogf.write('## error - cannot open %s to read map - no offsets will be available for output files') + #rerla@beast galaxy]$ head test-data/tinywga.bim + #22 rs2283802 0 21784722 4 2 + #22 rs2267000 0 21785366 4 2 rgbin = os.path.split(rexe)[0] # get our rg bin path #plinktasks = [' --freq',' --missing',' --mendel',' --hardy',' --check-sex'] # plink v1 fixes that bug! # if we could, do all at once? Nope. Probably never. @@ -1276,7 +1292,7 @@ if __name__ == "__main__": subjects,subjectTops = subjectRep(froot=repout,outfname=asubjf,newfpath=newfpath, logf=alogf) # writes the subject_froot.xls file markers,markerTops = markerRep(froot=repout,outfname=amarkf,newfpath=newfpath, - logf=alogf) # marker_froot.xls + logf=alogf,maplist=maplist) # marker_froot.xls nbreaks = 100 s = '## starting plotpage, newfpath=%s,m=%s,s=%s/n' % (newfpath,markers[:2],subjects[:2]) alogf.write(s)