本稿では、発現量の正規化（FPKM/RPKM/TPM）を行ったRNA-Seqデータに対して、低発現遺伝子の除外・Fold change・統計計算を行う手法およびPythonによる実装を解説する。

RNA-Seqデータ解析のワークフロー

RNA-Seq (バルク) のデータ解析は以下の流れで進む (参考 https://sc.ddbj.nig.ac.jp/advanced_guides/Rhelixa_RNAseq/ 2022年6月18日アクセス) (ステップ6以降は研究者により異なる場合が多い)。

生データのクオリティーチェック (FastQC)→クオリティ情報に基づいたデータのトリミング (Trimimomatic)
リファレンスゲノムへのマッピング (Hisat2やSTARなど) および bam/sam変換 (Samtool)
リード数のカウント (featureCount)
発現量の正規化 (自力計算)
低発現遺伝子の除外・Fold changeの算出・統計計算 (自力計算)
PCA・UMAPや階層的クラスタリング・ヒートマップによるデータ概要の確認 (PythonのUMAP, clustermapモジュール)
GO解析・パスウェイ解析などの高次解析 (自由, GSEAやDAVID, metascapeなどが有名)

筆者はステップ1-3については国立遺伝学研究所 (National Institute of Genetics) のスーパーコンピュータシステムにおける Rhelixa RNAseqパイプラン (https://sc.ddbj.nig.ac.jp/advanced_guides/Rhelixa_RNAseq/ 2022年6月18日アクセス) を使用させていただき、以降はPythonを用い、自力で解析を進めた。本稿では、ステップ5についてPythonを用いて実行する方法を解説する。

使用データ

本稿では、先行研究の公開データを用いてデモ解析を行う。
生後１ヶ月、２ヶ月、４ヶ月齢の雌雄両方のマウス由来の海馬を用いたRNA-Seq解析データを利用する (Bundy et al. BMC Genet. 2017, #GSE83931)。

Age	Sex	Sample Number
1 month	Male	5
1 month	Female	5
2 month	Male	5
2 month	Female	5
4 month	Male	5
4 month	Male	5

こちらの公開データのXlsxファイル（GSE83931_Counts_and_FPKM.xlsx）の中の「FPKM」シートのデータを解析に使用しているが、"Associated.Gene.Name"と"Description"のカラムは削除してからUMAP解析に使用した。
データ例:

GeneID	1071D	1071E	・・・	1056G	1056H
Gene_A	1.2	1.1	・・・	2.2	3.2
Gene_B	1.0	0.9	・・・	1.1	1.1
・・・	・・・	・・・	・・・	・・・	・・・
Gene_Y	321.4	123.4	・・・	234.5	345.6
Gene_Z	0.0	0.0	・・・	0.0	0.0

低発現遺伝子の除外

元データはリファレンスゲノム情報に載っている全ての遺伝子の発現量情報を有しているので、37,315遺伝子分の発現量情報が載っている。
実は低発現遺伝子の除外は実行しなくても良いが、実行するメリットは:

データが軽くなるので、以降の解析の際に取り扱いがしやすくなる。
ヒートマップを作成する際に、発現量が低い遺伝子に引っ張られてしまうことがなくなる。
Benjamini-Hochberg法など、P値を遺伝子数で補正してQ値を算出する場合に、余計な遺伝子が残っていると厳しくなりすぎてしまう。
絶対的な差は小さいにもかかわらず、fold changeが大きく見える現象（下表）を未然に防げる。

GeneID	Average of GroupA (TPM)	Average of GroupB (TPM)	Log2 Fold change
Gene X	0.2	0.1	1.0
Gene Y	250	125	1.0

Gene Yの方がグループA-Bの間で表現系に影響を与えるような差がありそうだが、Fold changeだけを見ると、どちらも同じである。

しかし、先述のように低発現遺伝子の除外は必ずしも実行する必要があるわけではなく、本稿で使用する公開データにおいても、実は低発現遺伝子が残った形で公開されている。そのため、低発現遺伝子の除外方法（基準）は割と研究室ごとの「自家製」であることが多いらしい。本稿では最も簡便な方法の一つである、「すべてのサンプルにおいて発現量が一定以下の場合に除外」という基準を用いる。今回は、FPKM < 1.0の遺伝子を除外することにする (適当に設定した値なので、各自自分の基準を設定されてほしい)。

Pythonでの実装

使用モジュールおよびデータの読み込み。本稿では上記公開データのFPKM値データを使用している。

RNASeq_Analysis.py

import numpy as np
import pandas as pd
import math
from scipy import stats
import scipy as sp
from statsmodels.stats.multitest import multipletests

#データ読み込み
file_path = 'データパス'
rsdf = pd.read_excel(file_path, header=0, index_col=None)
print('Number of the genes >', len(rsdf))

まず下記辞書配列を定義して、生物学的グループ、サンプル名の関係を定義する。下記は各自のデータによって可変なので留意してほしい。

RNASeq_Analysis.py

# 生物学的グループの設定 (手動でやる)
bio_groups_FPKM = \
{'1month_male':['1071D','1071E','1071G','1071H','1071I'], \
'1month_female':['1071P','1071Q','1071X','1072B','1071T'], \
#省略
}

下記コードでは、ある遺伝子が除外の基準を満たすか判定する関数を定義している。実際に解析を行う際には各自の基準に沿った関数を定義してほしい。

RNASeq_Analysis.py

# あるひとつのサンプルでもFPKMが1を超えていたら、その遺伝子は除外されない。
# 除外される時は、Trueを返す。
def drop_judge(gene_index):
    judge = True
    for g in bio_groups_FPKM:
        for s in bio_groups_FPKM[g]:
            if rsdf.at[gene_index, s] > 1.0:
                judge = False
                break
        else:
            continue
        break
    return judge

下記コードでは実際にデータの除外を実行し、除外された遺伝子名をテキストファイルにして出力している。

RNASeq_Analysis.py

# 除外される遺伝子のインデックスを取得し、drop_list配列に格納
drop_list = []
for i in rsdf.index:
    if drop_judge(i): drop_list.append(i)
print('Number of the genes eliminated >', len(drop_list))

# drop_listの遺伝子を元データから除外。
rsdf_drop = rsdf.copy()
rsdf_drop.drop(drop_list, inplace=True)

#残った遺伝子をエクセルファイルにして出力
rsdf_drop.to_excel('List_of_Genes_remained.xlsx')
print('Number of the genes remained >', len(rsdf_drop))

# 除外された遺伝子のリストをテキストファイル形式で出力。
with open('Eliminated_Gene_List_Old_RNASeq.txt', 'w') as f:
    for i in drop_list:
        f.write(rsdf.at[i, 'Name']+'\n')

Log2倍率変化 (fold change) の計算

続いて下記コードでは、Log2 (fold change) の計算の準備を進めていく。まず、生物学的グループごとの発現量の平均を算出する。

RNASeq_Analysis.py

# 生物学的グループごとの発現量の平均の算出を行い、エクセルファイルとして出力
for i in rsdf_drop.index:
    for g in bio_groups_FPKM:
        sum = 0
        for s in bio_groups_FPKM[g]: sum += rsdf_drop.at[i, s]
        rsdf_drop.at[i, g+'_av'] = sum / len(bio_groups_FPKM[g])
rsdf_drop.to_excel('Eliminated_Average.xlsx')

Log2 (fold change)の計算を行なっていく。今回はすべての生物学的グループ間の Log2 (fold change) の計算を行う。本稿で使用している公開データの例では、6生物学的グループあるので、その組み合わせ、₆C₂ = 15 とおりの計算を行う。

RNASeq_Analysis.py

# 組み合わせの取得
import itertools
groups = bio_groups_FPKM.keys()
comparison = list(itertools.combinations(groups, 2))

# Log2 fold changeの計算
for c in comparison:
    for i in rsdf_drop.index:
        rsdf_drop.at[i, 'Log2FC_'+c[0]+'_'+c[1]] = math.log2(rsdf_drop.at[i, c[0]+'_av']/rsdf_drop.at[i, c[1]+'_av'])
rsdf_drop.to_excel('Eliminated_Average_Log2FoldChange.xlsx')

統計計算 (Student's t-test・ Benjamini-Hochberg法)

続いてStudent's t-testのP valueの計算および、Benjamini-Hochberg法による多重比較の補正を行なっていく。Benjamini-Hochberg法によるQ valueの算出はRNA-Seq解析においてメジャーな手法の一つであるが、詳細については別記事を参照されてほしい。

RNASeq_Analysis.py

# Student's t-testの実行
for c in comparison:
    for i in rsdf_drop.index:
        # 2つの生物学的グループの発現量のリストを作成
        group_a = []; group_b = []
        for s in bio_groups_FPKM[c[0]]:
            group_a.append(rsdf_drop.at[i,s])
        for s in bio_groups_FPKM[c[1]]:
            group_b.append(rsdf_drop.at[i,s])

        # 実際にStudent's t-testを行う部分
        st, stp = stats.ttest_ind(group_a, group_b)
        rsdf_drop.at[i, 'StuP_'+c[0]+'_'+c[1]] = stp
rsdf_drop.to_excel('Eliminated_Average_Log2FoldChange_StuP.xlsx')

最後に Q valueの計算を行う。その前にP valueの列に欠損値があると計算が進まないので、"1" で穴埋めする。

RNASeq_Analysis.py

# 欠損値の穴埋め
for c in comparison:
    rsdf_drop.fillna({'StuP_'+c[0]+'_'+c[1]:1}, inplace=True)

# Q valueの算出
for c in comparison:
    pvalues = rsdf_drop['StuP_'+c[0]+'_'+c[1]].to_numpy()
    # Benjamini-Hochberg Q valueの算出を実際に行う部分
    rjarray, corrp, aS, aB = multipletests(pvals=pvalues, alpha=0.05, method='fdr_bh', is_sorted=False, returnsorted=False)
    rsdf_drop['BHQ_'+c[0]+'_'+c[1]] = corrp
rsdf_drop.to_excel('Eliminated_Average_Log2FoldChange_StuP_BHQ.xlsx')

おわりに

本稿で公開したコードはほとんどそのままコピペして利用可能であるが、(1) RNA-Seqデータのファイルパス、(2) 生物学的グループ、サンプル名の関係を定義した辞書配列の設定、の2点は各自でやらないといけないことに留意してほしい。