annotate tools/rgenetics/rgClean.py @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
1 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 # galaxy tool xml files can define a galaxy supplied output filename
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 # that must be passed to the tool and used to return output
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4 # here, the plink log file is copied to that file and removed
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 # took a while to figure this out!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 # use exec_before_job to give files sensible names
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 # ross april 14 2007
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 # plink cleanup script
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 # ross lazarus March 2007 for camp illumina whole genome data
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11 # note problems with multiple commands being ignored - eg --freq --missing --mendel
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 # only the first seems to get done...
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 ##Summary statistics versus inclusion criteria
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 ##
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 ##Feature As summary statistic As inclusion criteria
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 ##Missingness per individual --missing --mind N
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 ##Missingness per marker --missing --geno N
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19 ##Allele frequency --freq --maf N
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 ##Hardy-Weinberg equilibrium --hardy --hwe N
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 ##Mendel error rates --mendel --me N M
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 # call as plinkClean.py $i $o $mind $geno $hwe $maf $mef $mei $outfile
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 # note plinkClean_code.py does some renaming before the job starts
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 <command interpreter="python2.4">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 rgClean.py '$input_file.extra_files_path' '$input_file.metadata.base_name' '$title' '$mind' '$geno' '$hwe' '$maf'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29 '$mef' '$mei' '$out_file1' '$out_file1.files_path' '$userId'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33 import sys,shutil,os,subprocess, glob, string, tempfile, time
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 from rgutils import galhtmlprefix, timenow, plinke
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 prog = os.path.split(sys.argv[0])[-1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 myversion = 'January 4 2010'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 verbose=False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40 def fixoutaff(outpath='',newaff='1'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 """ quick way to create test data sets - set all aff to 1 or 2 for
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 some hapmap data and then merge
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 [rerla@beast galaxy]$ head tool-data/rg/library/pbed/affyHM_CEU.fam
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 1341 14 0 0 2 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 1341 2 13 14 2 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 1341 13 0 0 1 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 1340 9 0 0 1 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 1340 10 0 0 2 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 nchanged = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 fam = '%s.fam' % outpath
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 famf = open(fam,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53 fl = famf.readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 famf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 for i,row in enumerate(fl):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 lrow = row.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 if lrow[-1] <> newaff:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 lrow[-1] = newaff
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 fl[i] = ' '.join(lrow)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 fl[i] += '\n'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 nchanged += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 fo = open(fam,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 fo.write(''.join(fl))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 fo.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65 return nchanged
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 def clean():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 if len(sys.argv) < 16:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 print >> sys.stdout, '## %s expected 12 params in sys.argv, got %d - %s' % (prog,len(sys.argv),sys.argv)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 print >> sys.stdout, """this script will filter a linkage format ped
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 and map file containing genotypes. It takes 14 parameters - the plink --f parameter and"
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 a new filename root for the output clean data followed by the mind,geno,hwe,maf, mef and mei"
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 documented in the plink docs plus the file to be returned to Galaxy
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 called as:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 <command interpreter="python">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 rgClean.py '$input_file.extra_files_path' '$input_file.metadata.base_name' '$title' '$mind'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 '$geno' '$hwe' '$maf' '$mef' '$mei' '$out_file1' '$out_file1.files_path'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 '$relfilter' '$afffilter' '$sexfilter' '$fixaff'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 </command>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87 plog = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 inpath = sys.argv[1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 inbase = sys.argv[2]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 killme = string.punctuation + string.whitespace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 trantab = string.maketrans(killme,'_'*len(killme))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 title = sys.argv[3].translate(trantab)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 mind = sys.argv[4]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 geno = sys.argv[5]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95 hwe = sys.argv[6]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 maf = sys.argv[7]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 me1 = sys.argv[8]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98 me2 = sys.argv[9]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99 outfname = sys.argv[10]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 outfpath = sys.argv[11]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 relf = sys.argv[12]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 afff = sys.argv[13]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 sexf = sys.argv[14]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 fixaff = sys.argv[15]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 output = os.path.join(outfpath,outfname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 outpath = os.path.join(outfpath,title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 outprunepath = os.path.join(outfpath,'ldprune_%s' % title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 os.makedirs(outfpath)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112 bfile = os.path.join(inpath,inbase)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 outf = file(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114 vcl = [plinke,'--noweb','--bfile',bfile,'--make-bed','--out',
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115 outpath,'--set-hh-missing','--mind',mind,
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 '--geno',geno,'--maf',maf,'--hwe',hwe,'--me',me1,me2]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117 # yes - the --me parameter takes 2 values - mendels per snp and per family
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118 if relf == 'oo': # plink filters are what they leave...
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119 vcl.append('--filter-nonfounders') # leave only offspring
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 elif relf == 'fo':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 vcl.append('--filter-founders')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 if afff == 'affonly':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 vcl.append('--filter-controls')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124 elif relf == 'unaffonly':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125 vcl.append('--filter-cases')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 if sexf == 'fsex':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 vcl.append('--filter-females')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 elif relf == 'msex':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 vcl.append('--filter-males')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 p=subprocess.Popen(' '.join(vcl),shell=True,cwd=outfpath)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 retval = p.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132 plog.append('%s started, called as %s' % (prog,' '.join(sys.argv)))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 outf.write(galhtmlprefix % prog)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 outf.write('<ul>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 plogf = '%s.log' % os.path.join(outfpath,title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137 plogl = file(plogf,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 plog += [x.strip() for x in plogl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 plog += ['###Cannot open plink log file %s' % plogf,]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 # if fixaff, want to 'fix' the fam file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142 if fixaff <> '0':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 nchanged = fixoutaff(outpath=outpath,newaff=fixaff)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 plog += ['## fixaff was requested %d subjects affection status changed to %s' % (nchanged,fixaff)]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 pf = file(plogf,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 pf.write('\n'.join(plog))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 pf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 globme = os.path.join(outfpath,'*')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 flist = glob.glob(globme)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150 flist.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151 for i, data in enumerate( flist ):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152 outf.write('<li><a href="%s">%s</a></li>\n' % (os.path.split(data)[-1],os.path.split(data)[-1]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 outf.write('</ul>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 outf.write("</ul></br></div></body></html>")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158 if __name__ == "__main__":
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159 clean()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160