Mercurial > repos > artbio > small_rna_map

annotate small_rna_map.py @ 0:1ad5d040f85f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
author artbio
date Mon, 10 Jul 2017 10:31:49 -0400
parents
children 2e0dc6032a98
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
1 import os
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
2 from os.path import basename
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
3 import pysam
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
4 import argparse
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
5 import numpy
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
6 import collections
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
7
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
8
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
9 def Parser():
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
10 the_parser = argparse.ArgumentParser()
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
11 the_parser.add_argument('--input', dest="input", required=True, nargs='+', help='input BAM files')
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
12 the_parser.add_argument('--output', action="store", type=str, help="output tabular file")
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
13 args = the_parser.parse_args()
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
14 return args
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
15
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
16
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
17 def Ref_length(bamfile):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
18 """
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
19 bamfile.references: list of reference names, and bamfile.lengths is list of their lengths
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
20 """
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
21 ref_length={}
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
22 for e1, e2 in zip(bamfile.references, bamfile.lengths):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
23 ref_length[e1]=e2
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
24 return ref_length
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
25
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
26
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
27 def Ref_coord_pol(bamfile):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
28 ref_coord ={}
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
29 rnames =list(bamfile.references)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
30 for read in bamfile:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
31 if not read.is_unmapped:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
32 polarity = "F" if (not read.is_reverse) else "R"
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
33 read_rname = read.reference_name
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
34 read_position = read.pos + 1
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
35 the_key = (read_rname, read_position,polarity)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
36 read_length = read.qlen
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
37 try:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
38 ref_coord[the_key][0] = ref_coord[the_key][0]+1
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
39 ref_coord[the_key][1].append(read_length)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
40 except:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
41 ref_coord[the_key] = [1,[read_length]]
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
42 if read_rname in rnames:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
43 rnames.remove(read_rname)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
44 for rname in rnames:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
45 ref_coord[(rname, 0, "F")]=[0,[0]]
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
46 return collections.OrderedDict(sorted(ref_coord.items())) #dictionary {(read_rname,read_position,polarity):[nmbr_read,list_of_alignment_length]}
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
47
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
48
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
49 def calcul_stat(dictionary,ref):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
50 dico = {}
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
51 for keys in dictionary:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
52 key = keys[0] #reference name
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
53 value = dictionary[keys][0] # position
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
54 try:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
55 dico[key].append(value)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
56 except:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
57 dico[key] = []
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
58 dico[key].append(value)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
59 the_max = max(dico[ref])
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
60 return the_max
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
61
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
62
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
63 def filename(infile):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
64 filename = basename(infile.filename)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
65 the_name = filename.split('.')
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
66 return the_name[0]
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
67
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
68
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
69 def main(infile, output):
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
70 with open(output,"w") as outfile:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
71 outfile.write("Dataset" + "\t" + "Chromosome"+ "\t" + "Chrom_length" + "\t" + "Coordinate" + "\t" + "Nbr_reads" + "\t" + "Polarity" + "\t" + "Max" + "\t" + "Mean" + "\t" + "Median" + "\n")
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
72 for bamfile in infile:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
73 with pysam.AlignmentFile(bamfile,"rb") as bamfile:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
74 ref_coord = Ref_coord_pol(bamfile)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
75 ref_length = Ref_length(bamfile)
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
76 tabulation = "\t"
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
77 for read in ref_coord:
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
78 table = (str(filename(bamfile)), str(read[0]), str(ref_length[str(read[0])]), str(read[1]), str(ref_coord[read][0]), str(read[2]), str(calcul_stat(ref_coord,read[0])), str(round(numpy.mean(ref_coord[read][1]),1)), str(numpy.median(ref_coord[read][1])))
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
79 outfile.write(tabulation.join(table) + "\n")
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
80
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
81
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
82 if __name__ == "__main__":
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
83 args= Parser()
1ad5d040f85f planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/small_rna_map commit fa452d860cf550c7524df59e77b36fd39e3e2a45
artbio
parents:
diff changeset
84 main(args.input, args.output)