【pytximportの使い方】転写産物ごとの発現量データを遺伝子ごとの発現量に変換する

RNA-seq
ぴよこ
kallistoで解析したデータをその後の解析にどうやって使えばいいのかわからない。調べたらtximportが便利みたいだけど、R使ったことないからできなさそう。。。
まさる博士
そんな君に朗報！tximportのpythonバージョンであるpytximportがあるよ

RNAseq解析においてアライメントフリーの発現定量ツールであるkallistoの出力データは転写産物ごとの発現量であるため、遺伝子ごとの発現量を解析したい場合には、データを遺伝子ごとに変換する必要があります。

本記事では、pytximportというツールを用いて、kallistoの出力データを遺伝子ごとの発現量に変換する方法について解説します。

kallistoって何？という方は、以下のリンクをご覧ください。

関連記事
【kallistoの使い方】アラインメントフリーで遺伝子発現量を計算する

より高速な発現定量データ解析手法として、発現定量を主眼とした、ゲノムにアラインメントしない方法（アラインメントフリー）が広まりつつあります。 アラインメントフリーの手法では、これまでに知られたすべてのトランスクリプトーム配列を使って、[…]

 

目次

pytximportとは

pytximportの前に、tximportについて説明します。

tximportとは

tximportは、RNA-Seqデータの解析において、転写産物の定量化をサポートする R パッケージです。特に、転写産物ごとの定量化データを遺伝子ごとのカウントに変換するために使用されます。

https://bioconductor.org/packages/release/bioc/html/tximport.html

以下に tximport の主な機能と目的について説明します。

主な機能

  1. トランスクリプトから遺伝子へのマッピング: tximport は、転写産物（トランスクリプト）の定量化結果を遺伝子レベルにマッピングする機能を提供します。これにより、トランスクリプトのカウントデータを遺伝子のカウントデータに変換できます。
  2. 複数のフォーマットに対応: tximport は、異なるツールやフォーマットから出力された定量化データを読み込み、統一された形式で出力します。これには、Salmon、Kallistoなどのツールからのデータが含まれます。
  3. 統合的なデータ処理: tximport を使用することで、転写産物レベルから遺伝子レベルにデータを変換する過程を簡略化し、さらにそのデータを DESeq2 や edgeR などの下流解析ツールで使用できる形式に整えます。

以上のように、tximportはRNAseqの解析において重要な役割を果たしますが、使用するためにはR プログラミング言語が必要となります。

そのため、私のようなpythonは使えるけど、Rはほとんど使ったことない、という人には少し抵抗があり、pythonから使えるものがあればなー、と思っていました。

そんな人々の悩みを解決するものが、今回紹介するpytximportです。

pytximportとは

pytximport は、tximportのpythonバージョンであり、tximportの機能を python 、あるいはコマンドラインから実行できるようにしてくれています。

下記のリンクから、pytximportの公式ドキュメントにアクセスできます。

pytximportの公式ドキュメントはこちら

pytximportのインストール

pipでインストールできます。
pip install pytximport
インストールできれば、準備完了です。

pytximportの使い方

ここから、pytximportを用いて、kallistoの出力形式であるabundance.h5から、遺伝子ごとの発現量に変換する手順を説明していきます。

pytximportは、以下の2種類の方法で使用することができます。

  1. pythonでimportして使用する
  2. コマンドラインで実行する

まず、①の方法から解説していきます。

①pythonでimportして使用する

まず、必要なライブラリをインストールします。

import numpy as np
import pandas as pd
from pytximport import tximport

トランスクリプトと遺伝子の対応表を作成する

遺伝子ごとの発現量に変換するためには、転写産物（トランスクリプト）と遺伝子名を対応づける必要があります。その方法として、BioMartを用いる方法を下記の記事で解説しました。

関連記事
【BioMartの使い方】geneIDとtranscriptIDの対応表を作成する

ぴよこ RNA-seqのデータを解析したけど、結果が転写産物IDごとになってる。遺伝子ごとに変換したいんだけど、どうやってやればいいの？ まさる博士 それなら、EnsemblのBioMartを使うといいよ！[…]

なんと、pytximportではこの処理をコードで簡単に実行することができます。

以下のコードで、Ensemblのデータをもとに、ヒト転写産物とヒト遺伝子名の対応表をpandasのデータフレーム形式として作成できます。

from pytximport.utils import create_transcript_to_gene_map

transcript_gene_map_human = create_transcript_to_gene_map(species="human", field="external_gene_name")
transcript_gene_map_human.head(5)

出力結果

create_transcript_to_gene_mapの入力は以下のとおりです。

  • species: “human” or “mouse”。デフォルトはhuman。
  • field: “ensembl_gene_id”, “external_gene_name”, “external_transcript_name”。もし遺伝子IDや転写産物名との対応を作成したい場合は、この値を変更する。

ファイルとして出力する場合は、以下のコードでtsvファイル（タブ区切りファイル）として保存してください。

transcript_gene_map_human.to_csv("transcript_gene_map_human.tsv", index=False, sep=\t)

 

kallistoのデータを遺伝子ごとの発現量に変換する

つぎに、kallistoのデータを遺伝子ごとの発現量に変換します。使用するデータは、以下の2つのデータがあるとします。

  • ./data1/abundance.h5
  • ./data2/abundance.h5

kallistoの出力データがない場合は、下記の記事を参考に作成してみてください。

関連記事
【kallistoの使い方】アラインメントフリーで遺伝子発現量を計算する

より高速な発現定量データ解析手法として、発現定量を主眼とした、ゲノムにアラインメントしない方法（アラインメントフリー）が広まりつつあります。 アラインメントフリーの手法では、これまでに知られたすべてのトランスクリプトーム配列を使って、[…]

それと、先ほど作成したtranscriptIDと遺伝子名の対応データ(transcript_gene_map_human)を使用します。

下記のコードを実行します。

txi = tximport(
file_paths=["./data1/abundance.h5", "./data2/abundance.h5"],
data_type="kallisto",
transcript_gene_map=transcript_gene_map_human,
counts_from_abundance=None,
output_type="anndata",
save_path="gene_count.csv",
save_path_overwrite=True
)

この、コードを実行するとgene_count.csvに遺伝子ごとに2つのサンプルの発現量が変換された出力されます。

gene_count.csvの中身

tximportの入力は以下の通りです。

  • file_paths: 入力となる定量化ファイルへのパス。ここではkallistoの出力であるabundance.h5を指定している。他にもsalmonのquant.sfなどにも対応している。
  • data_type: kallisto, salmonなど、定量化ファイルの種類。デフォルトはsalmon。
  • transcript_gene_map: transcriptと遺伝子の対応。pd.DataFrameあるいはファイルパスを指定する。transcript_idとgene_idの2つの列が必要。
  • counts_from_abundance: abundanceに基づいてカウント推定値を計算するかどうか。デフォルトはNone。
  • output_type: 出力データの種類。anndataあるいはxarray。デフォルトはanndata。
  • save_path: 遺伝子ごとの発現量を保存するファイルパス。
  • save_path_overwrite: 保存パスがすでに存在する場合に上書きするかどうか。

 

②コマンドラインで実行する

コマンドラインから実行するには、ターミナルで以下のコマンドを実行します。

pytximport -i ./data1/abundance.h5 -i ./data2/abundance.h5 -t "kallisto" -m "transcript_gene_map_human.tsv" --output_type="anndata" --save_path_overwrite -o gene_count_cli.csv

*転写産物と遺伝子名の対応表のファイルは、csv形式だとエラーとなったので、tsvの形式で保存しておくようにしてください。

おわりに

この記事では、pytximportを用いてkallistoの出力データを遺伝子ごとの発現量に変換する方法を解析しました。

主にpythonを使っている人には、ありがたいツールです。

 

RNA-seqの最新記事8件