# 2022/10/28 進捗報告

## 研究概要
* [研究①](https://hackmd.io/XVKNoDqUR1q7u4tR4jcGKg)を参照すること。

## 過去の進捗報告(秋学期分)
* [第1回進捗報告(2022/4/21)](https://drive.google.com/file/d/1SiSPZDtQBaMe5kN2wp1nBQwWRELAXlOx/view)
* [第2回進捗報告(2022/5/10)](https://hackmd.io/Ekh91XAWQ9etnda4NvwuPg)
* [第3回進捗報告(2022/5/26)](https://hackmd.io/hjFDDlCZTxOD1IvmsHubVw)
* [第4回進捗報告(2022/6/7)](https://hackmd.io/R6j16PqcS_OuBFDGBmP_dQ)
* [第5回進捗報告(2022/6/21)](https://hackmd.io/9nYV1QgUTO2EV6Y8PGfQBQ)
* [第6回進捗報告(2022/7/5)](https://hackmd.io/QpKwmhhCQDqe1dk41CV_WA)
* [第7回進捗報告(2022/7/19)](https://hackmd.io/qlntToDITMKLbXSiPqFDCg)
* [第8回進捗報告(2022/10/14)](https://hackmd.io/WLDYfyviQ225stpXX6OM6w)

## 今回の進捗
### TBM参加
* 大きな失敗はなく発表できたと思う。
* 卒論発表会ではもっと長い原稿で発表しなくてはならないので大変そう。

### 3タンパク質分類タスクの成功（一応）
* ソースコードの共有の上手い方法がわからないので今回はここにコードを記載

import torch
from transformers import BertModel, BertTokenizer
from pyknp import Juman
import re
import numpy as np
import pandas as pd
import urllib.request 
import matplotlib.pyplot as plt
import sklearn #機械学習のライブラリ
from sklearn.decomposition import PCA #主成分分析器
import csv
from Bio import SeqIO
import numpy as np
import scipy.spatial.distance as distance

tokenizer = BertTokenizer.from_pretrained("prot_bert", do_lower_case=False )
model = BertModel.from_pretrained("prot_bert")
fastafile = "REall.fasta"
id_part = []
desc_part = []
seq= []

with open(fastafile, "r") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        id_part.append(record.id)
        desc_part.append(record.description)
        seq.append(record.seq)
        print("id:", id_part[-1])
        print("desc:", desc_part[-1])
        print("seq:", seq[-1].replace(""," "))
    
sequence_Example = []
for i in range(len(seq)):
    sequence_Example.append(re.sub(r"[UZOB]", "X", str(seq[i].replace(""," "))))
encoded_input = tokenizer(sequence_Example[:], return_tensors='pt',depth=8,max_length=413,padding=True)
output = model(**encoded_input)

data = output['pooler_output']
data = data.detach().cpu().numpy()

print(data)

from sklearn.decomposition import PCA

print(len(data))
print(len(data[0]))

pca = PCA()  # PCA を行ったり PCA の結果を格納したりするための変数を、pca として宣言
pca.fit(data) 

# 分析結果を元にデータセットを主成分に変換する
transformed = pca.fit_transform(data)

df = pd.DataFrame(transformed, columns=["PC{}".format(x + 1) for x in range(len(data))])
print(df)

# 主成分をプロットする
#plt.subplot(1, 2, 2)
fig, ax = plt.subplots()

for i in range(0,44):
    if "ACE2" in desc_part[i]: 
        plt.scatter(transformed[i, 0], transformed[i, 1], c='r',alpha=0.5)
    elif "TMPRSS2" in desc_part[i]: 
        plt.scatter(transformed[i, 0], transformed[i, 1], c='b',alpha=0.5)
    else:
        plt.scatter(transformed[i, 0], transformed[i, 1],c='g',alpha=0.5)

plt.title('principal component')
plt.xlabel('pc1')
plt.ylabel('pc2')


# グラフを表示する
plt.show()

# 主成分の次元ごとの寄与率を出力する
print("PC1:",pca.explained_variance_ratio_[0])

* ACE2、TMPRSS2、CaMK2を以下の図のように分類できた
* ACE2=赤、TMPRSS2=青、CaMK2=緑
* CaMK2は比較的良く分類できていそう。
* 他はまあまあか？
* 寄与率はPC1で0.99999934
* PC2とはオーダが異なる
* PC1 で ACE2 、PC2 でTMPRSS2 の違いが説明されている（？）
* →ACE2 の多様性は TMPRSS2 の多様性に比べて大きい

### 元の論文のタンパク質を利用したPCAタスク
* 元の論文は[こちら](https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1009149)
* 現在のソースコードは以下

from Bio import Entrez
from Bio import SeqIO
import datetime


def timestamp():
    # returns current time stamp
    dt = datetime.datetime.now()
    ts = dt.strftime('%d%b%Y-%H:%M:%S')
    return ts


def elapsed_time(dt_start, dt_end):
    # calculate delta between start and end
    ds = (dt_end - dt_start).total_seconds()
    h = int(ds // 3600)
    m = int((ds % 3600) // 60)
    s = int(ds % 60)
    # construct str line
    str_line = str(h).zfill(2) + 'h ' + str(m).zfill(2) + 'm ' + str(s).zfill(2) + 's'
    return str_line

print('process started: ', timestamp())
# mark start time
dt_start = datetime.datetime.now()

# declare who you are
Entrez.email =  '[e-mail]' 
Entrez.api_key = 'd262971e56d8dd17639201a57a05d9020607'

# get acc ids
acc_ids = pd.read_csv("./Results/ID.txt", encoding="UTF-8", header=None)
 
# Genbank ID数を取得
id_num = len(acc_ids)

# GWnbank ID の数だけForで情報を取得する
for i in tqdm(range(0, id_num), desc = Color.CYAN + 'Progress rate' + Color.END):
 
        handle = Entrez.efetch(db = 'Protein', id = acc_ids.iloc[i,0], rettype = 'gb', retmode = 'text')
 
        for record_info in SeqIO.parse(handle, "gb"): # gbを指定するとGenbank形式の情報を取得する
                
                # リスト格納
                record_list = []
                
                # 配列全体を格納
                parent_sequence = record_info.seq

 
                # 学名を格納
                organism = record_info.annotations['organism']
 
                # Accession numberを格納
                acc = record_info.id
                                
                # 各領域ごとに配列を切り分けて配分
                for feature in record_info.features:            
                        
                        
                        # 不明な遺伝子領域 misc_feature 以外の場合
                        if feature.type == 'source':
 
                            sequence_slice = feature.location.extract(parent_sequence)
                        
                            # 配列の末端取得
                            strand_int = feature.location.strand
                            if strand_int == 1:
                                strand = 'top' # 5'末端
                            elif strand_int == -1:
                                strand = 'btm' # 3'末端
                            else:
                                strand = 'cannotbe'
 
                            # 配列の開始位置指定
                            start = feature.location.parts[0].start.position + 1
                            # 配列の終了位置指定
                            end   = feature.location.parts[-1].end.position
                        
                        
                            # Fasta形式にした時の情報部分(> ~~~ の部分)の作成
                            id_line = acc + '_' + organism + '(' + str(start) + '..' + str(end) + ', strand=' + strand + ')'
 
                            # レコードから取得した情報をFasta形式にする
                            record_list = makeSeqRecord(parent_sequence, id_line)
 
                            # 一応0以上かチェック
                            type_count = len(record_list)
                            
                            # merge fasta file name
                            mfa = './Results/Fasta/mergefa.fasta'
                            with open(mfa, mode = 'a') as mfa:
                                    if type_count > 0:
                                            SeqIO.write(record_list, mfa, 'fasta')
                                    else:
                                        pass

# STDOUT
print()
print('process finished: ', timestamp())
et = elapsed_time(dt_start, datetime.datetime.now())
print('elapsed time: ', et)

* 現在できること
* 元論文の中からいくつかaccessionコードを抜き出してID.txtにコピーし、それをもとにDNA配列をfastaファイルにまとめる
* 実現させたいこと
* 元論文の個数、fastaファイルにまとめたい(配列数多い場合には他の方法があるかもしれない)
* DNAではProtBERTが使えないので、タンパク質配列でfastaにまとめたい
* 元論文に添付されている表にあるaccessionコードは、すべてDNAのもの
* タンパク質のコードはまた異なる
* DNAのコードからタンパク質のに変更する方法がわからない(法則は特になく、良いツールなども見当たらない)
* 配列に多く入っている改行を取り除きたい

## 次回の進捗予定
* 元タンパク質からのPCAの成功
* ProtBERTに関する勉強
* スライド作成＆論文執筆にも取り掛かり始めたい(今の状態で結果が書けるのか)