<commit_before>#!/usr/bin/python
# version 1 7-5-2012 unification of the SmRNAwindow class

import sys, subprocess
from collections import defaultdict
from numpy import mean, median, std
from scipy import stats

def get_fasta (index="/home/galaxy/galaxy-dist/bowtie/5.37_Dmel/5.37_Dmel"):
  '''This function will return a dictionary containing fasta identifiers as keys and the
  sequence as values. Index must be the path to a fasta file.'''
  p = subprocess.Popen(args=["bowtie-inspect","-a", "0", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
  outputlines = p.stdout.readlines()
  p.wait()
  item_dic = {}
  for line in outputlines:
    if (line[0] == ">"):
      try:
        item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
      except: pass
      current_item = line[1:].rstrip().split()[0] #take the first word before space because bowtie splits headers !
      item_dic[current_item] = ""
      stringlist=[]
    else:
      stringlist.append(line.rstrip() )
  item_dic[current_item] = "".join(stringlist) # for the last item
  return item_dic

def get_fasta_headers (index):
  p = subprocess.Popen(args=["bowtie-inspect","-n", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
  outputlines = p.stdout.readlines()
  p.wait()
  item_dic = {}
  for line in outputlines:
    header = line.rstrip().split()[0] #take the first word before space because bowtie splits headers !
    item_dic[header] = 1
  return item_dic


def get_file_sample (file, numberoflines):
  '''import random to use this function'''
  F=open(file)
  fullfile = F.read().splitlines()
  F.close()
  if len(fullfile) < numberoflines:
    return "sample size exceeds file size"
  return random.sample(fullfile, numberoflines)

def get_fasta_from_history (file):
  F = open (file, "r")
  item_dic = {}
  for line in F:
    if (line[0] == ">"):
      try:
        item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
      except: pass
      current_item = line[1:-1].split()[0] #take the first word before space because bowtie splits headers !
      item_dic[current_item] = ""
      stringlist=[]
    else:
      stringlist.append(line[:-1])
  item_dic[current_item] = "".join(stringlist) # for the last item
  return item_dic

def antipara (sequence):
    antidict = {"A":"T", "T":"A", "G":"C", "C":"G", "N":"N"}
    revseq = sequence[::-1]
    return "".join([antidict[i] for i in revseq])

def RNAtranslate (sequence):
    return "".join([i if i in "AGCN" else "U" for i in sequence])
def DNAtranslate (sequence):
    return "".join([i if i in "AGCN" else "T" for i in sequence])

def RNAfold (sequence_list):
  thestring= "\n".join(sequence_list)
  p = subprocess.Popen(args=["RNAfold","--noPS"], stdin= subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
  output=p.communicate(thestring)[0]
  p.wait()
  output=output.split("\n")
  if not output[-1]: output = output[:-1] # nasty patch to remove last empty line
  buffer=[]
  for line in output:
    if line[0] in ["N","A","T","U","G","C"]:
      buffer.append(DNAtranslate(line))
    if line[0] in ["(",".",")"]:
      fields=line.split("(")
      energy= fields[-1]
      energy = energy[:-1] # remove the ) parenthesis
      energy=float(energy)
      buffer.append(str(energy))
  return dict(zip(buffer[::2], buffer[1::2]))

def extractsubinstance (start, end, instance):
  ''' Testing whether this can be an function external to the class to save memory'''
  subinstance = SmRNAwindow (instance.gene, instance.sequence[start-1:end], start)
  subinstance.gene = "%s %s %s" % (subinstance.gene, subinstance.windowoffset, subinstance.windowoffset + subinstance.size - 1)
  upcoordinate = [i for i in range(start,end+1) if instance.readDict[i] ]
  downcoordinate = [-i for i in range(start,end+1) if instance.readDict[-i] ]
  for i in upcoordinate:
    subinstance.readDict[i]=instance.readDict[i]
  for i in downcoordinate:
    subinstance.readDict[i]=instance.readDict[i]
  return subinstance

class HandleSmRNAwindows:
  def __init__(self, alignmentFile="~", alignmentFileFormat="tabular", genomeRefFile="~", genomeRefFormat="bowtieIndex", biosample="undetermined"):
    self.biosample = biosample
    self.alignmentFile = alignmentFile
    self.alignmentFileFormat = alignmentFileFormat # can be "tabular" or "sam"
    self.genomeRefFile = genomeRefFile
    self.genomeRefFormat = genomeRefFormat # can be "bowtieIndex" or "fastaSource"
    self.instanceDict = {}
    if genomeRefFormat == "bowtieIndex":
      self.itemDict = get_fasta (genomeRefFile)
    elif genomeRefFormat == "fastaSource":
      self.itemDict = get_fasta_from_history (genomeRefFile)
    for item in self.itemDict:
      self.instanceDict[item] = SmRNAwindow(item, sequence=self.itemDict[item], windowoffset=1, biosample=self.biosample) # create as many instances as there is items
    self.readfile()

  def readfile (self) :
    if self.alignmentFileFormat == "tabular":
      F = open (self.alignmentFile, "r")
      for line in F:
        fields = line.split()
        polarity = fields[1]
        gene = fields[2]
        offset = int(fields[3])
        size = len (fields[4])
        self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
      F.close()
      return self.instanceDict
    elif self.alignmentFileFormat == "sam":
      F = open (self.alignmentFile, "r")
      dict = {"0":"+", "16":"-"}
      for line in F:
        if line[0]=='@':
            continue
        fields = line.split()
        if fields[2] == "*": continue
        polarity = dict[fields[1]]
        gene = fields[2]
        offset = int(fields[3])
        size = len (fields[9])
        self.instanceDict[gene].addread (polarity, offset, size) # sam format is already 1-based coordinates
      F.close()
    elif self.alignmentFileFormat == "bam":
      import pysam
      samfile = pysam.Samfile(self.alignmentFile)
      for read in samfile:
        if read.tid == -1:
          continue # filter out unaligned reads
        if read.is_reverse:
          polarity="-"
        else:
          polarity="+"
        gene = samfile.getrname(read.tid)
        offset = read.pos
        size = read.qlen
        self.instanceDict[gene].addread (polarity, offset+1, size) # pysam converts coordinates to 0-based (https://media.readthedocs.org/pdf/pysam/latest/pysam.pdf)
      return self.instanceDict

class SmRNAwindow:

  def __init__(self, gene, sequence="ATGC", windowoffset=1, biosample="Undetermined"):
    self.biosample = biosample
    self.sequence = sequence
    self.gene = gene
    self.windowoffset = windowoffset
    self.size = len(sequence)
    self.readDict = defaultdict(list) # with a {+/-offset:[size1, size2, ...], ...}
    
  def addread (self, polarity, offset, size):
    '''ATTENTION ATTENTION ATTENTION'''
    ''' We removed the conversion from 0 to 1 based offset, as we do this now during readparsing.'''
    if polarity == "+":
      self.readDict[offset].append(size)
    else:
      self.readDict[-(offset + size -1)].append(size)
    return

  def barycenter (self, upstream_coord=None, downstream_coord=None):
    '''refactored 24-12-2013 to save memory and introduce offset filtering see readcount method for further discussion on that
    In this version, attempt to replace the dictionary structure by a list of tupple to save memory too'''
    upstream_coord = upstream_coord or self.windowoffset
    downstream_coord = downstream_coord or self.windowoffset+self.size-1
    window_size = downstream_coord - upstream_coord +1
    def weigthAverage (TuppleList):
      weightSum = 0
      PonderWeightSum = 0
      for tuple in TuppleList:<commit_msg>Added the CountFeatures method to the HandleSmRNAwindows class<commit_after>#!/usr/bin/python
# version 1 7-5-2012 unification of the SmRNAwindow class

import sys, subprocess
from collections import defaultdict
from numpy import mean, median, std
from scipy import stats

def get_fasta (index="/home/galaxy/galaxy-dist/bowtie/5.37_Dmel/5.37_Dmel"):
  '''This function will return a dictionary containing fasta identifiers as keys and the
  sequence as values. Index must be the path to a fasta file.'''
  p = subprocess.Popen(args=["bowtie-inspect","-a", "0", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
  outputlines = p.stdout.readlines()
  p.wait()
  item_dic = {}
  for line in outputlines:
    if (line[0] == ">"):
      try:
        item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
      except: pass
      current_item = line[1:].rstrip().split()[0] #take the first word before space because bowtie splits headers !
      item_dic[current_item] = ""
      stringlist=[]
    else:
      stringlist.append(line.rstrip() )
  item_dic[current_item] = "".join(stringlist) # for the last item
  return item_dic

def get_fasta_headers (index):
  p = subprocess.Popen(args=["bowtie-inspect","-n", index], stdout=subprocess.PIPE, stderr=subprocess.STDOUT) # bowtie-inspect outputs sequences on single lines
  outputlines = p.stdout.readlines()
  p.wait()
  item_dic = {}
  for line in outputlines:
    header = line.rstrip().split()[0] #take the first word before space because bowtie splits headers !
    item_dic[header] = 1
  return item_dic


def get_file_sample (file, numberoflines):
  '''import random to use this function'''
  F=open(file)
  fullfile = F.read().splitlines()
  F.close()
  if len(fullfile) < numberoflines:
    return "sample size exceeds file size"
  return random.sample(fullfile, numberoflines)

def get_fasta_from_history (file):
  F = open (file, "r")
  item_dic = {}
  for line in F:
    if (line[0] == ">"):
      try:
        item_dic[current_item] = "".join(stringlist) # to dump the sequence of the previous item - try because of the keyerror of the first item
      except: pass
      current_item = line[1:-1].split()[0] #take the first word before space because bowtie splits headers !
      item_dic[current_item] = ""
      stringlist=[]
    else:
      stringlist.append(line[:-1])
  item_dic[current_item] = "".join(stringlist) # for the last item
  return item_dic

def antipara (sequence):
    antidict = {"A":"T", "T":"A", "G":"C", "C":"G", "N":"N"}
    revseq = sequence[::-1]
    return "".join([antidict[i] for i in revseq])

def RNAtranslate (sequence):
    return "".join([i if i in "AGCN" else "U" for i in sequence])
def DNAtranslate (sequence):
    return "".join([i if i in "AGCN" else "T" for i in sequence])

def RNAfold (sequence_list):
  thestring= "\n".join(sequence_list)
  p = subprocess.Popen(args=["RNAfold","--noPS"], stdin= subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
  output=p.communicate(thestring)[0]
  p.wait()
  output=output.split("\n")
  if not output[-1]: output = output[:-1] # nasty patch to remove last empty line
  buffer=[]
  for line in output:
    if line[0] in ["N","A","T","U","G","C"]:
      buffer.append(DNAtranslate(line))
    if line[0] in ["(",".",")"]:
      fields=line.split("(")
      energy= fields[-1]
      energy = energy[:-1] # remove the ) parenthesis
      energy=float(energy)
      buffer.append(str(energy))
  return dict(zip(buffer[::2], buffer[1::2]))

def extractsubinstance (start, end, instance):
  ''' Testing whether this can be an function external to the class to save memory'''
  subinstance = SmRNAwindow (instance.gene, instance.sequence[start-1:end], start)
  subinstance.gene = "%s %s %s" % (subinstance.gene, subinstance.windowoffset, subinstance.windowoffset + subinstance.size - 1)
  upcoordinate = [i for i in range(start,end+1) if instance.readDict[i] ]
  downcoordinate = [-i for i in range(start,end+1) if instance.readDict[-i] ]
  for i in upcoordinate:
    subinstance.readDict[i]=instance.readDict[i]
  for i in downcoordinate:
    subinstance.readDict[i]=instance.readDict[i]
  return subinstance

class HandleSmRNAwindows:
  def __init__(self, alignmentFile="~", alignmentFileFormat="tabular", genomeRefFile="~", genomeRefFormat="bowtieIndex", biosample="undetermined"):
    self.biosample = biosample
    self.alignmentFile = alignmentFile
    self.alignmentFileFormat = alignmentFileFormat # can be "tabular" or "sam"
    self.genomeRefFile = genomeRefFile
    self.genomeRefFormat = genomeRefFormat # can be "bowtieIndex" or "fastaSource"
    self.instanceDict = {}
    if genomeRefFormat == "bowtieIndex":
      self.itemDict = get_fasta (genomeRefFile)
    elif genomeRefFormat == "fastaSource":
      self.itemDict = get_fasta_from_history (genomeRefFile)
    for item in self.itemDict:
      self.instanceDict[item] = SmRNAwindow(item, sequence=self.itemDict[item], windowoffset=1, biosample=self.biosample) # create as many instances as there is items
    self.readfile()

  def readfile (self) :
    if self.alignmentFileFormat == "tabular":
      F = open (self.alignmentFile, "r")
      for line in F:
        fields = line.split()
        polarity = fields[1]
        gene = fields[2]
        offset = int(fields[3])
        size = len (fields[4])
        self.instanceDict[gene].addread (polarity, offset+1, size) # to correct to 1-based coordinates of SmRNAwindow
      F.close()
      return self.instanceDict
    elif self.alignmentFileFormat == "sam":
      F = open (self.alignmentFile, "r")
      dict = {"0":"+", "16":"-"}
      for line in F:
        if line[0]=='@':
            continue
        fields = line.split()
        if fields[2] == "*": continue
        polarity = dict[fields[1]]
        gene = fields[2]
        offset = int(fields[3])
        size = len (fields[9])
        self.instanceDict[gene].addread (polarity, offset, size) # sam format is already 1-based coordinates
      F.close()
    elif self.alignmentFileFormat == "bam":
      import pysam
      samfile = pysam.Samfile(self.alignmentFile)
      for read in samfile:
        if read.tid == -1:
          continue # filter out unaligned reads
        if read.is_reverse:
          polarity="-"
        else:
          polarity="+"
        gene = samfile.getrname(read.tid)
        offset = read.pos
        size = read.qlen
        self.instanceDict[gene].addread (polarity, offset+1, size) # pysam converts coordinates to 0-based (https://media.readthedocs.org/pdf/pysam/latest/pysam.pdf)
      return self.instanceDict

  def CountFeatures (self, GFF3="path/to/file"):
    featureDict = defaultdict(int)
    F  = open (GFF3, "r")
    for line in F:
      if line[0] ==  "#": continue
      fields = line[:-1].split()
      chrom, feature, leftcoord, rightcoord, polarity = fields[0], fields[2], fields[3], fields[4], fields[6]
      featureDict[feature] += self.instanceDict[chrom].readcount(upstream_coord=leftcoord, downstream_coord=rightcoord, polarity="both")
    F.close()
    return featureDict

class SmRNAwindow:

  def __init__(self, gene, sequence="ATGC", windowoffset=1, biosample="Undetermined"):
    self.biosample = biosample
    self.sequence = sequence
    self.gene = gene
    self.windowoffset = windowoffset
    self.size = len(sequence)
    self.readDict = defaultdict(list) # with a {+/-offset:[size1, size2, ...], ...}
    
  def addread (self, polarity, offset, size):
    '''ATTENTION ATTENTION ATTENTION'''
    ''' We removed the conversion from 0 to 1 based offset, as we do this now during readparsing.'''
    if polarity == "+":
      self.readDict[offset].append(size)
    else:
      self.readDict[-(offset + size -1)].append(size)
    return

  def barycenter (self, upstream_coord=None, downstream_coord=None):
    '''refactored 24-12-2013 to save memory and introduce offset filtering see readcount method for further discussion on that
    In this version, attempt to replace the dictionary structure by a list of tupple to save memory too'''
    upstream_coord = upstream_coord or self.windowoffset
    downstream_coord = downstream_coord or self.windowoffset+self.size-1
    window_size = downstream_coord - upstream_coord +1
    def weigthAverage (TuppleList):
      weightSum = 0
      PonderWeightSum = 0
      for tuple in TuppleList:
{'fraction_non_alphanumeric': 0.1266516470803142, 'fraction_numerical': 0.009925105035620775, 'mean_word_length': 3.5232718259432665, 'pattern_counts': {'":': 29, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2730331', 'n_tokens_mistral': 5372, 'n_tokens_neox': 5143, 'n_words': 1462}