ヒトゲノムの配列の特定されてないところを見てみたりする夏休みの絵日記的なやつだよ

## はじめに

この記事はすごいテキトーな記事なので真に受けるなよろし。

## 目的

ヒトゲノム計画では、人のDNAの配列が全部明らかにされたと言われていますが、実際にはわかっていない部分もかなりあるはずです。
どこの部分がわかっていないのか、グラフを描きます。

## ヒトゲノムのfastaのダウンロード

Gencodeからヒトゲノムの配列がダウンロードできます。
https://www.gencodegenes.org/human/



ここでは、贅沢にALLをダウンロードします。

```
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/GRCh38.p13.genome.fa.gz
```

解凍してもいいのですが、`zcat` コマンドを使えばある程度見ることができますね。

まずは `head` を表示してみます。

```sh
zcat GRCh38.p13.genome.fa.gz | head
```

```
>chr1 1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
```

はい、早速出てきました。Nの山です。つまり、配列が特定できていない部分がのっけの頭から来ているわけですね。
次に、一応染色体を確認しておきましょう。

```sh
zcat GRCh38.p13.genome.fa.gz | grep "^>"
```

![hoge.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/59a19a43-e6e6-58d3-7b73-9190622b4b0a.png)

まあ、こんな感じで大量に出てきます。

次に、本当に配列がTCGANだけの5種類の文字列から構成されているのか確認してみることにします。
きっとコマンドラインツールを使う方法もあると思うのですが、ここはCrystal言語を使って小さなプログラムを書いてしまうことにします。なぜCrystal言語かというと爆速だからですね。

```crystal:count.cr
counter = Hash(Char, Int32).new(0)

while l = gets
    l.chomp.each_char do |c|
        counter[c] += 1
    end
end

p counter
```

はい、Rubyわかる人には何の変哲もないコードですね。ポイントは、Hashの生成時に型を指定しているところだけです。ほかは全部Rubyと同じです。これがCrystalです。ビルドします。

```
crystal build counter.cr --release
```

このビルド時間を含めると実行速度でJulia等とあまり変わらないかも知れませんが気にしないことにします。
では、TCGANの数を数えてみましょう

```
zcat GRCh38.p13.genome.fa.gz | grep -v "^>" | ./count
```

`grep -v "^>"` で行頭の染色体などを定義している行を省いています。結果は

```
{'N' => 161331703, 'T' => 916993511, 'A' => 914265135, 'C' => 635937481, 'G' => 638590158}
```

となりました。確かに,TCGAおよびN以外の文字列が存在しないことが確認できました。次はNの連続の分布を数えてみましょう。つまりNが何文字連続するか、を調べてみてみることにします。

```crystal:nseq.cr

temp = 0 

while l = gets
    l.chomp.each_char do |c|
        if c == 'N'
             temp += 1
        elsif temp != 0
            puts temp
            temp = 0
        end
    end
end
puts temp if temp != 0
```

ビルドします。

```
crystal build nseq.cr --release
```

まずはNの連続する箇所がいくつあるか調べてみましょう。

```
zcat GRCh38.p13.genome.fa.gz | grep -v "^>" | ./nseq | wc -l
```

1234 

そんなに多くないことがわかります。これを長い順に並べて表示するとこんな感じになります。

![fuga.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/81b16db2-e9b0-661c-99cf-f984aee27fde.png)

1っていうのが不思議な気がするのと、100とか50000とかが多いなあといった印象ですかね。

## 本題

さて、適当に10000文字ごとに、N, AT, CGの割合を計算してくれるようなプログラムを書きます。

```crystal:n2
tcgan = Hash(Char, Int32).new(0)

target = ARGV[0]
chr = ""
loc = 1
flag = false

while l = gets
  if l.starts_with?(">")
    exit if flag
    if l == target
      puts "loc\tN\tAT\tCG"
      flag = true
    end
    next
  end
  if flag
    l.chomp.each_char do |c|
      tcgan[c] += 1
      loc += 1
      if loc % 10000 == 0
        total = tcgan.values.sum.to_f
        ta = (tcgan['A'] + tcgan['T']) / total
        cg = (tcgan['G'] + tcgan['C']) / total
        n  = tcgan['N'] / total
        puts "#{loc}\t#{n}\t#{ta}\t#{cg}"
        tcgan = Hash(Char, Int32).new(0)
      end
    end
  end
end

```

ためしに実行します。

```
zcat GRCh38.p13.genome.fa.gz | ./n2 ">chr1 1" | head
```

```
loc	N	AT	CG
10000	1.0	0.0	0.0
20000	0.0001	0.4076	0.5923
30000	0.0	0.4826	0.5174
40000	0.0	0.5288	0.4712
50000	0.0	0.6439	0.3561
60000	0.0	0.6346	0.3654
70000	0.0	0.6669	0.3331
80000	0.0	0.6199	0.3801
90000	0.0	0.6294	0.3706
```

だいたいうまくいってるようですね。Rubyで同じ動作のプログラムを作成すると、実行時間がすごくがかかってしまいますがCrystalは超高速です。これをコマンドラインでuplotに投げていきます。

[uplot](https://github.com/kojix2/uplot)というのは、UnicodePlots.rbを使ってターミナル上にグラフを表示できる私が個人的に作っているRuby製のツールです。

ここからターミナルにグラフを描いて、それをスクリーンショットを取ってQiitaに貼り付けるというかなり意味のわからないことをしています。

### chr1

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/8458d1d8-eacf-e37d-c507-5e96a041c7a6.png)

### chr2

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/0746ca87-450c-b866-489a-a2f5d74b5b7f.png)

### chr3

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/a8cf956d-b8fb-240b-c9e3-0d2ce36fa07d.png)

### chr4

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/3a351669-5178-d637-3b15-fb24244b0d72.png)

### chr5

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/b7c663ab-8cf2-affd-8092-0724a9131de1.png)

### chr6

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1252e02a-0f3f-b779-da1d-cb37a226b484.png)

### chr7

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1089ac02-8a83-85b8-dc82-4a96aed3db6a.png)

### chr8

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/853a0cbc-b419-b879-3d21-5b9f560b556e.png)

### chr9

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/96f92235-e6a6-f9a2-a163-9711d207e1e5.png)

### chr10

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/8d1af4fb-a6e6-b70b-e9fd-ae31175e7f00.png)

### chr11

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/38d27e64-e2d3-2bec-2f73-2df755181ed9.png)

### chr12

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/ac8952c3-538d-e705-1dcf-3b7d1016ed1d.png)

### chr13

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1794251b-0004-372b-3bb3-bcecb5f127c3.png)

### chr14

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/6c54e36d-1504-6612-a8fb-c708a518da9e.png)

### chr 15

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1cd832f6-188f-c1ac-67b8-76361ce7c502.png)

### chr 16

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/8c7744a0-b32a-a65b-4b2d-dab743128bfd.png)

### chr 17

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/cb8eeea2-26a2-a6fc-edba-b501a10b3756.png)

### chr 18

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/25a09986-9b79-661f-a5ec-fb9196b65f53.png)

### chr 19

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/f8ee331c-835e-1f64-c458-dca4e8c44bb0.png)

### chr 20

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/97d4cc2f-52f1-2587-410d-c0f06d71fcdc.png)

### chr 21

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/00627029-8701-4ee3-f2f2-17d391245252.png)

scatterも追加
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/523b58b6-cc2c-51e6-dff4-ee171e937f46.png)

### chr 22

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/1952e8ff-b233-e73e-9593-9a434a8bc18e.png)

scatterも追加
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/c0b06759-8f0f-2979-f89f-1eef2e4e28c5.png)

### chr X

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/f8dc7140-7f31-5dcc-2ccc-78502caab669.png)

### chr Y

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/554aceb6-473f-79b4-6164-63669350e225.png)

Wikipediaによると、
主にヘテロクロマチン領域が未解読であり、その領域は活性がないそうです。転写とかほとんどされないという意味だと思います。
つまり生物学的にはあまり意味がないということなのかも知れません。

しかしこうやってみると、配列のわかっていない部分はずいぶんあるんだな〜という感じですね。

詳しい方とかいらっしゃいましたらご自由に突っ込んでください。

この記事は以上です。