【研究】RNA-seq：前処理編まとめv20240307.1
こんにちは、あるいは、こんばんは。
RNA-seqの生データをどのように処理してから、データ解析に移れば良いのか悩まれる方も多いと思います。
これまでに試行してきました内容を整理した記事をメモとして残しておきたいと思います。
いつの日か、どこかのだれかの参考にでもなれば幸いです。
注意
● わたしは、プログラミングに疎いため、無駄な処理や美しくないコード
も多く含まれているかと思います。
● もし、改善点がございましたら、そっと、やさしく、コメントしていた
だけますと幸いです。
● すぐコメント確認と修正対応ができないかも知れませんが、真摯に受け
止め、精進してまいります。
PCスペック
Rosetta 2でターミナル.appを動かします。
設定方法は以前の記事を参考にしていただけますと幸いです。
検索すると、もっと分かりやすい情報があると思います。
環境構築
Rosetta 2でターミナル.appを使用するよ！
ターミナルアプリを起動し、以下のコードを入力して、「i386」と返ってくれば、Intel CPUを使用していることと同じ挙動で実行できます。
archi386フォルダの準備
処理過程で必要となるファイルを格納するためのフォルダを作成しておきます。
わたしのフォルダの階層構造は以下のとおりです。
## Homebrewでtreeをインストール
brew install tree
## treeで階層構造を表示
tree /Volumes/4T/RNAseq
/Volumes/4T/RNAseq
└── rawdata
├── Sample01_1.fastq.gz
├── Sample01_2.fastq.gz
...省略しました...
├── Sample18_1.fastq.gz
├── Sample18.fastq.gz
├── JN########_18samples_md5sum_DownloadLink.txt
├── md5sum.txt
└── sample_list.txt
2 directories, 39 files
「/Volumes/4T/RNAseq」に以下の4つのフォルダを作成します。
● downloads
● trimmed
● bamfiles
● featureCounts
以下のコードを実行します。
# ディレクトリのリスト
directories=("downloads" "trimmed" "bamfiles" "featureCounts")
# ディレクトリの存在をチェックして再帰的に作成する関数
check_and_create_directory() {
if [ ! -d "$1" ]; then
echo "Creating directory: $1"
mkdir -p "$1"
else
echo "Directory already exists: $1. Nothing to do."
fi
}
# リスト内の各ディレクトリに対して処理を行う
for dir in "${directories[@]}"; do
check_and_create_directory "$dir"
done
# ディレクトリの内容を再帰的に表示する
tree -L 1 /VolumeCreating directory: downloads
Creating directory: trimmed
Creating directory: bamfiles
Creating directory: featureCounts
/Volumes/4T/RNAseq
├── bamfiles
├── downloads
├── featureCounts
├── rawdata
└── trimmed
6 directories, 0 files「rawdata」と同じ階層にフォルダを作成することができました。
ちなみに。。。
もう一度、フォルダを作成しようとすると以下のエラーコメントが返ってきます。
Directory already exists: downloads. Nothing to do.
Directory already exists: trimmed. Nothing to do.
Directory already exists: bamfiles. Nothing to do.
Directory already exists: featureCounts. Nothing to do.conda（mamba）をインストール
今回、外付けのSSD（4Tと名付けました！）にデータを格納しています。
以前のSDXC 1TBよりも安定して高速に処理できることを願っています。
というのも、途中で処理が止まってしまったり、明らかに処理が遅くなっていたりすることがありました。
それが原因かはっきりしませんが、出来上がった「featureCount.txt」が他の研究者と一致していませんでした。
以下は前回の記事と同様なので読み飛ばしていただいて大丈夫です。
さっそく、「Mambaforge3」を導入します。
過去の環境が残っているようでしたら、使用してもしなくても大丈夫です。
調子が悪いようでしたら、いっそのことデリりましょう。
Rosetta2から起動したターミナルに以下のコードを入力します。
mambaforgeの再構築の準備（まず削除）
# デリる（削除する）
rm -rf ~/mambaforge/
vi ~/.bashrc表示される画面にて「dd」「dd」…と押していき、記入された情報を「行」単位で削除していきます。
全て消し去ったら、「:wq」（コロン→w→qの順）です。
最後に以下のコードを実行してリセットが完了するっぽいです。
source ~/.bashrcmambaforgeの再導入
ワーキングディレクトリは、好きな場所で大丈夫です。
今回は、外付けSSD 4T「/Volumes/4TB/RNAseq/downloads」というフォルダに、Mambaforge3.shというファイルを保存し、実行します。
# ワーキングディレクトリの変更: ここにダウンロードされます
cd /Volumes/4T/RNAseq/downloads
# 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を導入していただければと思います。
（オプション）仮想環境の構築
Anacondaの仮想環境を構築しておくと、不具合があったときに削除したり、別の仮想環境で実行することができると思います。
以下のコードを実行します。
# 仮想環境の構築
conda create -n newenv
# アクティベートすることで、仮想環境で実行できます
conda activate rnaseqプロンプトの冒頭が
（base）→（newenv）になっていれば
仮想環境がアクティベートされています。
次回から、「conda activate newenv」を実行しましょう！
仮想環境名（newenv）は、独自に設定できますので、プロジェクト名にしておくと良いかも知れません（例：rnaseq）
chnannel設定
Bioconda推奨の設定をしておきます。以下のコードを順番に実行します。
デフォルトで設定されているようですので、省略可能と思われます。
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --set channel_priority strictツールのインストール
一連の作業で使用するツールをインストールします。
# mambaのアップデート
conda install -c conda-forge mamba -y
# インストール（subreadのみバージョンを指定しておきます）
mamba install pigz -y
mamba install -c bioconda seqkit trimmomatic hisat2 samtools subread=1.5.2 -ysubreadというのは、featureCountによるカウントを行うために必要なものです。最新バージョンではエラーが起こり、解決が難しかったため「1.5.2」を指定しています。
解析前に確認！（20240307追記）
作業中断後に再解析する場合
上述の環境構築を行い、後述の解析を進めるうちに、作業の中断やPCの再起動を余儀なくされた場合には、以下のコードを実行してから作業を再開しましょう。
source ~/.bashrc
conda activate newenvコマンド（パッケージ）がないとエラーが表示される場合
処理中に必要となるツールがインストールされていない場合には、上述の仮想環境をアクティベートしてから、インストールし直しましょう。
どうしても作業の再開ができない場合
そんなこともあります！
原因解明に時間を取られても良くないと思いますので、思い切ってデリ（削除）ってください（上述）。
もうちょっとソフトにデリる（削除する）のであれば、conda仮想環境（newenv）を削除して、再構築しましょう。
その場合は、以下のコードをお試し下さい。
# アクティベート中のconda仮想環境の中止
conda deactivate
# 仮想環境（newenvの場合）を削除
mamba remove -n newenv --allconda仮想環境の構築（上述）から再開しましょう。
FASTAの簡易解析
seqkitを用いて、FASTAの簡易的な解析を行います。
以下のコードを実行します。
# ワーキングディレクトリ
cd /Volumes/4T/RNAseq/rawdata
# seqkitの実行
seqkit stats *.fastq.gzオプションも様々あり、SeqKitでは色々な操作ができるようです。詳細はSeqKitの公式ページをご確認ください。
-a、--all：すべての統計情報を出力
-b、--basename：塩基名のみを出力
-E、--fq-encoding string：fastqの品質エンコーディングを指定
-G、--gap-letters string：ギャップを表す文字を指定
-h、--help：統計情報に関するヘルプを表示
-e、--skip-err：エラーをスキップし、警告メッセージのみ表示
-S、--skip-file-check：指定ファイルのチェックをスキップ
-T、--tabular：タブ形式で出力 など
以下が表示されます。36データの解析に約30分かかりました。
念の為、コピーしてテキストデータとして保存しておきます。
トリミング
Trimmomatic
Trimmomaticを用いてアダプター配列のトリミングを行います。
アダプター配列を指定して処理を進めることになります。
ますは必要となるファイルのダウンロードです。
以下のコードを実行します。
# ワーキングディレクトリの変更
cd /Volumes/4T/RNAseq/rawdata/downloads
# ダウンロードと解凍
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unpigz Trimmomatic-0.39.zipArchive: Trimmomatic-0.39.zip
creating: Trimmomatic-0.39/
inflating: Trimmomatic-0.39/LICENSE
inflating: Trimmomatic-0.39/trimmomatic-0.39.jar
creating: Trimmomatic-0.39/adapters/
inflating: Trimmomatic-0.39/adapters/NexteraPE-PE.fa
inflating: Trimmomatic-0.39/adapters/TruSeq2-PE.fa
inflating: Trimmomatic-0.39/adapters/TruSeq2-SE.fa
inflating: Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa
inflating: Trimmomatic-0.39/adapters/TruSeq3-PE.fa
inflating: Trimmomatic-0.39/adapters/TruSeq3-SE.fa 解凍後に幾つかファイルが格納されていることが確認できます。
Trimmomaticを実行します（約5時間／36データ）。
# ワーキングディレクトリの変更
cd /Volumes/4T/RNAseq/rawdata
# サンプルリストとアダプター配列情報のパス
SAMPLE_LIST_FILE=sample_list.txt
ILLUMINACLIP_FILE=../downloads/Trimmomatic-0.39/adapters/TruSeq2-PE.fa:2:30:10
## trimmomaticの実行
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
echo "Processing sample: $SAMPLE"
trimmomatic PE -threads 12 \
-trimlog ../trimmed/${SAMPLE}.txt \
${SAMPLE}_1.fastq.gz ${SAMPLE}_2.fastq.gz \
../trimmed/paired_${SAMPLE}_1.trim.fastq.gz \
../trimmed/unpaired_${SAMPLE}_1.trim.fastq.gz \
../trimmed/paired_${SAMPLE}_2.trim.fastq.gz \
../trimmed/unpaired_${SAMPLE}_2.trim.fastq.gz \
ILLUMINACLIP:ILLUMINACLIP_FILE
echo "Sample $SAMPLE processing complete."
doneマッピング
HISAT2
HISAT2を用いてリードデータをマッピングします。
まずリファレンスゲノム配列のインデックスを準備しておきます。
作成するディレクトリは、リードデータと同じ場所でも大丈夫です。
今回は「/Volumes/4T/RNAseq/downloads」に格納しています。
以下のコードを実行します（約40分）。
# ワーキングディレクトリの変更
cd "/Volumes/4T/RNAseq/downloads"
# リファレンスゲノム配列のダウンロードと解凍
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.fna.gz
unpigz -d GCF_000001635.27_GRCm39_genomic.fna.gz
## インデックスの作成
hisat2-build GCF_000001635.27_GRCm39_genomic.fna GCF_000001635.27_GRCm39_genomic_index -p 12以下の8つのインデックスが出力されます。
# インデックス
index_GRCm39_genomic.1.ht2
index_GRCm39_genomic.2.ht2
index_GRCm39_genomic.3.ht2
index_GRCm39_genomic.4.ht2
index_GRCm39_genomic.5.ht2
index_GRCm39_genomic.6.ht2
index_GRCm39_genomic.7.ht2
index_GRCm39_genomic.8.ht2
作成したインデックスを使用して、リードデータをマッピングします。
以下のコードを実行します（時間がかかるためO/Nか日を改めて）。
# ワーキングディレクトリの変更
cd /Volumes/4T/RNAseq/rawdata
# サンプルリストとインデックス情報のパス
SAMPLE_LIST_FILE=sample_list.txt
INDEX=../downloads/GCF_000001635.27_GRCm39_genomic_index
# HISAT2（マッピング）の実行
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
echo "Processing sample: $SAMPLE"
FILE1=../trimmed/paired_${SAMPLE}_1.trim.fastq.gz
FILE2=../trimmed/paired_${SAMPLE}_2.trim.fastq.gz
hisat2 -p 1 -x ${INDEX} -1 ${FILE1} -2 ${FILE2} \
| samtools sort -@ 1 > ../bamfiles/${SAMPLE}_sorted.bam
samtools index -@ 1 ../bamfiles/${SAMPLE}_sorted.bam
echo "Sample $SAMPLE processing complete."
doneカウント
featureCounts（subread=1.5.2）
featureCountsを用いてマッピングとソートされたBAMファイルから、リードカウントの情報を取得していきます。
まずアノテーション情報のデータを準備しておきます。
格納するディレクトリは、好きな場所で大丈夫です。
今回は「/Volumes/4T/RNAseq/downloads」に格納しています。
以下のコードを実行します。
# ワーキングディレクトリの変更
cd /Volumes/4T/RNAseq/downloads
# アノテーション情報のデータのダウンロードと解凍
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/4T/RNAseq/rawdata
# アノテーション情報と出力ファイルのパス
ANNOTATION=../downloads/GCF_000001635.27_GRCm39_genomic.gtf
OUTPUT_FILE=../featureCounts/featureCounts.txt
# featureCountの実行
featureCounts -F GTF -t CDS -g gene_id -a ${ANNOTATION} -o ${OUTPUT_FILE} ../bamfiles/*_sorted.bam「ANNOTATION=…」：上述のダウンロード・解凍したGTFファイルです。
「*_sorted.bam」：ワイルドカードにより、ワーキングディレクトリに含まれるファイルを読み込むようにしました。
以下のデータが出力されます。
featureCount.txt
featureCount.txtを解析ソフト（エクセル、R、Python、SPSSなど）
に読み込み、処理することになります。
まとめ
いかがでしょうか。
コードを実行中に修正を加えながら記載しておりましたので、noteもデータ前処理もミスが起こっている可能性がございます。
記事公開後に、最終産物「featureCount.txt」を確認して、動作状況と解析再現性などを評価したいと思います。
初版執筆時点（2024/3/7 16:22:22）では解析途中です。。。
更新
20240307.0：初版
20240307.1：「解析前に確認！」追記
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？