Getting started with BaMORC

knitr::opts_chunk$set(comment = "#>", collapse = TRUE)

BaMORC quickstart guide

The goal of this document is to get you up and running with BaMORC as quickly as possible. BaMORC stands for the Bayesian Model Optimized Reference Correction (BaMORC) Method for Assigned and Unassigned Protein NMR Spectra. For a detailed explanation of the algorithm, please refer to "Automatic 13C chemical shift reference correction for unassigned protein NMR spectra".

There are two important parts to BaMORC: the bamorc() reference correction function for assigned 13C protein NMR spectra, and the unassigned_bamorc() reference correction function for unassigned 13C protein NMR spectra. In the first section, you'll learn about the basics of running both of these functions. In the second section, you will dive into the basics of input data process, and in third section, you'll learn a little more about the functions used behind-the-scenes by the BaMORC algorithm.

BaMORC basics

To make a correction, first load BaMORC, then call bamorc() or unassigned_bamorc() with the correct arguments:


For assigned protein NMR spectra:

Here we will using the built-in data to demonstrate the arguments that will be passed into bamorc().

## Arguments:
sequence = paste(RefDB_data$carbonDat[[1]]$AA, collapse = "")
secondary_structure = paste(RefDB_data$carbonDat[[1]]$SS, collapse = "")
chemical_shifts_input = RefDB_data$carbonDat[[1]][, c(4,5)]
from= -5
to = 5

## Running bamorc() function:
bamorc(sequence = sequence, secondary_structure = secondary_structure, 
       chemical_shifts_input = chemical_shifts_input, from = -5, to = 5)

The length of the sequence and secondary_structure should be the same, if not please assign secondary_structure = NULL, however, the peaklist groups number could be more or less than the length of the sequence. Often peak lists will have fewer CA/CB pairs than the length of the protein sequence.

Printing an argument object gives you some useful information: the actual format, the size, the object type, and if it's a string.

For unassigned protein NMR spectra:

Next, we will use the built-in data to demonstrate the arguments that will be passed in the unassigned_bamorc(). The output will be slightly different each time it runs, due to the randomness of the optimization.

## Arguments:enerate a temperary sample NMR spectra file and later will be 
## removed.
sample_file_path = system.file("extdata", "bpti_HNcoCACB.txt", package = "BaMORC")

## Running unassigned_bamorc() function:
# unassigned_bamorc will take a while to run
unassigned_bamorc(peakList_file_loc = sample_file_path, sequence = sequence,
                 secondary_structure = NULL, from = -5, to = 5)

The data processing:

The data passed in both above functions should be pre-processed. Luckly, we provide a variety of helper functions. See the following on how to process the data.

The sample generating functions:

There are three file parsing functions within the BaMORC package: read_raw_file(), read_nmrstar_file(), and read_db_file(). These functions handle a wide range of the input files.

The reading functions will return sequence and chemical shifts:

## Download a BMRB file
bmrb_download(4020, output_dir = "./")

## Read in BMRB file and procec
file_path = "bmr4020.str"
head(read_nmrstar_file(file_path = file_path))

## Delete downloaded BMRB file
id = 4022
output <- read_db_file(id = id)

Estimating econdary structure from sequence:

Secondary structure information predicted based on the protein sequence is required for the optimization. This is predicted through the JPred and jpredapi.

secondary_structure <- jpred_fetcher(protein_sequence = protein_sequence)


Behind-the-scenes functions

The BaMORC algorithm minimizes the difference between the actual relative cummulative frequency (RCF) of the protein sequence and the estimated RCF predicted from the chemical shifts information. The statistical functions used by the algorithm are shown below with examples. For detailed function descriptions, please see the reference:

calculate_aa_prob(chi_squared_stat = c(0.05, 0.1, 0.5), df = 2)
calculate_chi_squared_stat(cacb_pair = c(54,45))
## Arguments:
sequence = paste(RefDB_data$carbonDat[[1]]$AA, collapse = "")
secondary_structure = paste(RefDB_data$carbonDat[[1]]$SS, collapse = "")

## Function:
calculate_rcf(sequence = sequence, secondary_structure = secondary_structure)
## chemicalShifts and aaFreq are predefined sample variables for demo purpose 
## within the BaMORC Package.
calculate_mse(step_ca = 1, step_cb = 1, dat_cacb = chemicalShifts[, c(3,4)], 
              aa_Freq = aaFreq)

Try the BaMORC package in your browser

Any scripts or data that you put into this service are public.

BaMORC documentation built on May 1, 2019, 6:35 p.m.