Synthetic TOCSY Macro

File contents

from ccpnmr.analysis.core.Util import getAnalysisSpectrum

def syntheticTocsyMacro(argServer):

  answer = argServer.askYesNo('Use all HSQC experiments? (Otherwise you can select just one)')
  
  analysis = argServer.parent
  project = argServer.getProject()
  nmrProject = project.currentNmrProject
  
  experiments = []
  for experiment in nmrProject.sortedExperiments():
    refExperiment = experiment.refExperiment
 
    if refExperiment and (refExperiment.name == 'H[C]'):
      experiments.append(experiment)
  
  if not answer:
    experiment = argServer.chooseObject(experiments)
    
    if experiment:
      experiments = [experiment,]
  
    else:
      argServer.showWarning('No experiment')
      return
  
  for experiment in experiments:
    peakList = makeTocsyFromHsqc(experiment)
    
    if peakList:
      spectrum = peakList.dataSource
      analysis.initAnalysisSpectrum(spectrum)
      analysis.initSpectrumPeakList(spectrum)
      analysis.initPeakList(peakList, initialising=False)
      analysis.checkCreateSpectrumViews(spectrum)
  
from ccpnmr.analysis.core.ExperimentBasic import setRefExperiment
from ccpnmr.analysis.core.ExperimentBasic import getPrimaryExpDimRef
from ccpnmr.analysis.core.ExperimentBasic import initExpTransfers
from ccpnmr.analysis.core.PeakBasic import makePeakListFromShifts
from ccpnmr.analysis.core.PeakBasic import setManualPeakIntensity

def makeTocsyFromHsqc(hsqc, suffix='HH'):

  nmrProject = hsqc.nmrProject
  project = hsqc.root
  hsqcSpectrum = hsqc.findFirstDataSource()

  shiftList = hsqc.shiftList
  molSystem = hsqc.findFirstMolSystem()
  
  if not molSystem:
    print "WARNING: Cannot use experiment %s - no molecular system set" % hsqc.name
    return
  
  print 'Working on %s' % hsqc.name
  
  # Make experiment
  
  name = hsqc.name + suffix
  
  if nmrProject.findFirstExperiment(name=name):
    name = name + '1'
  
    i = 1
    while nmrProject.findFirstExperiment(name=name):
      i += 1
      name = name[:1] + '%d' % i
  
  tocsy = nmrProject.newExperiment(name=name, numDim=2, shiftList=shiftList)
  
  for expDim in hsqc.expDims:
    expDimRef = getPrimaryExpDimRef(expDim)
    if '1H' in expDimRef.isotopeCodes:
      sf = expDimRef.sf
      dataDimRef = expDimRef.findFirstDataDimRef()
      break
    
  expDims = tocsy.sortedExpDims()
  expDims[0].newExpDimRef(sf=sf, isotopeCodes=("1H",), unit='ppm')
  expDims[1].newExpDimRef(sf=sf, isotopeCodes=("1H",), unit='ppm')
  
  
  nmrExpPrototype = project.findFirstNmrExpPrototype(name='H_H.TOCSY')
  refExperiment = nmrExpPrototype.findFirstRefExperiment(name='H_H.TOCSY')
  setRefExperiment(tocsy, refExperiment)
  initExpTransfers(tocsy)
  
  # Make spectrum
  
  spectrum = tocsy.newDataSource(name="1", numDim=2, dataType="processed")
 
  dataDim = dataDimRef.dataDim
  numPoints = dataDim.numPoints
  numPointsOrig = dataDim.numPoints
  valuePerPoint = dataDim.valuePerPoint
  refPoint = dataDimRef.refPoint
  refValue = dataDimRef.refValue

  for i, expDim in enumerate(tocsy.sortedExpDims()):
    dim = expDim.dim
    freqDataDim = spectrum.newFreqDataDim(dim=dim, expDim=expDim,
                                          isComplex=False,
                                          numPoints=numPoints,
                                          numPointsOrig=numPointsOrig,
                                          valuePerPoint=valuePerPoint)
                                          
    expDimRef = expDim.findFirstExpDimRef()
    freqDataDim.newDataDimRef(expDimRef=expDimRef,
                              refPoint=refPoint,
                              refValue=refValue)
                            
  # Make data store

  hsqcMatrix =  hsqcSpectrum.dataStore
  dataLocationStore = hsqcMatrix.dataLocationStore

  kwArgs = {'path':'uknown',
            'dataUrl':hsqcMatrix.dataUrl,
            'blockSizes':hsqcMatrix.blockSizes,
            'headerSize':hsqcMatrix.headerSize,
            'isBigEndian':hsqcMatrix.isBigEndian,
            'isComplex':hsqcMatrix.isComplex,
            'nByte':hsqcMatrix.nByte,
            'numPoints':[numPoints] * 2,
            'numberType':hsqcMatrix.numberType}
 
  dataStore = dataLocationStore.newBlockedBinaryMatrix(**kwArgs)
  spectrum.dataStore = dataStore
  

  # Make peak list

  peakList = makePeakListFromShifts(spectrum,
                                    useUnassigned=True,
                                    shiftList=shiftList,
                                    molSystem=molSystem,
                                    interResidueTransfer=False,
                                    labellingScheme=None,
                                    labellingThreshold=0.0)
 
  if peakList:
    for peak in peakList.peaks:
      setManualPeakIntensity(peak, 1.0, intensityType='height')
      setManualPeakIntensity(peak, 1.0, intensityType='volume')
 
    print '    Made %d peaks' % len(peakList.peaks)
  else:
    print '    Peak synthesis failed'
  
  return peakList