SDS PAGEの定量をできるだけ自動でしたい（2）

目的

　ポリアクリルアミド電気泳動（SDS-PAGE）はタンパク質を分子量で分離することができ、液中のタンパク質の分布を半定量的に確認できます。生命科学では多用され、私自身も多く活用しています。しかし、その解析は泳動後のゲルを視覚的に観察するか、目的のバンドに対してimage J等で濃度を換算し、相対的に定量するようなことしかしておりませんでした。

　今回は、複数のSDS-PAGEのゲルに対して、汎用的な方法でレーン毎を比較できるような手法を開発することを目的としました。具体的なゴールはまだできていませんが、同じようなことをしている方や誰かの参考になれば幸いです。

本投稿はSDS PAGEの定量をできるだけ自動でしたい（1）の後半となります。

方法

まだsmilingやsmearingに対してはあまり対処しておりませんが、大きく分けると

画像からレーンを切り取る
↓
レーンからクロマトグラムを作成する
↓
バックグラウンドを除去する
↓
タンパク質スタンダードラダーのレーンから分子量の検量線を作成する
↓
各クロマトグラムに多変量正規分布の当てはめを行う
↓
一つの正規分布が一つのバンドとみなし、面積を定量
↓
正規分布の位置を検量線によって分子量に変換
↓
バンドの分子量と濃度が算出できる

となっています。ここからさらに同一高さのバンドを同じものとみなしたりしたかったのですが、まだアイデアが湧いていません。良い意見があれば是非お願いいたします。

本記事では太字になっている部分までの流れを記載します。
前半については前記事で説明しています。

各クロマトグラムに多変量正規分布の当てはめを行う

得られたクロマトグラムに対してひとつのピークにひとつの正規分布を当てはめ、バンドサイズを定量します。

まず、ピーク位置とその高さを入手します。scipy.sginalのfind_peakを利用します。

from lmfit import models
import random

def get_peak(lane,th=0.001):
  x =np.arange(0,len(lane))

  #Find peaks
  peaks = signal.find_peaks(lane,height=th)
  height = peaks[1]['peak_heights'] #list of the heights of the peaks
  peak_pos = x[peaks[0]] #list of the peaks positions

  #Plotting
  fig = plt.figure()
  ax = fig.subplots()
  ax.plot(x,lane)
  ax.scatter(peak_pos, height, color = 'r', s = 15, marker = 'D', label = 'Maxima')
  ax.legend()
  ax.grid()
  plt.show()

  return peak_pos,height

peak_pos, height = get_peak(lane)

正規分布の当てはめにはlmfitを使用します。こちらの記事を参考に（というか丸パクリ）させていただきました。ありがとうございます。

pip install lmfit

#モデルを作るための関数
def generate_model(spec):
    composite_model = None
    params = None
    x = spec['x']
    y = spec['y']
    x_min = np.min(x)
    x_max = np.max(x)
    x_range = x_max - x_min
    y_max = np.max(y)
    '''
    モデルの数だけfor分を回す
    enumerateはlistの番号を返す
    '''
    #このforもよく使う
    #i にはインデックスが格納
    #basis_funcには要素が格納
    #specオブジェクトのモデルリスト内の要素が入ってくるはず
    for i, basis_func in enumerate(spec['model']):
        #prefixを指定　(ex) fm1,fm2,fm3....
        prefix = f'm{i}_' 
        
        #getattrメソッドにより属性の値を取得する
        #getattr(object, name[, default])のような使い方
        #modelにはmodelとprefixを含んだものがリストとして返る
        '''
        print(model)を実行した場合
        <lmfit.Model: Model(gaussian, prefix='m0_')>
        <lmfit.Model: Model(gaussian, prefix='m1_')>.....
        '''
        model = getattr(models, basis_func['type'])(prefix=prefix)
        
        
        if basis_func['type'] in ['GaussianModel', 'LorentzianModel', 'VoigtModel']: # for now VoigtModel has gamma constrained to sigma
            '''
            set_param_hintメソッドについて
            set_param_hint(name, **kwargs)の形をとるメソッド
            
            optinal引数には
            value:パラメータの値
            vary :フィット中にパラメータが変化するかどうか（デフォルトではTrue）
            min  :パラメータの最小値
            max  :パラメータの最大値
            expr :パラメータの拘束を行うための条件式（評価関数みたいなもの
            
            '''
            model.set_param_hint('sigma', min=1e-6, max=x_range)
            model.set_param_hint('center', min=x_min, max=x_max)
            model.set_param_hint('height', min=1e-6, max=1.1*y_max)
            model.set_param_hint('amplitude', min=1e-6)
            
            
            # default guess is horrible!! do not use guess()
            #デフォでは推測メソッドがあるけどくそなので使うなとのこと
            default_params = {
                prefix+'center': x_min + x_range * random.random(),
                prefix+'height': y_max * random.random(),
                prefix+'sigma': x_range * random.random()
            }
            
        #例外処理
        else:
            raise NotImplemented(f'model {basis_func["type"]} not implemented yet')
        if 'help' in basis_func:  # allow override of settings in parameter
            #itemメソッドは辞書オブジェクトのタプル型を返すメソッド
            for param, options in basis_func['help'].items():
                model.set_param_hint(param, **options)
        model_params = model.make_params(**default_params, **basis_func.get('params', {}))
        if params is None:
            params = model_params
            #make_paramsメソッドで作成されたデフォ値が入る
        else:
            params.update(model_params)
            #拘束paramが指定されてたらそれでオブジェクト（構造体）のアップデートをかける
        if composite_model is None:
            composite_model = model
        else:
            composite_model = composite_model + model
    #返り値は複数
    return composite_model, params

generate_model関数に渡す関数として、get_modelではクロマトグラムとピーク位置・高さ、初期分散を指定しています。これによりフィッティングを行います。

def get_model(lane,peak_pos,height,sigma=0.05):
  x=np.arange(0,len(lane))
  y=lane

  spec = {
    'x': x,#ここは各自変える
    'y': y,
    'model': [
        {'type': 'GaussianModel','params': {'center': 0, 'height': 0, 'sigma': sigma}} for i in range(len(peak_pos))
        ]
  }

  for i,(p,h) in enumerate(zip(peak_pos,height)):
    spec['model'][i]['params']['center'] = p
    spec['model'][i]['params']['height'] = h

  model, params = generate_model(spec)
  output = model.fit(spec['y'], params, x=spec['x'])
  output.plot(data_kws={'markersize':  1},title="Model fit")
  plt.show()

  fig, ax = plt.subplots(figsize=(10,4))
  ax.scatter(spec['x'], spec['y'], s=4)
  components = output.eval_components(x=spec['x'])
  print(len(spec['model']))
  for i, model in enumerate(spec['model']):
      ax.plot(spec['x'], components[f'm{i}_'])
  plt.show()
   

  df = best_values(spec, output)
  

  return df

def best_values(spec, output):
    model_params = {
        'GaussianModel':   ['amplitude', 'sigma'],
        'LorentzianModel': ['amplitude', 'sigma'],
        'VoigtModel':      ['amplitude', 'sigma', 'gamma']
    }
    best_values = output.best_values
    df = pd.DataFrame(index=np.arange(0,len(spec['model'])),columns=['center','amplitude','sigma'])
    for i, model in enumerate(spec['model']):
        prefix = f'm{i}_'
        for col in ["center",'amplitude', 'sigma']:
          df.loc[i,col] = best_values[prefix + col]
    return df

df = get_model(lane,peak_pos,height,sigma=1)

一つの正規分布が一つのバンドとみなし、面積を定量

ただし、これにより得られたモデルではほとんど高さのない正規分布やノイズのような正規分布が現れたりします。それらは排除し、さらに amplitude によって面積を定義して、全体を1とします。

df = df[(0.01<df["sigma"])&(df["sigma"]<50)]
df = df[(0.01<df["amplitude"])&(df["amplitude"]<1)]
df["amplitude"] = df["amplitude"]/df["amplitude"].sum()

この操作をすべてのレーンに対して行います。

lane_df = pd.DataFrame()
for i,lane in enumerate(lanes):
  peak_pos, height = get_peak(lane)
  df = get_model(lane,peak_pos,height,sigma=1)
  df["lane"] = i

  df = df[(0.01<df["sigma"])&(df["sigma"]<50)]
  df = df[(0.01<df["amplitude"])&(df["amplitude"]<1)]
  df["amplitude"] = df["amplitude"]/df["amplitude"].sum()
  lane_df = pd.concat([lane_df,df])

正規分布の位置を検量線によって分子量に変換

先に作った分子量の検量線を用いて、分布の中央にあたる分子量を計算します。ただし、分布なのでsigmaから分子量の上と下も計算しておきます。

lane_df["mw"] = lane_df["center"].map(lambda x: ladder[int(x)])
lane_df["mw_high"] = (lane_df["center"]+(lane_df["sigma"])).map(lambda x: ladder[int(x)])
lane_df["mw_low"] = (lane_df["center"]-(lane_df["sigma"])).map(lambda x: ladder[int(x)])

バンドの分子量と濃度が算出できる

これによりバンドの分子量と相対濃度が決定できました。ただし、バンドの位置が同じかどうか等は分子量が完全には一致しないと思いますので、修正すべき点かと思います。また、正規分布が必ずしも当てはまらない可能性もありますし、サチレーションしている場合は正確な値にならないかと思います。また、最後の図示はあえてしませんでしたが、どのような見せ方をすべきかも迷うところです。一例を↓に置きますが、もう一度ゲルを作っているような感じになってしまいました。
どなたかがより良いモノを作ってくださるのを待っています。

参考文献

Python:lmfitを使ったカーブフィッティングその2(https://gsmcustomeffects.hatenablog.com/entry/2019/01/01/030623)