bcftoolsでvcfファイルの編集(個人メモ・追記予定)

※2023年7月19日 一部修正

VCF形式

詳しい説明は省略。GATK公式による解説がおそらく一番分かりやすい。 VCFファイルを整形したり、ほしい多型情報を抽出するツールとしてはbcftoolsとVCFtoolsが定番となっている。個人的にはbcftoolの方が、

  • biocondaで配布している
  • bamファイルからの多型抽出コマンドも備えていて、多型抽出からフィルタリングをシームレスに行える

という点で扱いやすいと感じた。

bcftoolsのインストール方法

anaconda環境があればbiocondaからインストールできる。

conda install -c bioconda bcftools

もしくはgithubからインストールする。

VCFの連結

bcftools concat を使う。引数-aをつけると最後のファイルので1行目を上書き。サンプル名、フィールド等が同じであればこれをつける。

# make list
echo -e "file-1.vcf.gz\nfile-2.vcf.gz\nfile-3.vcf.gz" > vcflist

# concatenate tables
bcftools concat -a -f vcflist -o vcf_conc.vcf.gz

# add index
bcftools index -f vcf_conc.vcf.gz

VCFの情報を確認

# list samples
bcftools query -l file.bcf

# number of samples
bcftools query -l file.bcf | wc -l

# list of positions
bcftools query -f '%POS\n' file.bcf

VCFのフィルタリング

-i '条件'を引数に取ることで、条件に合う多型サイトを抽出できる。条件に合うものを除外するには-e '条件'とする。 複数サンプルのmulti-vcfファイルの場合、一つでも条件を満たすサンプルがあれば、そのサイトは抽出される。全サンプルが満たす条件で抽出する場合には、-eを使う。 パイプで繋げる時は-Ou (output uncompressed BCF)をつけると早い。

(例) 1. FILTERフィールドがPASSであり(GTAKなどで事前にフィルターをかけていた場合)、かつQAUL (quality)が100以上であり、 2. 1~4番目の個体のうちgenotypeが0/0である個体が3個体以上の多型を抽出

bcftools view -O u -i 'FILTER="PASS" & QUAL>100' file.vcf.gz |
bcftools view -i 'COUNT(FMT/GT[0-3]="0/0")>2' -o file_filtered.vcf

※フィールドのINFOとFMTの違い INFOは全サンプルを総合した数値。FMT(FORMAT)は個別サンプルにつけられる値を表す。たとえばINFO/DPは全サンプルのdepth合計だが、FMT/DPとすると1サンプルごとのdepthを示す。 ※&と&&の違い &で条件を繋げると、全条件が一つのサンプル内で満たされなければいけない。&&で繋げた場合、各条件がいずれかのサンプルで満たされていればよい。

VCFをテーブルに変換

bcftools query で必要な情報を取り出してテーブルに保存できる。-f '条件' の形式。-H:ヘッダーをつけて保存

(例) タブ区切りでCHROM, POS, REF, ALT, QUAL, AC, AFの値、その右にGT, AD, DPの値を個体ごとにタブで区切って出力、テキスト形式で保存。

bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%AC\t%AF[\t%GT %AD %DP]\n' file_filered.vcf.gz > file_filtered_GTADDP.txt

ジェノタイプ変換

setGTプライグインが提供されている。ただしデータの変換はvcfにあった情報を破壊することになるので慎重に行うべき。一度保存してから外部のツールを使うべきかもしれない。 (例) GQ (Genotype Quality: 個体ごとのGTの確からしさ)が20未満の場合、そのGTを./.に変換

bcftools +setGT input.vcf.gz -o transformed.vcf.gz -- -t q -i 'GQ<20' -n "./."

参考サイト

http://samtools.github.io/bcftools/howtos/filtering.html

。これに関しては情報が少ない。

TrimmomaticでBGIシーケンサーのアダプターを指定する

ここ数年、中国のゲノム関連企業であるBGI次世代シーケンサーがゲノム解読に使われることが増えてきた。配列の増幅にローリングサークル複製機構を使っているなど、illumina社のシーケンサーとはだいぶシステムが異なるらしい。1当然、使われているアダプター配列も異なる。

アダプタートリミングツールであるTrimmomaticはデフォルトでillumina社の主要なシーケンサーに対応しているが、BGIのBGISEQやMGISEQを使用した場合は自分でアダプター配列を指定する必要がある。

検索するとアダプター配列の一覧がヒットし2、その中にアダプターを取り除くには以下の配列を使いなさいとあった。

>Forward_filter
AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA
>Reverse_filter
AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG

これのfastaファイルを作り、Trimmomaticのアダプター格納フォルダに保存する(conda installした場合、shareディレクトリ以下にある)。あとはTrimmomaticのILLUMINACLIPで指定するファイルを変更すれば良い。

#paired endの実行例
adapter="/Users/username/miniconda3/share/trimmomatic-0.39-1/adapters/BGI_adapters.fa"

java -jar /Users/username/miniconda3/share/trimmomatic-0.39-1/trimmomatic.jar \
PE \
-threads 16 \
-phred33 \
-trimlog log.txt \
R1.fastq.gz R2.fastq.gz \ #input files: paired end reads
R1_paired_output.fastq.gz R1_unpaired_output.fastq.gz \ #paired and unpaired trimmed reads from input file 1
R2_paired_output.fastq.gz R2_unpaired_output.fastq.gz \ #paired and unpaired trimmed reads from input file 2
ILLUMINACLIP:${adapter}:2:30:10 \
LEADING:20 \
TRAILING:20 \
SLIDINGWINDOW:4:15 \
MINLEN:36

参考サイト様

  1. [bioinfo]DNBSEQ-G400のSEデータtrimming設定(trimmomatic) – Takeshi Igawa, Ph.D.

  2. Trimmomatic | FASTQ クリーニングツール

集団構造解析ソフトAdmixtureの使い方

形態などのデータからそれなりに多様性があることがわかっている個体群について、各個体の全ゲノム的な多型データが得られているとする。

このとき、 1. その集団を遺伝的に同じ起源を持つグループ(=祖先集団)に分けた場合、どのように分けられるか? 2. 祖先集団の交雑により成立した個体はあるか?

を推定する解析を 集団構造解析 と呼ぶ(あくまで自分の分野の場合)。

集団構造解析に用いられるツールとしては、ひと昔前にはSTRUCTUREがもっぱら使われていたが、近年はより計算が早いとされるAdmixtureによる解析が主流になっている。

自分の研究でも活用させてもらったが、データを読み込む際につまずいたので備忘録的に使い方(の注意点)を書いておく。

インストール

公式サイトからDLして展開する。残念ながらWindows版はない模様。admixtureのアイコンをクリックするとターミナルが開くので、あとはマニュアル通りで基本的には問題ない。

ファイルの準備

ただし、読み込みファイルの形式には注意が必要。マニュアルには、.bed, .ped, .genoの3種類の形式が読めると書いているが、.pedもしくは .genoファイルを使いたい場合は、同じディレクトリに対応する.mapファイルを置いておく必要がある。また、.pedファイルは見本通りにエクセルとテキストエディタで作成しても認識されず、必ずPLINKというソフトフェアで保存したものでなければならない。

解析ファイルの作り方

解析ファイルの作成にはいったんTASSELを挟むとうまくいった。

1.PLINKをインストールして任意のフォルダに展開。

2.手持ちのジェノタイプデータをエクセル等で作成してhapmap形式(.txt)で保存。形式はこんな感じ

AlleleID alleles chrom pos strand assembly# center protLSID assayLSID panelLSID QCcode ind.1 ind.2 ind.3
marker1 G/A 1 1 NA NA NA NA NA NA NA GG GG GG
marker2 C/T 1 2 NA NA NA NA NA NA NA CC CC CC
marker3 A/G 1 3 NA NA NA NA NA NA NA AA AA AA
marker4 G/C 1 4 NA NA NA NA NA NA NA GG GG GG

3.TASSELで保存したファイルを開き、「Save as...」からPLINK形式を選んで保存。.plkファイルができていることを確認。

4.PLINKを起動し、以下のコマンドを打つ

plink --file 入力ファイル名.plk --recode12 --out 出力ファイル名

5.出力したmapファイルと.pedファイルを、Admixtureを展開したのと同じフォルダに入れれば準備完了

解析

解析の流れとしては、

  1. 適切な祖先集団数(K)の推定:マニュアル2.1章の通り

  2. 指定したKの値でAdmixture解析、Q値出力

  3. 各個体のQ値をプロット:Rのbarplotが一般的

という感じ。具体的なコマンドは以下のブログに詳しく解説されている。 https://indo-european.eu/human-ancestry/admixture-ancestry-components-r-statistics-plink-convertf-bed-and-ped-files/

ついでながら以下にggplot2でbarplotを書く場合のスクリプトを置いておく。ご参考まで。

#K=2
q.k2 <- read_csv(file="acc_name.csv") %>>% # read accessions name
  cbind(read_table(file="acc.plk_fmt.2.Q", col_names = c("q1", "q2"))) %>>% # add Q value for each accession
  dplyr::arrange(q1, q2, species)

#K=3
q.k3 <- read_csv(file="acc_name.csv") %>>%
    cbind(read_table(file="acc.plk_fmt.3.Q", col_names = c("q1", "q2", "q3"))) %>>%
    #note: round the values to properly order the accessions
    dplyr::mutate_if(is.numeric, funs(round(., digits=5))) %>>%
    dplyr::arrange(species, q1, q2, q3)  

#barplot by ggplot2
#K=2
    qg.k2 <- q.k2 %>>%
      tidyr::gather(key = "Q_number", value = "Q_portion", -c("species", "accession")) %>>%
      dplyr::mutate(K_value = "K = 2")

    #K=3
    qg.k3 <-  q.k3 %>>%
      tidyr::gather(key = "Q_number", value = "Q_portion", -c("species", "accession")) %>>%
      dplyr::mutate(K_value = "K = 3")

    qg <- rbind(qg.k2, qg.k3)

    qg %>>%
      ggplot(aes(y = Q_portion, x = accession, fill = Q_number)) +
      geom_bar(stat="identity") +
      labs(x="", y="Q value")+
      theme(
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1),
        legend.position = "left") +
      scale_x_discrete(limits = lab.order) +
      scale_fill_manual(values = c("blue3", "forestgreen", "darkorange3"),
                        name = "")+
      facet_wrap(~K_value, nrow = 2, strip.position = "left")

プロダクションマッチフェス走り方考察

(2020年3月22日:誤字を色々修正)

担当の白菊ほたる上位フェスやら副業やらで忙しくて気がついたら総選挙結果発表も終わっていた。 ちゃんみお、加蓮、りあむちゃん、おめでとうございます。

さて、前回のフェスではかつての雪辱を晴らすため、自分なりに本気で走り方を考えて取り組んだ。初めて傭兵にも出た。 結果、奇跡的に2000位に食い込むことができた。その記念というか記録として、1枚取りを想定した走り方の考察を書き残しておく。

攻編成の組み方

フェス用編成の組み方については、すでに秀逸な記事がいくつもあるので詳しくは書かない。 要点だけまとめると、

  • 攻フロントメンバーの属性を揃える(1色か2色)
  • 揃えた色の属性の攻をアップする特技持ちメンバーをできるだけ加える
  • 攻ステータスは高いに越したことはない。守は気にしない

というぐらい。フロント5人にはカードのコストでいうと24以上を揃えるのが望ましい(上位報酬狙うなら必須)。

走り方について

傭兵に出るべきか

傭兵とは、フェス期間のみ一時的に別のプロダクションで戦うこと。 発揮値1%の差が勝敗を分けることもあるので、フェストレーニングルームはLv.7、参加者10人以上(コンボボーナス確保のため)のプロダクションに所属しておくべき。 自分の今いるプロダクションがこれを満たしていないなら、Twitterなどで検索して募集広告が出ているプロダクションにお願いしてみよう。 2000位を目指すつもりであれば断られることはないように思う。

どのラウンドを走るか

基本的に昼間に開催される試合のほうが参加人数が少ないので、一戦ランキング上位報酬が得やすい。 そのかわりに、プロダクションの人数が少ないとコンボやハイテンションモードにもなりにくいので、その辺はバランスを見極める必要がある。

投入コストについて

フェスでは通常の対戦とは異なり、20%, 40%, 50%, 60%, 80%, 100%の中から選ぶことができる。 通説では40%で叩き続けるのが最も効率(投入コストあたりの獲得ポイント)が高くセオリーとされている

とはいえ、自分で実際に調べたことはなかったので、この機会に集計をとってみることにした。 20%は明らかに効率が悪いので除外し、40%, 50%, 60%で集計すると以下のようになった (自攻コスト480、攻発揮値の合計約61万、ルームLv.7、コンボ・ハイテンションモードMax、対戦相手は1桁台の上位プロ)

投入コスト割合 試合数 勝利数 平均獲得ポイント コスト1%あたりの平均獲得ポイント
40% 33 31 3737 93.4
50% 22 22 4575 91.5
60% 7 7 4806 80.1

というわけで、通説を覆すことはできず、やはり40%が一番効率が良いようだ。 よって、40%選択→ブラウザバックの繰り返しをひたすら続けていれば問題ない。

ところで、この数字は負け試合を含めた平均のため、バトルに敗北しなければ更に効率は上がる。 もし時間に余裕があるならば、相手の守コストを逐一確認して、負けそうならば多めに攻コストを投入することで更に最適化できると思われる。

コンボ数とポイント加算

上位プロダクション同士の対戦ならばコンボ数は数万に達することもあるが、流石にこれに比例してボーナスが加算されるわけではなく、 コンボボーナスはコンボ60回でMaxとなる。 コンボ数2000ぐらいのときに私が計算した結果、獲得ポイントに占めるコンボボーナスの割合は平均約52%だった。

消費したアイテムと時間

イベント開始前にエナドリ換算でおよそ3050本持っていた。(マイエナ800本、ハーフ2000本、テンチャ10本、生エナジードリンク150本) それがイベント終了時に残ったのは生エナジードリンク100本ほど...2000位ボーダーであった2150万ポイントのためにエナドリ1900本以上使ったことになる。 他の方がどの程度使うのか分からないが、上位常連でもなければここまで貯めるのは本当に苦労する。エナドリは他のイベントで使ってしまいがちなので、節約に気をつけなくてはいけない(戒め)。

かかった時間は、5日間合わせて5時間程度。ボーダーが平均の2倍程度に上がった今回でもこの程度なので、他のイベントに比べれば時間効率はかなりいい方だと思う。

まとめ

以上を踏まえて走り方の流れをまとめると、 1. 攻フロントメンバー強化 2. 今のプロダクションの戦力確認。不安ならば傭兵に出る。1枚取りを確約すれば受け入れ先はある可能性は高い 3. ボーダーを予想してスケジュールを立てる 4. 40%選択→ブラウザバックの繰り返し

上では書かなかったが、何気にボーダー予測が重要でこれを誤るとすべてが水泡に帰すことになる。 イベントの開催時期や報酬によって相当に変動するので、とある弓兵さんのブログimcgborderで下調べは欠かせない。 ツイッターなどでボーダー予測している方もいるが盲信は禁物だ。今回のボーダーも誰が予想したか...

発足から7年以上たっても何が起こるかわからない、それがモバマス

以上、わずかでも参考になれば幸いです。

Rのデータフレームで検索と置換をする方法

(2019年12月26日追記) 普段からエクセルを使い慣れている人にとっては、「検索と置換」はよく使う機能だと思う。 Rのデータフレームで同様の処理をする機会があり、素人なりに方法を調べたのでまとめてみた。

データフレームの準備

適当なデータでデータフレームを作成。ベクトルを合成してデータフレームを作ると、デフォルトでは文字列がfactor型になってしまうので、引数でstringsAsFactors = FALSEを指定する。

x <- c("abc","cde","efg","ghi","ijk","klm")
x2 <- c("Alice","Bob","Charles","Deven","Eve","Feldman")
y <- c(TRUE,FALSE,TRUE,TRUE,TRUE,FALSE)
z <- c(1,3,5,2,4,7)
df<- data.frame(X=x, X2=x2, Y=y, Z=z, stringsAsFactors = FALSE)
df

    X      X2     Y Z
1 abc   Alice  TRUE 1
2 cde     Bob FALSE 3
3 efg Charles  TRUE 5
4 ghi   Deven  TRUE 2
5 ijk     Eve  TRUE 4
6 klm Feldman FALSE 7

列を指定して検索&置換

まずは、データフレーム中の特定の列だけ検索対象にしたいときの方法。試しにX2列を対象に、「e」を「***」に置換してみる。

#Rの基本パッケージを使う方法
df1 <- df
df1["X2"] <- lapply(df1["X2"], gsub, pattern="e", replacement = "_")
df1

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cde     Bob FALSE 3
3 efg Charl_s  TRUE 5
4 ghi   D_v_n  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7

#dplyrを使う方法
library(dplyr)
library(stringr)
df2 <- dplyr::mutate(df,X2=gsub(X2,pattern="e",replacement = "_", ignore.case = TRUE))
df2
 
    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cde     Bob FALSE 3
3 efg Charl_s  TRUE 5
4 ghi   D_v_n  TRUE 2
5 ijk     _v_  TRUE 4
6 klm F_ldman FALSE 7


df3 <- dplyr::mutate(df,X2=gsub(X2,pattern="e",replacement = "_", ignore.case = FALSE))
df3

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cde     Bob FALSE 3
3 efg Charl_s  TRUE 5
4 ghi   D_v_n  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7


df4 <- dplyr::mutate(df,X2=str_replace(X2,pattern="e",replacement = "_"))
df4

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cde     Bob FALSE 3
3 efg Charl_s  TRUE 5
4 ghi   D_ven  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7

gsubでは、ignore.case = TRUEを指定することで大文字小文字区別なく置換される。 ここで、gsubstr_replaceでは置換のされ方が違うことに注意(df3とdf4)。前者では要素内でヒットするすべての文字列が置換されるが、後者では最初のヒットのみ置換される。str_replace_allを使うとgsubと同様に全てのヒットが置換される。

また、名前付きベクトルをpatternに指定すれば複数の変換を一度にできる

df4.5 <- dplyr::mutate(df, X2 = str_replace_all(X2, pattern = c("e" = "_", "a" = "_2")))
df4.5

    X       X2     Y Z
1 abc    Alic_  TRUE 1
2 cde      Bob FALSE 3
3 efg Ch_2rl_s  TRUE 5
4 ghi    D_v_n  TRUE 2
5 ijk      Ev_  TRUE 4
6 klm F_ldm_2n FALSE 7

データフレーム全体を検索&置換

次に、データフレーム全体を検索対象にしたいときの方法。

#Rの基本パッケージを使う方法
df5 <- data.frame(lapply(df, function(x){
gsub(pattern="e", replacement = "_", x)
}),stringsAsFactors = FALSE)
df5

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cd_     Bob FALSE 3
3 _fg Charl_s  TRUE 5
4 ghi   D_v_n  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7


#dplyrを使う方法
df6 <- dplyr::mutate_all(df, ~gsub(.,pattern="e",replacement = "_"))
df6

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cd_     Bob FALSE 3
3 _fg Charl_s  TRUE 5
4 ghi   D_v_n  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7

df7 <- dplyr::mutate_all(df, ~str_replace(.,pattern="e",replacement = "_"))
df7

    X      X2     Y Z
1 abc   Alic_  TRUE 1
2 cd_     Bob FALSE 3
3 _fg Charl_s  TRUE 5
4 ghi   D_ven  TRUE 2
5 ijk     Ev_  TRUE 4
6 klm F_ldman FALSE 7

Rの基本パッケージを使う方法でdata.frameをかぶせているのは、lapplyで返されるのがリスト型なので、データフレームに直している。

注意点

データフレーム全体を検索対象にすると、factor型やdouble型の列もすべてcharacter型にされてしまう。

> mode(df$Y)
[1] "logical"
> mode(df$Z)
[1] "numeric"
> mode(df5$Y)
[1] "character"
> mode(df6$Z)
[1] "character"
> mode(df5$Y)
[1] "character"
> mode(df6$Z)
[1] "character"

よって必要に応じて型の変換を行う。

df3$Y = as.logical(df6$Y)
df3$Z = as.numeric(df6$Z)

参考サイト様

r - Replace all occurrences of a string in a data frame - Stack Overflow

dplyrのmutate_if()とかについて - Technically, technophobic.

{dplyr}のパイプの中でNAを0に置換する - Note of Pediatric Surgery

stringr: Rの文字列をまともな方法で処理する - Heavy Watal

https://stringr.tidyverse.org/reference/str_replace.html

アイドルチャレンジ大運動会編(2017年10月上位報酬白菊ほたる)参加記録

アイドルチャレンジ白菊ほたるイベのまとめ

アイチャレ史に残る激しい戦いだった。最終2000位ボーダーは34,169,844pt。一番最近の8R制アイチャレ(ベーカリー編)が2200万ptぐらいだったことを考えると、1.5倍以上に跳ね上がった。

理由は…白菊ほたる担当Pとしてはほたる人気のおかげと思いたいけどどうだろう。リミテッドガチャや火曜シンデレラのキャンペーンでCPブレッドがたくさん供給されたことが一因かもしれない。

ちなみに自分は約2900万ptで2172位。見積もりが甘かったなーと悔しい思いをしたので反省も込めて少し記録を書いておく。

 

効率の良い走り方について

効率良く走るためには、とにかくスコアの高いレッスンチャレンジレッスンだけ行うということに尽きる。下にスコアあたりの獲得ポイントの比をまとめてみた。

f:id:Fujinitaka:20171106011942p:plain

スコア3300万のチャレンジレッスンは350万の時よりも3倍近く効率が良い。レッスンを行うなら3300万だけにするか、せめて青で網掛けしたものだけに絞るべきかと。あとはCP1//6の3の倍数でちょうど倒せるようなレッスンを優先して行うようにしたい。

あと、時間とスタドリに余裕があるなら本気ゲージMAXのときだけレッスンする方が良い。インターバル期間はストレスなく本気ゲージが上がるので、この期間に3人全員マックスにしておくと結構楽だった。

ところで960万のチャレンジレッスンが1000万のよりも効率がいいようだ。確かこれはほたるのレッスンだったと思うので、上位報酬キャラのレッスンが一番効率良いのかもしれない。

 

消費アイテム数

微課金で弱フロント(開始時攻604063、守557127)だったので、アイテムは大量に必要になった。

開始時がCP換算で50本程度、ストアで購入したCPが105本、ガチャで手に入れたのがCP換算190本程度。ほぼ全て使い切ったので、メダルで手に入れたものを除いて345本も消費していた。

マイエナハーフは600個、マイスタは500個ぐらい消費していた。

 

メダルチャレンジやガチャの結果

メダルチャレンジでは、海さん、むつみちゃんのSR2枚とも揃えることが出来た。

リミテッドガチャは、ガチャチケ300枚使ったおかげで、ジュエル30000個使ってSR6枚全員取ることが出来た。さらにイベント後に7枚イベント専用アイドルチケが貯まったので、これで運良くSR幸子を引くことが出来た。

正直今までほぼ無課金でガチャの仕組みとか分かってなかったので、リミテでガチャチケ使えばSRがこれだけ取れるのは衝撃だった。

 

感想と反省

今回始めて走ってみましたが、最後のボーダーの伸びが想像以上で、モバPを甘く見すぎていた…というのが一番の反省点。あとリミテッドガチャを引いたのがイベント最後の方だったので、もっと早く引いてればパワー持ちSRをフリトレで高値取引できたはず。

それでも、イベント終了後には、リミテで引いたSRを元手にフリトレで上位ほたるをお迎えできたので(スタドリ190程度だったのでSR幸子の半分くらい)、この辺一種のお金かけた人に対する救済措置なのかなーとも思った。他のイベントではどうなのか、これからはもっと気を配るようにしたい。

 

このイベントをきっかけにほたる人気がもっと上がると良いな…今まで記憶にある限り総選挙50位以内入ってないし…最近は、関裕美、森久保乃々とワンステップスというユニットを組んでるので、このユニットきっかけに声がつくストーリーを妄想している。近ごろ森久保も貫禄が出てきたし、新人の2人を引っ張る役を任されて…とか考えるとすごい良い話になりそう。

igraphでシンデレラガールズのアイドル紹介相関図

アイドルマスターシンデレラガールズ第6回総選挙が終わりましたね。楓さんおめでとうございます。

さて、今回の総選挙では、登場アイドルがお題に沿って他のアイドルを紹介する特設ページ開設され、話題になった。

www.nicovideo.jp

全アイドルを網羅する内容の丁寧さもさることながら、自分としては『おっ、これはいい有向重みなしグラフ…』と思ったので、Rのigraphパッケージ*1を使って相関図を作成してみた。

まずは、以下のように1列目に紹介する側、2列めに紹介される側のアイドルをタブ区切りで記入したファイルを作成し、エッジリストを作成。

卯月 まゆ
卯月 加蓮
卯月 美嘉
卯月 みく
卯月 雅
有香 亜季
有香 拓海
有香 早苗
有香 アヤ
有香 あやめ
ゆかり 星花
ゆかり 久美子

(以下略)

これをigraphに取り込んで、以下のスクリプトで相関図が描ける。

library(igraph)#igraphパッケージの呼び出し
cg <- read.table(file="cg_idol_edgelist.txt", header=FALSE)#エッジリストを読み込み
g <- graph.data.frame(cg[1:2], directed=T)#グラフオブジェクトに渡す
V(g)$color <- "pink"#ノードの色をcute仕様に
V(g)$size <- 6#ノードサイズの決定
V(g)$color[63:125] <- "cornflowerblue"#ノードの色をcool仕様に
V(g)$color[126:183] <- "yellow"#ノードの色をpassion仕様に
V(g)$label.cex <- 1.2#文字サイズ変更
png("cg_idol.png",width=2000,height=2000)pngで保存する準備
plot(g,layout=layout.kamada.kawai)#グラフをプロット
dev.off()#デバイスを閉じる

f:id:Fujinitaka:20170524023745j:plain

…いやー非常にごちゃごちゃしていて分かりにくい。

ただ、クラスタを形成しているように見えないことから、少なくとも全アイドルまんべんなく紹介されていることが示唆されている。また、よく見ると、左の方にアヤさん、渚さん、櫂くん、洋子さんが配置されていて、スポーツ系が集まっていることや、右端には蘭子、飛鳥くん、聖ちゃん、音羽さん、詩織さん、小梅ちゃんが位置してるのは、魂の共鳴を感じるCool不思議組のクラスタっぽいかな-とは思える。

全体的な傾向を見ると、例えば平均アイドル紹介数は以下で求まる。

gsize(g)/gorder(g)#エッジ数/ノード数
[1] 4.945355

任意のノード間の平均パス数を求めてみると、

mean_distance(g, directed = TRUE)

[1] 3.344262

となる。つまり、ランダムに1組のペアを選んだとき、一方のアイドルからもう一方のアイドルまでは平均約3.3回のクリックでたどり着ける計算になる。おもったより短くて済むという印象。

他のアイドルへの距離(最低必要クリック数)を求めて、テキストに一覧を出力するには、

dis <- distances(g, mode="out")#outでノードから出ていくエッジを数える
write.table(dis,"cg_dis_table.txt",append=T,quote=F,col.names=T)

とすればよい。

 

次に、マイ担当アイドルである「白菊ほたる」に注目して、周囲の関係性を調べてみた。まず、ほたるが紹介する・ほたるを紹介するアイドルを求める。

incident_edges(g,58,mode="all")
$ほたる
+ 9/905 edges (vertex names):
[1] ほたる->ネネ ほたる->朋 ほたる->晴 ほたる->海
[5] ほたる->巴 加奈 ->ほたる 裕美 ->ほたる 惠 ->ほたる
[9] 智香 ->ほたる

これで、ほたるの周囲の関係性は把握できた。惠さんと智香から矢印出てるのは意外な気がする。

つづいて、ほたるが所属するユニットであるGBNS(ガールズビーネクストステップ)のメンバーには、どのように到達できるか調べてみた。

shortest_paths(g,"ほたる","千鶴",mode="out")
$vpath
$vpath[[1]
+ 5/183 vertices, named:
[1] ほたる 海 晶葉 亜子 千鶴

shortest_paths(g,"ほたる","裕美",mode="out")
$vpath
$vpath[[1]
+ 5/183 vertices, named:
[1] ほたる 巴 夏樹 涼 裕美

shortest_paths(g,"ほたる","泰葉",mode="out")
$vpath
$vpath[[1]
+ 4/183 vertices, named:
[1] ほたる ネネ 薫 泰葉

という感じで、最短ルートを出すこともできる。こう見ると、アイドル全員ちゃんと繋がってるんだな…という温かみのようなものを感じる。

 

このグラフ、更に掘り下げて次数中心性など調べてみると面白そうだけど、今は時間と知識がちょっと足りないので別の機会に。とはいえこのエッジリストがあれば、誰でも同様の解析できるので、どこかのサイトにアップしようかと思う(果たして需要はあるのか?)。

 

8/19追記:

Dropboxにてエッジリストを公開しました。

Dropbox - cg_idol_edgelist.txt

 

参考サイト様:

第6回シンデレラガール総選挙のアイドル紹介をjson化 · isaisstillalive/imas_cg_hash@384e2a2 · GitHub

注:メアリーと美世さんだけ抜けている

igraph R manual pages

農薬の散布間隔図をRのigraphパッケージでプロットしてみた - もうカツ丼でいいよな

R seminar on igraph

先行文献:

【モバマス】第六回シンデレラガールズ総選挙「アイドル紹介」 アイドル毎に選ばれた人数まとめ:もばます!