msp_sr_bowtie_parser: smRtools.py annotate

author	drosofff
date	Wed, 27 May 2015 17:19:15 -0400
parents
children	ca3845fb0b31

rev	line source
0 b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	1 #!/usr/bin/python
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	2 # version 1 7-5-2012 unification of the SmRNAwindow class
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	3
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	4 import sys, subprocess
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	5 from collections import defaultdict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	6 from numpy import mean, median, std
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	7 from scipy import stats
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	8
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	9 def get_fasta (index="/home/galaxy/galaxy-dist/bowtie/5.37_Dmel/5.37_Dmel"):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	10 '''This function will return a dictionary containing fasta identifiers as keys and the
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	11 sequence as values. Index must be the path to a fasta file.'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	12 p = subprocess.Popen(args=["bowtie-inspect","-a", "0", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	13 outputlines = p.stdout.readlines()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	14 p.wait()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	15 item_dic = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	16 for line in outputlines:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	17 if (line[0] == ">"):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	18 try:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	19 item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	20 except: pass
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	21 current_item = line[1:].rstrip().split()[0] #take the first word before space because bowtie splits headers !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	22 item_dic[current_item] = ""
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	23 stringlist=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	24 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	25 stringlist.append(line.rstrip() )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	26 item_dic[current_item] = "".join(stringlist) # for the last item
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	27 return item_dic
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	28
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	29 def get_fasta_headers (index):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	30 p = subprocess.Popen(args=["bowtie-inspect","-n", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	31 outputlines = p.stdout.readlines()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	32 p.wait()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	33 item_dic = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	34 for line in outputlines:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	35 header = line.rstrip().split()[0] #take the first word before space because bowtie splits headers !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	36 item_dic[header] = 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	37 return item_dic
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	38
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	39
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	40 def get_file_sample (file, numberoflines):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	41 '''import random to use this function'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	42 F=open(file)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	43 fullfile = F.read().splitlines()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	44 F.close()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	45 if len(fullfile) < numberoflines:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	46 return "sample size exceeds file size"
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	47 return random.sample(fullfile, numberoflines)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	48
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	49 def get_fasta_from_history (file):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	50 F = open (file, "r")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	51 item_dic = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	52 for line in F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	53 if (line[0] == ">"):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	54 try:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	55 item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	56 except: pass
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	57 current_item = line[1:-1].split()[0] #take the first word before space because bowtie splits headers !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	58 item_dic[current_item] = ""
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	59 stringlist=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	60 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	61 stringlist.append(line[:-1])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	62 item_dic[current_item] = "".join(stringlist) # for the last item
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	63 return item_dic
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	64
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	65 def antipara (sequence):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	66 antidict = {"A":"T", "T":"A", "G":"C", "C":"G", "N":"N"}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	67 revseq = sequence[::-1]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	68 return "".join([antidict[i] for i in revseq])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	69
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	70 def RNAtranslate (sequence):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	71 return "".join([i if i in "AGCN" else "U" for i in sequence])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	72 def DNAtranslate (sequence):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	73 return "".join([i if i in "AGCN" else "T" for i in sequence])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	74
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	75 def RNAfold (sequence_list):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	76 thestring= "\n".join(sequence_list)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	77 p = subprocess.Popen(args=["RNAfold","--noPS"], stdin= subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	78 output=p.communicate(thestring)[0]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	79 p.wait()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	80 output=output.split("\n")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	81 if not output[-1]: output = output[:-1] # nasty patch to remove last empty line
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	82 buffer=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	83 for line in output:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	84 if line[0] in ["N","A","T","U","G","C"]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	85 buffer.append(DNAtranslate(line))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	86 if line[0] in ["(",".",")"]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	87 fields=line.split("(")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	88 energy= fields[-1]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	89 energy = energy[:-1] # remove the ) parenthesis
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	90 energy=float(energy)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	91 buffer.append(str(energy))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	92 return dict(zip(buffer[::2], buffer[1::2]))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	93
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	94 def extractsubinstance (start, end, instance):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	95 ''' Testing whether this can be an function external to the class to save memory'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	96 subinstance = SmRNAwindow (instance.gene, instance.sequence[start-1:end], start)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	97 subinstance.gene = "%s %s %s" % (subinstance.gene, subinstance.windowoffset, subinstance.windowoffset + subinstance.size - 1)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	98 upcoordinate = [i for i in range(start,end+1) if instance.readDict.has_key(i) ]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	99 downcoordinate = [-i for i in range(start,end+1) if instance.readDict.has_key(-i) ]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	100 for i in upcoordinate:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	101 subinstance.readDict[i]=instance.readDict[i]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	102 for i in downcoordinate:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	103 subinstance.readDict[i]=instance.readDict[i]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	104 return subinstance
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	105
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	106 class HandleSmRNAwindows:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	107 def __init__(self, alignmentFile="~", alignmentFileFormat="tabular", genomeRefFile="~", genomeRefFormat="bowtieIndex", biosample="undetermined", size_inf=None, size_sup=1000, norm=1.0):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	108 self.biosample = biosample
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	109 self.alignmentFile = alignmentFile
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	110 self.alignmentFileFormat = alignmentFileFormat # can be "tabular" or "sam"
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	111 self.genomeRefFile = genomeRefFile
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	112 self.genomeRefFormat = genomeRefFormat # can be "bowtieIndex" or "fastaSource"
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	113 self.alignedReads = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	114 self.instanceDict = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	115 self.size_inf=size_inf
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	116 self.size_sup=size_sup
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	117 self.norm=norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	118 if genomeRefFormat == "bowtieIndex":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	119 self.itemDict = get_fasta (genomeRefFile)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	120 elif genomeRefFormat == "fastaSource":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	121 self.itemDict = get_fasta_from_history (genomeRefFile)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	122 for item in self.itemDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	123 self.instanceDict[item] = SmRNAwindow(item, sequence=self.itemDict[item], windowoffset=1, biosample=self.biosample, norm=self.norm) # create as many instances as there is items
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	124 self.readfile()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	125
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	126 def readfile (self) :
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	127 if self.alignmentFileFormat == "tabular":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	128 F = open (self.alignmentFile, "r")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	129 for line in F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	130 fields = line.split()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	131 polarity = fields[1]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	132 gene = fields[2]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	133 offset = int(fields[3])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	134 size = len (fields[4])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	135 if self.size_inf:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	136 if (size>=self.size_inf and size<= self.size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	137 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	138 self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	139 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	140 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	141 self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	142 F.close()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	143 return self.instanceDict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	144 # elif self.alignmentFileFormat == "sam":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	145 # F = open (self.alignmentFile, "r")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	146 # dict = {"0":"+", "16":"-"}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	147 # for line in F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	148 # if line[0]=='@':
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	149 # continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	150 # fields = line.split()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	151 # if fields[2] == "*": continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	152 # polarity = dict[fields[1]]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	153 # gene = fields[2]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	154 # offset = int(fields[3])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	155 # size = len (fields[9])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	156 # if self.size_inf:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	157 # if (size>=self.size_inf and size<= self.size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	158 # self.instanceDict[gene].addread (polarity, offset, size)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	159 # self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	160 # else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	161 # self.instanceDict[gene].addread (polarity, offset, size)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	162 # self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	163 # F.close()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	164 elif self.alignmentFileFormat == "bam" or self.alignmentFileFormat == "sam":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	165 import pysam
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	166 samfile = pysam.Samfile(self.alignmentFile)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	167 for read in samfile:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	168 if read.tid == -1:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	169 continue # filter out unaligned reads
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	170 if read.is_reverse:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	171 polarity="-"
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	172 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	173 polarity="+"
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	174 gene = samfile.getrname(read.tid)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	175 offset = read.pos
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	176 size = read.qlen
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	177 if self.size_inf:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	178 if (size>=self.size_inf and size<= self.size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	179 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	180 self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	181 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	182 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	183 self.alignedReads += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	184 return self.instanceDict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	185
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	186 # def size_histogram (self):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	187 # size_dict={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	188 # size_dict['F']= defaultdict (int)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	189 # size_dict['R']= defaultdict (int)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	190 # size_dict['both'] = defaultdict (int)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	191 # for item in self.instanceDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	192 # buffer_dict_F = self.instanceDict[item].size_histogram()['F']
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	193 # buffer_dict_R = self.instanceDict[item].size_histogram()['R']
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	194 # for size in buffer_dict_F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	195 # size_dict['F'][size] += buffer_dict_F[size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	196 # for size in buffer_dict_R:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	197 # size_dict['R'][size] -= buffer_dict_R[size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	198 # allSizeKeys = list (set (size_dict['F'].keys() + size_dict['R'].keys() ) )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	199 # for size in allSizeKeys:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	200 # size_dict['both'][size] = size_dict['F'][size] + size_dict['R'][size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	201 # return size_dict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	202 def size_histogram (self): # in HandleSmRNAwindows
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	203 '''refactored on 7-9-2014 to debug size_histogram tool'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	204 size_dict={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	205 size_dict['F']= defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	206 size_dict['R']= defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	207 size_dict['both'] = defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	208 for item in self.instanceDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	209 buffer_dict = self.instanceDict[item].size_histogram()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	210 for polarity in ["F", "R"]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	211 for size in buffer_dict[polarity]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	212 size_dict[polarity][size] += buffer_dict[polarity][size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	213 for size in buffer_dict["both"]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	214 size_dict["both"][size] += buffer_dict["F"][size] - buffer_dict["R"][size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	215 return size_dict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	216
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	217 def CountFeatures (self, GFF3="path/to/file"):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	218 featureDict = defaultdict(int)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	219 F = open (GFF3, "r")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	220 for line in F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	221 if line[0] == "#": continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	222 fields = line[:-1].split()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	223 chrom, feature, leftcoord, rightcoord, polarity = fields[0], fields[2], fields[3], fields[4], fields[6]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	224 featureDict[feature] += self.instanceDict[chrom].readcount(upstream_coord=int(leftcoord), downstream_coord=int(rightcoord), polarity="both", method="destructive")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	225 F.close()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	226 return featureDict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	227
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	228 class SmRNAwindow:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	229
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	230 def __init__(self, gene, sequence="ATGC", windowoffset=1, biosample="Undetermined", norm=1.0):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	231 self.biosample = biosample
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	232 self.sequence = sequence
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	233 self.gene = gene
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	234 self.windowoffset = windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	235 self.size = len(sequence)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	236 self.readDict = defaultdict(list) # with a {+/-offset:[size1, size2, ...], ...}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	237 self.matchedreadsUp = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	238 self.matchedreadsDown = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	239 self.norm=norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	240
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	241 def addread (self, polarity, offset, size):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	242 '''ATTENTION ATTENTION ATTENTION'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	243 ''' We removed the conversion from 0 to 1 based offset, as we do this now during readparsing.'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	244 if polarity == "+":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	245 self.readDict[offset].append(size)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	246 self.matchedreadsUp += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	247 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	248 self.readDict[-(offset + size -1)].append(size)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	249 self.matchedreadsDown += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	250 return
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	251
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	252 def barycenter (self, upstream_coord=None, downstream_coord=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	253 '''refactored 24-12-2013 to save memory and introduce offset filtering see readcount method for further discussion on that
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	254 In this version, attempt to replace the dictionary structure by a list of tupple to save memory too'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	255 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	256 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	257 window_size = downstream_coord - upstream_coord +1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	258 def weigthAverage (TuppleList):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	259 weightSum = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	260 PonderWeightSum = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	261 for tuple in TuppleList:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	262 PonderWeightSum += tuple[0] * tuple[1]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	263 weightSum += tuple[1]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	264 if weightSum > 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	265 return PonderWeightSum / float(weightSum)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	266 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	267 return 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	268 forwardTuppleList = [(k, len(self.readDict[k])) for k in self.readDict.keys() if (k > 0 and abs(k) >= upstream_coord and abs(k) <= downstream_coord)] # both forward and in the proper offset window
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	269 reverseTuppleList = [(-k, len(self.readDict[k])) for k in self.readDict.keys() if (k < 0 and abs(k) >= upstream_coord and abs(k) <= downstream_coord)] # both reverse and in the proper offset window
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	270 Fbarycenter = (weigthAverage (forwardTuppleList) - upstream_coord) / window_size
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	271 Rbarycenter = (weigthAverage (reverseTuppleList) - upstream_coord) / window_size
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	272 return Fbarycenter, Rbarycenter
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	273
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	274 def correlation_mapper (self, reference, window_size):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	275 '''to map correlation with a sliding window 26-2-2013'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	276 if window_size > self.size:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	277 return []
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	278 F=open(reference, "r")
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	279 reference_forward = []
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	280 reference_reverse = []
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	281 for line in F:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	282 fields=line.split()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	283 reference_forward.append(int(float(fields[1])))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	284 reference_reverse.append(int(float(fields[2])))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	285 F.close()
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	286 local_object_forward=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	287 local_object_reverse=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	288 ## Dict to list for the local object
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	289 for i in range(1, self.size+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	290 local_object_forward.append(len(self.readDict[i]))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	291 local_object_reverse.append(len(self.readDict[-i]))
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	292 ## start compiling results by slides
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	293 results=[]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	294 for coordinate in range(self.size - window_size):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	295 local_forward=local_object_forward[coordinate:coordinate + window_size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	296 local_reverse=local_object_reverse[coordinate:coordinate + window_size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	297 if sum(local_forward) == 0 or sum(local_reverse) == 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	298 continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	299 try:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	300 reference_to_local_cor_forward = stats.spearmanr(local_forward, reference_forward)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	301 reference_to_local_cor_reverse = stats.spearmanr(local_reverse, reference_reverse)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	302 if (reference_to_local_cor_forward[0] > 0.2 or reference_to_local_cor_reverse[0]>0.2):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	303 results.append([coordinate+1, reference_to_local_cor_forward[0], reference_to_local_cor_reverse[0]])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	304 except:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	305 pass
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	306 return results
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	307
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	308 def readcount (self, size_inf=0, size_sup=1000, upstream_coord=None, downstream_coord=None, polarity="both", method="conservative"):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	309 '''refactored 24-12-2013 to save memory and introduce offset filtering
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	310 take a look at the defaut parameters that cannot be defined relatively to the instance are they are defined before instanciation
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	311 the trick is to pass None and then test
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	312 polarity parameter can take "both", "forward" or "reverse" as value'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	313 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	314 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	315 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "both":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	316 return self.matchedreadsUp + self.matchedreadsDown
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	317 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "forward":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	318 return self.matchedreadsUp
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	319 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "reverse":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	320 return self.matchedreadsDown
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	321 n=0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	322 if polarity == "both":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	323 for offset in xrange(upstream_coord, downstream_coord+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	324 if self.readDict.has_key(offset):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	325 for read in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	326 if (read>=size_inf and read<= size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	327 n += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	328 if method != "conservative":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	329 del self.readDict[offset] ## Carefull ! precludes re-use on the self.readDict dictionary !!!!!! TEST
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	330 if self.readDict.has_key(-offset):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	331 for read in self.readDict[-offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	332 if (read>=size_inf and read<= size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	333 n += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	334 if method != "conservative":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	335 del self.readDict[-offset]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	336 return n
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	337 elif polarity == "forward":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	338 for offset in xrange(upstream_coord, downstream_coord+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	339 if self.readDict.has_key(offset):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	340 for read in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	341 if (read>=size_inf and read<= size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	342 n += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	343 return n
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	344 elif polarity == "reverse":
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	345 for offset in xrange(upstream_coord, downstream_coord+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	346 if self.readDict.has_key(-offset):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	347 for read in self.readDict[-offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	348 if (read>=size_inf and read<= size_sup):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	349 n += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	350 return n
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	351
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	352 def readsizes (self):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	353 '''return a dictionary of number of reads by size (the keys)'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	354 dicsize = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	355 for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	356 for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	357 dicsize[size] = dicsize.get(size, 0) + 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	358 for offset in range (min(dicsize.keys()), max(dicsize.keys())+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	359 dicsize[size] = dicsize.get(size, 0) # to fill offsets with null values
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	360 return dicsize
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	361
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	362 # def size_histogram(self):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	363 # norm=self.norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	364 # hist_dict={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	365 # hist_dict['F']={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	366 # hist_dict['R']={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	367 # for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	368 # for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	369 # if offset < 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	370 # hist_dict['R'][size] = hist_dict['R'].get(size, 0) - 1*norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	371 # else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	372 # hist_dict['F'][size] = hist_dict['F'].get(size, 0) + 1*norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	373 # ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	374 # if not (hist_dict['F']) and (not hist_dict['R']):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	375 # hist_dict['F'][21] = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	376 # hist_dict['R'][21] = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	377 # ##
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	378 # return hist_dict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	379
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	380 def size_histogram(self, minquery=None, maxquery=None): # in SmRNAwindow
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	381 '''refactored on 7-9-2014 to debug size_histogram tool'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	382 norm=self.norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	383 size_dict={}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	384 size_dict['F']= defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	385 size_dict['R']= defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	386 size_dict['both']= defaultdict (float)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	387 for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	388 for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	389 if offset < 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	390 size_dict['R'][size] = size_dict['R'][size] - 1*norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	391 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	392 size_dict['F'][size] = size_dict['F'][size] + 1*norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	393 ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	394 if not (size_dict['F']) and (not size_dict['R']):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	395 size_dict['F'][21] = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	396 size_dict['R'][21] = 0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	397 ##
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	398 allSizeKeys = list (set (size_dict['F'].keys() + size_dict['R'].keys() ) )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	399 for size in allSizeKeys:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	400 size_dict['both'][size] = size_dict['F'][size] - size_dict['R'][size]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	401 if minquery:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	402 for polarity in size_dict.keys():
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	403 for size in xrange(minquery, maxquery+1):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	404 if not size in size_dict[polarity].keys():
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	405 size_dict[polarity][size]=0
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	406 return size_dict
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	407
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	408 def statsizes (self, upstream_coord=None, downstream_coord=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	409 ''' migration to memory saving by specifying possible subcoordinates
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	410 see the readcount method for further discussion'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	411 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	412 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	413 L = []
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	414 for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	415 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	416 for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	417 L.append(size)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	418 meansize = mean(L)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	419 stdv = std(L)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	420 mediansize = median(L)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	421 return meansize, mediansize, stdv
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	422
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	423 def foldEnergy (self, upstream_coord=None, downstream_coord=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	424 ''' migration to memory saving by specifying possible subcoordinates
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	425 see the readcount method for further discussion'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	426 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	427 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	428 Energy = RNAfold ([self.sequence[upstream_coord-1:downstream_coord] ])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	429 return float(Energy[self.sequence[upstream_coord-1:downstream_coord]])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	430
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	431 def Ufreq (self, size_scope, upstream_coord=None, downstream_coord=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	432 ''' migration to memory saving by specifying possible subcoordinates
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	433 see the readcount method for further discussion. size_scope must be an interable'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	434 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	435 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	436 freqDic = {"A":0,"T":0,"G":0,"C":0, "N":0}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	437 convertDic = {"A":"T","T":"A","G":"C","C":"G","N":"N"}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	438 for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	439 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	440 for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	441 if size in size_scope:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	442 startbase = self.sequence[abs(offset)-self.windowoffset]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	443 if offset < 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	444 startbase = convertDic[startbase]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	445 freqDic[startbase] += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	446 base_sum = float ( sum( freqDic.values()) )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	447 if base_sum == 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	448 return "."
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	449 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	450 return freqDic["T"] / base_sum * 100
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	451
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	452 def Ufreq_stranded (self, size_scope, upstream_coord=None, downstream_coord=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	453 ''' migration to memory saving by specifying possible subcoordinates
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	454 see the readcount method for further discussion. size_scope must be an interable
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	455 This method is similar to the Ufreq method but take strandness into account'''
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	456 upstream_coord = upstream_coord or self.windowoffset
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	457 downstream_coord = downstream_coord or self.windowoffset+self.size-1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	458 freqDic = {"Afor":0,"Tfor":0,"Gfor":0,"Cfor":0, "Nfor":0,"Arev":0,"Trev":0,"Grev":0,"Crev":0, "Nrev":0}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	459 convertDic = {"A":"T","T":"A","G":"C","C":"G","N":"N"}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	460 for offset in self.readDict:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	461 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	462 for size in self.readDict[offset]:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	463 if size in size_scope:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	464 startbase = self.sequence[abs(offset)-self.windowoffset]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	465 if offset < 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	466 startbase = convertDic[startbase]
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	467 freqDic[startbase+"rev"] += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	468 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	469 freqDic[startbase+"for"] += 1
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	470 forward_sum = float ( freqDic["Afor"]+freqDic["Tfor"]+freqDic["Gfor"]+freqDic["Cfor"]+freqDic["Nfor"])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	471 reverse_sum = float ( freqDic["Arev"]+freqDic["Trev"]+freqDic["Grev"]+freqDic["Crev"]+freqDic["Nrev"])
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	472 if forward_sum == 0 and reverse_sum == 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	473 return ". \| ."
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	474 elif reverse_sum == 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	475 return "%s \| ." % (freqDic["Tfor"] / forward_sum * 100)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	476 elif forward_sum == 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	477 return ". \| %s" % (freqDic["Trev"] / reverse_sum * 100)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	478 else:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	479 return "%s \| %s" % (freqDic["Tfor"] / forward_sum * 100, freqDic["Trev"] / reverse_sum * 100)
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	480
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	481
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	482 def readplot (self):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	483 norm=self.norm
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	484 readmap = {}
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	485 for offset in self.readDict.keys():
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	486 readmap[abs(offset)] = ( len(self.readDict.get(-abs(offset),[]))norm , len(self.readDict.get(abs(offset),[]))norm )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	487 mylist = []
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	488 for offset in sorted(readmap):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	489 if readmap[offset][1] != 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	490 mylist.append("%s\t%s\t%s\t%s" % (self.gene, offset, readmap[offset][1], "F") )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	491 if readmap[offset][0] != 0:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	492 mylist.append("%s\t%s\t%s\t%s" % (self.gene, offset, -readmap[offset][0], "R") )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	493 ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	494 if not mylist:
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	495 mylist.append("%s\t%s\t%s\t%s" % (self.gene, 1, 0, "F") )
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	496 ###
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	497 return mylist
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	498
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	499 def readcoverage (self, upstream_coord=None, downstream_coord=None, windowName=None):
b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/ drosofff parents: diff changeset	500 '''Use by MirParser tool'''

0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

1 #!/usr/bin/python

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

2 # version 1 7-5-2012 unification of the SmRNAwindow class

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

3

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

4 import sys, subprocess

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

5 from collections import defaultdict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

6 from numpy import mean, median, std

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

7 from scipy import stats

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

8

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

9 def get_fasta (index="/home/galaxy/galaxy-dist/bowtie/5.37_Dmel/5.37_Dmel"):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

10 '''This function will return a dictionary containing fasta identifiers as keys and the

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

11 sequence as values. Index must be the path to a fasta file.'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

12 p = subprocess.Popen(args=["bowtie-inspect","-a", "0", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

13 outputlines = p.stdout.readlines()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

14 p.wait()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

15 item_dic = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

16 for line in outputlines:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

17 if (line[0] == ">"):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

18 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

19 item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

20 except: pass

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

21 current_item = line[1:].rstrip().split()[0] #take the first word before space because bowtie splits headers !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

22 item_dic[current_item] = ""

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

23 stringlist=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

24 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

25 stringlist.append(line.rstrip() )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

26 item_dic[current_item] = "".join(stringlist) # for the last item

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

27 return item_dic

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

28

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

29 def get_fasta_headers (index):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

30 p = subprocess.Popen(args=["bowtie-inspect","-n", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

31 outputlines = p.stdout.readlines()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

32 p.wait()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

33 item_dic = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

34 for line in outputlines:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

35 header = line.rstrip().split()[0] #take the first word before space because bowtie splits headers !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

36 item_dic[header] = 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

37 return item_dic

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

38

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

39

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

40 def get_file_sample (file, numberoflines):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

41 '''import random to use this function'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

42 F=open(file)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

43 fullfile = F.read().splitlines()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

44 F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

45 if len(fullfile) < numberoflines:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

46 return "sample size exceeds file size"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

47 return random.sample(fullfile, numberoflines)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

48

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

49 def get_fasta_from_history (file):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

50 F = open (file, "r")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

51 item_dic = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

52 for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

53 if (line[0] == ">"):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

54 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

55 item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

56 except: pass

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

57 current_item = line[1:-1].split()[0] #take the first word before space because bowtie splits headers !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

58 item_dic[current_item] = ""

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

59 stringlist=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

60 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

61 stringlist.append(line[:-1])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

62 item_dic[current_item] = "".join(stringlist) # for the last item

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

63 return item_dic

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

64

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

65 def antipara (sequence):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

66 antidict = {"A":"T", "T":"A", "G":"C", "C":"G", "N":"N"}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

67 revseq = sequence[::-1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

68 return "".join([antidict[i] for i in revseq])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

69

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

70 def RNAtranslate (sequence):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

71 return "".join([i if i in "AGCN" else "U" for i in sequence])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

72 def DNAtranslate (sequence):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

73 return "".join([i if i in "AGCN" else "T" for i in sequence])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

74

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

75 def RNAfold (sequence_list):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

76 thestring= "\n".join(sequence_list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

77 p = subprocess.Popen(args=["RNAfold","--noPS"], stdin= subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

78 output=p.communicate(thestring)[0]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

79 p.wait()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

80 output=output.split("\n")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

81 if not output[-1]: output = output[:-1] # nasty patch to remove last empty line

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

82 buffer=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

83 for line in output:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

84 if line[0] in ["N","A","T","U","G","C"]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

85 buffer.append(DNAtranslate(line))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

86 if line[0] in ["(",".",")"]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

87 fields=line.split("(")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

88 energy= fields[-1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

89 energy = energy[:-1] # remove the ) parenthesis

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

90 energy=float(energy)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

91 buffer.append(str(energy))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

92 return dict(zip(buffer[::2], buffer[1::2]))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

93

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

94 def extractsubinstance (start, end, instance):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

95 ''' Testing whether this can be an function external to the class to save memory'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

96 subinstance = SmRNAwindow (instance.gene, instance.sequence[start-1:end], start)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

97 subinstance.gene = "%s %s %s" % (subinstance.gene, subinstance.windowoffset, subinstance.windowoffset + subinstance.size - 1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

98 upcoordinate = [i for i in range(start,end+1) if instance.readDict.has_key(i) ]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

99 downcoordinate = [-i for i in range(start,end+1) if instance.readDict.has_key(-i) ]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

100 for i in upcoordinate:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

101 subinstance.readDict[i]=instance.readDict[i]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

102 for i in downcoordinate:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

103 subinstance.readDict[i]=instance.readDict[i]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

104 return subinstance

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

105

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

106 class HandleSmRNAwindows:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

107 def __init__(self, alignmentFile="~", alignmentFileFormat="tabular", genomeRefFile="~", genomeRefFormat="bowtieIndex", biosample="undetermined", size_inf=None, size_sup=1000, norm=1.0):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

108 self.biosample = biosample

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

109 self.alignmentFile = alignmentFile

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

110 self.alignmentFileFormat = alignmentFileFormat # can be "tabular" or "sam"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

111 self.genomeRefFile = genomeRefFile

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

112 self.genomeRefFormat = genomeRefFormat # can be "bowtieIndex" or "fastaSource"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

113 self.alignedReads = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

114 self.instanceDict = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

115 self.size_inf=size_inf

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

116 self.size_sup=size_sup

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

117 self.norm=norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

118 if genomeRefFormat == "bowtieIndex":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

119 self.itemDict = get_fasta (genomeRefFile)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

120 elif genomeRefFormat == "fastaSource":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

121 self.itemDict = get_fasta_from_history (genomeRefFile)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

122 for item in self.itemDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

123 self.instanceDict[item] = SmRNAwindow(item, sequence=self.itemDict[item], windowoffset=1, biosample=self.biosample, norm=self.norm) # create as many instances as there is items

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

124 self.readfile()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

125

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

126 def readfile (self) :

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

127 if self.alignmentFileFormat == "tabular":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

128 F = open (self.alignmentFile, "r")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

129 for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

130 fields = line.split()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

131 polarity = fields[1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

132 gene = fields[2]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

133 offset = int(fields[3])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

134 size = len (fields[4])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

135 if self.size_inf:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

136 if (size>=self.size_inf and size<= self.size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

137 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

138 self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

139 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

140 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

141 self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

142 F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

143 return self.instanceDict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

144 # elif self.alignmentFileFormat == "sam":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

145 # F = open (self.alignmentFile, "r")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

146 # dict = {"0":"+", "16":"-"}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

147 # for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

148 # if line[0]=='@':

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

149 # continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

150 # fields = line.split()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

151 # if fields[2] == "*": continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

152 # polarity = dict[fields[1]]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

153 # gene = fields[2]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

154 # offset = int(fields[3])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

155 # size = len (fields[9])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

156 # if self.size_inf:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

157 # if (size>=self.size_inf and size<= self.size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

158 # self.instanceDict[gene].addread (polarity, offset, size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

159 # self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

160 # else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

161 # self.instanceDict[gene].addread (polarity, offset, size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

162 # self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

163 # F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

164 elif self.alignmentFileFormat == "bam" or self.alignmentFileFormat == "sam":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

165 import pysam

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

166 samfile = pysam.Samfile(self.alignmentFile)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

167 for read in samfile:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

168 if read.tid == -1:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

169 continue # filter out unaligned reads

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

170 if read.is_reverse:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

171 polarity="-"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

172 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

173 polarity="+"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

174 gene = samfile.getrname(read.tid)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

175 offset = read.pos

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

176 size = read.qlen

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

177 if self.size_inf:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

178 if (size>=self.size_inf and size<= self.size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

179 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

180 self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

181 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

182 self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

183 self.alignedReads += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

184 return self.instanceDict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

185

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

186 # def size_histogram (self):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

187 # size_dict={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

188 # size_dict['F']= defaultdict (int)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

189 # size_dict['R']= defaultdict (int)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

190 # size_dict['both'] = defaultdict (int)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

191 # for item in self.instanceDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

192 # buffer_dict_F = self.instanceDict[item].size_histogram()['F']

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

193 # buffer_dict_R = self.instanceDict[item].size_histogram()['R']

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

194 # for size in buffer_dict_F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

195 # size_dict['F'][size] += buffer_dict_F[size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

196 # for size in buffer_dict_R:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

197 # size_dict['R'][size] -= buffer_dict_R[size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

198 # allSizeKeys = list (set (size_dict['F'].keys() + size_dict['R'].keys() ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

199 # for size in allSizeKeys:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

200 # size_dict['both'][size] = size_dict['F'][size] + size_dict['R'][size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

201 # return size_dict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

202 def size_histogram (self): # in HandleSmRNAwindows

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

203 '''refactored on 7-9-2014 to debug size_histogram tool'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

204 size_dict={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

205 size_dict['F']= defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

206 size_dict['R']= defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

207 size_dict['both'] = defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

208 for item in self.instanceDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

209 buffer_dict = self.instanceDict[item].size_histogram()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

210 for polarity in ["F", "R"]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

211 for size in buffer_dict[polarity]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

212 size_dict[polarity][size] += buffer_dict[polarity][size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

213 for size in buffer_dict["both"]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

214 size_dict["both"][size] += buffer_dict["F"][size] - buffer_dict["R"][size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

215 return size_dict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

216

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

217 def CountFeatures (self, GFF3="path/to/file"):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

218 featureDict = defaultdict(int)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

219 F = open (GFF3, "r")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

220 for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

221 if line[0] == "#": continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

222 fields = line[:-1].split()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

223 chrom, feature, leftcoord, rightcoord, polarity = fields[0], fields[2], fields[3], fields[4], fields[6]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

224 featureDict[feature] += self.instanceDict[chrom].readcount(upstream_coord=int(leftcoord), downstream_coord=int(rightcoord), polarity="both", method="destructive")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

225 F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

226 return featureDict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

227

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

228 class SmRNAwindow:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

229

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

230 def __init__(self, gene, sequence="ATGC", windowoffset=1, biosample="Undetermined", norm=1.0):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

231 self.biosample = biosample

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

232 self.sequence = sequence

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

233 self.gene = gene

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

234 self.windowoffset = windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

235 self.size = len(sequence)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

236 self.readDict = defaultdict(list) # with a {+/-offset:[size1, size2, ...], ...}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

237 self.matchedreadsUp = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

238 self.matchedreadsDown = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

239 self.norm=norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

240

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

241 def addread (self, polarity, offset, size):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

242 '''ATTENTION ATTENTION ATTENTION'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

243 ''' We removed the conversion from 0 to 1 based offset, as we do this now during readparsing.'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

244 if polarity == "+":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

245 self.readDict[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

246 self.matchedreadsUp += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

247 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

248 self.readDict[-(offset + size -1)].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

249 self.matchedreadsDown += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

250 return

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

251

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

252 def barycenter (self, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

253 '''refactored 24-12-2013 to save memory and introduce offset filtering see readcount method for further discussion on that

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

254 In this version, attempt to replace the dictionary structure by a list of tupple to save memory too'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

255 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

256 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

257 window_size = downstream_coord - upstream_coord +1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

258 def weigthAverage (TuppleList):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

259 weightSum = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

260 PonderWeightSum = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

261 for tuple in TuppleList:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

262 PonderWeightSum += tuple[0] * tuple[1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

263 weightSum += tuple[1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

264 if weightSum > 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

265 return PonderWeightSum / float(weightSum)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

266 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

267 return 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

268 forwardTuppleList = [(k, len(self.readDict[k])) for k in self.readDict.keys() if (k > 0 and abs(k) >= upstream_coord and abs(k) <= downstream_coord)] # both forward and in the proper offset window

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

269 reverseTuppleList = [(-k, len(self.readDict[k])) for k in self.readDict.keys() if (k < 0 and abs(k) >= upstream_coord and abs(k) <= downstream_coord)] # both reverse and in the proper offset window

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

270 Fbarycenter = (weigthAverage (forwardTuppleList) - upstream_coord) / window_size

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

271 Rbarycenter = (weigthAverage (reverseTuppleList) - upstream_coord) / window_size

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

272 return Fbarycenter, Rbarycenter

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

273

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

274 def correlation_mapper (self, reference, window_size):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

275 '''to map correlation with a sliding window 26-2-2013'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

276 if window_size > self.size:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

277 return []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

278 F=open(reference, "r")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

279 reference_forward = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

280 reference_reverse = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

281 for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

282 fields=line.split()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

283 reference_forward.append(int(float(fields[1])))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

284 reference_reverse.append(int(float(fields[2])))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

285 F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

286 local_object_forward=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

287 local_object_reverse=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

288 ## Dict to list for the local object

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

289 for i in range(1, self.size+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

290 local_object_forward.append(len(self.readDict[i]))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

291 local_object_reverse.append(len(self.readDict[-i]))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

292 ## start compiling results by slides

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

293 results=[]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

294 for coordinate in range(self.size - window_size):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

295 local_forward=local_object_forward[coordinate:coordinate + window_size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

296 local_reverse=local_object_reverse[coordinate:coordinate + window_size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

297 if sum(local_forward) == 0 or sum(local_reverse) == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

298 continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

299 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

300 reference_to_local_cor_forward = stats.spearmanr(local_forward, reference_forward)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

301 reference_to_local_cor_reverse = stats.spearmanr(local_reverse, reference_reverse)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

302 if (reference_to_local_cor_forward[0] > 0.2 or reference_to_local_cor_reverse[0]>0.2):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

303 results.append([coordinate+1, reference_to_local_cor_forward[0], reference_to_local_cor_reverse[0]])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

304 except:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

305 pass

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

306 return results

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

307

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

308 def readcount (self, size_inf=0, size_sup=1000, upstream_coord=None, downstream_coord=None, polarity="both", method="conservative"):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

309 '''refactored 24-12-2013 to save memory and introduce offset filtering

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

310 take a look at the defaut parameters that cannot be defined relatively to the instance are they are defined before instanciation

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

311 the trick is to pass None and then test

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

312 polarity parameter can take "both", "forward" or "reverse" as value'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

313 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

314 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

315 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "both":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

316 return self.matchedreadsUp + self.matchedreadsDown

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

317 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "forward":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

318 return self.matchedreadsUp

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

319 if upstream_coord == 1 and downstream_coord == self.windowoffset+self.size-1 and polarity == "reverse":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

320 return self.matchedreadsDown

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

321 n=0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

322 if polarity == "both":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

323 for offset in xrange(upstream_coord, downstream_coord+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

324 if self.readDict.has_key(offset):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

325 for read in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

326 if (read>=size_inf and read<= size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

327 n += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

328 if method != "conservative":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

329 del self.readDict[offset] ## Carefull ! precludes re-use on the self.readDict dictionary !!!!!! TEST

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

330 if self.readDict.has_key(-offset):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

331 for read in self.readDict[-offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

332 if (read>=size_inf and read<= size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

333 n += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

334 if method != "conservative":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

335 del self.readDict[-offset]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

336 return n

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

337 elif polarity == "forward":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

338 for offset in xrange(upstream_coord, downstream_coord+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

339 if self.readDict.has_key(offset):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

340 for read in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

341 if (read>=size_inf and read<= size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

342 n += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

343 return n

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

344 elif polarity == "reverse":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

345 for offset in xrange(upstream_coord, downstream_coord+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

346 if self.readDict.has_key(-offset):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

347 for read in self.readDict[-offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

348 if (read>=size_inf and read<= size_sup):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

349 n += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

350 return n

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

351

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

352 def readsizes (self):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

353 '''return a dictionary of number of reads by size (the keys)'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

354 dicsize = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

355 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

356 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

357 dicsize[size] = dicsize.get(size, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

358 for offset in range (min(dicsize.keys()), max(dicsize.keys())+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

359 dicsize[size] = dicsize.get(size, 0) # to fill offsets with null values

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

360 return dicsize

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

361

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

362 # def size_histogram(self):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

363 # norm=self.norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

364 # hist_dict={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

365 # hist_dict['F']={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

366 # hist_dict['R']={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

367 # for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

368 # for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

369 # if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

370 # hist_dict['R'][size] = hist_dict['R'].get(size, 0) - 1*norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

371 # else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

372 # hist_dict['F'][size] = hist_dict['F'].get(size, 0) + 1*norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

373 # ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

374 # if not (hist_dict['F']) and (not hist_dict['R']):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

375 # hist_dict['F'][21] = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

376 # hist_dict['R'][21] = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

377 # ##

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

378 # return hist_dict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

379

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

380 def size_histogram(self, minquery=None, maxquery=None): # in SmRNAwindow

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

381 '''refactored on 7-9-2014 to debug size_histogram tool'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

382 norm=self.norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

383 size_dict={}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

384 size_dict['F']= defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

385 size_dict['R']= defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

386 size_dict['both']= defaultdict (float)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

387 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

388 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

389 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

390 size_dict['R'][size] = size_dict['R'][size] - 1*norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

391 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

392 size_dict['F'][size] = size_dict['F'][size] + 1*norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

393 ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

394 if not (size_dict['F']) and (not size_dict['R']):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

395 size_dict['F'][21] = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

396 size_dict['R'][21] = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

397 ##

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

398 allSizeKeys = list (set (size_dict['F'].keys() + size_dict['R'].keys() ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

399 for size in allSizeKeys:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

400 size_dict['both'][size] = size_dict['F'][size] - size_dict['R'][size]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

401 if minquery:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

402 for polarity in size_dict.keys():

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

403 for size in xrange(minquery, maxquery+1):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

404 if not size in size_dict[polarity].keys():

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

405 size_dict[polarity][size]=0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

406 return size_dict

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

407

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

408 def statsizes (self, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

409 ''' migration to memory saving by specifying possible subcoordinates

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

410 see the readcount method for further discussion'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

411 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

412 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

413 L = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

414 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

415 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

416 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

417 L.append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

418 meansize = mean(L)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

419 stdv = std(L)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

420 mediansize = median(L)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

421 return meansize, mediansize, stdv

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

422

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

423 def foldEnergy (self, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

424 ''' migration to memory saving by specifying possible subcoordinates

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

425 see the readcount method for further discussion'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

426 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

427 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

428 Energy = RNAfold ([self.sequence[upstream_coord-1:downstream_coord] ])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

429 return float(Energy[self.sequence[upstream_coord-1:downstream_coord]])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

430

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

431 def Ufreq (self, size_scope, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

432 ''' migration to memory saving by specifying possible subcoordinates

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

433 see the readcount method for further discussion. size_scope must be an interable'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

434 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

435 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

436 freqDic = {"A":0,"T":0,"G":0,"C":0, "N":0}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

437 convertDic = {"A":"T","T":"A","G":"C","C":"G","N":"N"}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

438 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

439 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

440 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

441 if size in size_scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

442 startbase = self.sequence[abs(offset)-self.windowoffset]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

443 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

444 startbase = convertDic[startbase]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

445 freqDic[startbase] += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

446 base_sum = float ( sum( freqDic.values()) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

447 if base_sum == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

448 return "."

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

449 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

450 return freqDic["T"] / base_sum * 100

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

451

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

452 def Ufreq_stranded (self, size_scope, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

453 ''' migration to memory saving by specifying possible subcoordinates

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

454 see the readcount method for further discussion. size_scope must be an interable

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

455 This method is similar to the Ufreq method but take strandness into account'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

456 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

457 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

458 freqDic = {"Afor":0,"Tfor":0,"Gfor":0,"Cfor":0, "Nfor":0,"Arev":0,"Trev":0,"Grev":0,"Crev":0, "Nrev":0}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

459 convertDic = {"A":"T","T":"A","G":"C","C":"G","N":"N"}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

460 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

461 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

462 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

463 if size in size_scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

464 startbase = self.sequence[abs(offset)-self.windowoffset]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

465 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

466 startbase = convertDic[startbase]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

467 freqDic[startbase+"rev"] += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

468 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

469 freqDic[startbase+"for"] += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

470 forward_sum = float ( freqDic["Afor"]+freqDic["Tfor"]+freqDic["Gfor"]+freqDic["Cfor"]+freqDic["Nfor"])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

471 reverse_sum = float ( freqDic["Arev"]+freqDic["Trev"]+freqDic["Grev"]+freqDic["Crev"]+freqDic["Nrev"])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

472 if forward_sum == 0 and reverse_sum == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

473 return ". | ."

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

474 elif reverse_sum == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

475 return "%s | ." % (freqDic["Tfor"] / forward_sum * 100)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

476 elif forward_sum == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

477 return ". | %s" % (freqDic["Trev"] / reverse_sum * 100)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

478 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

479 return "%s | %s" % (freqDic["Tfor"] / forward_sum * 100, freqDic["Trev"] / reverse_sum * 100)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

480

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

481

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

482 def readplot (self):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

483 norm=self.norm

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

484 readmap = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

485 for offset in self.readDict.keys():

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

486 readmap[abs(offset)] = ( len(self.readDict.get(-abs(offset),[]))*norm , len(self.readDict.get(abs(offset),[]))*norm )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

487 mylist = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

488 for offset in sorted(readmap):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

489 if readmap[offset][1] != 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

490 mylist.append("%s\t%s\t%s\t%s" % (self.gene, offset, readmap[offset][1], "F") )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

491 if readmap[offset][0] != 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

492 mylist.append("%s\t%s\t%s\t%s" % (self.gene, offset, -readmap[offset][0], "R") )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

493 ## patch to avoid missing graphs when parsed by R-lattice. 27-08-2014. Test and validate !

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

494 if not mylist:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

495 mylist.append("%s\t%s\t%s\t%s" % (self.gene, 1, 0, "F") )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

496 ###

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

497 return mylist

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

498

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

499 def readcoverage (self, upstream_coord=None, downstream_coord=None, windowName=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

500 '''Use by MirParser tool'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

501 upstream_coord = upstream_coord or 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

502 downstream_coord = downstream_coord or self.size

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

503 windowName = windowName or "%s_%s_%s" % (self.gene, upstream_coord, downstream_coord)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

504 forORrev_coverage = dict ([(i,0) for i in xrange(1, downstream_coord-upstream_coord+1)])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

505 totalforward = self.readcount(upstream_coord=upstream_coord, downstream_coord=downstream_coord, polarity="forward")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

506 totalreverse = self.readcount(upstream_coord=upstream_coord, downstream_coord=downstream_coord, polarity="reverse")

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

507 if totalforward > totalreverse:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

508 majorcoverage = "forward"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

509 for offset in self.readDict.keys():

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

510 if (offset > 0) and ((offset-upstream_coord+1) in forORrev_coverage.keys() ):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

511 for read in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

512 for i in xrange(read):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

513 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

514 forORrev_coverage[offset-upstream_coord+1+i] += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

515 except KeyError:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

516 continue # a sense read may span over the downstream limit

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

517 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

518 majorcoverage = "reverse"

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

519 for offset in self.readDict.keys():

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

520 if (offset < 0) and (-offset-upstream_coord+1 in forORrev_coverage.keys() ):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

521 for read in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

522 for i in xrange(read):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

523 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

524 forORrev_coverage[-offset-upstream_coord-i] += 1 ## positive coordinates in the instance, with + for forward coverage and - for reverse coverage

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

525 except KeyError:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

526 continue # an antisense read may span over the upstream limit

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

527 output_list = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

528 maximum = max (forORrev_coverage.values()) or 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

529 for n in sorted (forORrev_coverage):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

530 output_list.append("%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.biosample, windowName, n, float(n)/(downstream_coord-upstream_coord+1), forORrev_coverage[n], float(forORrev_coverage[n])/maximum, majorcoverage))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

531 return "\n".join(output_list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

532

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

533

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

534 def signature (self, minquery, maxquery, mintarget, maxtarget, scope, zscore="no", upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

535 ''' migration to memory saving by specifying possible subcoordinates

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

536 see the readcount method for further discussion

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

537 scope must be a python iterable; scope define the *relative* offset range to be computed'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

538 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

539 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

540 query_range = range (minquery, maxquery+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

541 target_range = range (mintarget, maxtarget+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

542 Query_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

543 Target_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

544 frequency_table = dict ([(i, 0) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

545 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

546 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

547 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

548 if size in query_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

549 Query_table[offset] = Query_table.get(offset, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

550 if size in target_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

551 Target_table[offset] = Target_table.get(offset, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

552 for offset in Query_table:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

553 for i in scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

554 frequency_table[i] += min(Query_table[offset], Target_table.get(-offset -i +1, 0))

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

555 if minquery==mintarget and maxquery==maxtarget: ## added to incorporate the division by 2 in the method (26/11/2013), see signature_options.py and lattice_signature.py

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

556 frequency_table = dict([(i,frequency_table[i]/2) for i in frequency_table])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

557 if zscore == "yes":

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

558 z_mean = mean(frequency_table.values() )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

559 z_std = std(frequency_table.values() )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

560 if z_std == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

561 frequency_table = dict([(i,0) for i in frequency_table] )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

562 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

563 frequency_table = dict([(i, (frequency_table[i]- z_mean)/z_std) for i in frequency_table] )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

564 return frequency_table

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

565

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

566 def hannon_signature (self, minquery, maxquery, mintarget, maxtarget, scope, upstream_coord=None, downstream_coord=None):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

567 ''' migration to memory saving by specifying possible subcoordinates see the readcount method for further discussion

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

568 note that scope must be an iterable (a list or a tuple), which specifies the relative offsets that will be computed'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

569 upstream_coord = upstream_coord or self.windowoffset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

570 downstream_coord = downstream_coord or self.windowoffset+self.size-1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

571 query_range = range (minquery, maxquery+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

572 target_range = range (mintarget, maxtarget+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

573 Query_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

574 Target_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

575 Total_Query_Numb = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

576 general_frequency_table = dict ([(i,0) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

577 ## filtering the appropriate reads for the study

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

578 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

579 if (abs(offset) < upstream_coord or abs(offset) > downstream_coord): continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

580 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

581 if size in query_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

582 Query_table[offset] = Query_table.get(offset, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

583 Total_Query_Numb += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

584 if size in target_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

585 Target_table[offset] = Target_table.get(offset, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

586 for offset in Query_table:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

587 frequency_table = dict ([(i,0) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

588 number_of_targets = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

589 for i in scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

590 frequency_table[i] += Query_table[offset] * Target_table.get(-offset -i +1, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

591 number_of_targets += Target_table.get(-offset -i +1, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

592 for i in scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

593 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

594 general_frequency_table[i] += (1. / number_of_targets / Total_Query_Numb) * frequency_table[i]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

595 except ZeroDivisionError :

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

596 continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

597 return general_frequency_table

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

598

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

599 def phasing (self, size_range, scope):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

600 ''' to calculate autocorelation like signal - scope must be an python iterable'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

601 read_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

602 total_read_number = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

603 general_frequency_table = dict ([(i, 0) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

604 ## read input filtering

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

605 for offset in self.readDict:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

606 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

607 if size in size_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

608 read_table[offset] = read_table.get(offset, 0) + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

609 total_read_number += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

610 ## per offset read phasing computing

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

611 for offset in read_table:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

612 frequency_table = dict ([(i, 0) for i in scope]) # local frequency table

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

613 number_of_targets = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

614 for i in scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

615 if offset > 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

616 frequency_table[i] += read_table[offset] * read_table.get(offset + i, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

617 number_of_targets += read_table.get(offset + i, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

618 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

619 frequency_table[i] += read_table[offset] * read_table.get(offset - i, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

620 number_of_targets += read_table.get(offset - i, 0)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

621 ## inclusion of local frequency table in the general frequency table (all offsets average)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

622 for i in scope:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

623 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

624 general_frequency_table[i] += (1. / number_of_targets / total_read_number) * frequency_table[i]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

625 except ZeroDivisionError :

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

626 continue

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

627 return general_frequency_table

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

628

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

629

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

630

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

631 def z_signature (self, minquery, maxquery, mintarget, maxtarget, scope):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

632 '''Must do: from numpy import mean, std, to use this method; scope must be a python iterable and defines the relative offsets to compute'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

633 frequency_table = self.signature (minquery, maxquery, mintarget, maxtarget, scope)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

634 z_table = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

635 frequency_list = [frequency_table[i] for i in sorted (frequency_table)]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

636 if std(frequency_list):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

637 meanlist = mean(frequency_list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

638 stdlist = std(frequency_list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

639 z_list = [(i-meanlist)/stdlist for i in frequency_list]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

640 return dict (zip (sorted(frequency_table), z_list) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

641 else:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

642 return dict (zip (sorted(frequency_table), [0 for i in frequency_table]) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

643

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

644 def percent_signature (self, minquery, maxquery, mintarget, maxtarget, scope):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

645 frequency_table = self.signature (minquery, maxquery, mintarget, maxtarget, scope)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

646 total = float(sum ([self.readsizes().get(i,0) for i in set(range(minquery,maxquery)+range(mintarget,maxtarget))]) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

647 if total == 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

648 return dict( [(i,0) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

649 return dict( [(i, frequency_table[i]/total*100) for i in scope])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

650

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

651 def pairer (self, overlap, minquery, maxquery, mintarget, maxtarget):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

652 queryhash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

653 targethash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

654 query_range = range (int(minquery), int(maxquery)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

655 target_range = range (int(mintarget), int(maxtarget)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

656 paired_sequences = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

657 for offset in self.readDict: # selection of data

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

658 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

659 if size in query_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

660 queryhash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

661 if size in target_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

662 targethash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

663 for offset in queryhash:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

664 if offset >= 0: matched_offset = -offset - overlap + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

665 else: matched_offset = -offset - overlap + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

666 if targethash[matched_offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

667 paired = min ( len(queryhash[offset]), len(targethash[matched_offset]) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

668 if offset >= 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

669 for i in range (paired):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

670 paired_sequences.append("+%s" % RNAtranslate ( self.sequence[offset:offset+queryhash[offset][i]]) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

671 paired_sequences.append("-%s" % RNAtranslate (antipara (self.sequence[-matched_offset-targethash[matched_offset][i]+1:-matched_offset+1]) ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

672 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

673 for i in range (paired):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

674 paired_sequences.append("-%s" % RNAtranslate (antipara (self.sequence[-offset-queryhash[offset][i]+1:-offset+1]) ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

675 paired_sequences.append("+%s" % RNAtranslate (self.sequence[matched_offset:matched_offset+targethash[matched_offset][i]] ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

676 return paired_sequences

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

677

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

678 def pairable (self, overlap, minquery, maxquery, mintarget, maxtarget):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

679 queryhash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

680 targethash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

681 query_range = range (int(minquery), int(maxquery)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

682 target_range = range (int(mintarget), int(maxtarget)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

683 paired_sequences = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

684

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

685 for offset in self.readDict: # selection of data

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

686 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

687 if size in query_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

688 queryhash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

689 if size in target_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

690 targethash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

691

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

692 for offset in queryhash:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

693 matched_offset = -offset - overlap + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

694 if targethash[matched_offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

695 if offset >= 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

696 for i in queryhash[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

697 paired_sequences.append("+%s" % RNAtranslate (self.sequence[offset:offset+i]) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

698 for i in targethash[matched_offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

699 paired_sequences.append( "-%s" % RNAtranslate (antipara (self.sequence[-matched_offset-i+1:-matched_offset+1]) ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

700 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

701 for i in queryhash[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

702 paired_sequences.append("-%s" % RNAtranslate (antipara (self.sequence[-offset-i+1:-offset+1]) ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

703 for i in targethash[matched_offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

704 paired_sequences.append("+%s" % RNAtranslate (self.sequence[matched_offset:matched_offset+i] ) )

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

705 return paired_sequences

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

706

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

707 def newpairable_bowtie (self, overlap, minquery, maxquery, mintarget, maxtarget):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

708 ''' revision of pairable on 3-12-2012, with focus on the offset shift problem (bowtie is 1-based cooordinates whereas python strings are 0-based coordinates'''

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

709 queryhash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

710 targethash = defaultdict(list)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

711 query_range = range (int(minquery), int(maxquery)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

712 target_range = range (int(mintarget), int(maxtarget)+1)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

713 bowtie_output = []

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

714

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

715 for offset in self.readDict: # selection of data

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

716 for size in self.readDict[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

717 if size in query_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

718 queryhash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

719 if size in target_range:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

720 targethash[offset].append(size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

721 counter = 0

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

722 for offset in queryhash:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

723 matched_offset = -offset - overlap + 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

724 if targethash[matched_offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

725 if offset >= 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

726 for i in queryhash[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

727 counter += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

728 bowtie_output.append("%s\t%s\t%s\t%s\t%s" % (counter, "+", self.gene, offset-1, self.sequence[offset-1:offset-1+i]) ) # attention a la base 1-0 de l'offset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

729 if offset < 0:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

730 for i in queryhash[offset]:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

731 counter += 1

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

732 bowtie_output.append("%s\t%s\t%s\t%s\t%s" % (counter, "-", self.gene, -offset-i, self.sequence[-offset-i:-offset])) # attention a la base 1-0 de l'offset

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

733 return bowtie_output

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

734

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

735

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

736 def __main__(bowtie_index_path, bowtie_output_path):

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

737 sequenceDic = get_fasta (bowtie_index_path)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

738 objDic = {}

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

739 F = open (bowtie_output_path, "r") # F is the bowtie output taken as input

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

740 for line in F:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

741 fields = line.split()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

742 polarity = fields[1]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

743 gene = fields[2]

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

744 offset = int(fields[3])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

745 size = len (fields[4])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

746 try:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

747 objDic[gene].addread (polarity, offset, size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

748 except KeyError:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

749 objDic[gene] = SmRNAwindow(gene, sequenceDic[gene])

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

750 objDic[gene].addread (polarity, offset, size)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

751 F.close()

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

752 for gene in objDic:

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

753 print gene, objDic[gene].pairer(19,19,23,19,23)

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

754

b996480cd604 planemo upload for repository https://bitbucket.org/drosofff/gedtools/

drosofff

parents:

diff changeset

755 if __name__ == "__main__" : __main__(sys.argv[1], sys.argv[2])

Mercurial > repos > drosofff > msp_sr_bowtie_parser

annotate smRtools.py @ 0:b996480cd604 draft