ヒト22番染色体を Burrows-Wheeler 変換してランの分布を見る

前回はバイナリファイルの中のランの分布を見るプログラムを書いた

• バイナリデータ中のランの度数分布をとる

これを使ってヒト染色体データを BWT したときのランの分布を見てみる．
ちなみに 22番染色体は https://github.com/zfy0701/Parallel-LZ77/tree/master/testData でゲット

• かなり長いランが出現している
• べき分布らしき形状
• これは繰り返し文字列が多い兆候．試しに DNA 文字列をシャッフルしたあとに BWT してランの分布を見てみると以下の通り

BWT のコード (sais.hxx を利用)

/* bwt.cpp */

#include <fstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <random>

#include "sais.hxx"

using namespace std;

template <class UINT_T>
vector<UINT_T> bwt(string filename){
  vector<UINT_T> S;
  ifstream fin(filename);
  UINT_T buf;

  bool do_shuffle = true;
  //bool do_shuffle = false;

  if(!fin){
    cerr << filename << "cannot open" << endl;
    exit(-1);
  }

  // Read file
  while(!fin.eof()){
    fin.read((char*)&buf, sizeof(buf));
    S.push_back(buf);
  }
  fin.close();

  if(do_shuffle){
    mt19937 rand(0);
    shuffle(S.begin(), S.end(), rand);
  }

  vector<int32_t> SA(S.size(),0);
  vector<UINT_T> BWT(S.size(),0);

  int32_t n = S.size();
  saisxx(S.begin(), SA.begin(), n);

  for(int i=0;i<S.size();++i){
    if(SA[i]==0){
      BWT[i] = S[S.size()-1];
    } else {
      BWT[i] = S[SA[i]-1];
    }      
  }

  return BWT;
}

int main(){
  string filename = "chr22.dna";
  auto BWT = bwt<uint8_t>(filename);

  ofstream out(filename + ".shuffle.bwt");
  for(auto c : BWT){
    out.write((char*)&c, sizeof(uint8_t));
  }
}