Mercurial > repos > xuebing > sharplabtool

comparison tools/vcf_tools/tools.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 #!/usr/bin/python
2
3 import os.path
4 import sys
5 import vcfPytools
6 from vcfPytools import __version__
7
8 # Determine whether to output to a file or stdout.
9 def setOutput(output):
10 if output == None:
11 outputFile = sys.stdout
12 writeOut = False
13 else:
14 output = os.path.abspath(output)
15 outputFile = open(output, 'w')
16 writeOut = True
17
18 return outputFile, writeOut
19
20 # Determine which file has priority for writing out records.
21 def setVcfPriority(priorityFile, vcfFiles):
22 if priorityFile == None: priority = 0
23 elif priorityFile == vcfFiles[0]: priority = 1
24 elif priorityFile == vcfFiles[1]: priority = 2
25 elif priorityFile.lower() == "merge": priority = 3
26 else:
27 print >> sys.stderr, "vcf file give priority must be one of the two input vcf files or merge."
28 exit(1)
29
30 return priority
31
32 # If the union or intersection of two vcf files is being performed
33 # and the output vcf file is to contain the information from both
34 # files, the headers need to be merged to ensure that all info and
35 # format entries have an explanation.
36 def mergeHeaders(v1, v2, v3):
37
38 # If either file does not have a header, terminate the program.
39 # In order to merge the headers, the different fields must be
40 # checked to ensure the files are compatible.
41 if not v1.hasHeader or not v2.hasHeader:
42 print >> sys.stderr, "Both vcf files must have a header in order to merge data sets."
43 exit(1)
44
45 v3.infoHeaderTags = v1.infoHeaderTags.copy()
46 v3.formatHeaderTags = v1.formatHeaderTags.copy()
47 v3.numberDataSets = v1.numberDataSets
48 v3.includedDataSets = v1.includedDataSets.copy()
49 v3.headerText = v1.headerText
50 v3.headerTitles = v1.headerTitles
51 v3.infoHeaderString = v1.infoHeaderString.copy()
52 v3.formatHeaderString = v1.formatHeaderString.copy()
53
54 # Merge the info field descriptions.
55 for tag in v2.infoHeaderTags:
56 if v1.infoHeaderTags.has_key(tag):
57 if v1.infoHeaderTags[tag][0] != v2.infoHeaderTags[tag][0] or \
58 v1.infoHeaderTags[tag][1] != v2.infoHeaderTags[tag][1]:
59 print v1.infoHeaderTags[tag][0]
60 print v1.infoHeaderTags[tag][1]
61 print v1.infoHeaderTags[tag][2]
62 print >> sys.stderr, "Input vcf files have different definitions for " + tag + " field."
63 exit(1)
64 else: v3.infoHeaderTags[tag] = v2.infoHeaderTags[tag]
65
66 # Merge the format field descriptions.
67 for tag in v2.formatHeaderTags:
68 if v1.formatHeaderTags.has_key(tag):
69 if v1.formatHeaderTags[tag][0] != v2.formatHeaderTags[tag][0] or \
70 v1.formatHeaderTags[tag][1] != v2.formatHeaderTags[tag][1]:
71 print >> sys.stderr, "Input vcf files have different definitions for " + tag + " field."
72 exit(1)
73 else: v3.formatHeaderTags[tag] = v2.formatHeaderTags[tag]
74
75 # Now check to see if the vcf files contain information from multiple
76 # records themselves and create an ordered list in which the data
77 # will appear in the file. For instance, of the first file has
78 # already got two sets of data and is being intersected with a file
79 # with one set of data, the order of data in the new vcf file will be
80 # the two sets from the first file followed by the second, e.g.
81 # AB=3/2/4, where the 3 and 2 are from the first file and the 4 is the
82 # value of AC from the second vcf. The header will have a ##FILE for
83 # each of the three files, so the origin if the data can be recovered.
84 if v1.numberDataSets == 0:
85 v3.includedDataSets[v3.numberDataSets + 1] = v1.filename
86 v3.numberDataSets += 1
87 if v2.numberDataSets == 0:
88 v3.includedDataSets[v3.numberDataSets + 1] = v2.filename
89 v3.numberDataSets += 1
90 else:
91 for i in range(1, v2.numberDataSets + 1):
92 v3.includedDataSets[v3.numberDataSets + 1] = v2.includedDataSets[i]
93 v3.numberDataSets += 1
94
95 # If either of the input files contain multiple data sets (e.g. multiple
96 # vcf files have undergone intersection or union calculations and all
97 # information has been retained) and the priority isn't set to 'merge',
98 # terminate the program. This is to ensure that the origin of the data
99 # doesn't get confused.
100 def checkDataSets(v1, v2):
101 if v1.numberDataSets + v2.numberDataSets != 0:
102 print >> sys.stderr, "\nERROR:"
103 print >> sys.stderr, "input vcf file(s) contain data sets from multiple vcf files."
104 print >> sys.stderr, "Further intersection or union operations must include --priority-file merge"
105 print >> sys.stderr, "Other tools may be incompatible with this format."
106 exit(1)
107
108 # Write the header to file.
109 def writeHeader (outputFile, v, removeGenotypes, taskDescriptor):
110 if not v.hasHeader:
111 v.headerText = "##fileformat=VCFv4.0\n##source=vcfPytools " + __version__ + "\n"
112 v.headerTitles = "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"
113 outputFile.write(v.headerText) if v.headerText != "" else None
114 print >> outputFile, taskDescriptor
115 for tag in v.infoHeaderString: print >> outputFile, v.infoHeaderString[tag]
116 for tag in v.formatHeaderString: print >> outputFile, v.formatHeaderString[tag]
117
118 # Write out a list of files indicating which data set belongs to which file.
119 if v.numberDataSets != 0:
120 for i in range(1, v.numberDataSets + 1):
121 print >> outputFile, "##FILE=<ID=" + str(i) + ",\"" + v.includedDataSets[i] + "\">"
122
123 if removeGenotypes:
124 line = v.headerTitles.rstrip("\n").split("\t")
125 newHeaderTitles = line[0]
126 for i in range(1,8):
127 newHeaderTitles = newHeaderTitles + "\t" + line[i]
128 newHeaderTitles = newHeaderTitles + "\n"
129 outputFile.write( newHeaderTitles )
130 else:
131 outputFile.write( v.headerTitles )
132
133 # Check that the two reference sequence lists are identical.
134 # If there are a different number or order, the results may
135 # not be as expected.
136 def checkReferenceSequenceLists(list1, list2):
137 errorMessage = False
138 if len(list1) != len(list2):
139 print >> sys.stderr, "WARNING: Input files contain a different number of reference sequences."
140 errorMessage = True
141 elif list1 != list2:
142 print >> sys.stderr, "WARNING: Input files contain different or differently ordered reference sequences."
143 errorMessage = True
144 if errorMessage:
145 print >> sys.stderr, "Results may not be as expected."
146 print >> sys.stderr, "Ensure that input files have the same reference sequences in the same order."
147 print >> sys.stderr, "Reference sequence lists observed were:\n\t", list1, "\n\t", list2
148
149 # Write out a vcf record to file. The record written depends on the
150 # value of 'priority' and could therefore be the record from either
151 # of the vcf files, or a combination of them.
152
153 def writeVcfRecord(priority, v1, v2, outputFile):
154 if priority == 0:
155 if v1.quality >= v2.quality: outputFile.write(v1.record)
156 else: outputFile.write(v2.record)
157 elif priority == 1: outputFile.write(v1.record)
158 elif priority == 2: outputFile.write(v2.record)
159 elif priority == 3:
160
161 # Define the missing entry values (depends on the number of data sets
162 # in the file).
163 info = ""
164 missingEntry1 = missingEntry2 = "."
165 for i in range(1, v1.numberDataSets): missingEntry1 += "/."
166 for i in range(1, v2.numberDataSets): missingEntry2 += "/."
167 secondList = v2.infoTags.copy()
168
169 # Build up the info field.
170 for tag in v1.infoTags:
171 if secondList.has_key(tag):
172 if v1.infoHeaderTags[tag][1].lower() != "flag": info += tag + "=" + v1.infoTags[tag] + "/" + v2.infoTags[tag] + ";"
173 del secondList[tag]
174 else:
175 if v1.infoHeaderTags[tag][1].lower() != "flag": info += tag + "=" + v1.infoTags[tag] + "/" + missingEntry2 + ";"
176
177 # Now include the info tags that are not populated in the first vcf file.
178 for tag in secondList:
179 if v2.infoHeaderTags[tag][1].lower() != "flag": info += tag + "=" + missingEntry1 + "/" + v2.infoTags[tag] + ";"
180
181 # Build the complete record.
182 info = info.rstrip(";")
183 record = v1.referenceSequence + "\t" + str(v1.position) + "\t" + v1.rsid + "\t" + v1.ref + "\t" + \
184 v1.alt + "/" + v2.alt + "\t" + v1.quality + "/" + v2.quality + "\t.\t" + info
185 print >> outputFile, record
186 else:
187 print >> sys.sterr, "Unknown file priority."
188 exit(1)