<jupyter_start><jupyter_text>マイクロアレイデータ解析 \~2021 Spring Semester\~- 程研HP: http://ui-tei.rnai.jp/- 実習wiki: http://ui-tei.rnai.jp/microarray/doku.php?id=2021今回の実習では、**「マイクロアレイを用いた網羅的遺伝子発現解析」**というタイトルで、網羅的に遺伝子の発現を解析する手法であるマイクロアレイの「原理」「操作手順」「解析方法」について、- wet(実験)：siRNAを導入した細胞からRNAを抽出し、マイクロアレイを行う。- dry(解析)：全mRNAの変動量を、マイクロアレイデータの特徴を踏まえて解析する。の両側面から学んでいただきますが、このNotebookでは、**dry(解析)** パートに関して、プログラミング言語: [Python](https://www.python.org/)を使って実際に手を動かしながら学んでいきます。  初めに今回dryパートの実習を担当させていただく、程研M1の岩崎（生情11期）です。よろしくお願いします。この実習が終わった時には、**「マイクロアレイ解析結構Pythonでできるし、あんなことやこんなことも自分でプログラムを書いて調べてみたい！！」** と思っていただけることを目指して頑張りますので、わからないことがあればなんでも遠慮なく質問してください！！また、**「解析手法に絶対はない」** ので、是非積極的に色々な解析を試みて、面白い結果を出してきてください！！！（ぜひ一緒に議論しましょう！）（※ このNotebookでは、**「全員が自分一人で解析できること」**を目指しているため、プログラミングが得意な人や、解析経験のある方からすると簡単すぎるかもしれません。その場合は説明を聞かなくて構いませんので、ご自身でどしどし進めてください！！） ***  0. 環境構築　それでは解析を始めていきましょう！！解析に必要なデータや、各種ツールを揃えていきます。<jupyter_code>! pip install "git+https://github.com/iwasakishuto/TeiLab-BasicLaboratoryWork-in-LifeScienceExperiments.git" --ignore-requires-python<jupyter_output><empty_output><jupyter_text>ここで```ERROR: XXX has requirement YYY== but you'll have YYY  which is incompatible.WARNING: The following packages were previously imported in this runtime```のようなエラーが出る分には問題ありません。以下のコマンドでエラーが出なければ準備はOKです！<jupyter_code>from teilab.question import ask
ret = ask(
    text="好きな言葉を入力してください。", 
    username="あなたの名前", 
    icon_emoji=":grinning:", 
    icon_url=None, 
    webhook_url=None,
)<jupyter_output><empty_output><jupyter_text>1. データの準備　続いて、データの準備に取り掛かります。先ほどインストールしたパッケージを用いてデータのダウンロードを行います。<jupyter_code>import numpy as np
import pandas as pd
password1 = ""
password2 = ""
from teilab.datasets import TeiLabDataSets
dataset = TeiLabDataSets()
path1 = dataset.get_data(password=password1)
path2 = dataset.get_data(password=password2)
print("Data1", path1)
print("Data2", path2)<jupyter_output><empty_output><jupyter_text>2. データの読み込み ```python 本来はこんな感じで読み込みをする。import osimport pandas as pdfrom teilab.utils._path import DATA_DIR Data (Password1)dirname1 = os.path.join(DATA_DIR, password1, "実習解析用データ")print("Data1 is @", dirname1)print(sorted(os.listdir(dirname1)))df1_1_1 = pd.read_csv(os.path.join(dirname1, 'US91503671_253949442637_S01_GE1_105_Dec08_1_1.txt'), sep="\t", header=9) Data (Password2)dirname2 = os.path.join(DATA_DIR, password2)print("Data2 is @", dirname2)print(sorted(os.listdir(dirname2)))df2_1_1 = pd.read_csv(os.path.join(dirname2, 'SG19378659_257236339458_S001_GE1_1200_Jun14_1_1.txt'), sep="\t", header=9)```` ```sh データの中身が知りたい方はこのコマンド! head -n20 {dataset.filePaths[0]}```<jupyter_code>dataset.samples.show_groups()
# 以下のコードで、楽に読み込みできる
df1_1_1 = dataset.read_data(no=5)
df2_1_1 = dataset.read_data(no=0)
# このノートブックで扱うgroupのnumberリストを取得しておきます。
group_numbers = dataset.samples.get_group_numbers(group_no=1)
conditions = dataset.samples.Condition[group_numbers]
for gn,cnd in zip(group_numbers, conditions):
  print(gn, cnd)<jupyter_output><empty_output><jupyter_text>3. データの統合 　それでは、昨日のExcelを使った解析と同様に、全サンプルのデータを一つにまとめていきましょう。　行番号とプローブ番号の対応関係は（同じタイミングの実験であれば）どのサンプルも同じであるので、サンプル $X$ のデータの隣にサンプル $Y$ のデータを concatenate すれば1枚のワークシート（DataFrame）にまとまります。　そこで、これを繰り返して、全サンプルの `gProcessedSignal` の値を1つのテーブルにまとめていきましょう。  3.1 アノテーションデータの読み込み　まず、各プローブのアノテーションデータ（どういうプローブかの説明。以下のカラムの情報が大事）を取り出します。なお、このデータはサンプルに寄らないので、8つのサンプルのうち一つから一度だけ取り出します。<jupyter_code>USE_COLS_ANNO  = [
    "FeatureNum", "ControlType", "ProbeName", "GeneName", "SystematicName"
]<jupyter_output><empty_output><jupyter_text>|column name|description||:-:|:-||`FeatureNum`|スポットの番号||`ControlType`|positive controlは `1` negative controlは `-1`それ以外（解析で用いる）は `0`||`ProbeName`|プローブ名||`GeneName`|遺伝子名||`SystematicName`|遺伝子名|<jupyter_code>df_anno = dataset.read_data(no=5, usecols=USE_COLS_ANNO)
df_anno.head(5)<jupyter_output><empty_output><jupyter_text>3.2 シグナル値の読み込み 続いて、各サンプルのシグナル強度( `gProcessedSignal` )のデータを取得します。なお、この時 `gIsWellAboveBG` が `0` のものは「（真の）シグナルがバックグラウンドのシグナルよりも低く、信頼できないデータである」ということを意味するため、取り除きます。<jupyter_code>USE_COLS_SYGNAL = [
    "gProcessedSignal", "gIsWellAboveBG"
]<jupyter_output><empty_output><jupyter_text>|column name|description||:-:|:-||`gProcessedSignal`|green(Cy-3)のシグナル強度（＝発現量）||`gIsWellAboveBG`|（真の）シグナルがバックグラウンドのシグナルより十分高いか？（＝信頼できるデータか）|<jupyter_code>df_combined = df_anno.copy(deep=True)
index = set(df_combined.index)
print(f"データ数(before): {len(df_combined)}")

for gn,cnd in zip(group_numbers, conditions):
  df_signal = dataset.read_data(no=gn, usecols=USE_COLS_SYGNAL)
  index = index & set(df_signal[(df_signal.gIsWellAboveBG==1)].index)
  df_combined = pd.concat([df_combined, df_signal[["gProcessedSignal"]].rename(columns={"gProcessedSignal" : cnd})], axis=1)

df_filtered_1 = df_combined.loc[index, :]
print(f"データ数(after) : {len(df_filtered_1)}")<jupyter_output><empty_output><jupyter_text>また、`ControlType` の値が $\pm1$ のものはコントロールであるため、`0` のもののみ取り出します。<jupyter_code>print("データ数(before):", len(df_filtered_1))
df_filtered_2 = df_filtered_1[df_filtered_1.ControlType == 0]
print("データ数(after) :", len(df_filtered_2))
# インデックスを振り直す。
df_filtered = df_filtered_2.reset_index(drop=True)
df_filtered.head(5)<jupyter_output><empty_output><jupyter_text>```python データをGoogleDriveに保存したい場合は、以下のコードを走らせてください。from google.colab import drivedrive.mount('/content/drive')df_filtered.to_excel("microarray_filtered.xlsx", index=False)```  4. データの前処理　無事にデータがダウンロードできたので、実験上のバイアス等を取り除くためにデータの前処理を行います。（ここでは、Summarizationのみを行います。）<jupyter_code>df_filtered.columns
# それぞれ単純に平均化を行う。
# df_filtered = df_filtered.groupby("GeneName").mean().reset_index()
df_filtered["mock"] = df_filtered[["mock(1)", "mock(2)"]].mean(axis=1)<jupyter_output><empty_output><jupyter_text>5. 解析 & 可視化ここでは、XYプロットとMAプロットを図示し、シグナル強度の分布を調べます。<jupyter_code>import plotly.express as px
import plotly.graph_objects as go
from plotly import offline
from plotly.subplots import make_subplots<jupyter_output><empty_output><jupyter_text>5.1 X-Yプロット - サンプル $X$ の `gProcessedSignal` の値を横軸- サンプル $Y$ の `gProcessedSignal` の値を縦軸にプロットしたものを **X-Yプロット** と呼びます。<jupyter_code>def XYplot(df, x, y, hover_name="", hover_data=[]):
    fig = px.scatter(
        df, x=x, y=y, hover_name=hover_name, hover_data=hover_data, 
        title=f"XY plot ({x} vs {y})"
    )
    return fig
fig = XYplot(df=df_filtered, x="mock", y="siVIM-270", hover_name="GeneName")
fig.show()<jupyter_output><empty_output><jupyter_text>> "vimentin" はどこにあるでしょうか？？  5.2 MAプロット- $log_2(Y/X)$ を縦軸 (Minus)- $log_{10}(XY)$ を横軸 (Average)にプロットしたものを **M-Aプロット** と呼びます。<jupyter_code>def MAplot(df, control, target, hover_data=[]):
  control_signals = df[control].values
  target_signals  = df[target].values
  X = np.log10(control_signals*target_signals)
  Y = np.log2(target_signals/control_signals)
  
  fig = go.Figure(data=go.Scatter(
      x=X, y=Y, hovertext=df[hover_data].values, mode="markers", marker_size=3
  ))
  fig.update_layout(
      title=f"MA plot ({control} vs {target})",
      xaxis_title = "$log_{10}" + f"({target} * {control})$", 
      yaxis_title = f"$log_2{target}/{control}$",
      
  )
  return fig
fig = MAplot(df=df_filtered, control="mock", target="siVIM-270", hover_data="GeneName")
fig.show()<jupyter_output><empty_output><jupyter_text>> - VIM はどこにあると予想できますか？？実際に確認してみてください！> - 縦軸が $0$ / $1$ / $-1$  とは何を意味しますか？？> - 横軸は何を意味しますか？？  6. 発展　ここからは、より良いデータ解析のためのいくつかの手法を紹介します。  6.1 データの正規化<jupyter_code>from scipy.stats.mstats import gmean
# 累積度数をプロットする。
def plotDensities(data, names, col=1, fig=None, title=""):
    fig = fig or make_subplots(rows=1, cols=1)
    for ith_data,name in zip(np.log2(data),names):
        hist, bin_edges = np.histogram(a=ith_data, bins=100, density=True)
        fig.add_trace(trace=go.Scatter(x=bin_edges[1:], y=hist, name=name, mode="lines"), row=1, col=col)
    fig.update_layout(
        title=title,
        xaxis_title = "$log_2(\\text{gProcessedSignal})$", 
        yaxis_title = "Density",
        
    )
    return fig
raw_data = df_filtered[conditions].T.values
print(f"raw_data.shape = {raw_data.shape}")
fig = plotDensities(raw_data, names=conditions, title="raw data.")
fig.show()<jupyter_output><empty_output><jupyter_text>6.1.1 75%tile最も単純な正規化手法1. 各サンプルごとに、発現量の小さい方から数えて順位75%に位置するものの値を求める。2. この75%tileの値は通常サンプルごとに異なるが、それらの（相乗）平均 `a` を求める。3. 各サンプルごとに、全プローブの値に 「 $a$ /そのサンプルにおける75%tileの値」をかける（つまり、全サンプルで75%tile値を $a$ に揃える。）<jupyter_code>def tile75_normalization(data):
    percentiles = np.percentile(a=data, q=75, axis=1)
    a = gmean(percentiles)
    return data*(a/percentiles)[:,np.newaxis]
tile75_data = tile75_normalization(raw_data)
print(f"tile75_data.shape = {tile75_data.shape}")
fig = plotDensities(tile75_data, names=conditions, title="75%tile")
fig.show()<jupyter_output><empty_output><jupyter_text>6.1.2 quantile法1. 各サンプルごとに発現量の値を順番に並べ替え、各順位の値をそれぞれ同順位のシグナル値の相乗平均で置き換える。2. その結果、全サンプルで分布が同一になる。<jupyter_code>def quantile_normalization(data):
    return gmean(a=np.sort(a=data, axis=1), axis=0)[np.argsort(np.argsort(data, axis=1), axis=1)]
quantiled_data = quantile_normalization(raw_data)
print(f"quantiled_data.shape = {quantiled_data.shape}")
fig = plotDensities(quantiled_data, names=conditions, title="quantile")
fig.show()<jupyter_output><empty_output><jupyter_text>6.1.3 正規化手法の比較各正規化処理後のデータを比較して見ましょう。<jupyter_code>fig = make_subplots(rows=1, cols=3)
plotDensities(raw_data,       names=sample_names, fig=fig, col=1)
plotDensities(tile75_data,    names=sample_names, fig=fig, col=2)
plotDensities(quantiled_data, names=sample_names, title="Comparison of normalization methods", fig=fig, col=3)
fig.show()<jupyter_output><empty_output><jupyter_text>6.2 seedマッチする遺伝子群の累積度数を調べる入力した配列を3'UTRにもつ遺伝子（アクセッション番号）のリストを表示するページ（[seedmatch](http://atlas.RNAi.jp/seedmatch/)）を用いて、siRNAのガイド鎖のseed（`UGAACUC`）と相補的な配列（`GAGTTCA`）が3'UTRに存在する遺伝子を検索する。  6.2.1 seed領域で調べてみる。まずは、siRNAのガイド鎖のseed（UGAACUC）と相補的な配列（GAGTTCA）が3'UTRに存在する遺伝子の発現量が本当に下がっているのか調べて見ましょう。つまり、オフターゲット効果を検証する、ということです。<jupyter_code>! gdown --id "0B1jjJ8m0VChCY28xRmhVcTRoamc" -O "seedmatch.txt"
df_matched_mRNAs = pd.read_csv("seedmatch.txt", sep="\t")
df_matched_mRNAs.columns = ["SystematicName", "NumHits"]
# seedmatchで検索したデータと紐付ける。
df_is_matched = pd.merge(df_filtered, df_matched_mRNAs, on="SystematicName", how="left").fillna(0)
# 累積度数曲線を描くために、ソートする。
df_is_matched["log2(RNA/mock)"] = np.log2(df_is_matched["siVIM-270"]/df_is_matched["mock"])
df_is_matched = df_is_matched.sort_values(by="log2(RNA/mock)").reset_index(drop=False)
df_is_matched.head(5)
def CFC_trace_create(data, name=""):
    """Create a CFC(Cumulative Frequency Curve)"""
    num_data = len(data)
    trace = go.Scatter(x=data, y=[(i+1)/num_data for i in range(num_data)], mode="lines", name=name)
    return trace
fig = {
    "data":[
        CFC_trace_create(data=df_is_matched["log2(RNA/mock)"].values, name="all"),
        CFC_trace_create(data=df_is_matched[df_is_matched["NumHits"]!=0]["log2(RNA/mock)"].values, name="seed matched mRNAs")
    ],
    "layout": go.Layout(title="The expression level(All vs seed matched mRNAs)", xaxis_title="$log_2(RNA/mock)$", yaxis_title="Cumulative frequency")
}
offline.iplot(fig)<jupyter_output><empty_output><jupyter_text>6.2.2 seedマッチで色々遊んでみる。　ここから先は、みなさんが興味を持った点について、思う存分データで遊んでいただく時間です。- マッチするシードの数（`NumHits`）って、多い方が抑制されてる？？- そもそもなんで2-8の7merでオフターゲット効果が起きるの？？  - 1-7や3-9は？？  - 6merや9merは？？  - Argonauteタンパク質と結合（loading）し、RISC(RNA induced silencing complex)を形成するが、構造的に…- 統計的に有意だと言える？？面白そうなことについては、積極的に調べて見てください！！※ 以下に、使えそうなツールを用意しておきました。ぜひ解析に役立ててください！！<jupyter_code>from teilab.seedmatch import get_matched_mRNAs
df_matched_mRNAs = get_matched_mRNAs("GAGTTCA")
df_matched_mRNAs.head(5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14800726203855796, 'fraction_numerical': 0.024379700873173683, 'mean_word_length': 8.938144329896907, 'pattern_counts': {'":': 3, '<': 77, '<?xml version=': 0, '>': 81, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '598945', 'license': 'MIT', 'n_tokens_mistral': 6473, 'n_tokens_neox': 6158, 'n_words': 548, 'path': 'notebook/Colaboratory/microarray2021S.ipynb', 'repo_name': 'iwasakishuto/TeiLab-BasicLaboratoryWork-in-LifeScienceExperiments', 'size': '11567'}