【研究】RNA-seq03：カウント編v20240305.0

2024年3月5日 20:01

こんにちは、あるいは、こんばんは。
執筆作業を終え、会議のための資料作成、会議の案内、次の執筆の企画？のようなことをしながら、ターミナルと夜な夜な格闘した成果をご報告します。

投稿した記事の順番が前後しておりますが（記事名でソートして表示させる方法がわからないです）、今回はマッピングしたデータのカウント作業を行いました。
正直、あっているのかどうか既にわからないですが、どこかの誰かの何かの参考にでもなれば幸いです。

PCスペック

MacBook Pro 14-inch, 2021
チップ：Apple M1 Max
コア：10（パフォーマンス 8, 効率性： 2）
メモリ：64 GB
macOS：macOS Sonoma 14.3

このMacについて

Rosetta 2でターミナル.appを動かしております。

今回の目標

前回までに、解析に使用するデータのマッピング処理まで済ませました。
　●　【研究】RNA-seq01：準備編
　●　【研究】RNA-seq01.5：納品データ確認編
　●　【研究】RNA-seq02：マッピング編：　データ使用します

今回はいよいよ統計解析に適用させるための処理を行います！

解析環境については、試行錯誤の結果、（繰り返しになるかもしれませんが）再構築の方法から掲載します。
もちろん、読み飛ばしていただいても大丈夫です！！

環境構築（再）

冒頭で、Rosetta 2でターミナル.appを動かしております。と記載しておりましたが、ナンノコッチャ？　という場合には、過去の記事をご確認下さい。
（えっ！？　そこは　飛ばすんかい！！）
お怒りはごもっともですが、ご容赦下さい。

Rosetta 2でターミナル.appを使用する

ターミナルアプリを起動し、以下のコードを入力して、「i386」と返ってくれば、Intel CPUを使用していることと同じ挙動で実行できます。

arch
# > i386

condaをインストール

最近、外付けのSDXC（1Tと名付けました！）にデータを格納しています。
ただ、読み書きする速度が遅いためか、処理に時間がかかります。
いずれはSSDに移行する予定です。

今回も、「Mambaforge3」を再導入します。
過去の環境が残っているようでしたら、使用してもしなくても大丈夫です。
調子が悪いようでしたら、いっそのことデリりましょう。
Rosetta2から起動したターミナルに以下のコードを入力します。

mambaforgeの再構築の準備（まず削除）

# デリる（削除する）
rm -rf ~/mambaforge/
vi ~/.bashrc

表示される画面にて「dd」「dd」…と押していき、記入された情報を「行」単位で削除していきます。
全て消し去ったら、「:wq」（コロン→w→qの順）です。
最後に以下のコードを実行してリセットが完了するっぽいです。

source ~/.bashrc

mambaforgeの再導入

ワーキングディレクトリは、好きな場所で大丈夫です。
わたしは、外付けSDXC 1T「/Volumes/1TB/RNAseq/03」というフォルダに、Mambaforge3.shというファイルを保存し、実行します。

# ワーキングディレクトリの変更: ここにダウンロードされます
cd /Volumes/1TB/RNAseq/03
# GitHubから「Mambaforge3.sh」としてダウンロードします
wget -O Mambaforge3.sh https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-MacOSX-x86_64.sh
# インストールします
bash Mambaforge3.sh -b -u
# パスを通します
echo ". ~/mambaforge/etc/profile.d/conda.sh" >> ~/.bashrc && \
echo "conda activate base" >> ~/.bashrc
source ~/.bashrc

もし、「wget」がないなどのエラーがありましたら、以前の記事を参考にしてHomebrewを導入していただければと思います。

chnannel設定

Bioconda推奨の設定をしておきます。以下のコードを順番に実行します。
デフォルトで設定されているようですので、省略可能と思われます。

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strict

ツールのインストール

一連の作業で使用していたツールをインストールします。
おそらく、「subread」のみで良いのかと思います。バージョンによりコードが機能しない不具合もございました。
執筆時点（2024/3/5）の内容で記載させていただきます。

## update mamba
conda install -c conda-forge mamba -y
## install
mamba install -c bioconda seqkit -y
mamba install -c bioconda trimmomatic -y
mamba install -c bioconda hisat2 -y
mamba install -c bioconda samtools -y
mamba install -c bioconda subread=1.5.2 -y   # バージョン指定しています

subreadというのは、featureCountによるカウントを行うために必要なものです。最新バージョンではエラーが起こり、解決が難しかったため「1.5.2」を指定しています。

カウント

featureCounts（subread=1.5.2）

featureCountsを用いてマッピングとソートされたBAMファイルから、リードカウントの情報を取得していきます。
まずアノテーション情報のデータを準備しておきます。
格納するディレクトリは、好きな場所で大丈夫です。
今回は「cd /Volumes/1TB/RNAseq/03」に格納しています。
以下のコードを実行します。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/03
# アノテーション情報のデータのダウンロード
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.gtf.gz
# 解凍
unpigz -d GCF_000001635.27_GRCm39_genomic.gtf.gz

アノテーション情報と前回マッピング・ソートしたデータを使用してカウント情報を整理します。
以下のコードを実行します（約20分／18データ）。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata
# アノテーション情報のデータのパス
ANNOTATION=/Volumes/1TB/RNAseq/03/GCF_000001635.27_GRCm39_genomic.gtf
# 出力ファイル
OUTPUT_FILE=featureCounts.txt
# featureCountの実行
featureCounts -F GTF -t CDS -g gene_id -a ${ANNOTATION} -o ${OUTPUT_FILE} *.sorted.bam

「ANNOTATION=…」：上述のダウンロード・解凍したGTFファイルです。
「*.sorted.bam」：ワイルドカードにより、ワーキングディレクトリに含まれるファイルを読み込むようにしました。ファイルが増えて判別しづらくなっているようでしたら、別のフォルダ「bamfiles」などを作成して整理しておくと良いかもしれません。

以下のデータが出力されます。

featureCounts.txt
featureCounts.txt.summary

featureCount.txt

featureCount.txtを解析ソフトに読み込み、処理することになります。
テキストデータになってしまえば、ナントカナリソウですね？！
前処理の段階で不具合がないか心配なところもありますが、まず、やってみるのが重要と信じています！

更新

20240305.0：初版

いいなと思ったら応援しよう！

この記事が参加している募集

#仕事について話そう

115,889件