Source Code for Module core.Mem

  1  """
 
  2  Created  2012
 
  3  @author: GieseS
 
  4  
 
  5  
 
  6  Little plotting script which is called in the analysis of different mappings to an artificial reference genome.
 
  7  It produces the following plots:
 
  8  
 
  9  
 
 10  
 
 11  
 
 12  1) ROC Curve
 
 13  2) Overview histograms for FP / TP.
 
 14  
 
 15  """ 
 16  
 
 17  
 
 18 -def ReadReferenceSAMfile(ref, compareList, readdic, entries): 
 19      """
 
 20      Function for reading a complete ReferenceFile into Memory.
 
 21      @type  ref: string
 
 22      @param ref: reference file.ArtRead
 
 23      @type  entries: int
 
 24      @param entries: size for np.array
 
 25      @type  compareList: list
 
 26      @param compareList: list containt 1s, where there is a difference in the genomes and 0s where the nucleotides are equal.
 
 27      @type  readdic: dictionary
 
 28      @param readdic: dictionary containing read ids and read qualities.
 
 29      @rtype:   array
 
 30      @return:  aligned read objects in an array.
 
 31      """ 
 32      # create array with max size = |reads|
 
 33      AlignedReadRefArray = np.zeros(entries, dtype=object) 
 34      tp = 0 
 35      start = time.time() 
 36      start100k = time.time() 
 37      #print ("\tREFERENCE:")
 
 38      fobj = open(ref, "r") 
 39      refdic = {} 
 40      k = 0 
 41      #print ("\tprocessed alignments in Million:")
 
 42      #print("\t\t"),
 
 43      for alignment in fobj: 
 44          k += 1 
 45          if k == 100000: 
 46              end100k = time.time() 
 47              print ("%f.." % (end100k - start100k)), 
 48  
 
 49          read, readname = isSaneAlignment(alignment, "ref", compareList, readdic) 
 50          if read == 0: 
 51              pass 
 52          
 
 53          else: 
 54              # Check if entry in AlignedReadRef == 0, 0 means no read with that ID
 
 55              if AlignedReadRefArray[readdic[readname].internalID] != 0: 
 56                  # read already in dic? Add the non-redundant values to the object
 
 57                  AlignedReadRefArray[readdic[readname].internalID].toObjself(read) 
 58                  #refdic[readname].toObjself(read)
 
 59              else: 
 60                  # new read? Just add it to the dictionary
 
 61                  AlignedReadRefArray[readdic[readname].internalID] = read 
 62                  tp += 1 
 63                  #refdic[readname] = read
 
 64  
 
 65      fobj.close() 
 66      end = time.time() 
 67      #print("\r\n")
 
 68      print ("\t %f " % (end - start)), 
 69      return(AlignedReadRefArray, tp) 
 70  
 
 71  
 
 72 -def ReadArtificialSAMfile(art, compareList, RefArray, readdic): 
 73      """
 
 74      Function for reading a complete ReferenceFile into Memory.
 
 75      @type  art: string
 
 76      @param art: artificial file.
 
 77      @type  RefArray: array
 
 78      @param RefArray: Results from reading the reference SAM file.
 
 79      @type  compareList: list
 
 80      @param compareList: list containt 1s, where there is a difference in the genomes and 0s where the nucleotides are equal.
 
 81      @type  readdic: dictionary
 
 82      @param readdic: dictionary containing read ids and read qualities.
 
 83      @rtype:   dictionary
 
 84      @return:  aligned read objects from the artificial reference, which where uniquely aligned.
 
 85      """ 
 86      
 
 87      start = time.time() 
 88      #print ("\tARTIFICIAL:")
 
 89      fobj = open(art, "r") 
 90      artdic = {} 
 91      k = 0 
 92      
 
 93      read = SkipHeader(fobj,compareList,readdic) 
 94      
 
 95      
 
 96      for alignment in fobj: 
 97          k += 1 
 98          if k % 1000000 == 0: 
 99              print ("%d.." %(k/1000000)), 
100          read, readname = isSaneAlignment(alignment, "art", compareList, readdic) 
101          if read == 0: 
102              pass 
103          else: 
104              # Check by internal naming if spot in array is taken (!= 0) by a read
 
105              index = returnIndex(readdic, readname) 
106              if RefArray[index] != 0: 
107                  # read already in dic? Check if the alignments are the same
 
108                  if (RefArray[index].isContained(read) == 0): 
109                      if readname in artdic: 
110                          artdic[readname].toObjself(read) 
111                      
 
112                      else: 
113                          artdic[readname] = read 
114                  else: 
115                      pass 
116                  
 
117              else: 
118                  # new read? Just add it to the dictionary
 
119                  artdic[readname] = read 
120             
 
121   
 
122      fobj.close() 
123      end = time.time() 
124      #print("\r\n")
 
125      print ("\t %f " % (end - start)), 
126      #print ("\tdone in %d seconds" % (end-start))
 
127      return(artdic) 
128