遺伝子発現情報の取得

前回は、GTFファイルから遺伝子ID情報を集約する流れまでした。
今回はいよいよ本丸の遺伝子発現情報を取り扱っていきます。
まずは前回ダウンロードしたファイル（1）を読み込み、とりあえず今回の解析では不要な「transcript_ids」列を削ったデータフレームを作成します。

# TPM data after RSEM(count data after mapping) of CCLE(1019 cell lines)
df3 = pd.read_table("CCLE_RNAseq_rsem_genes_tpm_20180929.txt")
df3 = df3.drop("transcript_ids", axis=1) # remove transcript_ids
print(df3.shape)
df3.head()

# print(df3.shape)
(57820, 1020)

カラム名に関して、一列目が「gene_id」、いわゆる「Ensembl ID」で、二列目以降のカラム名は「癌細胞の名前＋癌種」です。そして二列目以降の数値が遺伝子の発現情報（TPMカウントデータ）となっています。
これを見てわかる通り、今回のデータは57820遺伝子、1019細胞のTPMカウントデータです。

TPMデータでは、各細胞の相対的な発現量を比較しやすいように、各細胞の遺伝子発現量の総カウント数が100万カウントに統一されています。一応、最初の10細胞に関してそれを確認してみましょう。

df0 = pd.DataFrame(columns=["TPM"])
for i in range(10):
    df0_2 = pd.DataFrame([df3.iloc[:,(i+1)].sum().astype(int)], index= [df3.columns[i+1]], columns=["TPM"])
    df0 = pd.concat([df0, df0_2])
df0

多少のブレはありますが、およそどの細胞も100万カウントになっていることがわかりますね。

癌細胞名に関する情報の取得

次に、癌細胞名のアノテーションを行い、辞書を作ります。
まずは遺伝子発現情報データdf3から細胞名を抽出します。

df3_2= pd.DataFrame(df3.columns[1:], columns = ["CCLE_ID"])
df3_2.head()

(オプショナル：CCLEのアノテーションファイルから癌細胞の情報を取得)

この項はオプショナルですので次の本番まで読み飛ばしてもらっても一応大丈夫です。
まず前回ダウンロードしたファイル（2）を読み込みます。このファイルには癌細胞の由来等の情報が記載されています。

df_n = pd.read_table("Cell_lines_annotations_20181226.txt")
print(df_n.shape)
print(df_n.columns)
df_n.head()

(1461, 33)
Index(['CCLE_ID', 'depMapID', 'Name', 'Pathology', 'Site_Primary',
       'Site_Subtype1', 'Site_Subtype2', 'Site_Subtype3', 'Histology',
       'Hist_Subtype1', 'Hist_Subtype2', 'Hist_Subtype3', 'Gender',
       'Life_Stage', 'Age', 'Race', 'Geo_Loc', 'inferred_ethnicity',
       'Site_Of_Finding', 'Disease', 'Annotation_Source',
       'Original.Source.of.Cell.Line', 'Characteristics', 'Growth.Medium',
       'Supplements', 'Freezing.Medium', 'Doubling.Time.from.Vendor',
       'Doubling.Time.Calculated.hrs', 'type', 'type_refined',
       'PATHOLOGIST_ANNOTATION', 'mutRate', 'tcga_code'],
      dtype='object')

ここには1461細胞、33種類の情報が記載されています。お気づきの方はいると思いますが、細胞の数が1019ではありません。ただしそれは間違っているわけではなく、「遺伝子発現」情報のデータが記載されている癌細胞が1019個だということです。前回のブログで、CCLEのダウンロードページには遺伝子発現情報の他にも、miRNAやRPPAといったデータも記載されていたことを覚えてらっしゃるでしょうか？

CCLEはヒト癌細胞を網羅的に解析しようとする試みであるため、それぞれの解析で多少重複しない細胞があったとしてもおかしくありません。では実際に、今回使用する細胞の情報をここから抽出してみましょう。

df3_3 = pd.merge(df3_2, df_n, on = "CCLE_ID", how ="inner")
print(df3_3.shape)
# (1019, 33)
df3_3.describe(exclude='number').iloc[[0,1],[0,2,4,26,27,28]]

データをマージすると、ちゃんと細胞数が1019個になったことが確認できます。ちなみについでに抽出しているのは、今後の解析に使用しそうな情報です。欲しいのは癌種情報なのですが、下記の例を見てもらえればわかるように、それぞれの視点から情報が記載されているため、実は一意に決まるわけではありません。
(*細かいことですがdescribeで指定した列と下記で指定した列が異なるのは、describeで「非数要素」の結果を対象に列指定しているためです)

df3_3.iloc[:,[0,2,4,28,29,30]].head(10)

当初はなんとなく、'PATHOLOGIST_ANNOTATION'の方が良いと思っていたのですが、例えばName行の重複要素を5個ほど見てみると、'PATHOLOGIST_ANNOTATION'どころか、'Site_Primary'や'type'といった項でも要素が'NaN'になっています！これでは後で癌種分類できない！

df3_3[df3_3["Name"].duplicated()].iloc[:,[0,2,4,28,29,30]].head()

ということで、方針変更して'CCLE_ID'から癌種を取り出すことにします。
というか実は当初からそうしてたのですが、今回はちゃんと検証するために一応ファイル（2）の中身も見てみました(^-^)。前回のブログで 「一応」 ファイル（２）をダウンロードした、というのはこういう意味です（つまり結局は使わなかったので・・・）。

本番：CCLEの遺伝子発現データから直接癌細胞の情報を取得

ここからが本番となります。
df3_2 をよく見ると、'CCLE_ID'は、'細胞名' ＋ ' _ ' ＋ '癌種'という形で文字が結合されていることがわかります。しかし、'癌種'の中身も' _ 'で更に区切られてたりするからややこしい（例：CENTRAL_NERVOUS_SYSTEM）。区切り文字別にしてくれたほうが良いのに。。。

ということで、最初の' _ 'だけ区切って後は残すような形で細胞名と癌種名を分けてみます。

# Separation of cell line name and cancer type
df3_4 = df3_2.copy()
df3_4 = df3_4.iloc[:,0].str.split("_", n=1, expand=True)
df3_4 = df3_4.rename(columns={0:"cell_name", 1:"cancer_type"})
print(df3_4.shape)
print(df3_4["cancer_type"].nunique())
df3_4.head()

# df3_4.shape
(1019, 2)
# cancer_typeの種類
26

ちゃんと細胞名と癌種が分かれましたね。全てを合算すると漏れも無さそうです。これを見るとCCLEでは上位2種類の癌がかなり多いですね。

print(df3_4.iloc[:,1].value_counts().sum())
# 1019
pd.DataFrame(df3_4.iloc[:,1].value_counts())

それではこれから癌種分類するための辞書を作ります。Scikit-learnのLabelEncoderを使ってもいいんですが、後のことを考えると辞書対応表が欲しかったので、ゴリゴリつくっちゃいました。ちゃんと26種類の癌にそれぞれ分類番号(0〜25)が付いていることがわかります。

# 癌の種類名
ct_name= list(df3_4["cancer_type"].unique())
# 分類番号
ct_num = list(range(df3_4.nunique()[1]))
#辞書の作成(making dictionary)
dic = dict(zip(ct_name, ct_num))
dic

# 癌分類対応表
{'PROSTATE': 0,
 'STOMACH': 1,
 'URINARY_TRACT': 2,
 'CENTRAL_NERVOUS_SYSTEM': 3,
 'OVARY': 4,
 'HAEMATOPOIETIC_AND_LYMPHOID_TISSUE': 5,
 'KIDNEY': 6,
 'THYROID': 7,
 'SKIN': 8,
 'SOFT_TISSUE': 9,
 'SALIVARY_GLAND': 10,
 'LUNG': 11,
 'BONE': 12,
 'PLEURA': 13,
 'ENDOMETRIUM': 14,
 'PANCREAS': 15,
 'BREAST': 16,
 'UPPER_AERODIGESTIVE_TRACT': 17,
 'LARGE_INTESTINE': 18,
 'AUTONOMIC_GANGLIA': 19,
 'OESOPHAGUS': 20,
 'FIBROBLAST': 21,
 'CERVIX': 22,
 'LIVER': 23,
 'BILIARY_TRACT': 24,
 'SMALL_INTESTINE': 25}

続いて、上記辞書を用いて1019個の癌細胞を癌細胞ごとにラベル分けしたデータを作成します。最初の１０行のデータを見ると、辞書で指定したように、それぞれの癌細胞がcategoryの列で癌種ごとに数値で分類されていることがわかりますね。

df3_5 = df3_4.copy()
df3_5["category"] = df3_5.iloc[:,1]
df3_5["category"] = df3_5["category"].replace(dic)
df3_5 = pd.concat([df3_2, df3_5], axis =1)
print(df3_5.shape)
# (1019, 4)
df3_5.to_csv("cancer_type_20190527.csv")
df3_5.head(10)

とりあえずデータフレームが増えてごちゃごちゃしてきたので、「cancer_type_20190527.csv」という形で保存しちゃいます。jupyter notebookも重くなってきたのでちょっと仕切りなおしです。
ということで、次回のブログに続きます・・・（たぶん近日公開予定）。