🍀 TCGAのがんゲノム変異シグネチャーをRubyで抽出してみたった


![signeture.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/61ff1b09-76af-f5ea-d4e6-fe6955118e86.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/402f4c93-ccbc-d447-5c86-434834da98c8.png)

こんにちは。
バイオインフォマティクスネタです。TCGAの公開している突然変異のデータにNMF(非負値行列因子分解)をかけたら、変異のシグネチャーがとれるかな？と単純に考えてやってみたところ、割とうまくいきましたので記事にします。

<font color="DarkGray">気づいたことはコメントや編集リクエストで指摘して頂けるとうれしいです。</font>

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/ca788e72-91ea-9501-1eed-479aa277c648.png) ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/d1c781e8-d39a-c3b4-be2e-4230342157c8.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/52a8f91a-c24c-42be-5664-c0ad931a357e.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/c1e218c3-c3fb-dd0e-7f16-9a619ea04e25.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/00f82420-d021-6aad-5cf8-1bf4ece19371.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/33dbbbf5-f0c8-4f77-2cae-8486775de385.png) ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9c2f44e4-79cf-2053-abab-79f3713210f1.png)
（何のロゴかわかりますか？ 正解は一番最後）

# 🐣 〜基礎知識編〜
## 🌽 変異シグネチャーってなに？

### がん細胞の遺伝子変異にはパターンがある
　がんは遺伝子の病気であると言われています。がん細胞には多数の遺伝子変異があり、増殖、免疫の回避、独自の代謝など、がん細胞としての性質を獲得しているとされます。発がんにいたるまでに細胞はさまざまな遺伝子変異を蓄積します。遺伝子変異の原因としてタバコ・紫外線などの発がん物質が有名ですが、DNAの損傷を修復する機構の破綻も影響しているとされます。
　正常細胞からがん細胞ができるまでには長い時間をかけて環境に適した細胞が選択されていく必要があります。
　蓄積された変異の頻度を調べて、タバコにはタバコ、紫外線には紫外線の特徴的なパターンを探すのがシグネチャーの考え方です。

### 遺伝子変異のパターンを分解する
　がん細胞の突然変異は複数の原因が組み合わさって発生しています。これを分解するためにNMF（非負値行列因子分解）という手法を使います。

## Single Base Substitution (SBS) Signatures
　変異シグネチャーのなかで、[Single Base Substitution (SBS) Signatures](https://cancer.sanger.ac.uk/cosmic/signatures/SBS/)が最も有名です。その他にもシグネチャーはいろいろあるようですが、私は詳しくありません。
　このシグネチャーは点変異と、その左右の塩基をカウントします。

* 塩基の変異のパターン
* 5'末端側の塩基
* 3'末端側の塩基

の場合の数は全部で48通りあります。この48通りについて頻度をカウントして行列を作ってシグネチャーとして定義します。
（Wikipedia https://en.wikipedia.org/wiki/Mutational_signatures より）
![signature](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9b1df3b4-7bbf-2546-b505-4d1c51c9cc40.png)

これだけではよくわからないと思うので、次に具体例を示します。

## 🌌 COSMICの変異シグネチャーをながめる
https://cancer.sanger.ac.uk/cosmic/signatures/SBS/SBS4.tt
下の図に示したのがSignature4です。４番はタバコの有害物質の一つ、ベンゾピレンの変異シグネチャーです。

![Signature-4](https://cancer.sanger.ac.uk/signatures/SBS_vignettes/image014s.jpg)

　それぞれのColumn(列)のバーが、変異の頻度を示しています。左側の水色のバーはC→Aへの変異の頻度を示しています。黒はC→Gの変異、赤はC→Tの変異、灰色はT→Aの変異、黄緑色はT→Cの変異、ピンクがT→Gの変異を示しています。点変異パターンが6種類しかないことに気が付かれましたでしょうか。たとえばG→Tへの変異というパターンがありません。なぜかというと、それは相補的な塩基から見てG→T変異は、C→A変異と同じものであるからです。よってG→Tは、C→Aと同一のものとしてカウントします。

> この２つは同じこと
> A G>T C →
> | | | |
> T C>A G ←

　C→Aの部分を拡大します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/5757cabc-2d66-7970-6080-0432c1c520b2.png)

5'末端側が、A,C,G,Tの4通り、3'末端側もA,C,G,Tの4通りで、合計16通りの変異パターンがあります。CCAがCAAになる頻度が高く、TCGがTAGになる頻度はあまり高くないようですね。

Signatureの全体像を把握するためにCosmicに掲載されている30個のシグネチャー（Version2)を示します。
（最新版Version3はホームページを参照してください。）
![Signature](https://cancer.sanger.ac.uk/signatures_v2/Signature_patterns.png)

## 🍅 〜目標〜

　さて、今回はTCGAの突然変異データを使って塩基置換の行列式をつくり、[Rumale](https://github.com/yoshoku/rumale)という機械学習ライブラリを使ってNMF(非負値行列因子分解)を実行します。これでCOSMICのシグネチャーが一部だけでも再現できたら成功ということにします。

　突然変異のデータはcBioPortalで一般公開されているExonのものを使用します。TCGAの変異のデータはcBioPortalのAPI叩くだけで取得できます。（ただしWhole Genomeではないので精度は低くなりバイアスも入ると思います）

## 🍝 〜方法などについて〜
個別の方法について、過去のQiitaのエントリーへのリンクなど貼っておきます。

### RubyでTCGAのデータベース cBioPortal から情報を取る方法
daru-apiclientという超薄い自作ツールでBioPortalからデータを取得しています。
httpartyでデータを取得しDaruにデータを渡しています。より詳しい情報は下記のページをご覧ください。
[Rubyでがんのゲノム情報を機械学習してみた](https://qiita.com/kojix2/items/f8587c7b9166b72366b6)
https://github.com/kojix2/daru-apiclient

### NMFについて
数学まったくわかりません。アルゴリズムの内容はわからないです。
他の方の書かれたわかりすい記事やネットの記事を参照して雰囲気でやっています。
[非負値行列因子分解（NMF）をふわっと理解する](https://qiita.com/nozma/items/d8dafe4e938c43fb7ad1)

### RubyでNMFを行う方法
RubyでNMFを実行する方法についてはこちらの記事を参照ください。
[Rubyで顔のデータセットをNMF(非負値行列因子分解)してみた](https://qiita.com/kojix2/items/f252eebe84fd37437c41)

### RubyでT-SNEを行う方法
[【Ruby】t-SNEでMNISTを次元圧縮して可視化してみた](https://qiita.com/kojix2/items/3cf79aeae2f21decf4ee)

### Rumaleを実行するときはnumo-linalgをインストールしよう
RumaleでNMFやt-SNEを実行する場合は、numo-linalgの導入をオススメします。numo-linalgをインストールしなくてもRumaleは問題なく実行できるのですが、圧倒的な計算時間の差が発生します。特にコア数が多いCPUではその差が顕著です。numo-linalgは意外とインストールが難しいライブラリです。ソースコードからOpenBLASをインストールして、`gem install -- --with-blas-dir` オプションを指定するのが正攻法だと思います。しかし環境によってはOpenBLASをパッケージから入れて、単に`gem install` するだけでもインストールできるかも知れません。
[numo-linalgをmacで動かしてみた](https://qiita.com/daddygongon/items/73a83565b7969ffc5d1a)

### RubyでPlotlyを使う方法
超マイナーな方法になりますが、iruby-plotlyを使用しています。
控えめに言ってもめちゃくちゃ便利です。流行ってなくても便利なものは便利。
[Rubyでバイオリンプロットを作る【Plotly】](https://qiita.com/kojix2/items/94f9b58b11b9edfbb139)
https://github.com/zach-capalbo/iruby-plotly

### ゲノムの任意の位置の塩基配列を取得する方法
@percipere さまの下記の記事をもとにsamtoolsを使用しました
[ゲノムの任意の位置の塩基配列を取得したい](https://qiita.com/percipere/items/a72e99ab031b2b5a60ca)

以下、「データ準備編」と「クラスタリング＆可視化編」の2つのJupyter Notebookに分けて実行しています。

# 🏭 〜データ準備編〜

## 🐠 目次

### 準備
0. 🔖 セーブポイント用関数を準備する
1. 📖 ライブラリを読み込む

### データ取得編
2. 🌏 TCGAに接続する準備をする
3. 📡 TCGAのスタディのIDを取得する
4. 📡 Mutationのデータを取得する

### 前処理編
5. 🔧 Mutationのデータを一つのDaru🍺ファイルにまとめる
6. 💾 Mutationのデータを保存する
7. 🔧 SNP かつ Chromosome 特定できるMutationを抽出する
8. 💾 SNPのデータを保存する

### Samtools編
9. ✅ samtools の動作確認をする
10. ✅ 両隣の塩基を確認する前に、ncbiBuildのバージョンが37であることを確認しておく
11. 🔧 左右の塩基をsamtoolsで特定する
12. 📚 mutations_snp_trinucleotide のデータ構造の確認
13. 💾 トリヌクレオチドのデータを保存する

### 変異頻度カタログ作成編
14. 🔧 シグネチャーの変異パターンを生成する
15. 🔧 逆相補的配列を求める関数
16. 🔧 シグネチャーの変異パターン（反対側）を生成する
17. 🔧 インデックスラベルを生成する
18. 🔧 変異の頻度をカウントしてカタログを作る
19. 🔧 カタログをNArrayに変換する
20. 💾 catalogのデータを保存する

### NMF実行編
21. 📚 RumaleのNMFの使い方を確認する
22. 🔧 RumaleでNMFを実行して 💾 モデルを保存する

## 0. 🔖 セーブポイント用関数を準備する

```ruby
## 今回は Signature ディレクトリにオブジェクトを保存していく
require 'fileutils'
FileUtils.mkdir_p("Signature")

def sync(str, obj)
    begin
      if obj.nil?
          Marshal.load(File.binread("Signature/#{str}.dat"))
      else
        File.binwrite("Signature/#{str}.dat", Marshal.dump(obj))
        obj
      end
    rescue
      Marshal.load(File.binread("Signature/#{str}.dat"))
    end
end
```



## 1. 📖 ライブラリを読み込む


```ruby
require 'daru/apiclient'  # TCGAのデータを取得します
require 'numo/linalg'     # OpenBLAS
require 'rumale'          # 機械学習
require 'parallel'        # 並列計算
require 'iruby-plotly'    # 可視化
require 'awesome_print'   # 出力をきれいに
require 'rainbow'         # 出力をきれいに

# Paralell で使用するコア数の設定
# Ryzen7のマシンで実行したので15で設定しました
PCORE = 15
```
15



## 2. 🌏 TCGAに接続する準備をする


```ruby
c = Daru::APIClient.new "https://www.cbioportal.org/api"
```

\<Daru::APIClient:0x00005588fd84ff40 @c=#<#<Class:0x00005588fd84ff18>:0x00005588fd84ee38>>



## 3. 📡 TCGAのスタディのIDを取得する


```ruby
study_ids = c.get("/molecular-profiles").filter_rows{|r| r["molecularAlterationType"] == "MUTATION_EXTENDED"}
                            .filter_rows{|r| r["molecularProfileId"].include? "tcga"}["molecularProfileId"]
                            .to_a.find_all{|i| i.include? "_tcga_mutations"}
                            .map{|i| i.gsub("_mutations","")}
```




    ["acc_tcga", "blca_tcga", "brca_tcga", "cesc_tcga", "chol_tcga", "coadread_tcga", "dlbc_tcga", "esca_tcga", "gbm_tcga", "hnsc_tcga", "kich_tcga", "kirc_tcga", "kirp_tcga", "laml_tcga", "lgg_tcga", "lihc_tcga", "luad_tcga", "lusc_tcga", "meso_tcga", "ov_tcga", "paad_tcga", "pcpg_tcga", "prad_tcga", "sarc_tcga", "skcm_tcga", "stad_tcga", "tgct_tcga", "thca_tcga", "thym_tcga", "ucec_tcga", "ucs_tcga", "uvm_tcga"]



## 4. 📡 Mutationのデータを取得する

ここでは、TCGAのスタディごとにMutationのデータを取得します。

```ruby
mutations = study_ids.map.with_index do |id, index|
  puts [id, index]
  df = c.get("/molecular-profiles/#{id}_mutations/mutations", query: {sampleListId: "#{id}_sequenced", projection: "DETAILED"})
    
  df2 = df["gene","startPosition","endPosition","referenceAllele","variantAllele","sampleId", "variantType", "ncbiBuild", "studyId"]
end
0
```


## 5. 🔧 Mutationのデータを一つのDaru🍺ファイルにまとめる
Daruは速度や信頼性の面で課題があるのですが、処理の流れのわかりやすさを考えてDaru::DataFrameにまとめました。

```ruby
mutations2 = mutations.inject(&:concat)

# Daruの挙動を一応確認しておく
if mutations.map(&:size).sum == mutations2.size
    puts "OK"
    mutations = mutations2
    mutations2 = nil # ほそぼそとしたメモリ対策
else
    raise
end

# 先頭のデータを確認する
mutations.head(5)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/de2da97f-03b0-3345-7ab4-6d68c5b08bfe.png)



```ruby
# 総数
mutations.size
```

1118839



## 6. 💾 Mutationのデータを保存する


```ruby
mutations = sync("mutations", mutations)
```



## 7. 🔧 SNP かつ Chromosome 特定できるMutationを抽出する


```ruby
# SNPのみを選ぶ
mutations_snp = mutations.where(mutations.variantType.eq "SNP")

# Chromosome 特定できるMutationのみを選ぶ
# (keep_row_if メソッドは遅いのでfilter_rowsを使用)
mutations_snp = mutations_snp.filter_rows{|r| r["gene"]["chromosome"] != nil}

# 先頭のデータを確認する
mutations_snp.head(5)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9ae47cec-8035-bd74-9e15-d0f044579cd7.png)



## 8. 💾 SNPのデータを保存する


```ruby
mutations_snp = sync("mutations_snp", mutations_snp)
```


#### 💻 mutations変数をクリア


```ruby
mutations = nil
```

## 🚀 Samtoolを使って、変異の両隣の塩基を特定する

## 9. ✅ samtools の動作確認をする

* @percipere さまの記事をもとにsamtoolsをセットアップしておく
* [ゲノムの任意の位置の塩基配列を取得したい](https://qiita.com/percipere/items/a72e99ab031b2b5a60ca)


```ruby
puts `samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz 7:150696110-150696118`
```

>7:150696110-150696118
ATCCCCCAG
    


## 10. ✅ 両隣の塩基を確認する前に、ncbiBuildのバージョンが37であることを確認しておく

* GRCh37だけでなく、単なる37というのが含まれているが許容範囲でしょうかね


```ruby
mutations_snp["ncbiBuild"].value_counts
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/faa6047f-4326-4f07-8364-618897fe0b94.png)


## 11. 🔧 左右の塩基をsamtoolsで特定する


```ruby
# IRubyでは Parallelで並列計算の最中に出力を行うとエラーが発生するので logger を準備する
require 'logger'
File.open("Signature/samtools.log",'w'){|file| file = nil}
logger = Logger.new("Signature/samtools.log")
```

```ruby
start_time = Time.now

sample_ids = mutations_snp["sampleId"].to_a.uniq

mutations_snp_trinucleotide = Parallel.map(sample_ids, in_processes: PCORE, progress: "Samtools") do |sample_id|
    
    # ここの効率よくないかも
    mutations_per_sample_df = mutations_snp.where(mutations_snp["sampleId"].eq sample_id)
    
    # Parallelでマルチコアを生かした計算をするため、Rubyの配列に変換する. 最後の to_h はなくても動く
    # その場合Hashではなく、Daru::Vectorになる
    mutations_per_sample = mutations_per_sample_df.each_row.to_a.map(&:to_h)
   
    mutations_per_sample.map do |mutation|
        chromosome       = mutation["gene"]["chromosome"]
        start_position   = mutation["startPosition"]
        end_position     = mutation["endPosition"]
        raise "start_position != end_position" if start_position != end_position
        reference_allele = mutation["referenceAllele"]
        variant_allele   = mutation["variantAllele"]
        trinucleotide = `samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz \\
                         #{chromosome}:#{start_position-1}-#{end_position+1}`
                         .split("\n")[1]
        
        # GRChの検索結果がTCGAの結果と一致していない場合ログを残す
        # IRubyでは Parallelで並列計算の最中に出力を行うとエラーが発生するので logger を利用
        unless trinucleotide.nil?
            if trinucleotide[1] != reference_allele
                logger.info "#{mutation["studyId"]}　tri[1]:#{trinucleotide[1]} != ra:#{reference_allele}"
            end
        end
        
        [trinucleotide, reference_allele, variant_allele]
    end
end

puts time: (Time.now - start_time)
```
    
Samtools: |====================================================================|
{:time=>1287.954323263}


## 12. 📚 mutations_snp_trinucleotide のデータ構造の確認

```ruby
puts mutations_snp_trinucleotide.class
puts "%d samples" % mutations_snp_trinucleotide.size 
puts mutations_snp_trinucleotide[0]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/557acd58-8374-b270-4871-4fd76d06c15a.png)



## 13. 💾 トリヌクレオチドのデータを保存する


```ruby
mutations_snp_trinucleotide = sync("mutations_snp_trinucleotide", mutations_snp_trinucleotide)
```


## 14. 🔧 シグネチャーの変異パターンを生成する


```ruby
substitution_class = [['C' , 'A'],
 ['C' , 'G'],
 ['C' , 'T'],
 ['T' , 'A'],
 ['T' , 'C'],
 ['T' , 'G']]

ACGT = %w(A C G T)

mutation_types = []

substitution_class.each do |u,v|
  ACGT.each do |j|
    ACGT.each do |n|
      mutation_types << ["" << j << u << n, v]
    end
  end
end

mutation_types = mutation_types.flatten(0)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/55a66aeb-c5dd-edab-a2f5-2844f8ed8043.png)




## 15. 🔧 逆相補的配列を求める関数


```ruby
def reverse_complement(str)
  return nil if str == nil
  complement = {"A" => "T", "T" => "A", "C" => "G", "G" => "C"}
  str.reverse.chars.map do |s|
      complement[s]
  end.join
end
```


```ruby
## 動作確認
reverse_complement("AAACTG")
```

"CAGTTT"



## 16. 🔧 シグネチャーの変異パターン（反対側）を生成する


```ruby
mutation_types_reverse = mutation_types.map do |pattern|
    pattern.map(&method(:reverse_complement))
end
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1e8625b4-685f-22e4-fdfd-761443918150.png)


## 17. 🔧 インデックスラベルを生成する
これは📊グラフの軸ラベルなどに使用するものです

```ruby
signature_index = mutation_types.map{|i| i.join " > "}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/e979e75c-6e99-8485-8280-6ae2a2c9d2b2.png)

```ruby
File.binwrite("Signature/signature_index.dat", Marshal.dump(signature_index))
```

## 18. 🔧 変異の頻度をカウントしてカタログを作る
* カウントできずにエラーになったものが赤い文字で表示されます


```ruby
mutations_snp_trinucleotide.sum(&:size)
```

```ruby
catalog = mutations_snp_trinucleotide.map do |mutations_per_sample|
  mutation_count = Array.new(96, 0)
  mutations_per_sample.each do |mutation|
      m = mutation.values_at(0,2)
      index = mutation_types.index(m) || mutation_types_reverse.index(m)
      if index.nil?
          # 処理に失敗したものを表示する
          print m.inspect.red
          next
      end
      mutation_count[index] += 1
  end
  mutation_count
end
0
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/d58dedc1-d3dd-08bf-f84a-9f184b2962f9.png)



## 19. 🔧 NArrayに変換する


```ruby
catalog = Numo::DFloat[*catalog]
puts catalog.inspect
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/5ec90e7a-0d3e-abc7-665e-c25804bf27eb.png)



## 20. 💾 catalogのデータを保存する


```ruby
catalog = sync("catalog", catalog)
```



## 21. 📚 RumaleのNMFの使い方を確認する


```ruby
ri Rumale::Decomposition::NMF.new
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/6597e0ec-e66c-bb6d-ee7c-a28864479e1a.png)

    
    


## 22. 🔧 RumaleでNMFを実行して 💾 モデルを保存する

* NMFでは局所解しか求まらないので、NMF500回試行してみることにします。
* コンポーネントの数は30とします。
* modelを保存します
* ここでは他のツールや他の解析のことも考えて500回分試行していますが、このあと実は100回分しか使いません。


```ruby
FileUtils.mkdir_p("Signature/nmf")
Parallel.each(1..500, in_processes: PCORE, progress: "NMF") do |seed|
  decomposer = Rumale::Decomposition::NMF.new(n_components: 30, max_iter: 500, random_seed: seed)
  decomposer.fit_transform(catalog)
  File.binwrite("Signature/nmf/decomposer_#{seed}.dat", Marshal.dump(decomposer))
end
```

    
NMF: |=========================================================================|
1..500

これでデータ準備編はおしまいです。つぎに〜クラスタリング＆可視化編〜にうつります。


# 🎨 〜クラスタリング＆可視化編〜


## 🐳 目次

### 準備
1. 🔖 セーブポイント用関数を準備する
2. 📖 ライブラリを読み込む

### データローディング
3. 📂 シグネチャーのインデックスラベルを読み込む
4. 📂 NArrayのNMF結果を読み込む

### 可視化してみる
5. 📊 NMF結果のヒートマップを表示する
6. 🔧 標準化してみる
7. GaussianMixtureというのでクラスタリングする
8. 📊 ヒートマップを描きなおす

### t-SNE編
9. 📚 Rumaleのt-SNEの使い方を確認する
11. 🔧 t-SNEを実行する
12. 💾 t-SNEの結果を保存する
13. 📊 t-sneの結果を可視化する
15. 📊 クラスタリングもふくめて可視化する

### シグネチャー表示
14. 📊 シグネチャー候補の可視化

### そのほか
15. 💾 CSVファイルを書き出す

## 1. 🔖 セーブポイント用関数を準備する


```ruby
## 今回は Signature ディレクトリにオブジェクトを保存していく
require 'fileutils'
FileUtils.mkdir_p("Signature")

def sync(str, obj)
    begin
      if obj.nil?
          Marshal.load(File.binread("Signature/#{str}.dat"))
      else
        File.binwrite("Signature/#{str}.dat", Marshal.dump(obj))
        obj
      end
    rescue
      Marshal.load(File.binread("Signature/#{str}.dat"))
    end
end
```



## 2. 📖 ライブラリを読み込む


```ruby
require 'numo/narray'     # 行列計算
require 'numo/linalg'     # OpenBLAS t-SNE高速化に重要
require 'iruby-plotly'    # 可視化
require 'rumale'          # 機械学習
```



## 3. 📂 シグネチャーのインデックスラベルを読み込む


```ruby
signature_index = Marshal.load(File.binread("Signature/signature_index.dat"))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/e979e75c-6e99-8485-8280-6ae2a2c9d2b2.png)


## 4. 📂 NArrayのNMF結果を読み込む


```ruby
decomposer = (1..100).map do |i|
    Marshal.load(File.binread("Signature/nmf/decomposer_#{i}.dat"))
end
components = Numo::DFloat.vstack(decomposer.map(&:components))
puts components.inspect
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/c346baac-425e-27e7-590d-d52ee1b2e56d.png)


## 5. 📊 NMF結果のヒートマップを表示する


```ruby
IRuby.plot(
    {x: signature_index,
    z: components.to_a,
    type: "heatmap"},
    
    {title: "TCGA シグネチャー候補",
        xaxis: {tickangle: 300,
            tickfont: {size: 6}},
        height: 600,
        width: 1000})
```

![output_11_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/3081da76-2e56-6ddf-85ac-7558018d17f8.png)



## 6. 🔧 標準化してみる


```ruby
normalizer = Rumale::Preprocessing::StandardScaler.new
samples_std = normalizer.fit_transform(components.transpose)
samples_std = samples_std.transpose
puts samples_std.inspect
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/6839d2ca-a135-a3ee-bab5-27c374ffd17f.png)



```ruby
puts samples_std.var(axis: 1).minmax
```
[0.9999999999999953, 1.0000000000000049]


## 7. 🔧 GaussianMixtureというのでクラスタリングする
* KMeansはサンプルを等分しようとするので実データでは失敗します。
    * 複数のシグネチャーが同じぐらい出現するという条件でシミュレーションを実行するとKMeansでもうまくいきます
* こちらのアルゴリズムの内容もよくわかりません…


```ruby
gmm = Rumale::Clustering::GaussianMixture.new(n_clusters: 30, random_seed: 5)
pca = Rumale::Decomposition::PCA.new(n_components: 50)
samples_std_pca = pca.fit_transform(samples_std)
cluster_ids = gmm.fit_predict(samples_std_pca)
puts cluster_ids.inspect
```

Numo::Int32#shape=[3000]
[5, 9, 21, 29, 12, 7, 14, 24, 20, 27, 2, 11, 20, 25, 18, 29, 3, 18, 22, 7, ...]


## 8. 📊 ヒートマップを描きなおす


```ruby
IRuby.plot(
    {x: signature_index,
    z: samples_std[cluster_ids.sort_index, true].to_a,
    type: "heatmap"},
    {title: "TCGA シグネチャー候補 クラスター＆ソート後",
        xaxis: {tickangle: 300,
            tickfont: {size: 6}},
        height: 600,
        width: 1000}
    )
```

![output_18_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/d45a5568-0c48-7b16-4332-62c0bc88b38f.png)




## 9. 📚 Rumaleのt-SNEの使い方を確認する


```ruby
ri Rumale::Manifold::TSNE.new
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/ccee94b9-cf96-8bf1-9822-004a4805abd7.png)

    
    


⏰ t-SNEは時間がかかるので、何回か実行した時に時間を記録に残せるようにしておく


```ruby
times = []
```



## 10. 🔧 t-SNEを実行する


```ruby
start_time = Time.now
tsne = Rumale::Manifold::TSNE.new(n_components: 2, perplexity: 40.0, max_iter: 500, random_seed: 1)
tsne_representations = tsne.fit_transform(samples_std_pca)
times << Time.now - start_time
```

## 11. 💾 t-SNEの結果を保存する


```ruby
tsne_representations = sync("tsne_representations", tsne_representations)
0
```

## 12. 📊 t-sneの結果を可視化する


```ruby
x = tsne_representations[true, 0]
y = tsne_representations[true, 1]
IRuby.plot(x: x.to_a,
    y: y.to_a,
    mode: 'markers',
    type: "scatter",
    title: "〜NMFで得られたシグネチャ候補をt-SNEで可視化〜",
    width: 800,
    height: 600)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/ebc96310-cd2c-dd89-2d54-852fb23b6e7a.png)


## 13. 📊 クラスタリングもふくめて可視化する


```ruby
data = cluster_ids.to_a.uniq.sort.map do |l|
    {
        x: x[cluster_ids.eq(l)].to_a,
        y: y[cluster_ids.eq(l)].to_a,
        type: "scatter",
        mode: 'markers',
        name: "sig#{l}"
    }
end
IRuby.plot(data, {title: "〜NMFで得られたシグネチャ候補をt-SNEで可視化〜", height: 600, width: 800})
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/76eba517-a508-189f-1ac3-321e727df474.png)


## 14. 🎨 シグネチャー候補の可視化


```ruby
data = 30.times.map do |l|
    s2 = components[cluster_ids.eq(l).where, true]
    s = s2.mean(0)
    color = %w(red orange green blue purple brown).map{|i| [i] * 16}.flatten
    {y: s.to_a,
     x: signature_index,
    xaxis: "x#{l+1}",
    yaxis: "y#{l+1}",
    marker: {color: color},
    type: "bar"}
end

layout = {
      grid: {rows: 10, columns: 3,
      pattern: 'independent'}}
(1..30).each do |i|
    layout["xaxis#{i}".to_sym] = {showticklabels: false}
end

IRuby.plot(data, layout)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/d375e7f2-ec06-6272-8f51-0bd10e1e6db6.png)


## 15. 💾 CSVファイルを書き出す
* 今回はRubyで頑張りましたが、🍊 Orange などのツールを使うとGUIで便利に解析することもできます。
* 他のツールでも何かできるように、CSVファイルで書き出します。

```ruby
require 'csv'
CSV.open("Signature/nmf/signatures.csv", "w") do |csv|
    csv << signature_index
    samples_std.to_a.each do |ar|
        csv << ar
    end
end
0
```

## 結果
　ここまでスクロールしてくれてありがとうございます。
　今回は、GaussianMixtureというアルゴリズムをつかってクラスタリングを行うという単純な方法を採用しています。結果ですが、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9ec66eb5-2f7d-43ce-dbd8-b1aba9bed3ef.png)

Signature1, Signature2, Signature6, Signature7, Signature10, Signature13, Signature14, Signature17, Signature21, Signature28 あたりは検出できた気がします。Signature4 は少し怪しいですが、検出できたことにします。


　Sig1やSig10が何個も検出されています。これらの重複して検出されたSignatureはt-SNE で似たような位置にあることがわかります。候補8や候補13は、サンプル数が少ないのでNMFで誤検出された局所解かもしれません。

　候補12や候補21は、少なくない回数検出されているにもかかわらずCosmicに似たようなSignatureが登録されていないものです。おそらく何らかのミスかアーチファクトだと思いますが、よくわかりません。エクソンのみの変異で本当に結果が出るか心配でしたが、思った以上にたくさんのシグネチャーが検出できました。

　今回ある程度検出に成功したシグネチャー

|COSMIC signature|関連する発がん要因|主ながん種|
|-----------------|------------------|------------------|
Signature1|5-methylcytosineの脱メチル化, 加齢|広く認められる|
Signature2|APOBEC酵素活性|腺がんを中心に広くみられる|
Signature4|喫煙|肺がん・頭頚部がん・肝臓がんなど|
Signature6|DNAミスマッチ修復異常|悪性リンパ腫など|
Signature7|紫外線|メラノーマ|
Signature10|POLE酵素異常|大腸がん・子宮がんなど|
Signature13|APOBEC酵素活性|腺がんを中心に広くみられる|
Signature14|POLE変異とミスマッチ修復異常|悪性リンパ腫など|
Signature17|不明|食道がん・胃がんなど|
Signature21|不明|胃がんなど|
Signature28|不明|大腸がん・子宮がんなど|

## おまけ 🍊 〜Orangeでもう少し詳しく結果をみてみよう〜
　Orangeはスロベニアのリュブリャナ大学が開発しているGUIでデータ解析や機械学習ができるツールです。Jupyterや生Pythonを覚えなくてもクリック操作で解析ができてしまうという、本当に素晴らしいツールです。NMFも[拡張機能](https://blog.biolab.si/2013/02/06/orange-nmf-add-on/)を使用すれば実行できるようで、APIを叩くところを別にすれば、このページで紹介した処理はOrangeでみんな実行できてしまいそうです。ここでは後編の「クラスタリング＆可視化編」をOrangeでやってみることにします。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/f194975c-5781-d16e-b62b-9a170792c6b1.png)

スクリーンショットのようにワークフローを定義して実行します。はい、便利すぎですよ (ToT) 💦
右のt-SNEのウィンドウで観察したいサンプルを選択すると、自動的に下のウィンドウでシグネチャーを描いています。
今描かれているのはSignature4だと思われます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1cb2f100-8624-fc39-efe9-4c8fe0107a88.png)

さらに細かく見ると、ここの小さい集団はSignature8に近いみたいだということがわかりました。

全体的にSignature10の影響が非常に強いです。Signature6や、Signature23は、それに近い候補が存在するが、Signature10番が強烈すぎて混ざってしまって正しく検出できてないのでは？という印象を受けました。

ごく小数ながらSignature20番に近いパターンが検出されていたりします。しかしせいぜい2個〜4個の候補が20番にちょっと似ている風だっただけなので偶然にすぎないかもしれません。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/a09c1c3f-b601-06d6-a0d9-c1d2239ea677.png)

また候補13番はサンプル数が少ないと思っていましたが、単なる雑音として無視できなぐらい再現性を持って検出されるようです。うーん、これは一体何でしょう。どこか私のミスが隠れているのでしょうか。

　<font color="Gray">Orangeはあまり注目されていませんがなぜでしょうか。それは解析APIを提供することでデータセンターに情報を集積するという大手ITプラットフォーマーの企業戦略と合致しないためかもしれません。ツールの良し悪しを判断するのは機械学習エンジニアです。しかし、付加価値を高める競争に日々晒されているであろうエンジニアからみると、Orangeは中途半端な機能しか持たない使えないツールに見えるかもしれません。あるいは逆に、それまでエンジニアが職人技でコードを書いて実現してきたプロセスの一部を、ボタンのクリックを数回やるだけで実行できるようにしているという意味で、機械学習職人の価値を毀損するツールのようにも見えて、反発したくなる心理が生じるかもしれません。ご心配なく。そういう心理は、エンジニアの世界に限らずどこのでもよくあることです。
　しかし、これから「非エンジニア」の人々が使うデータ解析ツールは、PythonやJuliaではなく、Excel+Orangeのようなものになるんじゃないかなあと。個人的にはそんな気がちょっとしたのでした。</font>

## 🎨 おまけ 〜扉絵の作り方〜

ペーパータオルを用意します [^1]
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/5b5e4cf4-e9a9-768f-97f3-71d7cf097193.png)
<font color="Lavender">ぶっちゃけちゃんとした画用紙がある場合はそっちの方がいいです。</font>

ボールペンで下絵を書いて、色鉛筆で塗ります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/8b0038c2-4fb7-b280-60a8-1e0344226242.png)
スキャナーで画像を取り込みます。
<font color="Lavender">↓スキャナーだとシワが写り込んでしまうので、カメラの方がいいかも</font>
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/674838aa-8e76-d10e-34b3-c4564f82b9bf.png)
Gimpで画像を開いて、レベル補正など適当に一通りやったら、`linear invert` して `Hue-Chroma` で色合いを調整します。これでほとんど完成です。あとはお好みで `symmetric nearest neighbor` など気に入ったフィルターをいくつか追加します。

扉絵のできあがり！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/f110d625-e03f-422c-7a0d-cc38df002a69.png)


### 何のロゴかわかるかな？
この記事を書くにあたってお世話になったツールやサイトのロゴです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/ca788e72-91ea-9501-1eed-479aa277c648.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/d1c781e8-d39a-c3b4-be2e-4230342157c8.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/52a8f91a-c24c-42be-5664-c0ad931a357e.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/c1e218c3-c3fb-dd0e-7f16-9a619ea04e25.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/00f82420-d021-6aad-5cf8-1bf4ece19371.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/33dbbbf5-f0c8-4f77-2cae-8486775de385.png)　![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/9c2f44e4-79cf-2053-abab-79f3713210f1.png)
左から順番に、[Ruby](https://www.ruby-lang.org/ja/), [Rumale](https://github.com/yoshoku/rumale), [Plotly](https://plot.ly/), [cBioPortal](http://www.cbioportal.org/), [COSMIC](https://cancer.sanger.ac.uk/cosmic), [Orange](https://orange.biolab.si/), [GIMP](https://www.gimp.org/) でした。

この記事は以上です。

[^1]: コンフォート サービスタオル200　https://pro.crecia.co.jp/product/search/index.php/item?id=11007