[[ノート>ノート/ノート]]~
訪問者数 &counter();      最終更新 &lastmod();

*RNAseq〜岸本データの後続分析 [#gcc5a764]

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

***データ [#sfe06552]
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を取出してみると

||	gene|	Start|	Length|	Anc|
|gene_3114|	arpB|	1786475|	474|	7.176203468|
|gene_3115|	arpB|	1786959|	1416|	1.441322222|
|gene_3782|	gatR|	2153495|	447|	6.087732338|
|gene_3786|	gatR|	2155197|	339|	163.2194667|
|gene_2015|	icd|	1172515|	1251|	1782.60336|
|gene_2058|	icd|	1190066|	165|	12.36916525|
|gene_6818|	ilvG|	3929665|	984|	1955.874255|
|gene_6819|	ilvG|	3930728|	582|	942.1393394|
|gene_0495|	insA|	 290296|	276|	0|
|gene_1774|	insA|	1040161|	276|	6.572986365|
|gene_6204|	insA|	3561389|	276|	0.821623296|
|gene_7778|	insA|	4498297|	276|	37.7946716|
|gene_0035|	insB|	  19812|	504|	3.149555967|
|gene_0475|	insB|	 278403|	504|	0|
|gene_0494|	insB|	 289874|	504|	2.6996194|
|gene_1775|	insB|	1040355|	504|	64.34092903|
|gene_7637|	insB|	4413646|	504|	3.599492533|
|gene_7779|	insB|	4498491|	294|	30.85279314|
|gene_7781|	insB|	4498785|	210|	78.82888648|
|gene_0647|	insC|	 381730|	411|	0|
|gene_1727|	insC|	1016493|	411|	0|
|gene_2221|	insC|	1280354|	411|	0|
|gene_5485|	insC|	3164049|	411|	1.103494061|
|gene_7744|	insC|	4477997|	411|	0|
|gene_0648|	insD|	 382080|	924|	0|
|gene_1728|	insD|	1016861|	906|	0|
|gene_2222|	insD|	1280722|	906|	0|
|gene_2529|	insD|	1450103|	906|	0|
|gene_3625|	insD|	2052755|	906|	0|
|gene_5486|	insD|	3164417|	906|	0|
|gene_7745|	insD|	4478347|	924|	2.6996194|
|gene_0535|	insE|	 315706|	309|	0|
|gene_0536|	insE|	 315715|	300|	0|
|gene_2021|	insE|	1177320|	309|	0|
|gene_3785|	insE|	2154038|	300|	0|
|gene_0467|	insH|	 273326|	981|	0|
|gene_0505|	insH|	 294457|	981|	0|
|gene_0973|	insH|	 566399|	981|	0|
|gene_1163|	insH|	 678257|	981|	0|
|gene_1843|	insH|	1081651|	1017|	0|
|gene_2218|	insH|	1278220|	981|	0|
|gene_2394|	insH|	1377059|	1017|	0.222977414|
|gene_2401|	insH|	1380609|	1017|	0.222977414|
|gene_2466|	insH|	1409928|	981|	16.18120497|
|gene_3011|	insH|	1725750|	981|	0|
|gene_3316|	insH|	1893378|	981|	0|
|gene_3623|	insH|	2050108|	1017|	28.31813152|
|gene_3690|	insH|	2085698|	1017|	0|
|gene_3792|	insH|	2158195|	1017|	5.351457925|
|gene_3983|	insH|	2274059|	1017|	0|
|gene_5382|	insH|	3108121|	981|	0|
|gene_5806|	insH|	3343608|	981|	0|
|gene_6318|	insH|	3630088|	1017|	0|
|gene_6713|	insH|	3868728|	1017|	0|
|gene_0463|	insI|	 269828|	1152|	16.53516882|
|gene_2530|	insI|	1451540|	1152|	1.181083487|
|gene_7761|	insI|	4487236|	1152|	11.22029313|
|gene_0028|	insL|	  15440|	1119|	0|
|gene_1028|	insL|	 598319|	1119|	0|
|gene_4345|	insL|	2499206|	1119|0|
|gene_0462|	insN|	 269467|	405|	9.518658032|
|gene_7760|	insN|	4486967|	267|	40.76728622|
|gene_0464|	insO|	 271055|	426|	117.11025|
|gene_7764|	insO|	4488728|	597|	20.89152701|
|gene_6038|	kefG|	3458512|	555|	0|
|gene_6039|	kefG|	3458512|	552|	0|
|gene_2164|	ldrB|	1247821|	135|	0|
|gene_2167|	ldrB|	1248356|	135|	68.87029047|
|gene_2169|	ldrB|	1248891|	135|	0|
|gene_2465|	lomR|	1409640|	156|	15.99005337|
|gene_2469|	lomR|	1410996|	171|	1.326128828|
|gene_3836|	molR|	2181468|	825|	3.573314406|
|gene_3838|	molR|	2182404|	1938|	8.658841171|
|gene_3840|	molR|	2184304|	960|	4.251900555|
|gene_7426|	phnE|	4301437|	369|	4.916380045|
|gene_7427|	phnE|	4301655|	621|	5.842654547|
|gene_1997|	potA|	1161851|	1137|	0|
|gene_1998|	potA|	1161851|	1119|	0|
|gene_6273|	rhsB|	3597098|	4236|	3.051411163|
|gene_6498|	rhsB|	3746869|	294|	13.11243709|
|gene_1713|	sulA|	1009402|	516|	0|
|gene_1714|	sulA|	1009402|	510|	0|
|gene_3689|	wbbL|	2085199|	282|	99.71360166|
|gene_3691|	wbbL|	2086719|	474|	189.4517716|
|gene_1830|	ycdN|	1071684|	111|	4.085910443|
|gene_1832|	ycdN|	1071794|	720|	11.65335708|
|gene_2077|	ycgH|	1198254|	1521|	32.05465244|
|gene_2078|	ycgH|	1199859|	1017|	105.0223618|
|gene_2212|	ychG|	1273079|	591|	44.50946097|
|gene_2213|	ychG|	1273621|	231|	54.97406778|
|gene_2284|	yciX|	1316143|	129|	54.49464277|
|gene_2285|	yciX|	1316253|	189|	278.3607559|
|gene_2528|	ydbA|	1447574|	2559|	13.20376569|
|gene_2532|	ydbA|	1452872|	3324|	120.0697148|
|gene_3511|	yedN|	1995151|	192|	1.181083487|
|gene_3512|	yedN|	1995352|	321|	0.70644246|
|gene_3570|	yedS|	2017854|	486|	5.132609723|
|gene_3572|	yedS|	2018349|	210|	16.1977164|
|gene_3574|	yedS|	2018642|	405|	19.03731606|
|gene_3597|	yeeL|	2036079|	327|	3.467401064|
|gene_3599|	yeeL|	2036427|	705|	0.964970339|
|gene_4039|	yfaS|	2314679|	486|	15.39782917|
|gene_4041|	yfaS|	2315180|	4104|	11.65888261|
|gene_4169|	yfcC|	2401962|	1521|	0.298182813|
|gene_4171|	yfcC|	2402004|	1479|	0|
|gene_7747|	yjgX|	4479813|	432|	58.26678538|
|gene_7748|	yjgX|	4480202|	360|	30.23573728|
|gene_0530|	ykgM|	312798|	141|	3697.444681|
|gene_0531|	ykgM|	312938|	267|	2342.420321|
|gene_0940|	ylbE|	548781|	1260|	0.179974627|
|gene_0941|	ylbE|	549039|	1002|	0|
|gene_2630|	yncI|	1512768|	747|	21.55358782|
|gene_2632|	yncI|	1513558|	201|	15.79478813|
|gene_2579|	yncK|	1485325|	168|	31.0456231|
|gene_2580|	yncK|	1485544|	288|	37.7946716|
|gene_2084|	ypjA|	1202348|	213|	5.323193183|
|gene_4790|	ypjA|	2756055|	4500|	23.98701824|

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

内容を見ると(岸本先生による)、同じ配列がコピーされている場合と、1つの配列が途中に余分な夾雑物が入って切れているように見える場合がある。さらに移動の途中で短くなったり長くなったりしているケースがある。これらの判断は難しいが、
-他の資料(EcoCycなど)から「元」が分かるとき、部分Xと部分Yがそれぞれの一部であり、結合するとほぼ「元」と同じになる場合、
--更に、部分XとYが数塩基〜数アミノ酸離れているだけ、つまり数塩基の夾雑物が入っているだけであるとみられる場合、
--間隔はかなり長いので、議論の信頼度は落ちるが、配列で見ると繋がる場合
-同じ配列がコピーされている場合 〜 長さがほぼ同じで配列もほぼ同じ
-1つは他資料からほぼフルセットと思われるが、もう一方は短い断片だと思われる場合
-よく分からない場合
に分けられる。

それぞれ、発現量について次のような対応をすることにする。
-同じ配列がコピーされている場合
--一方の部分配列の発現値が0の場合 ⇒ 発現値0の部分配列は削除してよいだろう。幸い、サンプルに依らずこのパターンになっているので、発現値0の部分配列は機能していないと思われる。
--それぞれの部分配列が発現している場合 ⇒ それぞれの部分配列からアミノ酸〜タンパク質が作られていると思われるので、発現量(TPM補正済み)は加えてしまう。
--それぞれの部分配列が発現している場合 ⇒ それぞれの部分配列からアミノ酸〜タンパク質が作られていると思われるので、発現量(TPM補正後の)は加えてしまう。
-遺伝子が近い距離で分割されている場合
--もともと一連の遺伝子と思われ、この形でも(発現量があるので)タンパク質が作られると思われる。
--もともと一連の遺伝子と思われ、この形でも(発現量があるので)タンパク質が作られると思われる。両者の発現量(TPM補正後の)の平均を取ってみる。(もしかすると、小さい側/最小値を取るべきかもしれない)部分配列の長さによる発現量への依存性はTPM補正により正規化されるので、単純に平均してよかろう。
-一方がフルセットで他方が断片と思われるとき
--断片の側を削除し(結果を捨て)、フルセット側だけを残してみる。

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

***やりたいこと: [#r40d87d3]

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

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

***例外的パターンについて [#l5fc7858]

0 - ... - 0 - 0 - 0 つまり一切発現無し、があり得る。これはおそらく分析対象から除外してよいだろう。
 
0 - 途中にnon-0がある - ... これらのケースは、例外視する必要はなく、1つのパターンと考えてよかろう。Ancが0であるケースは、「なかったものが出てくるようになった」という意味で少し考える必要があるかも知れないが、とりあえず1つのパターンとして考えることにする。

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

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


***変動の数値指標について [#s91230ce]

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

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

***途中の道草 [#t554de2c]
今はlog(tpm/Anc)を取った表。

|	|gene	|Start	|Length	|Anc	|43B	|45a_minus	|45A_minus	|45L	|1_2-1	|2_5-1|
|gene_0001	|thrL	|190	|66	|0	|0.518801371	|0.120822012	|0.130553092	|0.662278199	|0.736981706	|0.782980429|
|gene_0003	|thrA	|338	|2463	|0	|0.363545812	|0.217293869	|0.207090028	|0.462360214	|0.457629592	|0.265097062|
|gene_0005	|thrB	|2802	|933	|0	|0.375643354	|0.338447083	|0.313037592	|0.533665319	|0.35669746	|0.17628472|
|gene_0007	|thrC	|3735	|1287	|0	|0.351946181	|0.178209151	|0.160556929	|0.449396654	|0.309469379	|0.092155761|
|gene_0008	|yaaX	|5235	|297	|0	|0.361804821	|0.497406189	|0.565856219	|0.483361196	|0.254979922	|0.141397537|
|gene_0009	|yaaA	|5684	|777	|0	|0.127790066	|-0.16041128	|-0.046404433	|0.015105221	|0.04435137	|-0.083863423|
|gene_0011	|yaaJ	|6530	|1431	|0	|-0.061913911	|0.120802403	|0.215477533	|0.132366925	|0.108788741	|-0.028453563|
|gene_0014	|talB	|8239	|954	|0	|0.025617922	|-0.296533334	|-0.408331467	|0.077236367	|0.241795445	|0.203333867|
|gene_0016	|mog	|9307	|588	|0	|-0.18900972	|-0.26035833	|-0.255735167	|-0.114590533	|-0.140329364	|-0.159416034|
|gene_0017	|yaaH	|9929	|567	|0	|-0.151212522	|-0.346290801	|-0.274411814	|-0.177810069	|-0.161112427	|-0.103483574|

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

&ref(CompareTPM.pdf);

**全体をクラスタ化 [#y7bd23d4]

 %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()

結論は、
|&ref(./cluster_1.png,250x180);|&ref(./cluster_2.png,250x180);|&ref(./cluster_3.png,250x180);|&ref(./cluster_4.png,250x180);|&ref(./cluster_5.png,250x180);|
|&ref(./cluster_6.png,250x180);|&ref(./cluster_7.png,250x180);|&ref(./cluster_8.png,250x180);|&ref(./cluster_9.png,250x180);|&ref(./cluster_10.png,250x180);|

-Clusterは作れる (時間はかかるが)
-距離はベクトル間のユークリッド距離を使った。
-できたClusterの妥当性がよくわからない。見た目はクラスター内の比較で全く違う。
-特に気になるのが0の扱いで、グラフ上はほぼ-6(log(0)は-∞だが代わりに0.000001を入れてあるため)として描かれている点である。これの動きが非常に気になる。

***おまけ [#n4bcb1c4]
必須遺伝子の発現パターンはどうなっているのか?

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS