! python --version

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

import pandas as pd

df = pd.read_excel("result.xlsx")

df.head(3)

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

import re
import requests
from tqdm import tqdm

guide_match_array = "gagttca"

front_patterns =r'<meta name="ncbi_uidlist" content="(.+)" /><meta name="ncbi_filter"'

back_request_format = "https://www.ncbi.nlm.nih.gov/sviewer/viewer.fcgi?id={}&db=nuccore&report=genbank&conwithfeat=on&hide-cdd=on&retmode=html&withmarkup=on&tool=portal&log$=seqview&maxdownloadsize=1000000"

back_patterns = ">([atgc ]*)</span>"

def gene2array(GeneName):
    try:
        #=== httpリクエストを飛ばしてプログラムで検索を行う ===
        front_response = requests.get('https://www.ncbi.nlm.nih.gov/nuccore/', params={'term': GeneName})
        genome_id = re.findall(front_patterns, front_response.text)[0]  # ex) 1519241814
        back_response = requests.get(back_request_format.format(genome_id))
        sequences = re.findall(back_patterns, back_response.text) # ex) ['atgcatgcat gcatgca...', 'at..', .., '..atgc']
        sequence = "".join(sequences).replace(" ","")
        
        #=== 配列の長さ ===
        length = len(sequence)
        
        #=== seed領域と一致するポジションの5'末端からの割合 ===
        position = re.search(guide_match_array, sequence)
        p = position.end()/length if position else 0
        
        return [sequence,length,p]
    except:
        return [None, None, None]

# GeneNameを全て調べる。
gene_names = df.SystematicName.values

sample_gene_name = gene_names[17]

print("sample gene name: {}".format(sample_gene_name))

gene2array(sample_gene_name)

# 階層構造は以下の通り。
! tree microarray

publicDNSs = [
    "3-113-33-209","3-112-149-206","3-112-213-203",
    "13-113-177-99","54-249-50-45","13-231-179-134",
    "13-115-166-8","3-112-73-40","52-198-222-27",
    "13-115-220-149","18-182-16-178"
]

begins = [  i  * 2000 for i in range(11)]
ends   = [(i+1)* 2000 for i in range(11)]

# 以下のファイルが作られた。
! tree AWS

import os

file_names = [name for name in sorted(os.listdir("AWS")) if name[-4:] == ".csv"]

# 確認
file_names

df_all = pd.DataFrame()
for filename in file_names:
    tmp = pd.read_csv(os.path.join('AWS',filename))
    df_all = pd.concat([df_all,tmp])

# 結合
df_cocated = pd.concat([df_all[['sequence','length','p']].reset_index(drop=True),df.reset_index(drop=True)], axis=1)

df_cocated.head(3)

# かなりNaN(配列データが取得できなかったもの)が多い
df_cocated.isnull().sum()

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

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

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

# 正しく配列データが取得できたものだけを利用する。
df_used = df_cocated[df_cocated.length>0].reset_index(drop=True)

print("取得できたmRNA: {}".format(len(df_used)))

# 取得したmRNAの配列の長さの分布
plt.hist(df_used.length)
plt.title("The histgram of the length of mRNA")
plt.xlabel("length")
plt.ylabel("N")
plt.show()

print("7merマッチ出現量の期待値: {}".format(df_used.length.sum()/(4**7)))

# 取得した配列の発現量の分布
plt.hist(df_used["log2(siRNA/mock)"], bins=20)
plt.title("The histgram of the expression level of mRNA")
plt.xlabel("expression level")
plt.ylabel("N")
plt.show()

# シード領域とマッチする領域の出てくる場所の割合
def match_rate(match_array, array, length, rate=False):
    position = re.search(match_array, array)
    if position:
        p = position.end()/length if rate else position.end()
    else:
        p = 0
    return p

# うまく取得できていなかったので、これらのデータを利用
df_used["position"] = df_used.apply(lambda x: match_rate('gagttca', x.sequence, x.length), axis=1)

df_used["position_rate"] = df_used.apply(lambda x: match_rate('gagttca', x.sequence, x.length, rate=True), axis=1)

df_used[["length","log2(siRNA/mock)"]].corr()

df_used_with_seed = df_used[df_used.position>0]

print("取得できたマッチ領域を含むmRNA: {}".format(len(df_used_with_seed)))

df_used_with_seed[["length","position","position_rate","log2(siRNA/mock)"]].corr()

sns.heatmap(
    df_used_with_seed[["length","position","position_rate","log2(siRNA/mock)"]].corr(),
    vmin=-1.0,
    vmax=1.0,
    center=0,
    annot=True, # True:格子の中に値を表示
    fmt='.3f',
)
plt.show()

def seedmatchNmer(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'

seedmatchNmer(df_used_with_seed, 7, g_match_array, arrayname="guide", return_data=False)

def seedNCBINmer(df, N, match_array, arrayname="", return_data=False):
    """
    関数の概要：色々な領域のNmerの発現量を調べる。
    @param df         ：元となるデータフレーム 
    @param N          ：領域の長さ
    @param match_array：ガイド鎖やパッセンジャー鎖のシード領域と相補的なmRNAの配列（5'-3'）
    """
    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)):
        array = match_array[p-N:p]
        colname = "{}-{}mer({}-{})".format(arrayname, N,LENGTH-p+1,LENGTH-p+N)
        COLNAMES.append(colname)
        
        df_Nmer[colname] = df_Nmer.apply(lambda x: match_rate(array, x.sequence, x.length), axis=1)

        # マッチする配列を持っているものを抜き出して
        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)

df_guide_7mer, guide_7mer_colnames = seedNCBINmer(df_used, 7, g_match_array.lower(), "guide", return_data=True)

feature_columns = guide_7mer_colnames+["length","position","position_rate"]

# 相関関係を調べる。
sns.heatmap(
    df_guide_7mer[feature_columns].corr(),
    vmin=-1.0,
    vmax=1.0,
    center=0,
    fmt='.3f',
)
plt.show()

from sklearn.ensemble import RandomForestRegressor
from sklearn import preprocessing

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

X = df_guide_7mer[feature_columns].values.astype(float)
y = df_guide_7mer["log2(siRNA/mock)"].values.reshape(-1,1).astype(float)

tree.fit(X, y)

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

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