annotate Tryp_T.py @ 10:320bdfa4d927 draft default tip

Deleted selected files
author johnheap
date Wed, 04 Jul 2018 11:37:20 -0400
parents 1e2f57c43854
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
1 """
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
2 * Copyright 2018 University of Liverpool
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
3 * Author: John Heap, Computational Biology Facility, UoL
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
4 * Based on original scripts of Sara Silva Pereira, Institute of Infection and Global Health, UoL
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
5 *
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
6 * Licensed under the Apache License, Version 2.0 (the "License");
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
7 * you may not use this file except in compliance with the License.
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
8 * You may obtain a copy of the License at
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
9 *
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
10 * http://www.apache.org/licenses/LICENSE-2.0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
11 *
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
12 * Unless required by applicable law or agreed to in writing, software
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
13 * distributed under the License is distributed on an "AS IS" BASIS,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
15 * See the License for the specific language governing permissions and
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
16 * limitations under the License.
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
17 *
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
18 """
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
19
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
20
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
21 import subprocess
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
22 import pandas as pd
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
23 import re
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
24 import os
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
25 import sys
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
26 import matplotlib as mpl
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
27 mpl.use('Agg')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
28 import matplotlib.pyplot as plt
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
29
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
30 pList = ['P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11', 'P12', 'P13', 'P14', 'P15']
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
31 quietString = "" #"">> Vap_log.txt 2>&1"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
32 def transcriptMapping(inputname, strain, forwardFN,reverseFN):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
33 #where is our Reference data -
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
34 dir_path = os.path.dirname(os.path.realpath(__file__))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
35 refName = dir_path+"/data/Reference/Tc148" #default
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
36 if strain == "Tc148":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
37 refName = dir_path+"/data/Reference/Tc148"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
38 if strain == "IL3000":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
39 refName = dir_path+"/data/Reference/IL3000"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
40 #argString = "bowtie2 -x Refe4rence/IL3000 -1 data/"+forwardFN+" -2 data/"+reverseFN+" -S "+inputname+".sam" #>log.txt
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
41 #argString = "bowtie2 -x Reference/Tc148 -1 data/"+forwardFN+" -2 data/"+reverseFN+" -S "+inputname+".sam" #>log.txt
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
42 argString = "bowtie2 -x "+refName+" -1 "+forwardFN+" -2 "+reverseFN+" -S "+inputname+".sam"+quietString #>log.txt
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
43 #print(argString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
44 returncode = subprocess.call(argString, shell=True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
45
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
46 def processSamFiles(inputname):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
47 #debug use a mapping sam file we have already found
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
48 #dir_path = os.path.dirname(os.path.realpath(__file__))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
49 #bugName = dir_path+"/data/T_Test" #defasult
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
50
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
51 cur_path = os.getcwd()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
52 samName = cur_path+"/"+inputname
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
53
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
54 #argString = "samtools view -bS "+bugName+" > "+inputname+".bam"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
55 argString = "samtools view -bS "+inputname+".sam > "+samName+".bam"+quietString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
56 #print(argString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
57 returncode = subprocess.call(argString, shell=True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
58
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
59
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
60 #argString = "samtools sort "+bugName+" -o "+inputname+".sorted"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
61 argString = "samtools sort "+samName+".bam -o "+samName+".sorted"+quietString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
62 #print("argstring = "+argString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
63 returncode = subprocess.call(argString, shell=True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
64
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
65 #argString = "samtools index "+bugName+".sorted "+inputname+".sorted.bai"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
66 argString = "samtools index "+samName+".sorted "+samName+".sorted.bai"+quietString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
67 #print("argstring = " + argString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
68 returncode = subprocess.call(argString, shell=True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
69
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
70
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
71
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
72
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
73 def transcriptAbundance(inputname, strain):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
74 dir_path = os.path.dirname(os.path.realpath(__file__))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
75 refName = dir_path + "/data/Reference/ORFAnnotation.gtf" # defasult
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
76 if strain == "Tc148":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
77 refName = dir_path + "/data/Reference/ORFAnnotation.gtf"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
78 if strain == "IL3000":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
79 refName = dir_path + "/data/Reference/IL3000.gtf"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
80 #argString = "cufflinks -G Reference/IL3000.gtf -o "+inputname+".cuff -u -p 8 "+inputname+".sorted"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
81 #argString = "cufflinks -G Reference/ORFAnnotation.gtf -o "+inputname+".cuff -u -p 8 "+inputname+".sorted"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
82 argString = "cufflinks -q -G "+refName+" -o "+inputname+".cuff -u -p 8 "+inputname+".sorted"+quietString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
83 returncode = subprocess.call(argString, shell = True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
84
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
85
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
86 def convertToFasta(inputName, strain): #equivalent to Sara's awk scripte
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
87 dir_path = os.path.dirname(os.path.realpath(__file__))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
88 refName = dir_path + "/data/Reference/ORFAnnotation.gtf" # default
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
89 if strain == "Tc148":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
90 refName = dir_path + "/data/Reference/148_prot.fasta"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
91 if strain == "IL3000":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
92 refName = dir_path + "data/Reference/IL3000_prot.fasta"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
93
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
94 cuff_df = pd.read_csv(inputName+".cuff/genes.fpkm_tracking", sep='\t')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
95 cuff_df = cuff_df[(cuff_df['FPKM'] > 0)]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
96 cuff_df.to_csv("cuffTest.csv")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
97 gene_id_List = cuff_df['gene_id'].tolist()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
98
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
99 #print(gene_id_List)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
100 #print ("Found from 8880="+str(found))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
101
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
102 # need to load in IL3000_prot.fasta
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
103 # for each line with >TcIL3000_1_1940
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
104 # search within cuff_df[gene_id] for match
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
105 # add it to the outfile. (need to save it as used by hmmer later
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
106 number = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
107 all = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
108 with open(inputName+"_6frame.fas", 'w') as outfile:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
109 ref = open(refName,'r')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
110 #ref = open(r"Reference/IL3000_prot.fasta",'r')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
111 n = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
112 line = ref.readline()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
113 while line:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
114 if line[0] == '>':
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
115 all = all+1
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
116 ln = line[1:] #remove >
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
117 ln = ln.rstrip() #remove /n /r etc
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
118 #print (ln)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
119 if ln in gene_id_List:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
120 number = number+1
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
121 outfile.write(line)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
122 line = ref.readline()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
123 if line:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
124 while line[0] != '>':
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
125 outfile.write(line)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
126 line=ref.readline()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
127 else:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
128 line = ref.readline()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
129 else:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
130 line =ref.readline()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
131 ref.close()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
132 print(str(len(gene_id_List))+":"+str(number)+" from "+str(all))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
133 return cuff_df
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
134
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
135 def HMMerMotifSearch(name, strain, cuff_df):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
136 motifs = ['1', '2a', '2b', '3', '4a', '4b', '4c', '5', '6', '7', '8a', '8b', '9a', '9b',
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
137 '9c', '10a', '10b', '11a', '11b', '12', '13a', '13b', '13c', '13d', '14', '15a', '15b', '15c']
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
138 dir_path = os.path.dirname(os.path.realpath(__file__))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
139 phylopath = dir_path + "/data/Motifs/Phylotype"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
140 lineCounts = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
141 compoundList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
142 for m in motifs:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
143 argString = "hmmsearch "+phylopath + m + ".hmm " + name + "_6frame.fas > Phy" + m + ".out"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
144 print(argString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
145 subprocess.call(argString, shell=True)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
146 hmmResult = open("Phy" + m + ".out", 'r')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
147 regex = r"Tc148[0-9]{1,8}"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
148 if strain == "Tc148":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
149 regex = r"Tc148[0-9]{1,8}"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
150 if strain == "IL3000":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
151 regex = r"TcIL3000_[0-9]{1,4}_[0-9]{1,5}"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
152 n = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
153 outList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
154 for line in hmmResult:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
155 m = re.search(regex, line)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
156 if m:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
157 outList.append(""+m.group())
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
158 n += 1
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
159 if re.search(r"inclusion", line):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
160 print("inclusion threshold reached")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
161 break
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
162 compoundList.append(outList)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
163 lineCounts.append(n)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
164 hmmResult.close()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
165 #print(lineCounts)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
166
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
167 #print(cuff_df)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
168 concatGroups = [1, 2, 1, 3, 1, 1, 1, 2, 3, 2, 2, 1, 4, 1, 3]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
169 countList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
170 weightList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
171 countIndex = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
172 totalCount = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
173 totalWeigth = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
174 for c in concatGroups:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
175 a = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
176 weight = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
177 for n in range(0, c):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
178 a = a + compoundList.pop(0)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
179 t = set(a)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
180 countList.append(len(t))
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
181 wa = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
182 for w in t:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
183 wt = cuff_df.loc[cuff_df['gene_id'] == w, 'FPKM'].iloc[0]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
184 #print(w)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
185 #print(wt)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
186 wa = wa+wt
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
187 weightList.append(wa)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
188 totalWeigth+=wa
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
189 totalCount += len(t)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
190 countList.append(totalCount)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
191 weightList.append(totalWeigth)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
192 #print(countList)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
193 #print("--------")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
194 #print(weightList)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
195 #print("--------")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
196 return countList,weightList
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
197
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
198 def relativeFrequencyTable(countList, name, htmlresource):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
199 relFreqList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
200 c = float(countList[15])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
201 for i in range(0, 15):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
202 relFreqList.append(countList[i] / c)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
203
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
204 data = {'Phylotype': pList, 'Relative Frequency': relFreqList}
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
205 relFreq_df = pd.DataFrame(data)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
206 j_fname = htmlresource+ "/" + name + "_t_relative_frequency.csv"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
207 relFreq_df.to_csv(j_fname)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
208 return relFreqList # 0-14 = p1-p15 counts [15] = total counts
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
209
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
210
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
211 def weightedFrequencyTable(countList, name, htmlresource):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
212 relFreqList = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
213 c = float(countList[15])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
214 for i in range(0, 15):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
215 relFreqList.append(countList[i] / c)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
216
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
217 data = {'Phylotype': pList, 'Weighted Frequency': relFreqList}
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
218 relFreq_df = pd.DataFrame(data)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
219 j_fname = htmlresource+ "/" + name + "_t_weighted_frequency.csv"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
220 relFreq_df.to_csv(j_fname)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
221 return relFreqList # 0-14 = p1-p15 counts [15] = total counts
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
222
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
223
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
224
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
225 def createStackedBar(name,freqList,strain,pdf,html_resource):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
226 palette = ["#0000ff", "#6495ed", "#00ffff", "#caff70",
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
227 "#228b22", "#528b8b", "#00ff00", "#a52a2a",
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
228 "#ff0000", "#ffff00", "#ffa500", "#ff1493",
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
229 "#9400d3", "#bebebe", "#000000", "#ff00ff"]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
230
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
231 VAP_148 = [0.072, 0.032, 0.032, 0.004, 0.007,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
232 0.005, 0.202, 0.004, 0.006, 0.014,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
233 0.130, 0.133, 0.054, 0.039, 0.265]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
234
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
235 VAP_IL3000 = [0.073, 0.040, 0.049, 0.018, 0.060,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
236 0.055, 0.054, 0.025, 0.012, 0.060,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
237 0.142, 0.100, 0.061, 0.078, 0.172]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
238 cmap = plt.cm.get_cmap('tab20')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
239 palette = [cmap(i) for i in range(cmap.N)]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
240
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
241 if strain == "Tc148":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
242 VAPtable = VAP_148
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
243 VAPname='Tc148\nGenome VAP'
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
244 if strain == "IL3000":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
245 VAPtable = VAP_IL3000
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
246 VAPname= 'IL3000\nGenome VAP'
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
247 width = 0.35 # the width of the bars: can also be len(x) sequence
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
248 plots = []
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
249 fpos = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
250 vpos = 0
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
251 for p in range(0, 15):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
252 tp = plt.bar(0, freqList[p], width, color= palette[p], bottom = fpos)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
253 fpos +=freqList[p]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
254
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
255 tp = plt.bar(1, VAPtable[p], width, color= palette[p], bottom = vpos)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
256 vpos +=VAPtable[p]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
257
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
258 plots.append(tp)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
259 plt.xticks([0,1],[name,VAPname])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
260 plt.legend(plots[::-1],['p15','p14','p13','p12','p11','p10','p9','p8','p7','p6','p5','p4','p3','p2','p1'])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
261 title = "Figure Legend: The transcriptomic Variant Antigen Profile of $\itTrypanosoma$ $\itcongolense$ estimated as phylotype " \
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
262 "proportion adjusted for transcript abundance and the reference genomic Variant Antigen Profile. " \
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
263 "\nData was produced with the 'Variant Antigen Profiler' (Silva Pereira and Jackson, 2018)."
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
264 #plt.title(title, wrap="True")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
265 #plt.text(-0.2, -0.05, title, va="top", transform=ax.transAxes, wrap="True")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
266 plt.text(-0.3, -0.15, title, va="top", wrap="True")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
267 plt.tight_layout(pad=1.5)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
268 plt.subplots_adjust(bottom = 0.3,top=0.99,left=0.125,right=0.9,hspace=0.2,wspace=0.2)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
269
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
270 plt.savefig(html_resource + "/stackedbar.png")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
271 if pdf == 'PDF_Yes':
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
272 plt.savefig(html_resource + "/stackedbar.pdf")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
273 #plt.show()
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
274
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
275
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
276 def createHTML(name,htmlfn,htmlresource,freqList,weightList):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
277 #assumes imgs are heatmap.png, dheatmap.png, vapPCA.png and already in htmlresource
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
278 htmlString = r"<html><title>T.congolense VAP</title><body><div style='text-align:center'><h2><i>Trypanosoma congolense</i> Variant Antigen Profile</h2><h3>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
279 htmlString += name
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
280 htmlString += r"<br>Transcriptomic Analysis</h3></p>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
281 htmlString += "<p style = 'margin-left:20%; margin-right:20%'>Table Legend: Variant Antigen Profiles of a transcriptome of <i>Trypanosoma congolense</i> estimated as phylotype proportion. " \
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
282 "Weighted frequency refers to the phylotype proportion based transcript abundance. " \
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
283 "Data was produced with the 'Variant Antigen Profiler' (Silva Pereira and Jackson, 2018).</p> "
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
284 htmlString += r"<style> table, th, tr, td {border: 1px solid black; border-collapse: collapse;}</style>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
285
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
286 htmlString += r"<table style='width:50%;margin-left:25%;text-align:center'><tr><th>Phylotype</th><th>Relative Frequency</th><th>Weighted Frequency</th></tr>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
287 tabString = ""
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
288 # flush out table with correct values
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
289 for i in range(0, 15):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
290 f = format(freqList[i], '.4f')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
291 w = format(weightList[i], '.4f')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
292 tabString += "<tr><td>phy" + str(i + 1) + "</td><td>" + f + "</td><td>" + w + "</td></tr>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
293 htmlString += tabString + "</table><br><br><br><br><br>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
294 htmlString += r"<p> <h3>Stacked Bar chart of Phylotype Frequency</h3> The 'weighted' relative frequency of each phylotype alongside the VAP of selected strain.</p>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
295 imgString = r"<img src = 'stackedbar.png' alt='Stacked bar chart of phylotype variation' style='max-width:100%'><br><br>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
296 htmlString += imgString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
297
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
298 # htmlString += r"<p><h3>The Deviation Heat Map and Dendogram</h3>The phylotype variation expressed as the deviation from your sample mean compared to the model dataset</p>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
299 # imgString = r"<img src = 'dheatmap.png' alt='Deviation Heatmap' style='max-width:100%'><br><br>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
300 # htmlString += imgString
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
301
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
302 # htmlString += r"<p><h3>The Variation PCA plot</h3>PCA analysis corresponding to absolute variation. Colour coded according to location</p>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
303 # imgString = r"<img src = 'vapPCA.png' alt='PCA Analysis' style='max-width:100%'><br><br>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
304 # htmlString += imgString + r"</div></body></html>"
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
305
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
306 with open(htmlfn, "w") as htmlfile:
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
307 htmlfile.write(htmlString)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
308
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
309 #argdict = {'name':2, 'pdfexport': 3, 'strain': 4, 'forward': 5, 'reverse': 6, 'html_file': 7, 'html_resource': 8}
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
310 def transcriptomicProcess(args,dict):
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
311 transcriptMapping(args[dict['name']], args[dict['strain']], args[dict['forward']], args[dict['reverse']]) #uses bowtie
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
312 processSamFiles(args[dict['name']]) #uses samtools
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
313 transcriptAbundance(args[dict['name']],args[dict['strain']]) #uses cufflinks -> ?.cuff/*.*
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
314 cuff_df = convertToFasta(args[dict['name']],args[dict['strain']])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
315 countList, weightList = HMMerMotifSearch(args[dict['name']],args[dict['strain']], cuff_df)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
316 relFreqList = relativeFrequencyTable(countList,args[dict['name']],args[dict['html_resource']])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
317 relWeightList = weightedFrequencyTable(weightList,args[dict['name']],args[dict['html_resource']])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
318 createStackedBar(args[dict['name']],relWeightList, args[dict['strain']],args[dict['pdfexport']],args[dict['html_resource']])
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
319 createHTML(args[dict['name']],args[dict['html_file']],args[dict['html_resource']], relFreqList, relWeightList)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
320
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
321 if __name__ == "__main__":
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
322 #print("Commencing Transcript Mapping")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
323 #transcriptMapping("T_Test", "Transcripts.1","Transcripts.2")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
324 #print("Processimg Sam Files")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
325 #processSamFiles("T_Test")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
326 #print("Assessing Transcript Abundance")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
327 #transcriptAbundance("T_Test")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
328 #print ("Converting to Fasta Subset")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
329 #cuff_df = convertToFasta("T_Test")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
330 #print("Commencing HMMer search")
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
331 #countList, weightList = HMMerMotifSearch("T_Test",cuff_df)
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
332 #relativeFrequencyTable(countList,'T_Test')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
333 #weightedFrequencyTable(weightList,'T_Test')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
334 relFreqList = [0.111842105,0.059210526,0.026315789,0.013157895,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
335 0.006578947,0.013157895,0.032894737,0.019736842,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
336 0.039473684,0.046052632,0.217105263,0.065789474,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
337 0.151315789,0.059210526,0.138157895]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
338
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
339 relWeightList = [0.07532571,0.05900545,0.009601452,0.042357532,0.01236219,0.001675663,0.04109726,
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
340 0.097464248,0.057491666,0.05826875,0.279457473,0.070004772,0.065329007,0.085361298,0.045197529]
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
341
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
342 createStackedBar('T_Test',relWeightList, 'Tc148','PDF_Yes','results')
1e2f57c43854 Uploaded
johnheap
parents:
diff changeset
343 createHTML("t_test","results/t_test.html","results",relFreqList,relWeightList)