annotate structurefold/get_reads/get_read.py @ 113:aedb21527abd draft

Uploaded
author tyty
date Tue, 14 Apr 2015 14:09:42 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
113
aedb21527abd Uploaded
tyty
parents:
diff changeset
1 #!/usr/bin/env python
aedb21527abd Uploaded
tyty
parents:
diff changeset
2 # -*- coding: utf-8 -*-
aedb21527abd Uploaded
tyty
parents:
diff changeset
3
aedb21527abd Uploaded
tyty
parents:
diff changeset
4 import sys
aedb21527abd Uploaded
tyty
parents:
diff changeset
5 from Bio import SeqIO
aedb21527abd Uploaded
tyty
parents:
diff changeset
6 import os
aedb21527abd Uploaded
tyty
parents:
diff changeset
7 from read_file import *
aedb21527abd Uploaded
tyty
parents:
diff changeset
8 import random
aedb21527abd Uploaded
tyty
parents:
diff changeset
9 import string
aedb21527abd Uploaded
tyty
parents:
diff changeset
10
aedb21527abd Uploaded
tyty
parents:
diff changeset
11 fasta_file = sys.argv[1]
aedb21527abd Uploaded
tyty
parents:
diff changeset
12 map_file = sys.argv[2]
aedb21527abd Uploaded
tyty
parents:
diff changeset
13 result_file = sys.argv[3]
aedb21527abd Uploaded
tyty
parents:
diff changeset
14
aedb21527abd Uploaded
tyty
parents:
diff changeset
15 syspathrs = os.getcwd()
aedb21527abd Uploaded
tyty
parents:
diff changeset
16
aedb21527abd Uploaded
tyty
parents:
diff changeset
17 os.system("samtools view -F 0xfff "+map_file+"|cut -f 3,4 > "+syspathrs+"map_info.txt")
aedb21527abd Uploaded
tyty
parents:
diff changeset
18
aedb21527abd Uploaded
tyty
parents:
diff changeset
19 fasta_sequences = SeqIO.parse(open(fasta_file),'fasta');
aedb21527abd Uploaded
tyty
parents:
diff changeset
20 length_seq = {};
aedb21527abd Uploaded
tyty
parents:
diff changeset
21 for seq in fasta_sequences:
aedb21527abd Uploaded
tyty
parents:
diff changeset
22 nuc = seq.id;
aedb21527abd Uploaded
tyty
parents:
diff changeset
23 length_seq[nuc] = len(seq.seq.tostring());
aedb21527abd Uploaded
tyty
parents:
diff changeset
24
aedb21527abd Uploaded
tyty
parents:
diff changeset
25
aedb21527abd Uploaded
tyty
parents:
diff changeset
26
aedb21527abd Uploaded
tyty
parents:
diff changeset
27 mapping = {}
aedb21527abd Uploaded
tyty
parents:
diff changeset
28 transcripts = []
aedb21527abd Uploaded
tyty
parents:
diff changeset
29
aedb21527abd Uploaded
tyty
parents:
diff changeset
30 f = open(syspathrs+"map_info.txt");
aedb21527abd Uploaded
tyty
parents:
diff changeset
31 for aline in f.readlines():
aedb21527abd Uploaded
tyty
parents:
diff changeset
32 tline = aline.strip();
aedb21527abd Uploaded
tyty
parents:
diff changeset
33 tl = tline.split('\t');
aedb21527abd Uploaded
tyty
parents:
diff changeset
34 if tl[0].strip() not in transcripts:
aedb21527abd Uploaded
tyty
parents:
diff changeset
35 transcripts.append(tl[0].strip());
aedb21527abd Uploaded
tyty
parents:
diff changeset
36 mapping[tl[0].strip()] = [];
aedb21527abd Uploaded
tyty
parents:
diff changeset
37
aedb21527abd Uploaded
tyty
parents:
diff changeset
38 mapping[tl[0].strip()].append(tl[1].strip());
aedb21527abd Uploaded
tyty
parents:
diff changeset
39
aedb21527abd Uploaded
tyty
parents:
diff changeset
40 distribution = {};
aedb21527abd Uploaded
tyty
parents:
diff changeset
41 coverage = {};
aedb21527abd Uploaded
tyty
parents:
diff changeset
42 for transcript in length_seq:
aedb21527abd Uploaded
tyty
parents:
diff changeset
43 distribution[transcript] = [];
aedb21527abd Uploaded
tyty
parents:
diff changeset
44 for i in range(0, length_seq[transcript]):
aedb21527abd Uploaded
tyty
parents:
diff changeset
45 distribution[transcript].append(0);
aedb21527abd Uploaded
tyty
parents:
diff changeset
46 sum_count = float(0);
aedb21527abd Uploaded
tyty
parents:
diff changeset
47 if transcript in mapping:
aedb21527abd Uploaded
tyty
parents:
diff changeset
48 for j in range(0, len(mapping[transcript])):
aedb21527abd Uploaded
tyty
parents:
diff changeset
49 index = mapping[transcript][j];
aedb21527abd Uploaded
tyty
parents:
diff changeset
50 #count = reads[mapping[transcript][j][0]];
aedb21527abd Uploaded
tyty
parents:
diff changeset
51 sum_count = sum_count + 1;
aedb21527abd Uploaded
tyty
parents:
diff changeset
52 distribution[transcript][int(index)-1] = distribution[transcript][int(index)-1] + 1;
aedb21527abd Uploaded
tyty
parents:
diff changeset
53 coverage[transcript] = float(sum_count)/float(length_seq[transcript]);
aedb21527abd Uploaded
tyty
parents:
diff changeset
54 else:
aedb21527abd Uploaded
tyty
parents:
diff changeset
55 coverage[transcript] = 0
aedb21527abd Uploaded
tyty
parents:
diff changeset
56
aedb21527abd Uploaded
tyty
parents:
diff changeset
57
aedb21527abd Uploaded
tyty
parents:
diff changeset
58
aedb21527abd Uploaded
tyty
parents:
diff changeset
59
aedb21527abd Uploaded
tyty
parents:
diff changeset
60
aedb21527abd Uploaded
tyty
parents:
diff changeset
61 h = file(result_file, 'w')
aedb21527abd Uploaded
tyty
parents:
diff changeset
62 for transcript in length_seq:
aedb21527abd Uploaded
tyty
parents:
diff changeset
63 h.write(transcript);
aedb21527abd Uploaded
tyty
parents:
diff changeset
64 h.write('\n')
aedb21527abd Uploaded
tyty
parents:
diff changeset
65 for i in range(0, length_seq[transcript]):
aedb21527abd Uploaded
tyty
parents:
diff changeset
66 h.write(str(distribution[transcript][i]))
aedb21527abd Uploaded
tyty
parents:
diff changeset
67 h.write('\t')
aedb21527abd Uploaded
tyty
parents:
diff changeset
68 h.write('\n')
aedb21527abd Uploaded
tyty
parents:
diff changeset
69 h.write('\n')
aedb21527abd Uploaded
tyty
parents:
diff changeset
70
aedb21527abd Uploaded
tyty
parents:
diff changeset
71 #os.system("rm -r "+syspathrs)
aedb21527abd Uploaded
tyty
parents:
diff changeset
72
aedb21527abd Uploaded
tyty
parents:
diff changeset
73
aedb21527abd Uploaded
tyty
parents:
diff changeset
74
aedb21527abd Uploaded
tyty
parents:
diff changeset
75 f.close();
aedb21527abd Uploaded
tyty
parents:
diff changeset
76 h.close()
aedb21527abd Uploaded
tyty
parents:
diff changeset
77
aedb21527abd Uploaded
tyty
parents:
diff changeset
78
aedb21527abd Uploaded
tyty
parents:
diff changeset
79
aedb21527abd Uploaded
tyty
parents:
diff changeset
80