RustでOpenCV使ってみた

RustからOpenCVのロジックを呼び出す理由は特に無いのですが、Rustが好きなのでやってみました。

参考

事前準備

MacOSでは brew install opencvOpenCVをインストールしたのち、以下のように Cargo.toml を作成することで使用可能でした。

[package]
name = "ocvtest"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
opencv = "0.84.5"

メイン

今回はArUcoマーカーを検出し、検出結果を画像に描画し表示するというプログラムを書いてみました。感想としては、例えば画像のリサイズ結果を入力画像用の変数に上書きするようなところでRust特有の参照の制限で少し冗長になってしまう気がしています(C++では入力と出力が同じでも大丈夫だったはず)。また、デフォルト引数の機能がないのも少し使いづらいかなと思います。

use opencv::prelude::*;
use opencv::core::{Scalar, Size};
use opencv::types::{VectorOfMat};
use opencv::imgcodecs::{imread, IMREAD_COLOR};
use opencv::objdetect::{
    ArucoDetector,
    ArucoDetectorTraitConst,
    DetectorParameters,
    RefineParameters,
    get_predefined_dictionary,
    PredefinedDictionaryType,
    draw_detected_markers,
};
use opencv::imgproc::{resize, INTER_LINEAR};
use opencv::highgui::{named_window, imshow, wait_key, WINDOW_NORMAL};

fn main() {

    // read image
    let img = imread("test.jpeg", IMREAD_COLOR).unwrap();

    // resize image
    let mut dst: Mat = Default::default();
    resize(&img, &mut dst, Size::new(0, 0), 0.2, 0.2, INTER_LINEAR).unwrap();
    let mut img = dst;

    // ArUco predefined dictionary
    let dictionary = get_predefined_dictionary(PredefinedDictionaryType::DICT_6X6_50).unwrap();

    // ArUco detector parameters
    let detector_params = DetectorParameters::default().unwrap();
    
    // ArUco refine parameters
    let refine_params = RefineParameters::new(10., 3., true).unwrap();

    // construct ArUco detector
    let aruco_detector = ArucoDetector::new(&dictionary, &detector_params, refine_params).unwrap();

    // declare some variables for the ArUco detection result
    let mut corners: VectorOfMat = Default::default();
    let mut ids: Mat = Default::default();
    let mut rejected_img_points: VectorOfMat = Default::default();

    // ArUco detection 
    aruco_detector.detect_markers(&img, &mut corners, &mut ids, &mut rejected_img_points).unwrap();

    // draw detection result
    let color = Scalar::new(0., 255., 0., 0.);
    draw_detected_markers(&mut img, &corners, &ids, color).unwrap();

    // show detection result on a window
    named_window("main", WINDOW_NORMAL).unwrap();
    imshow("main", &img).unwrap();
    wait_key(0).unwrap();

}

実行結果

ShinyProxyとnginxで解析用サーバを構築する

nginxとShinyProxyを使ってShinyアプリをホストするサーバを構築することがあったのでその作業を忘れないようにメモします。

以下システムのイメージ図です。

環境

  • Ubuntu 22.04
  • Docker 24.0.1
  • nginx 1.18.0
  • ShinyProxy 3.0.0

NGINXをリバースプロキシとして動かす

ShinyProxyへのリバースプロキシのための以下の設定を追加します。今回は /etc/nginx/sites-availables/default に追記する形で設定しましたがちゃんとした作法かは分かりません。 nginx自体は apt でインストールしました。

map $http_upgrade $connection_upgrade {
  default upgrade;
  '' close;
}

server {
  listen 80;

  location /shiny/ {
    proxy_pass http://127.0.0.1:8080;
    proxy_http_version 1.1;

    proxy_set_header Host $host;
    proxy_set_header X-Real-IP $remote_addr;
    proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;
    proxy_set_header X-Forwarded-Protocol $scheme;

    # websocket headers
    proxy_set_header Upgrade $http_upgrade;
    proxy_set_header Connection $connection_upgrade;
  }
}

(https://support.openanalytics.eu/t/using-shinyproxy-with-nginx-and-a-location-directive-other-than/693/2 (2023/07/01訪問)より引用、一部改変)

location の後の /shiny/ の部分は任意のURIにすることができます。この部分は後述する application.ymlcontent-path の設定と一貫する必要があります。

nginxの有効化と起動を行います。ここでエラーが出る場合は編集した設定ファイルにミスがある可能性があります。

systemctl enable nginx
systemctl start nginx

ShinyProxyをDocker上で動かす

https://github.com/openanalytics/shinyproxy-config-examples/tree/master/02-containerized-docker-engine から Dockerfileapplication.yml をダウンロードして作業ディレクトリに保存します。 application.yml の方は以下の設定を proxy: の前に追加します。

server:
  servlet:
    context-path: /shiny

(https://support.openanalytics.eu/t/using-shinyproxy-with-nginx-and-a-location-directive-other-than/693/2 (2023/07/01訪問)より引用、一部改変)

ShinyProxyのDockerイメージをビルドする前に、ShinyProxyと各Shinyアプリが通信するためのネットワークを以下のコマンドで作っておきます。ネットワーク名は任意ですが、application.yml 内の container-network の設定と一致する必要があります。先述のGitHubリポジトリからダウンロードするとその値はsp-example-netになっていますが、今回はshinyproxy-netに変更し、application.ymlの方も変更してあります。

sudo docker network create shinyproxy-net

(https://github.com/openanalytics/shinyproxy-config-examples/tree/master/02-containerized-docker-engine (2023/07/01訪問)より引用、一部改変)

以下のコマンドでビルドを行います。Dockerイメージにつける名前も任意です。これはapplication.ymlの内容が変わるたびに実行します。(もしくはdocker cpdocker restartを利用する方法でもいけると思います。)

sudo docker build . -t my-shinyproxy

(https://github.com/openanalytics/shinyproxy-config-examples/tree/master/02-containerized-docker-engine (2023/07/01訪問)より引用、一部改変)

ビルドしたDockerイメージからコンテナを起動します。--net オプションは先ほど作成したネットワーク名、最後の引数には先ほどビルドしたDockerイメージの名前を指定します。

sudo docker run --restart=always -d --name shinyproxy -v /var/run/docker.sock:/var/run/docker.sock:ro --group-add $(getent group docker | cut -d: -f3) --net shinyproxy-net -p 8080:8080 my-shinyproxy

(https://github.com/openanalytics/shinyproxy-config-examples/tree/master/02-containerized-docker-engine (2023/07/01訪問)より引用、一部改変)

確認

ブラウザから http://localhost/shiny/ にアクセスして、ShinyProxyのログイン画面が見られたら成功です。

参考

react-well-platesを使ってみた

前回ちらっと紹介したReactコンポーネントライブラリreact-well-platesを使ってみました。公式のStorybookが参考になりました。

下準備

まず、create-react-appによるアプリの骨組み作成とnpmによるライブラリのインストールを行います。

npx create-react-app my-wellplate-app
cd my-wellplate-app
npm install react-well-plates

メインコードの作成

my-wellplate-app/src/App.jsの内容を編集し、以下のような内容で保存します。propsrowscolumnsは必須項目のようです。

import logo from './logo.svg';
import './App.css';
import { WellPlate } from 'react-well-plates';

function App() {
  const props = {
    rows: 8,
    columns: 12
  };
  return (
    <WellPlate {...props} />
  );
}

export default App;

アプリの起動

以下のコマンドでアプリを起動します。自動でWebブラウザが立ち上がり基本的なwellplateの図が表示されます。

npm start

様々な設定

上記で示したWellPlateコンポーネントにはpropsを編集することで様々な設定ができます。

// 設定例
const props = {
    rows: 8, // N行のwellplateを表示 96なら8、384なら16といった感じ
    columns: 12, // M行のwellplateを表示 96なら12、384なら24といった感じ
    wellSize: 30, // 1wellのサイズデフォルトは多分40(px)
    wellStyle: () => ({
        fontSize: 'x-small',
    }),  // wellのstyle情報を返す関数を設定
    renderText: ({index, label}) => {
        return (
            <div>
                <div>{index}</div>
                <div>{label}</div>
            </div>
        );
    },  // wellの中に表示する内容を返す関数を設定
    displayAsGrid: true  // グリッド表示(格子状表示)
};

Well選択可能なwellplate

クリックによるwellの選択を可能にするためには<WllPlate />の代わりに<MultiWellPicker />を使います。また、ReactのuseState関数を使用します。

設定必須の項目として、valueonChangeが増えます。

import logo from './logo.svg';
import './App.css';
import { useState } from 'react';
import { MultiWellPicker } from 'react-well-plates';

function App() {
    const initialValue = [3, 4, 5];  // 最初に選択されているwellのインデックス
    const [value, setValue] = useState(initialValue);
    const props = {
        rows: 8,
        columns: 2,
        value: value,
        onChange: setValue,
        disabled: [0, 1, 2],  // 選択できないwellを指定できる(任意)
    };
    const selected_list = value.map((v) => <li>{v}</li>); // 選択したwellのインデックスをリスト表示するための処理
    return (
        <div>
            <MultiWellPicker {...props} />
            <ul>{selected_list}</ul>
        </div>
    );
}

export default App;

React使って96well plate表現してみた

menseki.hatenablog.jp

この記事を書いたときに96 well plateのインターフェースをjavascriptで書いてみたいと思っていたのですが、なんとなくReactのドキュメント読んで適当に書いてみたら自分でもよくわからないまま、なんか動くものができたので一旦ここに記録を残しておきたいと思います。

今この記事を書いてる時点ですでに当コードを書いてから時間が経っていて自分でもどうやって開発したかわからなくなってしまったのでコードの添付だけですが、確かcreate-react-appを使ってどうのこうのした気がします。

正直すでに世の中にはreact-well-platesなるものが存在するので実際にはそれを使うのが良い気がしますが、使い方がわかりません。

index.css

body {
  font: 14px "Century Gothic", Futura, sans-serif;
  margin: 20px;
}

ol, ul {
  padding-left: 30px;
}

.plate-row:after {
  clear: both;
  content: "";
  display: table;
}

.well {
  background: #fff;
  border: 1px solid #999;
  float: left;
  font-size: 24px;
  font-weight: bold;
  line-height: 34px;
  height: 34px;
  margin: 2px;
  padding: 0;
  text-align: center;
  width: 34px;
  border-radius: 17px;
}
.selected {
  border: 3px solid #999;
}

.well:focus {
  outline: none;
}
.label {
  background: #fff;
  color: #444;
  border: 0px solid #999;
  float: left;
  font-size: 24px;
  font-weight: bold;
  line-height: 34px;
  height: 34px;
  margin: 2px;
  padding: 0;
  text-align: center;
  width: 34px;
  border-radius: 17px;
}

.label:focus {
  outline: none;
}

.kbd-navigation .well:focus {
  background: #ddd;
}

.wellplate {
  display: flex;
  flex-direction: row;
  user-select: none;
}

.wellplate-plate {
  border: 1px solid #999;
  padding: 30px;
}
.plate {
  border: 2px solid #999;
  padding: 10px;
}

.game-info {
  margin-left: 20px;
}

index.js

import React from 'react';
import ReactDOM from 'react-dom/client';
import './index.css';

function Well(props) {
  const rowLabels = ['\\', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'];
  let className;
  let value;
  if (props.row === 0 || props.col === 0) {
    className = 'label';
    value = props.col === 0 ? rowLabels[props.row] : props.col;
  } else {
    className = 'well' + (props.isConfirmed || props.isSelected ? " selected" : "");
  }
  return(
    <button
      className={className}
      onMouseDown={props.onMouseDown}
      onMouseUp={props.onMouseUp}
      onMouseOver={props.onMouseOver}
    >
      {value}
    </button>
  );
}

class Plate extends React.Component {
  renderWell(row, col) {
    const rowLabels = ['\\', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'];
    return (
      <Well
        key={rowLabels[row]+col.toString()}
        row={row}
        col={col}
        isConfirmed={this.props.isConfirmed[(row-1)*12 + (col-1)]}
        isSelected={this.props.isSelected[(row-1)*12 + (col-1)]}
        onMouseDown={(e) => this.props.onMouseDown(e, row, col)}
        onMouseOver={(e) => this.props.onMouseOver(e, row, col)}
        onClick={(e) => e.stopPropagation() }
      />
    );
  }
  
  renderRow(row) {
    const rowLabels = ['\\', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z']
    return(
      <div key={row} className="plate-row">
        { Array(13).fill(null).map((_, col) => this.renderWell(row, col)) }
      </div>
    );
  }

  render() {
    return (
      <div className="plate" onMouseLeave={this.props.onMouseLeave} onMouseUp={this.props.onMouseUp}>
        {Array(9).fill(0).map((_, row) => this.renderRow(row))}
      </div>
    );
  }
}

class WellPlate extends React.Component {

  constructor(props) {
    super(props);
    this.state = {
      isConfirmed: Array(96).fill(false),
      isSelected: Array(96).fill(false),
      wellMouseDown: [null, null],
    };
  }

  selectWells(row0, col0, row, col) {
    if (row0 > row) [row0, row] = [row, row0];
    if (col0 > col) [col0, col] = [col, col0];

    let isSelected = Array(96).fill(false);
    if (row0 === 0 && col0 === 0) {
      [row, col] = [8, 12];
    } else if (row0 === 0) {
      row0 = 1;
      row = (row === 0 ? 8 : -1);
    } else if (col0 === 0) {
      col0 = 1;
      col = (col === 0 ? 12 : -1);
    }
    for (let r = row0-1; r < row; r++) {
      for (let c = col0-1; c < col; c++) {
        isSelected[r * 12 + c] = true;
      }
    }
    return isSelected;
  }

  handleMouseDown(e, row, col) {
    console.log(['down', row, col, e.ctrlKey]);
    if (!e.ctrlKey) {
      this.setState({
        isConfirmed: Array(96).fill(false),
      });
    }
    this.setState({
      wellMouseDown: [row, col],
    });
    e.stopPropagation();
  }
  handleMouseDown2(e) {

    this.setState({
      isConfirmed: Array(96).fill(false),
    });
  }
  
  handleMouseUp() {
    const isSelected = this.state.isSelected.slice();
    let isConfirmed = this.state.isConfirmed.slice();
    for (let i = 0; i < 96; ++i) {
      isConfirmed[i] |= isSelected[i];
    }
    this.setState({
      isConfirmed: isConfirmed,
      wellMouseDown: [null, null],
    });
  }
  
  handleMouseOver(e, row, col) {
    console.log(['over', row, col]);
    let row0, col0;
    if (this.state.wellMouseDown[0] === null && this.state.wellMouseDown[1] === null) {
      [row0, col0] = [row, col];
    } else {
      [row0, col0] = this.state.wellMouseDown;
    }
    const isSelected = this.selectWells(row0, col0, row, col);
    this.setState({
      isSelected: isSelected,
    })
  }

  handleMouseLeave() {
    this.setState({
      isSelected: Array(96).fill(false),
    })
  }
  
  render() {
    const isConfirmed = this.state.isConfirmed.slice();
    const isSelected = this.state.isSelected.slice();
    return (
      <div className="wellplate" onContextMenu={ (e) => e.preventDefault() }>
        <div className='wellplate-plate' onMouseDown={() => this.handleMouseDown2()}>
          <Plate
            isConfirmed={isConfirmed}
            isSelected={isSelected}
            onMouseDown={(e, row, col) => this.handleMouseDown(e, row, col)}
            onMouseUp={() => this.handleMouseUp()}
            onMouseOver={(e, row, col) => this.handleMouseOver(e, row, col)}
            onMouseLeave={() => this.handleMouseLeave()}
            onMouseDown2={() => this.handleMouseDown2()}
          />
        </div>
      </div>
    );
  }
}

const root = ReactDOM.createRoot(document.getElementById("root"));
root.render(<WellPlate />);

実行結果

Rで多重forを回避するTips

多重forが発生するケース

例えば3個のパラメータについて、いくつかのパターンの総組合せを試したいということ があるかと思います。

f <- function(x, y, z) {
  # x, y, zを使って何かするコード
  sprintf("x: %d, y: %d, z: %d", x, y, z)
}

x_list <- c(1, 2)
y_list <- c(100, 50, 20)
z_list <- c(1, 0)

results<- list()
for (x in x_list) {
  for (y in y_list) {
    for (z in z_list) {
      results <- c(results, f(x, y, z))
    }
  }
}
head(results)
##> [[1]]
##> [1] "x: 1, y: 100, z: 1"
##> 
##> [[2]]
##> [1] "x: 1, y: 100, z: 0"
##> 
##> [[3]]
##> [1] "x: 1, y: 50, z: 1"
##> 
##> [[4]]
##> [1] "x: 1, y: 50, z: 0"
##> 
##> [[5]]
##> [1] "x: 1, y: 20, z: 1"
##> 
##> [[6]]
##> [1] "x: 1, y: 20, z: 0"

このコードの嫌なところは例えば関数fに新しいパラメータwを追加してw_list <- c("a", "b", "c")というベクトルをパラメータの組合せに加えるということをしようとしたとき、変更しなければならない箇所が多いことが挙げられます。

また、見栄えも悪い気がします。

f <- function(x, y, z, w) {
  # x, y, z, wを使って何かするコード
  sprintf("x: %d, y: %d, z: %d, w: %s", x, y, z, w)
}

x_list <- c(1, 2)
y_list <- c(100, 50, 20)
z_list <- c(1, 0)
w_list <- c("a", "b", "c")

results<- list()
for (x in x_list) {
  for (y in y_list) {
    for (z in z_list) {
      for (w in w_list) {
        results <- c(results, f(x, y, z, w))
      }
    }
  }
}
head(results)
##> [[1]]
##> [1] "x: 1, y: 100, z: 1, w: a"
##> 
##> [[2]]
##> [1] "x: 1, y: 100, z: 1, w: b"
##> 
##> [[3]]
##> [1] "x: 1, y: 100, z: 1, w: c"
##> 
##> [[4]]
##> [1] "x: 1, y: 100, z: 0, w: a"
##> 
##> [[5]]
##> [1] "x: 1, y: 100, z: 0, w: b"
##> 
##> [[6]]
##> [1] "x: 1, y: 100, z: 0, w: c"

そんなときに私ならどうするかを紹介しようと思います。

あくまで最近の私がどうしているかでこれが最善ではありません。

expand.gridとpurrrを使う

以下にコードを示します。もし関数f の引数が増えたらその対応はwrapper_fの内部とparam_listの要素の追加になります。新しいコードブロックが増えるようなことはないので多少はシンプルかなと思っています。

欠点としては全組合せがいったんデータフレームに展開されるので組合せ数が爆発するとメモリに乗らなくなるかもしれないという危険があるので結局のところ一長一短ではあります。

f <- function(x, y, z, w) {
  # x, y, z, wを使って何かするコード
  sprintf("x: %d, y: %d, z: %d, w: %s", x, y, z, w)
}

wrapper_f <- function(p) {
  f(p$x, p$y, p$z, p$w)  #### fの引数の数が変わったら変更すべき場所(1) #####
}

param_list <- list(
  x = c(1, 2),
  y = c(100, 50, 20),
  z = c(1, 0),
  w = c("a", "b", "c") #### fの引数の数が変わったら変更すべき場所(2) #####
)
param_list |>
  expand.grid(stringsAsFactors = FALSE) |> # 全ての組合せを網羅するデータフレーム
  purrr::transpose() |> # データフレームを行ごとにイテレートする(直感的でないのがマイナスポイント)
  purrr::map(wrapper_f) |>
  head()
##> [[1]]
##> [1] "x: 1, y: 100, z: 1, w: a"
##> 
##> [[2]]
##> [1] "x: 2, y: 100, z: 1, w: a"
##> 
##> [[3]]
##> [1] "x: 1, y: 50, z: 1, w: a"
##> 
##> [[4]]
##> [1] "x: 2, y: 50, z: 1, w: a"
##> 
##> [[5]]
##> [1] "x: 1, y: 20, z: 1, w: a"
##> 
##> [[6]]
##> [1] "x: 2, y: 20, z: 1, w: a"

Rustをスクリプトとして実行する

rust-script(https://rust-script.org/)を使うとRustのソースコードスクリプトとして実行できます。Rust 1.54以上が必要です。

今回は公式のマニュアル(上記URL)からいくつか抜粋して説明します。

以下の説明はUNIX系OSで実行することを前提としています。

インストール

cargoをインストールしていれば、1コマンドで済みます。

$ cargo install rust-script

スクリプトを書く

cargo newでプロジェクトを作ってRustを書くときは依存クレートをCargo.tomlに記載しますが、rust-scriptの場合は//!から始まる行に記載します以下にコード例を示します。 1ファイルで完結するので気軽です。

#!/usr/bin/env rust-script
//! 2つの塩基配列をコマンドライン引数をとして受け取りローカルアラインメントするプログラム
//! 
//! ```cargo
//! [dependencies]
//! bio = "1.1.0"
//! ```
//!
//!  
use bio::alignment::pairwise::Aligner;
fn main() {
    let args: Vec<String> = std::env::args().collect();
    if args.len() < 3 {
        panic!("Please specify two arguments");
    }
    let seq1 = args[1].as_bytes();
    let seq2 = args[2].as_bytes();
    let score = |a: u8, b: u8| if a == b { 1i32 } else { -1i32 };
    let mut aligner = Aligner::new(-5, -1, &score);
    let alignment = aligner.local(seq1, seq2);
    println!("{}", alignment.pretty(seq1, seq2));
}

実行する

あとは他のスクリプト言語スクリプトと同じようにファイルに実行権限を付与するまたはrust-scriptコマンドに渡してあげれば実行できます。

コードを編集して一回目の実行はコンパイルが発生するので多少時間がかかりますが、コンパイル結果はキャッシュされるので二回目以降ははやいです。

$ chmod +x main.rs
$ ./main.rs AATGGCAGGACCA ATGGACCAGGA # rust-script main.rs でもOK
AATGGCA  GGACCA   
         ||||||   
       ATGGACCAGGA

ワンライナー

ワンライナー-eオプションで書けます。依存関係は --dep で指定します。

$ rust-script -e '2_f32.sqrt()'
1.4142135
$ rust-script --dep rand -e 'rand::random::<i32>()'
610692028

フィルター

Perlで標準入力の正規表現マッチ行だけを出力するといったことをする場合、

$ cat main.rs | perl -ne 'print if /^\/\/!/'
//! ```cargo
//! [dependencies]
//! bio = "1.1.0"
//! ```
//! 

となると思います。

これをrust-scriptでやると、--loopオプションとクロージャを用いて、

$ cat main.rs | rust-script --dep regex --loop 'let re = regex::Regex::new(r"^//! ").unwrap(); move |l| if re.is_match(l) { print!("{l}") }'
//! ```cargo
//! [dependencies]
//! bio = "1.1.0"
//! ```
//! 

となります。

Perlの方が簡潔ですが場合によってはRustの方が便利といったことがあるかもしれません。ただしTemplate機能を使えば多少は簡潔に書けるかもしれません。詳しくは公式マニュアルを参照してください。

snakemakeとの連携

rust-scriptで実行できるファイルはsnakemakeスクリプトとしても使用できます。このとき、Rustコード側ではsnakemakeインスタンスからinputoutputparams等で指定した値を参照できます。詳しくはこちら

rule all:
    input:
        "some_input.txt",
    output:
        "some_input.txt",
    params:
        "some_params",
    script:
        "rust_code_for_some_analysis.rs"

FMIndexで高速に塩基配列を検索してみた

今回は以下のリポジトリにあるCライブラリを使用してみた記録です。

https://github.com/TravisWheelerLab/AvxWindowFmIndex

ライブラリのビルド

まずは以下の要領でライブラリのビルド〜インストールを行いました。パス等は各環境で適切な値に変えます。

(MacOSの場合)

$ git clone https://github.com/TravisWheelerLab/AvxWindowFmIndex.git
$ cd AvxWindowFmIndex/
$ git submodule init
$ git submodule update
$ brew install gcc
$ # /path/to/gcc, /prefix/path は適宜変える
$ cmake -DCMAKE_C_COMPILER=/path/to/gcc -DCMAKE_INSTALL_PREFIX=/prefix/path .
$ make
$ make install

プログラムの作成

READMEを参考にFMIndexによる配列検索を行うプログラムをCで書きました。以下、コード全体です。main.cとして保存しています。

fastas/dnaSequence.fasta というファイルから配列を読み込み、コマンドライン引数に与えられた配列を検索し、ヒットしたポジションをプリントします。

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <AwFmIndex.h>

int main(int argc, char **argv) {
    char *indexFileSrc = "indexFiles/index.awfmi";
    char *fastaInputFileSrc = "fastas/dnaSequence.fasta";
    struct AwFmIndex *index;
    struct AwFmIndexConfiguration config = {
        .suffixArrayCompressionRatio = 8,
        .kmerLengthInSeedTable       = 12,
        .alphabetType                = AwFmAlphabetNucleotide,
        .keepSuffixArrayInMemory     = false,
        .storeOriginalSequence       = true
    };
    
    //enum AwFmReturnCode returnCode = awFmCreateIndexFromFasta(&index, &config, sequence, sequenceLength, indexFileSrc, true); origin
    enum AwFmReturnCode returnCode = awFmCreateIndexFromFasta(
        &index, &config, fastaInputFileSrc, indexFileSrc);
    if (awFmReturnCodeIsFailure(returnCode)) {
        printf("create index failed with return code %u\n", returnCode);
        exit(1);
    }
    uint64_t numQueries = argc-1; //first argument is the name of the executable, so skip it.
    struct AwFmKmerSearchList *searchList = awFmCreateKmerSearchList(numQueries);
    if (searchList == NULL) {
        printf("could not allocate memory for search list\n");
        exit(2);
    }
    
    searchList->count = numQueries;
    //initialize the queries in the search list
    for (size_t i = 0; i < numQueries; i++) {
        searchList->kmerSearchData[i].kmerLength = strlen(argv[i+1]);
        searchList->kmerSearchData[i].kmerString = argv[i+1];
    }
    
    //search for the queries
    int numThreads = 3;
    returnCode = awFmParallelSearchLocate(index, searchList, numThreads);
    if (awFmReturnCodeIsFailure(returnCode)) {
        printf("parallel search failed with return code %u\n", returnCode);
        exit(3);
    }
    
    //print all the hits.
    for (size_t kmerIndex = 0; kmerIndex < searchList->count; kmerIndex++) {
        const uint32_t numHitPositions = searchList->kmerSearchData[kmerIndex].count;
        uint64_t *positionList = searchList->kmerSearchData[kmerIndex].positionList;
      
        for (size_t i = 0; i < numHitPositions; i++) {
            printf("kmer at index %zu found at database position %zu.\n", kmerIndex, positionList[i]);
        }
    }
    
    //code cleanup
    awFmDeallocKmerSearchList(searchList);
    awFmDeallocIndex(index);
    
    exit(0);
}

結果の比較用にRustによる自前のFMIndex配列検索プログラムも書きました。以下にコード全体を示します。

#!/usr/bin/env rust-script
//!
//! ```cargo
//! [dependencies]
//! bio = "*"
//! 
//! ```
use std::path::Path;
use bio::io::fasta::{self, FastaRead};
use bio::alphabets;
use bio::data_structures::bwt::{bwt, less, Occ};
use bio::data_structures::fmindex::{BackwardSearchResult, FMIndex, FMIndexable};
use bio::data_structures::suffix_array::{suffix_array};
use std::collections::HashSet;

fn exact_search<P: AsRef<Path> + std::fmt::Debug>(ref_: P, que: P) {
    
    // read reference sequences and concatenate them
    let mut reader = fasta::Reader::from_file(ref_)
        .expect("Failed to open the reference fasta file.");
    let mut record = fasta::Record::new();

    let mut ref_seq = Vec::new();
    let mut ref_len = Vec::new();
    let mut ref_ids = HashSet::new();
    reader.read(&mut record).expect("Failed to parse the reference fasta record.");
    while !record.is_empty() {
        if record.check().is_err() {
            panic!("Rubbish record in the reference fasta");
        }
        ref_seq.extend_from_slice(record.seq());
        ref_seq.push(b'$');
        if ref_ids.contains(record.id()) {
            panic!("Record IDs are duplicated");
        }
        ref_ids.insert(record.id().to_string());
        ref_len.push((record.id().to_string(), record.seq().len()));
        reader.read(&mut record).expect("Failed to parse the reference fasta record.");
    }
    ref_seq.make_ascii_uppercase();
    //ref_seq.push(b'$');

    //  fmindex construction
    let alphabet = alphabets::Alphabet::new(&ref_seq);
    let sa = suffix_array(&ref_seq);
    let bwt = bwt(&ref_seq, &sa);
    let less = less(&bwt, &alphabet);
    let occ = Occ::new(&bwt, 128, &alphabet);
    let fm = FMIndex::new(&bwt, &less, &occ);

    // query fasta read
    let mut reader = fasta::Reader::from_file(que)
        .expect("Failed to open the query fasta file.");
    let mut record = fasta::Record::new();
    reader.read(&mut record).expect("Failed to parse the query fasta record.");
    while !record.is_empty() {
        if record.check().is_err() {
            panic!("Rubbish record in the query fasta");
        }
        let mut seq = record.seq().to_vec();
        seq.make_ascii_uppercase();
        if alphabet.is_word(&seq) {
            let interval = fm.backward_search(seq.iter());
            let mut positions = match interval {
                BackwardSearchResult::Complete(saint) => saint.occ(&sa),
                _ => Vec::new()
            };
            positions.sort();
            for i in positions {
                println!("{} found ad database position {}", record.id(), i);
            }
        }
        reader.read(&mut record)
            .expect("Failed to parse the query fasta record.");
    }

}

fn main() {
    exact_search("fastas/dnaSequence.fasta", "fastas/query.fasta");
}

実行結果

結局コンパイルオプションの適切な与え方はわかりませんでしたがコンパイルできたので適当なFASTAに対して 適当なクエリを検索してみました。以下、実行結果になります。ポジションは正確に出力できているかと思います。

ひとまず動いたので今回はここで終了です。SIMDあたりはもう少し勉強してみたいと思います。

$ # 入力ファイルのサイズ
$ seqkit stat fastas/dnaSequence.fasta
file                      format  type  num_seqs  sum_len  min_len  avg_len  max_len
fastas/dnaSequence.fasta  FASTA   DNA          5   50,000   10,000   10,000   10,000
$ # コンパイル(SIMDの有効化がよくわからなかった)
$ gcc main.c -I prefix/path/include -I AvxWindowFmIndex/lib/FastaVector/src -L prefix/path/lib/ -lawfmindex  -o main
$ export DYLD_LIBRARY_PATH=prefix/path/lib
$ # 実行結果
$ time ./main GCTCGCTGTTCTCACACATGAAACCTGACTGAATACGAT TATTGGAACATTGCACAACCTATAGACTGCGCGAAGGGCA GTGTCAACTTTAAGTTAGCCCATC
kmer at index 0 found at database position 18829.
kmer at index 1 found at database position 28108.
kmer at index 2 found at database position 49980.
./main GCTCGCTGTTCTCACACATGAAACCTGACTGAATACGAT  GTGTCAACTTTAAGTTAGCCCATC  2.90s user 0.19s system 99% cpu 3.110 total
$ # 実行結果(Rust版比較用プログラム。Indexファイル出力がないので早いのは当たり前。検索結果はちゃんと一致している。)
$ time ./main.rs
query1 found ad database position 18829
query2 found ad database position 28108
query3 found ad database position 49980
./main.rs  0.08s user 0.03s system 46% cpu 0.230 total