Mercurial > repos > drosofff > msp_blastparser_and_hits

comparison BlastParser_and_hits.py @ 0:69ea2a13947f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://bitbucket.org/drosofff/gedtools/
author drosofff
date Sun, 21 Jun 2015 14:31:29 -0400
parents
children 1964514aabde
comparison
equal deleted inserted replaced
-1:000000000000 0:69ea2a13947f
1 #!/usr/bin/python
2 # blastn blastx parser revised debugged: 3-4-2015. Commit issue.
3 # drosofff@gmail.com
4
5 import sys
6 import argparse
7 from collections import defaultdict
8
9 def Parser():
10 the_parser = argparse.ArgumentParser()
11 the_parser.add_argument('--blast', action="store", type=str, help="Path to the blast output (tabular format, 12 column)")
12 the_parser.add_argument('--sequences', action="store", type=str, help="Path to the fasta file with blasted sequences")
13 the_parser.add_argument('--fastaOutput', action="store", type=str, help="fasta output file of blast hits")
14 the_parser.add_argument('--tabularOutput', action="store", type=str, help="tabular output file of blast analysis")
15 the_parser.add_argument('--flanking', action="store", type=int, help="number of flanking nucleotides added to the hit sequences")
16 the_parser.add_argument('--mode', action="store", choices=["verbose", "short"], type=str, help="reporting (verbose) or not reporting (short) oases contigs")
17 args = the_parser.parse_args()
18 if not all ( (args.sequences, args.blast, args.fastaOutput, args.tabularOutput) ):
19 the_parser.error('argument(s) missing, call the -h option of the script')
20 if not args.flanking:
21 args.flanking = 0
22 return args
23
24 def median(lst):
25 lst = sorted(lst)
26 if len(lst) < 1:
27 return None
28 if len(lst) %2 == 1:
29 return lst[((len(lst)+1)/2)-1]
30 if len(lst) %2 == 0:
31 return float(sum(lst[(len(lst)/2)-1:(len(lst)/2)+1]))/2.0
32
33 def mean(lst):
34 if len(lst) < 1:
35 return 0
36 return sum(lst) / float(len(lst))
37
38 def getfasta (fastafile):
39 fastadic = {}
40 for line in open (fastafile):
41 if line[0] == ">":
42 header = line[1:-1]
43 fastadic[header] = ""
44 else:
45 fastadic[header] += line
46 for header in fastadic:
47 fastadic[header] = "".join(fastadic[header].split("\n"))
48 return fastadic
49
50 def insert_newlines(string, every=60):
51 lines = []
52 for i in xrange(0, len(string), every):
53 lines.append(string[i:i+every])
54 return '\n'.join(lines)
55
56 def getblast (blastfile):
57 '''blastinfo [0] Percentage of identical matches
58 blastinfo [1] Alignment length
59 blastinfo [2] Number of mismatches
60 blastinfo [3] Number of gap openings
61 blastinfo [4] Start of alignment in query
62 blastinfo [5] End of alignment in query
63 blastinfo [6] Start of alignment in subject (database hit)
64 blastinfo [7] End of alignment in subject (database hit)
65 blastinfo [8] Expectation value (E-value)
66 blastinfo [9] Bit score
67 blastinfo [10] Subject length (NEED TO BE SPECIFIED WHEN RUNNING BLAST) '''
68 blastdic = defaultdict (dict)
69 for line in open (blastfile):
70 fields = line[:-1].split("\t")
71 transcript = fields[0]
72 subject = fields[1]
73 blastinfo = [float(fields[2]) ] # blastinfo[0]
74 blastinfo = blastinfo + [int(i) for i in fields[3:10] ] # blastinfo[1:8] insets 1 to 7
75 blastinfo.append(fields[10]) # blastinfo[8] E-value remains as a string type
76 blastinfo.append(float(fields[11])) # blastinfo[9] Bit score
77 blastinfo.append(int(fields[12])) # blastinfo[10] Subject length MUST BE RETRIEVED THROUGH A 13 COLUMN BLAST OUTPUT
78 try:
79 blastdic[subject][transcript].append(blastinfo)
80 except:
81 blastdic[subject][transcript] = [ blastinfo ]
82 return blastdic
83
84 def getseq (fastadict, transcript, up, down, orientation="direct"):
85 def reverse (seq):
86 revdict = {"A":"T","T":"A","G":"C","C":"G","N":"N"}
87 revseq = [revdict[i] for i in seq[::-1]]
88 return "".join(revseq)
89 pickseq = fastadict[transcript][up-1:down]
90 if orientation == "direct":
91 return pickseq
92 else:
93 return reverse(pickseq)
94
95 def subjectCoverage (fastadict, blastdict, subject, QueriesFlankingNucleotides=0):
96 SubjectCoverageList = []
97 HitDic = {}
98 bitScores = []
99 for transcript in blastdict[subject]:
100 prefix = "%s--%s_" % (subject, transcript)
101 hitNumber = 0
102 for hit in blastdict[subject][transcript]:
103 hitNumber += 1
104 suffix = "hit%s_IdMatch=%s,AligLength=%s,E-val=%s" % (hitNumber, hit[0], hit[1], hit[8])
105 HitDic[prefix+suffix] = GetHitSequence (fastadict, transcript, hit[4], hit[5], QueriesFlankingNucleotides) #query coverage by a hit is in hit[4:6]
106 SubjectCoverageList += range (min([hit[6], hit[7]]), max([hit[6], hit[7]]) + 1) # subject coverage by a hit is in hit[6:8]
107 bitScores.append(hit[9])
108 subjectLength = hit [10] # always the same value for a given subject. Stupid but simple
109 TotalSubjectCoverage = len ( set (SubjectCoverageList) )
110 RelativeSubjectCoverage = TotalSubjectCoverage/float(subjectLength)
111 return HitDic, subjectLength, TotalSubjectCoverage, RelativeSubjectCoverage, max(bitScores), mean(bitScores)
112
113 def GetHitSequence (fastadict, FastaHeader, leftCoordinate, rightCoordinate, FlankingValue):
114 if rightCoordinate > leftCoordinate:
115 polarity = "direct"
116 else:
117 polarity = "reverse"
118 leftCoordinate, rightCoordinate = rightCoordinate, leftCoordinate
119 if leftCoordinate - FlankingValue > 0:
120 leftCoordinate -= FlankingValue
121 else:
122 leftCoordinate = 1
123 return getseq (fastadict, FastaHeader, leftCoordinate, rightCoordinate, polarity)
124
125 def outputParsing (F, Fasta, results, Xblastdict, fastadict, mode="verbose"):
126 F= open(F, "w")
127 Fasta=open(Fasta, "w")
128 if mode == "verbose":
129 print >>F, "# SeqId\t%Identity\tAlignLength\tStartSubject\tEndSubject\t%QueryHitCov\tE-value\tBitScore\n"
130 for subject in sorted (results, key=lambda x: results[x]["meanBitScores"], reverse=True):
131 print >> F, "#\n# %s" % subject
132 print >> F, "# Suject Length: %s" % (results[subject]["subjectLength"])
133 print >> F, "# Total Subject Coverage: %s" % (results[subject]["TotalCoverage"])
134 print >> F, "# Relative Subject Coverage: %s" % (results[subject]["RelativeSubjectCoverage"])
135 print >> F, "# Maximum Bit Score: %s" % (results[subject]["maxBitScores"])
136 print >> F, "# Mean Bit Score: %s" % (results[subject]["meanBitScores"])
137 for header in results[subject]["HitDic"]:
138 print >> Fasta, ">%s\n%s" % (header, insert_newlines(results[subject]["HitDic"][header]) )
139 print >> Fasta, "" # final carriage return for the sequence
140 for transcript in Xblastdict[subject]:
141 transcriptSize = float(len(fastadict[transcript]))
142 for hit in Xblastdict[subject][transcript]:
143 percentIdentity, alignLenght, subjectStart, subjectEnd, queryCov = hit[0], hit[1], hit[6], hit[7], "%.1f" % (abs(hit[5]-hit[4])/transcriptSize*100)
144 Eval, BitScore = hit[8], hit[9]
145 info = [transcript] + [percentIdentity, alignLenght, subjectStart, subjectEnd, queryCov, Eval, BitScore]
146 info = [str(i) for i in info]
147 info = "\t".join(info)
148 print >> F, info
149 else:
150 print >>F, "# subject\tsubject length\tTotal Subject Coverage\tRelative Subject Coverage\tMaximum Bit Score\tMean Bit Score"
151 for subject in sorted (results, key=lambda x: results[x]["meanBitScores"], reverse=True):
152 line = []
153 line.append(subject)
154 line.append(results[subject]["subjectLength"])
155 line.append(results[subject]["TotalCoverage"])
156 line.append(results[subject]["RelativeSubjectCoverage"])
157 line.append(results[subject]["maxBitScores"])
158 line.append(results[subject]["meanBitScores"])
159 line = [str(i) for i in line]
160 print >> F, "\t".join(line)
161 for header in results[subject]["HitDic"]:
162 print >> Fasta, ">%s\n%s" % (header, insert_newlines(results[subject]["HitDic"][header]) )
163 print >> Fasta, "" # final carriage return for the sequence
164 F.close()
165 Fasta.close()
166
167
168
169 def __main__ ():
170 args = Parser()
171 fastadict = getfasta (args.sequences)
172 Xblastdict = getblast (args.blast)
173 results = defaultdict(dict)
174 for subject in Xblastdict:
175 results[subject]["HitDic"], results[subject]["subjectLength"], results[subject]["TotalCoverage"], results[subject]["RelativeSubjectCoverage"], results[subject]["maxBitScores"], results[subject]["meanBitScores"] = subjectCoverage(fastadict, Xblastdict, subject, args.flanking)
176 outputParsing (args.tabularOutput, args.fastaOutput, results, Xblastdict, fastadict, args.mode)
177 if __name__=="__main__": __main__()