! python --version

# 見栄えの問題
import warnings
warnings.filterwarnings('ignore')

import os
import numpy as np
import pandas as pd

# データの形式
pd.read_csv("US91503671_253949442637_S01_GE1_105_Dec08_1_1.txt",sep="\t").head(15)

# 実習ファイルの形式
! tree | grep .txt

USE_COLS1 = ['FeatureNum', 'ControlType', 'ProbeName', 'GeneName', 'SystematicName']
USE_COLS2 = ['gProcessedSignal', 'gIsWellAboveBG']

samples = ['mock(1)','RNA','2OMe3','2OMe5','mock(5)','2OMe7','LNA3','LNA7']

# ※注意：sorted をつける必要あり！！（∵ os.listdir() の順番はtreeの順番とは異なる。）
file_paths = sorted([path for path in os.listdir() if path[-4:] == ".txt"])

# 確かに tree と同じ順番になっている事を確認。
file_paths

# 基盤となるデータフレーム 
df_all = pd.read_csv("US91503671_253949442637_S01_GE1_105_Dec08_1_1.txt",sep="\t",header=9, engine='python')[USE_COLS1]
all_index = df_all.index

for path,sample in zip(file_paths,samples):
    df = pd.read_csv(path,sep="\t",header=9, engine='python')[USE_COLS2]
    
    # gIsWellAboveBG が 1 のもののみを抽出する。
    select_index = df[(df.gIsWellAboveBG==1)].index
    all_index = set(all_index) & set(select_index)
    
    # gIsWellAboveBG はもういらないのでシグナル強度のみを残す。
    df = df[['gProcessedSignal']]
    df.columns = [sample]
    
    # 連結する
    df_all = pd.concat([df_all, df], axis=1)

print("データ数: {}".format(len(df_all)))

selected_index = sorted(list(all_index))

df_output = df_all.loc[selected_index,:]

df_output = df_output[df_output.ControlType == 0]

print("データ数: {}".format(len(df_output)))

# 一度保存
df_output.reset_index(drop=True).to_excel("result.xlsx", index=False)

import seaborn as sns
import matplotlib.pyplot as plt

sns.pairplot(df_output[samples])
plt.show()

def MAplot(df,col_x,col_y,col_c=None,cdict=None,title=True,axis=True,ymin=None,ymax=None,vim=True):
    # MAプロットの変数を定義
    X=df[col_x].values
    Y=df[col_y].values
    row = np.log10(X*Y)
    col = np.log2(Y/X)
    
    # 色を分ける場合は、カラム(col_c)と対応する色(cdict)を指定
    if col_c:
        group=df[col_c].values
        
        for g in np.unique(group):
            ix = np.where(group==g)
            plt.scatter(row[ix], col[ix], s=1, c=cdict[g],label=g)
    else:
        plt.scatter(row, col, s=1)
    
    # y軸の範囲[ymin, ymax]を指定
    plt.ylim(ymin,ymax)
    
    # +1,0,-1に水平線を描画
    plt.hlines(y=0, xmin=min(row), xmax=max(row), colors='r', linewidths=2)
    plt.hlines(y=1, xmin=min(row), xmax=max(row), colors='r', linewidths=1)
    plt.hlines(y=-1,xmin=min(row), xmax=max(row), colors='r', linewidths=1)
    
    # タイトルを入れる
    if title: plt.title("MA plot $X$ = {}, $y$ = {}".format(col_x, col_y))
    
    # 軸の名前を入れる
    if axis: 
        plt.xlabel("$\log_{10}XY$")
        plt.ylabel("$\log_{2}(Y/X)$")
    
    # ビメンチンだけ別に描画
    if vim:        
        df_vim = df[df.GeneName == "VIM"]
        vim_x = df_vim[col_x]
        vim_y = df_vim[col_y]
        vim_row = np.log10(vim_x*vim_y)
        vim_col = np.log2(vim_y/vim_x)
        plt.scatter(vim_row, vim_col, s=5, c="black", marker="*", label="vimentin")

N = len(samples)
fig = plt.figure(figsize=(16, 16))

for j,col_y in enumerate(samples):
    for i,col_x in enumerate(samples):
        if j<i: continue # 修飾と修飾をMAプロットで比較する意味があるかはわからないが一応
        plt.subplot(N, N, j*N + i+1)
        MAplot(df_output,col_x,col_y,title=None,axis=None)
        if i==0: plt.ylabel(col_y)
        if j==len(samples)-1: plt.xlabel(col_x)
plt.show()

sns.pairplot(df_output[["mock(1)","mock(5)"]])
plt.show()

MAplot(df_output,"mock(1)","mock(5)")
plt.show()

# mock2つを並べてプロット
def mockplot(df,col_y,col_c=None,cdict=None,title=True,axis=True,ymin=None,ymax=None,vim=True):
    fig = plt.figure(figsize=(12,4))
    plt.subplot(1,2,1)
    MAplot(df,"mock(1)",col_y,col_c=col_c,cdict=cdict,title=title,axis=axis,ymin=ymin,ymax=ymax,vim=vim)
    plt.subplot(1,2,2)
    MAplot(df,"mock(5)",col_y,col_c=col_c,cdict=cdict,title=title,axis=axis,ymin=ymin,ymax=ymax,vim=vim)

# mock-RNA

# X-Y plot
sns.pairplot(df_output[["mock(1)","mock(5)","RNA"]])
plt.show()

# MA plot
mockplot(df_output, col_y="RNA")
plt.show()

raw_data = df_output[samples]

raw_data.shape

# 相乗平均を求める
from scipy.stats.mstats import gmean

tile_75s = [sorted(raw_data[col])[int(len(raw_data)*0.75)] for col in samples]

a = gmean(tile_75s)

tile_75_data = a * df_output[samples]/tile_75s

tile_75_data.shape

quantile_data = gmean(np.array([sorted(raw_data[col]) for col in samples]),axis=0)

quantile_data.shape

X = [i for i in range(len(samples))]
Y = df_output[df_output.GeneName == "VIM"][samples].values[0]
plt.xticks(X,samples)
plt.ylabel("gProcessedSignal")
plt.bar(X,Y)
plt.title('The relationship between \n"sRNA modification" and "gProcessedSignal"')
plt.show()

df_sm = pd.read_excel("seedmatch.xlsx")

# データ形式の確認
df_sm.head(3)

print("データ数: {}".format(len(df_sm)))

# マージするためにdf_outputとカラム名を揃える。（なお、N_seedはシード領域の配列をいくつ含むか）
df_sm.columns = ["SystematicName","N_seed"]

df_with_Nseed = pd.merge(df_output,df_sm,on="SystematicName", how="left").fillna(0)

df_with_Nseed.groupby("N_seed").size()

# 先にN_seedが小さいものからプロットすることで、グラフを見やすくする。
df_with_Nseed_sorted = df_with_Nseed.sort_values(by="N_seed")

color_dict = {0: "#EAECEE", 1:"#FBEEE6", 2:"#EB984E",3:"#873600",4:"#641E16"}

mockplot(df_with_Nseed_sorted,col_y="RNA",col_c="N_seed",cdict=color_dict,ymin=-1,ymax=1)
plt.legend()
plt.show()

mockplot(df_with_Nseed_sorted,"2OMe3",col_c="N_seed",cdict=color_dict,ymin=-1,ymax=1)
plt.legend()
plt.show()

mockplot(df_with_Nseed_sorted,"2OMe5",col_c="N_seed",cdict=color_dict,ymin=-1,ymax=1)
plt.legend()
plt.show()

mockplot(df_with_Nseed_sorted,"LNA3",col_c="N_seed",cdict=color_dict,ymin=-1,ymax=1)
plt.legend()
plt.show()

# mockは、平均値を利用する。
df_with_Nseed["mock"] = df_with_Nseed[["mock(1)","mock(5)"]].mean(axis=1)

# 発現量を調べる。
df_with_Nseed["log2(siRNA/mock)"] = np.log2(df_with_Nseed["RNA"]/df_with_Nseed["mock"])

df_signal = df_with_Nseed.sort_values(by="log2(siRNA/mock)")

signal_seed = df_signal[df_signal["N_seed"]!=0]["log2(siRNA/mock)"].values
signal_no_seed = df_signal[df_signal["N_seed"]==0]["log2(siRNA/mock)"].values

plt.plot(signal_no_seed,[(i+1)/len(signal_no_seed) for i in range(len(signal_no_seed))], label="others")
plt.plot(signal_seed,[(i+1)/len(signal_seed) for i in range(len(signal_seed))], label="have seed regions")
plt.title('The difference between\n"no seed regions" and "have seed resions"')
plt.xlabel("$\log_{2}(siRNA/mock)$")
plt.ylabel("Cumulative frequency")
plt.legend()
plt.grid()
plt.vlines(x=0, ymin=0, ymax=1, colors='black', linewidths=1)
plt.xlim(-1,1)
plt.show()

signal_no_seed = df_signal[df_signal["N_seed"]==0]["log2(siRNA/mock)"].values
signal_1_seed  = df_signal[df_signal["N_seed"]==1]["log2(siRNA/mock)"].values
signal_2_seed  = df_signal[df_signal["N_seed"]==2]["log2(siRNA/mock)"].values
signal_3_seed  = df_signal[df_signal["N_seed"]==3]["log2(siRNA/mock)"].values
signal_4_seed  = df_signal[df_signal["N_seed"]==4]["log2(siRNA/mock)"].values

plt.plot(signal_no_seed,[(i+1)/len(signal_no_seed) for i in range(len(signal_no_seed))], label="others")
plt.plot(signal_1_seed,[(i+1)/len(signal_1_seed) for i in range(len(signal_1_seed))], label="have 1 seed regions")
plt.plot(signal_2_seed,[(i+1)/len(signal_2_seed) for i in range(len(signal_2_seed))], label="have 2 seed regions")
plt.plot(signal_3_seed,[(i+1)/len(signal_3_seed) for i in range(len(signal_3_seed))], label="have 3 seed regions")
plt.plot(signal_4_seed,[(i+1)/len(signal_4_seed) for i in range(len(signal_4_seed))], label="have 4 seed regions")
plt.title('The difference between\n"no seed regions" and "have seed resions"')
plt.xlabel("$\log_{2}(siRNA/mock)$")
plt.ylabel("Cumulative frequency")
plt.legend()
plt.grid()
plt.vlines(x=0, ymin=0, ymax=1, colors='black', linewidths=1)
plt.xlim(-1,1)
plt.show()

df_output["mock"] = df_output[["mock(1)","mock(5)"]].mean(axis=1)

df_output["log2(siRNA/mock)"] = np.log2(df_output["RNA"]/df_output["mock"])

# 毎回 seedmatch のGUIで検索するのは面倒なので、リクエストを飛ばすことにする。
import requests
import re

def seedNmer(df, N, match_array, arrayname="", return_data=False):
    """
    関数の概要：色々な領域のNmerの発現量を調べる。
    @param df         ：元となるデータフレーム 
    @param N          ：領域の長さ
    @param match_array：ガイド鎖やパッセンジャー鎖のシード領域と相補的なmRNAの配列（5'-3'）
    """
    PATTERNS = r'\n(.+)\t([0-9]+)'
    GMEANS = []
    COLNAMES = []
    LENGTH = len(match_array)
    df_Nmer = df.sort_values(by="log2(siRNA/mock)")
    plt.figure(figsize=(12,8))
    
    for p in reversed(range(len(match_array)+1)):
        colname = "{}-{}mer({}-{})".format(arrayname, N,LENGTH-p+1,LENGTH-p+N)
        COLNAMES.append(colname)

        # "seedmatch"にリクエスト飛ばしてseed領域に対応する配列を検索
        array = match_array[p-N:p]
        r = requests.post('http://atlas.rnai.jp/seedmatch/seedmatch.cgi', data={"seedseq": array})
        df_mRNA = pd.DataFrame(re.findall(PATTERNS, r.text))
        df_mRNA.columns=["SystematicName", colname]

        # マージ
        df_Nmer = pd.merge(df_Nmer, df_mRNA, on="SystematicName", how="left").fillna(0)
        # マッチする配列を持っているものを抜き出して
        signal_7mer = df_Nmer[df_Nmer[colname]!=0]["log2(siRNA/mock)"].values
        GMEANS.append(np.mean(signal_7mer))
        length = len(signal_7mer)
        # 累積密度をプロット
        plt.plot(signal_7mer,[(i+1)/length for i in range(length)], label=colname + ": {}".format(length))

        if p-N==0: break
    
    #=== プロット後処理 ===
    plt.title('The difference between\n"seed match resions" and "Expression levels"')
    plt.xlabel("$\log_{2}(siRNA/mock)$")
    plt.ylabel("Cumulative frequency")
    plt.grid()
    plt.vlines(x=0, ymin=0, ymax=1, colors='black', linewidths=2)
    plt.legend()
    plt.xlim(-1,1)
    plt.show()
    
    #===  ===
    X = [i for i in range(len(GMEANS))]
    label_before = [i for i in range(len(GMEANS)) if i %2 == 0]
    label_after  = ["{}-{}".format(p+1,p+N) for p in range(len(GMEANS)) if p %2 == 0]

    plt.barh(X,GMEANS)
    plt.title('The difference between\n"seed match resions" and "Expression levels"')
    plt.xlabel("7mer position")
    plt.ylabel("average $\log_{2}(siRNA/mock)$")
    plt.yticks(label_before,label_after)
    plt.show()
    
    if return_data: 
        return (df_Nmer, COLNAMES)

guide_array   = "UUGAACUCGGUGUUGAUGGCG" # 5'-3'
g_match_array = "CGCCATCAACACCGAGTTCAA" # 5'-3'

df_guide_7mer, guide_7mer_colnames = seedNmer(df_output, 7, g_match_array, arrayname="guide", return_data=True)

df_guide_7mer[guide_7mer_colnames].head(3)

passenger_array = "CCAUCAACACCGAGUUCAAGA" # 5'-3'
p_match_array   = "TCTTGAACTCGGTGTTGATGG" # 5'-3'

seedNmer(df_output, 7, p_match_array, arrayname="passenger")

seedNmer(df_output, 8, g_match_array, arrayname="guide")

seedNmer(df_output, 5, g_match_array, arrayname="guide")

df_output.head(3)

from sklearn.ensemble import RandomForestRegressor
from sklearn import preprocessing

tree = RandomForestRegressor(random_state=0, max_depth=15)

x = df_guide_7mer[guide_7mer_colnames].values.astype(float)
y = df_guide_7mer["log2(siRNA/mock)"].values.astype(float)

tree.fit(x, y)

plt.scatter(y,tree.predict(x), s=1)
plt.title('predicted "$\log_2(siRNA/mock)$" and answer "$\log_2(siRNA/mock)$"')
plt.ylabel("answer")
plt.xlabel("predict")
plt.show()

n_features = len(guide_7mer_colnames)
plt.figure(figsize=(12, 8))
plt.barh(range(n_features), tree.feature_importances_ , align='center')
plt.yticks(np.arange(n_features), guide_7mer_colnames)
plt.title("importances for predict the Expression levels")
plt.show()