ノート
訪問者数　281　　　　　　最終更新　2020-03-02 (月) 11:56:09

RNAseq〜岸本データの後続分析 †

TPM補正が終わった実験データをどう分析するか考える。

↑

データ †

TPM補正済みデータ（発現値の補正済みデータ）は

	Chr	Start	End	Strand	Length	10B.sam	10D_minus.sam	1p2-1.sam	1p2-2.sam	2-10B_minus.sam
gene_0001_thrL	AP012030.1	190	255	+	66	1119.230093	2070.357664	5212.719939	4193.551146	1955.183064
gene_0003_thrA	AP012030.1	338	2800	+	2463	1872.865147	4391.896215	5773.99795	6352.945352	3971.157253
gene_0005_thrB	AP012030.1	2802	3734	+	933	1925.903487	4443.130625	3779.116205	4639.666182	3926.988291
gene_0007_thrC	AP012030.1	3735	5021	+	1287	2288.887037	3590.528484	3897.465673	4230.197519	3276.631145
gene_0008_yaaX	AP012030.1	5235	5531	+	297	186.9824878	553.2073517	293.9131871	321.6221759	585.3916812
gene_0009_yaaA	AP012030.1	5684	6460	-	777	42.27207978	45.02772533	73.69652755	49.63567721	49.63092301
gene_0011_yaaJ	AP012030.1	6530	7960	-	1431	29.68188822	18.17449013	19.74707008	17.02170195	21.92416
gene_0014_talB	AP012030.1	8239	9192	+	954	3197.624287	397.0260641	1580.352153	1586.773069	443.3769857
gene_0016_mog	AP012030.1	9307	9894	+	588	218.27842	87.23291626	137.3536608	138.6934099	94.32622154
gene_0017_yaaH	AP012030.1	9929	10495	-	567	146.5658456	50.05540317	89.97128851	78.75914474	61.42555542

のようになっている。

注意：　同じgene名で、複数の異なるCDSになっている場合がある。つまりgene名ではユニークではなく、CDS位置でユニークになる。これらはトランスポゾン（insA〜Nなど）などゲノム内を飛び回る遺伝子らしい。

gene名で見て重複のあるCDSを取出してみると

geneStartLengthAnc
gene_3114arpB17864754747.176203468
gene_3115arpB178695914161.441322222
gene_3782gatR21534954476.087732338
gene_3786gatR2155197339163.2194667
gene_2015icd117251512511782.60336
gene_2058icd119006616512.36916525
gene_6818ilvG39296659841955.874255
gene_6819ilvG3930728582942.1393394
gene_0495insA2902962760
gene_1774insA10401612766.572986365
gene_6204insA35613892760.821623296
gene_7778insA449829727637.7946716
gene_0035insB198125043.149555967
gene_0475insB2784035040
gene_0494insB2898745042.6996194
gene_1775insB104035550464.34092903
gene_7637insB44136465043.599492533
gene_7779insB449849129430.85279314
gene_7781insB449878521078.82888648
gene_0647insC3817304110
gene_1727insC10164934110
gene_2221insC12803544110
gene_5485insC31640494111.103494061
gene_7744insC44779974110
gene_0648insD3820809240
gene_1728insD10168619060
gene_2222insD12807229060
gene_2529insD14501039060
gene_3625insD20527559060
gene_5486insD31644179060
gene_7745insD44783479242.6996194
gene_0535insE3157063090
gene_0536insE3157153000
gene_2021insE11773203090
gene_3785insE21540383000
gene_0467insH2733269810
gene_0505insH2944579810
gene_0973insH5663999810
gene_1163insH6782579810
gene_1843insH108165110170
gene_2218insH12782209810
gene_2394insH137705910170.222977414
gene_2401insH138060910170.222977414
gene_2466insH140992898116.18120497
gene_3011insH17257509810
gene_3316insH18933789810
gene_3623insH2050108101728.31813152
gene_3690insH208569810170
gene_3792insH215819510175.351457925
gene_3983insH227405910170
gene_5382insH31081219810
gene_5806insH33436089810
gene_6318insH363008810170
gene_6713insH386872810170
gene_0463insI269828115216.53516882
gene_2530insI145154011521.181083487
gene_7761insI4487236115211.22029313
gene_0028insL1544011190
gene_1028insL59831911190
gene_4345insL249920611190
gene_0462insN2694674059.518658032
gene_7760insN448696726740.76728622
gene_0464insO271055426117.11025
gene_7764insO448872859720.89152701
gene_6038kefG34585125550
gene_6039kefG34585125520
gene_2164ldrB12478211350
gene_2167ldrB124835613568.87029047
gene_2169ldrB12488911350
gene_2465lomR140964015615.99005337
gene_2469lomR14109961711.326128828
gene_3836molR21814688253.573314406
gene_3838molR218240419388.658841171
gene_3840molR21843049604.251900555
gene_7426phnE43014373694.916380045
gene_7427phnE43016556215.842654547
gene_1997potA116185111370
gene_1998potA116185111190
gene_6273rhsB359709842363.051411163
gene_6498rhsB374686929413.11243709
gene_1713sulA10094025160
gene_1714sulA10094025100
gene_3689wbbL208519928299.71360166
gene_3691wbbL2086719474189.4517716
gene_1830ycdN10716841114.085910443
gene_1832ycdN107179472011.65335708
gene_2077ycgH1198254152132.05465244
gene_2078ycgH11998591017105.0223618
gene_2212ychG127307959144.50946097
gene_2213ychG127362123154.97406778
gene_2284yciX131614312954.49464277
gene_2285yciX1316253189278.3607559
gene_2528ydbA1447574255913.20376569
gene_2532ydbA14528723324120.0697148
gene_3511yedN19951511921.181083487
gene_3512yedN19953523210.70644246
gene_3570yedS20178544865.132609723
gene_3572yedS201834921016.1977164
gene_3574yedS201864240519.03731606
gene_3597yeeL20360793273.467401064
gene_3599yeeL20364277050.964970339
gene_4039yfaS231467948615.39782917
gene_4041yfaS2315180410411.65888261
gene_4169yfcC240196215210.298182813
gene_4171yfcC240200414790
gene_7747yjgX447981343258.26678538
gene_7748yjgX448020236030.23573728
gene_0530ykgM3127981413697.444681
gene_0531ykgM3129382672342.420321
gene_0940ylbE54878112600.179974627
gene_0941ylbE54903910020
gene_2630yncI151276874721.55358782
gene_2632yncI151355820115.79478813
gene_2579yncK148532516831.0456231
gene_2580yncK148554428837.7946716
gene_2084ypjA12023482135.323193183
gene_4790ypjA2756055450023.98701824

これらについては、発現量解析の当初では放置しておいた（CDS位置で区別して処理して いた）が、遺伝子の発現量としてみると何らかの判断をしなければならない。

内容を見ると（岸本先生による）、同じ配列がコピーされている場合と、１つの配列が途中に余分な夾雑物が入って切れているように見える場合がある。さらに移動の途中で短くなったり長くなったりしているケースがある。これらの判断は難しいが、

それぞれ、発現量について次のような対応をすることにする。

これらによって、gene名を（CDS位置に依らず）ユニークにした。

↑

やりたいこと： †

サンプル間の変動パターンをgeneごとに求めて、gene間で比較して似ているパターン（相補的なパターンを含めて）を探す。もしパターン間の距離を定義できれば、距離の近いものをグループ化する（クラスタ化）と同時に、距離を元にグラフ化して考えることでgene間の関連（反応等のつながり）を考えることができる。

変動の対象として、データから系列に従って、　Anc - ... - 1_2-1 - 2_5-1 と　Anc - ... - 2_2-1 - 2_6-2　を抽出する。

↑

例外的パターンについて †

0 - ... - 0 - 0 - 0　つまり一切発現無し、があり得る。これはおそらく分析対象から除外してよいだろう。

0 - 途中にnon-0がある - ...　これらのケースは、例外視する必要はなく、１つのパターンと考えてよかろう。Ancが0であるケースは、「なかったものが出てくるようになった」という意味で少し考える必要があるかも知れないが、とりあえず１つのパターンとして考えることにする。

なお、後述のように、変動として系列上の前の値との比を考えると、0で割ることになるので、微小値に置き換えるなどの工夫が必要になるし、微小値にするとlog10を取ると負の大きな値になって多少始末が悪い。

同じgene名が複数のCDSに現れる件について。？？

↑

変動の数値指標について †

�＼亀�化：
発現量の絶対値は、今仮に意味が無いと考える。欲しいのはサンプル間で見たときの増減のパターン（増減の方向と多寡、それが一連のサンプル間でどういう連鎖か）なので、数値はgene間で正規化する必要があるが、１つの方法として同一gene内での変動比率を取ってしまう。
２つの比率が考えられる：１番目は基準値（たとえばAncか、または全体の平均値、最頻値のようなもの）に対する比率、２番目はサンプル間の比率。

��値のlog（log10）による圧縮
これは、発現量の絶対値がgeneによって桁が違う点を考慮するために導入を考えられるが、もし�，波耄┐鮃佑┐襪箸修譴曚桧嫐�は無い。また比率は、対数化した後では差になることにも留意する。
対数化の注意点は、発現値が0の場合にlogが取れないことで、この場合微小値に置き換えることが考えられる（但しlogを取った結果の値は負の大きな数になる）。系列上すべて0のケースは解析対象から外すことが考えられるが、一部の発現値のみ0のケースは意味があり得るので、外すことは考えない。

↑

途中の道草 †

今はlog(tpm/Anc)を取った表。

geneStartLengthAnc43B45a_minus45A_minus45L1_2-12_5-1
gene_0001thrL1906600.5188013710.1208220120.1305530920.6622781990.7369817060.782980429
gene_0003thrA338246300.3635458120.2172938690.2070900280.4623602140.4576295920.265097062
gene_0005thrB280293300.3756433540.3384470830.3130375920.5336653190.356697460.17628472
gene_0007thrC3735128700.3519461810.1782091510.1605569290.4493966540.3094693790.092155761
gene_0008yaaX523529700.3618048210.4974061890.5658562190.4833611960.2549799220.141397537
gene_0009yaaA568477700.127790066-0.16041128-0.0464044330.0151052210.04435137-0.083863423
gene_0011yaaJ653014310-0.0619139110.1208024030.2154775330.1323669250.108788741-0.028453563
gene_0014talB823995400.025617922-0.296533334-0.4083314670.0772363670.2417954450.203333867
gene_0016mog93075880-0.18900972-0.26035833-0.255735167-0.114590533-0.140329364-0.159416034
gene_0017yaaH99295670-0.151212522-0.346290801-0.274411814-0.177810069-0.161112427-0.103483574

上記データに対してclusteringを行った結果（全部行うと距離計算に時間が非常にかかるので100個だけ）

fileCompareTPM.pdf

↑

全体をクラスタ化 †

%matplotlib inline
#
import pandas as pd
import numpy as np
from scipy.spatial import distance
from scipy.cluster.hierarchy import linkage, dendrogram
import matplotlib.pyplot as plt
import math
import os
import pickle

def dfnormalize(row):  # Ancが0ならオール0、そうでなければlog(u/Anc)
    Anc = row['Anc']
    rest = row.to_list()[3:]
    #print('rest\n', rest)
    if Anc==0:
        result = [0] * len(rest)
    else:
        result = [math.log10(u/Anc) if u>0.000001 else math.log10(0.000001/Anc) for u in rest] 
    #print('result\n', result)
    output = pd.Series(result, index=(row.index.to_list()[3:]))
    output['gene'] = row['gene']
    output['Anc'] = 0 if Anc==0 else 0
    output['Start'] = row['Start']
    output['Length'] = row['Length']
    #print('output\n', output)
    return(output)

#def myeuc(u):   # Euclidean距離を計算する関数〜mapするために用意
#    #print('u\n', u, '\ndfl.loc[target]\n', dfl.loc[target])
#    result = distance.euclidean(u, dfl[['Anc', '43B', '45a_minus', '45A_minus', '45L', '1_2-1', '2_5-1']].loc[target])
#    #print('result:', result)
#    return(result)

picklefname = 'DistanceTest.pickle'
slist = ['Anc', '43B', '45a_minus', '45A_minus', '45L', '1_2-1', '2_5-1']
    
if not os.path.exists(picklefname):
    fname = 'count_tpm.tsv'
    df = pd.read_csv(fname, sep='\t', index_col=0)
    #print(df.columns.to_list())
    df = df.rename(columns=
    {'10B.sam': '45a_2-10Bplus', 
     '10D_minus.sam': '45a_10D_minus', 
     '1p2-1.sam': '1_2-1', 
     '1p2-2.sam': '1_2-2', 
     '2-10B_minus.sam': '45a_2-10B_minus', 
     '2p5-1.sam': '2_5-1', 
     '2p6-1.sam': '2_6-1', 
     '43B.sam': '43B', 
     '45A_minus.sam': '45A_minus', 
     '45A_plus.sam': '45A_plus', 
     '45L.sam': '45L', 
     '45a10D_plus.sam': '45a_10D_plus', 
     '45aIII6c_plus.sam': '45a_III6c_plus', 
     '45a_minus.sam': '45a_minus', 
     '45a_plus.sam': '45a_plus', 
     '45alll6c_minus.sam': '45a_III6c_minus', 
     '45b_minus.sam': '45b_minus', 
     '45b_plus.sam': '45b_plus', 
     '45c_minus.sam': '45c_minus', 
     '45c_plus.sam': '45c_plus', 
     '45d7B_minus.sam': '45d_7B_minus', 
     '45d7B_plus.sam': '45d_7B_plus', 
     'Anc.sam': 'Anc', 
     'PwOw_minus.sam': 'PwOw_minus', 
     'PwOw_plus.sam': 'PwOw_plus'           
    })
    df['gene'] = [u[10:] for u in df.index.to_list()]
    df.index = [u[:9] for u in df.index.to_list()]
     
    dfdup = df[df.duplicated(subset='gene', keep=False)]\
    [['gene', 'Start', 'Length', 'Anc', '43B', '45A_minus', '45L', \
          '1_2-1', '2_5-1']].sort_values(['gene', 'Start'])
    dfdup.to_excel('DuplicatedCDS.xlsx')
    print(dfdup)
    
    df1 = df.copy()[['gene', 'Start', 'Length', \
          'Anc', '43B', '45a_minus', '45b_minus', '45c_minus', '45A_minus', '45L', \
          '1_2-1', '2_5-1']]

    df1 = df1[df1['Anc']!=0]   # Ancが0のものを除く（Ancで割るから）
    df1x = df1[:].apply(dfnormalize, axis=1)
    df1x = df1x[['gene', 'Start', 'Length', \
             'Anc', '43B', '45a_minus', '45A_minus', '45L', '1_2-1', '2_5-1']]
    df1x.to_excel('CompareTPM.xlsx')

############
# Line Graphs
############
    #df1g = df1x[['Anc', '43B', '45a_minus', '45A_minus', '45L', '1_2-1', '2_5-1']]
    #min = 0.4  # 絶対値が0.4以上のデータ点だけプロット
    #df1g = df1g[(abs(df1g['45a_minus'])>min) & (abs(df1g['45A_minus'])>min) &\
    #            (abs(df1g['45L'])>min) & (abs(df1g['1_2-1'])>min) & (abs(df1g['2_5-1'])>min) ]
    #
    #df1g.T.plot()
    #plt.show()

    df2 = df.copy()[slist]
    # df2 = df2 + 1  # 値0を避けるため ⇒　すべきではない。むしろ微小な正数にすべきだろう。
    df2 = df2 + 0.00000001
    dfl = np.log10(df2)   # 先にlog10を取る
    dfl_t = dfl.T

    # オール0の行は除かなければならない
    dfl_copy = dfl.replace(0.0, np.nan).dropna(how='all', axis=0)
    dfl = dfl.loc[dfl_copy.index]
    print('dfl\n', dfl.head()); print()
    dfl.to_pickle(picklefname)
else:
    dfl = pd.read_pickle(picklefname)

pickle2fname = 'CompareTPM2.pickle'
if not os.path.exists(pickle2fname):
    #dfl = dfl.head(100)

    # targetと他のgeneとの対の距離を計算する
    genenamelist = dfl.T.columns.to_list()
    print(genenamelist)
    #for target in genenamelist[:50]:
    for target in genenamelist:
        #print('dfl[slist].loc[target]\n', dfl[slist].loc[target])
        dfl['D_'+target] = dfl[slist].apply(lambda x: \
            distance.euclidean(x, dfl[slist].loc[target]), axis=1)

        print('target:', target)
        #dfx = dfl.sort_values('d', ascending=True)[:10]
        #print('dfx\n', dfx)
        #dfg = dfx.drop('d', axis=1).T
        #dfg.plot()
        #plt.ylabel('$log_{10}(TPM)$')
        #plt.title(target)
        #plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
        #plt.show()
       
    print(dfl.drop(columns=slist).head(), '\n')
    dfl.to_pickle('CompareTPM2.pickle')
else:
    dfl = pd.read_pickle(pickle2fname)

pickleLinkagefname = 'CompareTPMLinkage.pickle'
if not os.path.exists(pickleLinkagefname):
    dArray = distance.squareform(dfl.drop(columns=slist))
    result = linkage(dArray, method = 'average')
    node_labels = [u[2:] for u in dfl.drop(columns=slist).columns.to_list()]
    
    with open(pickleLinkagefname, 'wb') as pfw:
        pickle.dump([result, node_labels], pfw)
else:
    with open(pickleLinkagefname, 'rb') as pf:
        result, node_labels = pickle.load(pf)

plt.figure(figsize=(100,100), dpi=200, facecolor='w', edgecolor='k')
dendrogram(result, labels=node_labels)
plt.savefig('CompareTPM.pdf')
plt.show()
print('complete')

更に、出力のCompareTPMLinkage.pickleを読んで、fclusterでクラスタを出力

%matplotlib inline
# 最後の部分（Linkage計算より後ろの部分）だけ
import pandas as pd
import numpy as np
from scipy.spatial import distance
from scipy.cluster.hierarchy import linkage, dendrogram
import matplotlib.pyplot as plt
import os
import pickle

# gene_numberからgene_nameへの変換辞書
fname = 'count_tpm.tsv'
df = pd.read_csv(fname, sep='\t', index_col=0)
gene_name_dict = {u[:9]: u[10:] for u in df.index.to_list()}
#print(gene_name_dict)

pickleLinkagefname = 'CompareTPMLinkage.pickle'
with open(pickleLinkagefname, 'rb') as pf:
    result, node_labels = pickle.load(pf)
# print(result[:10])
NUM_CLUSTERS = 10
for num in range(10, NUM_CLUSTERS+1):
    labels = fcluster(result, t=num, criterion='maxclust')
    #fclusterは、入力がどのクラスタに属するか（クラスタ番号 labels）を返す
    #print(num, labels)
    # クラスタごとに、それに属する入力をリストとして表示
    clusters = []
    for cl_id in range(1, num+1):
        l = [gene_name_dict[ node_labels[n] ] for n in range(0,len(labels)) if labels[n]==cl_id]
        #print(' ', cl_id, l)
        clusters.append([cl_id, l])
    with open('clusters_'+str(num)+'.pickle', 'wb') as pwf:
        pickle.dump(clusters, pwf)
print('complete')

クラスタごとに、それぞれに属するgeneの発現変動をグラフ表示する。

%matplotlib inline
# クラスタに属するgeneのグラフを描く
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pickle

fname = 'count_tpm.tsv'
df = pd.read_csv(fname, sep='\t', index_col=0)
gene_name_dict = {u[:9]: u[10:] for u in df.index.to_list()}
#print(gene_name_dict)

NUM = 10
with open('clusters_'+str(NUM)+'.pickle', 'rb') as pf:
    clusters = pickle.load(pf)

picklefname = 'DistanceTest.pickle'
dfl = pd.read_pickle(picklefname)
#print(dfl.index.to_list())
for ucl_id, l in clusters:
    print(l)
    ############
    # Line Graphs
    dfl['gene'] = [gene_name_dict[u] for u in dfl.index]
    #print(l, dfl[dfl['gene'].isin(l)])
    dfg = dfl[dfl['gene'].isin(l)]
    dfg = dfg[['Anc', '43B', '45a_minus', '45A_minus', '45L', '1_2-1', '2_5-1']]
    dfg = dfg.iloc[0:10, :]
    dfg.T.plot()
    plt.show()

結論は、

cluster_1.pngcluster_2.pngcluster_3.pngcluster_4.pngcluster_5.png
cluster_6.pngcluster_7.pngcluster_8.pngcluster_9.pngcluster_10.png
↑

おまけ †

必須遺伝子の発現パターンはどうなっているのか？

添付ファイル: filecluster_9.png 487件 [詳細] filecluster_10.png 479件 [詳細] filecluster_8.png 460件 [詳細] filecluster_7.png 484件 [詳細] filecluster_6.png 459件 [詳細] filecluster_5.png 486件 [詳細] filecluster_4.png 441件 [詳細] filecluster_3.png 449件 [詳細] filecluster_2.png 413件 [詳細] filecluster_1.png 450件 [詳細] fileCompareTPM.pdf 419件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2020-03-02 (月) 11:56:09 (1122d)