BioinfoにおけるPandas,Matplotlibの基礎
この記事は、Pythonで実践 生命科学データの機械学習 (https://www.yodosha.co.jp/yodobook/book/9784758122634/) の内容を含んでいます。
私は現在バイオインフォマティクス研究室に所属する学生です。
勉強した事をアウトプットする場として用いていますため、何卒ご理解のほどよろしくお願いします
(＞人＜;)
機械学習の一連の流れ
機械学習を行う前には学習が効率に行われるようにデータを"いい感じに"処理する必要がある。
例を挙げると、欠損値の処理、Objectの対応、重複の処理などが挙げられる。
必要に応じて、データの可視化を行い、機械学習に持ち込むといった流れである。
ライブラリ
scikit-learnは機械学習ライブラリの一つである。アルゴリズムが予め実装されている。
PyTorchは深層学習ライブラリである。
ワークフロー
基本的には、以下の順番で機械学習モデルを実装する。
1.データを入手
2.データの概要理解 可視化を行う
3.データの前処理
4.機械学習モデルの選定、訓練
5.モデルの決定、微調整
6.テストデータにモデルを適応、評価をする
ここではデータの前処理に重点をおいて記事を書いていく
データの前処理
基本的にデータの前処理で扱うパッケージは以下の通りである。
GEOpraseはGoogle Colaboratoryにはデフォルトで入っていない為、DLが必要である
※ GEO（Gene Expression Omnibus）データベースからデータを取得・解析するためのPythonライブラリである。
# パッケージのインポートとDL
!pip install GEOparse
%matplotlib inline
import GEOparse
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
以下に重要だと思った操作を挙げていく
pd.set_option('display.max_columns', 50)
pandasの表示設定を変更し、データフレームを表示するときに最大50列まで表示可能にする。
デフォルトでは表示できる列数に制限があるため、列が多いデータフレームを確認する際に便利
data_table = gse.gpls['GEOデータセット名'].table
gse.gpls はGEOデータセット（データ名）のプラットフォームデータを取得する。
このデータにはプローブ情報（遺伝子ID、シーケンス情報など）が含まれる。
data_table.info()
#出力結果の例
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 48762 entries, 0 to 48761
Data columns (total 10 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 ID 48762 non-null object
1 Gene Symbol 47500 non-null object
2 Gene Title 47500 non-null object
3 Organism 48762 non-null object
4 Chromosome 46000 non-null object
5 Start 48762 non-null int64
6 End 48762 non-null int64
7 Strand 48762 non-null object
8 Sequence 48762 non-null object
9 Source 48762 non-null object
dtypes: int64(2), object(8)
memory usage: 3.7+ MB
info() の主な表示内容
データフレームのクラスとインデックス情報
データフレームの行数、インデックスの範囲が表示される。
各列の情報
列名
列に含まれる非欠損値の数（Non-Null Count）
各列のデータ型（dtype）
全体のメモリ使用量
データフレーム全体が使用しているメモリ量を表示
欠損値を確認するときは、以下のコードを使うと便利
#各列の欠損値の数を確認するには以下を使用
print(data_table.isnull().sum())
#データ型のみを確認する場合
print(data_table.dtypes)
例えば、ある行で特定の文字型のみを抽出したいとする。
data_table_hs = data_table[data_table['列の名前'] == '抽出したい文字列']
data_table_hs
このようにしたら抽出ができる
特定の文字型ではなく、NaNを削除するときは、
data_table_hs_sym = data_table_hs.dropna(subset=['列の名前'])
data_table_hs_sym
['列の名前']列の欠損値をdropnaで除去
データ処理した列は、_列の略称を末尾に付けるとおすすめ
文字列の最初の三文字の内訳をカウントしたいときは、
print(data_table_hs_sym['列の名前'].str[:3].value_counts())
この時に、ある最初の三文字のみを抽出したい場合は、
data_table_hs_sym_NM = data_table_hs_sym[data_table_hs_sym['列の名前'].str[:3] == '抽出したい三文字']
data_table_hs_sym_NM
これで可能になる。
ある列での文字列ごとのカウント数が欲しい場合は、
data_phenotype['列の名前'].value_counts()
これで、同じ文字列がいくつあるのか確認できる。
その際に、欲しい二つの文字型データのみ欲しいという場合は、
normal_ID = data_phenotype[
data_phenotype['列の名前'] == '抽出したい文字列'
].index
tumor_ID = data_phenotype[
gse_phenotype['列の名前'] == '抽出したい文字列'
].index
このようにそれぞれ分類することも可能である。
遺伝子発現量Dataの可視化
ここで重要なのは可視化である。データ数が膨大な上に直感的にわかりにくい為に予想外の欠損値や記述ミスに気づきにくい。そのような事態を防ぐために可視化を行い、データの概略を掴むことが重要である。
GEOparse ライブラリを使用して GEO データセット内の発現量データ（例: VALUE 列）をサンプルごとにピボット（再構築）し、1つのデータフレームにまとめる操作をう。
# pivot_samplesで発現量をまとめる
pivoted_samples = gse.pivot_samples('VALUE')
pivoted_samples
欠損値カウントする時のコードをもう一度まとめる。
#列の欠損値の数をカウントする
print(pivoted_samples.isnull().sum())
# value_countsでpivoted_samples.isnull().sum()にどのような値があったか確認する
print(pivoted_samples.isnull().sum().value_counts())
#行の欠損値の数をカウントする
pivoted_samples.isnull().sum(axis=1).value_counts()
#欠損値を含む行だけ抽出する
pivoted_samples[pivoted_samples.isnull().any(axis=1)]
このようにして欠損値の数を把握することができる。
欠損値処理をする前に、データの群の特性を理解する方法も記す。
必ずといっていいほど、最初には基本統計量を計算する。
pivoted_samples.describe()
箱ひげ図をさまざまな書き方をここでまとめて記す。
箱ひげ図一覧
(上に空行が必要)
# 全サンプルを図示すると見にくいため，5サンプルのみ
#boxplotを用いた。
pivoted_samples.iloc[:, 0:5].boxplot()
以下はSeabornライブラリで描いた箱ひげ図である。
# Seabornで描いた箱ひげ図
# Seabornはsnsという形でインポートすることが多い (import seaborn as sns)
sns.boxplot(data=pivoted_samples.iloc[:, 0:5])
データの全体傾向をより知りたい場合は、violinplotを使うといい。
# violinplotは箱ひげ図よりもデータの全体傾向を把握しやすいというメリットがある.
# 「ヴァイオリン」の縁の部分はカーネル密度推定で描かれている
sns.violinplot(data=pivoted_samples.iloc[:, 0:5])
大規模データセット用の箱ひげ図として用いられるのは、boxenplotである。
# boxenplotは大規模なデータセット用の箱ひげ図と言える.
# 四分位しかわからない箱ひげ図と比べ，細かい分位をも理解できる
sns.boxenplot(data=pivoted_samples.iloc[:, 0:5])
また、displotでヒストグラムとカーネル密度推定を同時に見ることができる。
# pivoted_samples.iloc[:, 0]は少し見にくいのでpivoted_samples.iloc[:, 1]を見る
current_palette = sns.color_palette(n_colors=24) # seabornでの色の指定を行うため
sns.displot(
pivoted_samples.iloc[:, 1],
kde=True,
facecolor=current_palette[0], # 青色がヒストグラム
color=current_palette[1],
) # オレンジの線がカーネル密度推定: ヴァイオリンプロットと同じもの
実行結果
また、散布図の行列の形式で出力する場合、pairplotがいい。
sns.pairplot(pivoted_samples.iloc[:, 0:5], height=2.5)
実行結果
各サンプルの相関関係を計算し、まとめたい時はヒートマップがおすすめである。
sns.heatmap(pivoted_samples.iloc[:, 0:5].corr(), square=True, annot=True)
実行結果
これらの結果は冒頭に記した書籍のデータを用いて作成したものである。是非とも購入して、実践的な解析の流れを掴んでほしい。
遺伝子情報の調べ方
PubMed で原著論文を検索し、遺伝子の役割や関連研究を確認。
多くの遺伝子は肝細胞がんのマーカーや予後因子として研究されているが、解明されていないものも多い。
NCBI Gene を活用して以下を確認：
遺伝子機能の要約
ゲノム情報
組織ごとの発現プロファイル