| 1 | #! /usr/bin/env python |
|---|
| 2 | import sys |
|---|
| 3 | import os |
|---|
| 4 | import json |
|---|
| 5 | import operator |
|---|
| 6 | import time |
|---|
| 7 | from ete2 import Tree, SeqGroup |
|---|
| 8 | from config import EpacConfig |
|---|
| 9 | from subprocess import call |
|---|
| 10 | |
|---|
| 11 | class hmmer: |
|---|
| 12 | def __init__(self, config, refalign = None, query = None, refprofile = None, discard = None, seqs = None, minp = 0.9): |
|---|
| 13 | self.cfg = config |
|---|
| 14 | self.refalign = refalign |
|---|
| 15 | self.query = query |
|---|
| 16 | self.refprofile = refprofile |
|---|
| 17 | self.hmmbuildpath = config.hmmer_home + "/hmmbuild" |
|---|
| 18 | self.hmmalignpath = config.hmmer_home + "/hmmalign" |
|---|
| 19 | if self.refprofile == None: |
|---|
| 20 | self.refprofile = config.tmp_fname("%NAME%.hmm") |
|---|
| 21 | self.stockname = config.tmp_fname("%NAME%.stock") |
|---|
| 22 | self.trimed = config.tmp_fname("%NAME%.trimed.afa") |
|---|
| 23 | self.output = config.tmp_fname("%NAME%.aligned.afa") |
|---|
| 24 | self.merged = config.tmp_fname("%NAME%.merged.afa") |
|---|
| 25 | self.discardpath = discard |
|---|
| 26 | self.seqs = seqs |
|---|
| 27 | self.minp = minp |
|---|
| 28 | self.minl = 100 |
|---|
| 29 | |
|---|
| 30 | def _remove(self, filename): |
|---|
| 31 | if os.path.exists(filename): |
|---|
| 32 | os.remove(filename) |
|---|
| 33 | |
|---|
| 34 | def __del__(self): |
|---|
| 35 | self._remove(self.stockname) |
|---|
| 36 | self._remove(self.trimed) |
|---|
| 37 | self._remove(self.output) |
|---|
| 38 | |
|---|
| 39 | def __processHMMseq(self, seqin): |
|---|
| 40 | newseq = "" |
|---|
| 41 | for s in seqin: |
|---|
| 42 | if s == ".": |
|---|
| 43 | pass |
|---|
| 44 | elif s == "-": |
|---|
| 45 | newseq = newseq + s |
|---|
| 46 | elif s.isupper(): |
|---|
| 47 | newseq = newseq + s |
|---|
| 48 | return newseq |
|---|
| 49 | |
|---|
| 50 | def build_hmm_profile(self, informat="afa"): |
|---|
| 51 | #hmmbuild --informat afa refotu.hmm ref_outs_547.fas |
|---|
| 52 | call_str = [self.hmmbuildpath, "--symfrac", "0.0", "--informat", informat, self.refprofile, self.refalign] |
|---|
| 53 | if self.cfg.verbose: |
|---|
| 54 | print "\n" + ' '.join(call_str) + "\n" |
|---|
| 55 | call(call_str) #, stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT) |
|---|
| 56 | return self.refprofile |
|---|
| 57 | |
|---|
| 58 | def hmm_align(self): |
|---|
| 59 | #hmmalign -o 454.stock refotu.hmm 454input.fna.min100.fasta |
|---|
| 60 | call_str = [self.hmmalignpath,"-o", self.stockname, self.refprofile, self.query] |
|---|
| 61 | if self.cfg.verbose: |
|---|
| 62 | print "\n" + ' '.join(call_str) + "\n" |
|---|
| 63 | call(call_str) #, stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT) |
|---|
| 64 | return self.stockname |
|---|
| 65 | |
|---|
| 66 | def get_hmm_refalignment(self): |
|---|
| 67 | sites = [] |
|---|
| 68 | hmp = open(self.refprofile) |
|---|
| 69 | l = hmp.readline() |
|---|
| 70 | start = False |
|---|
| 71 | while l!="": |
|---|
| 72 | if l.startswith("//"): |
|---|
| 73 | break |
|---|
| 74 | if start: |
|---|
| 75 | l = l.strip() |
|---|
| 76 | ll = l.split() |
|---|
| 77 | usedsite = int(ll[5]) |
|---|
| 78 | sites.append(usedsite) |
|---|
| 79 | l = hmp.readline() |
|---|
| 80 | l = hmp.readline() |
|---|
| 81 | else: |
|---|
| 82 | if l.startswith("HMM "): |
|---|
| 83 | start = True |
|---|
| 84 | l = hmp.readline() |
|---|
| 85 | l = hmp.readline() |
|---|
| 86 | l = hmp.readline() |
|---|
| 87 | l = hmp.readline() |
|---|
| 88 | l = hmp.readline() |
|---|
| 89 | hmp.close() |
|---|
| 90 | align = SeqGroup(self.refalign) |
|---|
| 91 | fout = open(self.trimed, "w") |
|---|
| 92 | for entr in align.get_entries(): |
|---|
| 93 | fout.write(">" + entr[0] + "\n") |
|---|
| 94 | for pos in sites: |
|---|
| 95 | fout.write(entr[1][pos-1]) |
|---|
| 96 | fout.write("\n") |
|---|
| 97 | fout.close() |
|---|
| 98 | return self.trimed, len(sites) |
|---|
| 99 | |
|---|
| 100 | def parse_HMM(self, l_ref): |
|---|
| 101 | """stock format""" |
|---|
| 102 | cnt = 0 |
|---|
| 103 | fin = open(self.stockname) |
|---|
| 104 | line = fin.readline() |
|---|
| 105 | seqs = {} |
|---|
| 106 | while line!="": |
|---|
| 107 | if line.startswith("//"): |
|---|
| 108 | break |
|---|
| 109 | elif line.startswith("#=GS"): |
|---|
| 110 | cnt = 0 |
|---|
| 111 | pass |
|---|
| 112 | elif line.startswith("#"): |
|---|
| 113 | pass |
|---|
| 114 | elif line.startswith("\n"): |
|---|
| 115 | cnt = cnt + 1 |
|---|
| 116 | else: |
|---|
| 117 | line = line.strip() |
|---|
| 118 | if cnt == 1: |
|---|
| 119 | l2 = line.split() |
|---|
| 120 | ss = self.__processHMMseq(l2[1]) |
|---|
| 121 | seqs[l2[0]] = ss |
|---|
| 122 | else: |
|---|
| 123 | l2 = line.split() |
|---|
| 124 | seq = seqs[l2[0]] |
|---|
| 125 | ss = self.__processHMMseq(l2[1]) |
|---|
| 126 | seqs[l2[0]] = seq + ss |
|---|
| 127 | line = fin.readline() |
|---|
| 128 | fin.close() |
|---|
| 129 | fout = open(self.output, "w") |
|---|
| 130 | foutdiscard = open(self.discardpath, "w") |
|---|
| 131 | for key in seqs.keys(): |
|---|
| 132 | |
|---|
| 133 | #key is the sequence name which is the id |
|---|
| 134 | numleft = count_non_gap(seqs[key]) |
|---|
| 135 | numall = len(self.seqs.get_seq(key)) |
|---|
| 136 | |
|---|
| 137 | if numall > 0: |
|---|
| 138 | pleft = float(numleft) / float(numall) |
|---|
| 139 | else: |
|---|
| 140 | pleft = 0 |
|---|
| 141 | |
|---|
| 142 | if numleft >= self.minl and pleft >= self.minp: |
|---|
| 143 | fout.write(">" + key + "\n") |
|---|
| 144 | seq = seqs[key] |
|---|
| 145 | lappd = l_ref - len(seq) |
|---|
| 146 | if lappd > 0: |
|---|
| 147 | appd = "-" * lappd |
|---|
| 148 | seq = seq + appd |
|---|
| 149 | elif lappd < 0: |
|---|
| 150 | print("Warning: query sequence > ref sequence") |
|---|
| 151 | |
|---|
| 152 | fout.write(seq + "\n") |
|---|
| 153 | else: |
|---|
| 154 | foutdiscard.write(">" + key + "\n") |
|---|
| 155 | fout.close() |
|---|
| 156 | return self.output |
|---|
| 157 | |
|---|
| 158 | def hmm_alignment(self, ref_align, query, outfolder, lmin = 100): |
|---|
| 159 | if not os.path.exists(self.refprofile): |
|---|
| 160 | self.build_hmm_profile() |
|---|
| 161 | self.hmm_align() |
|---|
| 162 | final_ref, ref_len = self.trim_refalign_hmm() |
|---|
| 163 | final_query = self.parse_HMM(l_ref = ref_len, minl = lmin) |
|---|
| 164 | |
|---|
| 165 | def align(self): |
|---|
| 166 | #aquire reference alignment that hmm would use |
|---|
| 167 | refaln, numsite = self.get_hmm_refalignment() |
|---|
| 168 | #alignment |
|---|
| 169 | self.hmm_align() |
|---|
| 170 | #process alignment |
|---|
| 171 | queryaln = self.parse_HMM(l_ref = numsite) |
|---|
| 172 | #merge refrence and query alignment |
|---|
| 173 | merge_alignment(aln1 = refaln, aln2 = queryaln, fout = self.merged, numsites = numsite) |
|---|
| 174 | return self.merged |
|---|
| 175 | |
|---|
| 176 | |
|---|
| 177 | class muscle: |
|---|
| 178 | def __init__(self, config): |
|---|
| 179 | self.cfg = config |
|---|
| 180 | self.musclepath = config.muscle_home + "/muscle" |
|---|
| 181 | self.outname = config.tmp_fname("%NAME%.afa") |
|---|
| 182 | |
|---|
| 183 | def merge(self, aln1, aln2): |
|---|
| 184 | #muscle -profile -in1 existing_msa.afa -in2 new_seq.fa -out combined.afa |
|---|
| 185 | call_str = [self.musclepath,"-profile", "-in1", aln1, "-in2", aln2, "-out", self.outname] |
|---|
| 186 | if self.cfg.debug: |
|---|
| 187 | print "\n" + ' '.join(call_str) + "\n" |
|---|
| 188 | call(call_str) |
|---|
| 189 | return self.outname |
|---|
| 190 | |
|---|
| 191 | |
|---|
| 192 | def merge_alignment(aln1, aln2, fout, numsites): |
|---|
| 193 | seqs1 = SeqGroup(aln1) |
|---|
| 194 | seqs2 = SeqGroup(aln2) |
|---|
| 195 | if len(seqs1) == 0 or len(seqs2) == 0: |
|---|
| 196 | print("No sequences aligned! ") |
|---|
| 197 | sys.exit() |
|---|
| 198 | with open(fout, "w") as fo: |
|---|
| 199 | for seq in seqs1.iter_entries(): |
|---|
| 200 | if len(seq[1].strip()) == numsites: |
|---|
| 201 | fo.write(">" + seq[0] + "\n" + seq[1] + "\n") |
|---|
| 202 | else: |
|---|
| 203 | print("Error in alignment ....") |
|---|
| 204 | sys.exit() |
|---|
| 205 | for seq in seqs2.iter_entries(): |
|---|
| 206 | if len(seq[1].strip()) == numsites: |
|---|
| 207 | fo.write(">" + seq[0] + "\n" + seq[1] + "\n") |
|---|
| 208 | else: |
|---|
| 209 | print("Error in alignment ....") |
|---|
| 210 | sys.exit() |
|---|
| 211 | |
|---|
| 212 | |
|---|
| 213 | def count_non_gap(seqin): |
|---|
| 214 | cnt = 0 |
|---|
| 215 | for s in seqin: |
|---|
| 216 | if s!="-": |
|---|
| 217 | cnt = cnt + 1 |
|---|
| 218 | return cnt |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | if __name__ == "__main__": |
|---|
| 222 | print("This is main") |
|---|
| 223 | cfg = EpacConfig() |
|---|
| 224 | hm = hmmer(config = cfg, refalign = "example/t1_trimed.fa") |
|---|
| 225 | #trimed = hm.process_ref_alignment() |
|---|
| 226 | pf = hm.build_hmm_profile() |
|---|
| 227 | print(pf) |
|---|
