こんにちは、あるいは、こんばんは。
執筆作業、会議資料作成、会議もろもろ、ようやく一時のオフィス籠もりタイムができましたので、格闘した格闘中の成果をご報告します。

（2024/2/9追記, 2/29追記）

PCスペック

Apple silicon M1でなんとか動かそうと奮闘しておりましたが、

でした。使用できないツールが多く、なんとか動かそうと努力しましたが、対応していないので無理なようです。
Rosetta 2でターミナル.appを動かしております。

今回の目標

前回までに、解析に使用するデータの確認まで済ませました。
　●　【研究】RNA-seq01：準備編　：　今回、環境リセットします
　●　【研究】RNA-seq01.5：納品データ確認編　：　データ使用します

今回は実データを用いて解析の前処理をし始めます！

解析環境については、今後、使用する可能性がありますが、基本的にリセットして一から再構築します。

ちなみに。。。
Intel CPUのMacBook Proを持っておりますが、そちらでも動くと思います。
（途中まで動作を確認しております）

環境構築

今までに構築した環境とは異なる環境を使用していきます。
もし、これまでにいろいろなプログラムをご使用の場合には、バックアップを取るなどしてから実施してください。

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

arch
# > i386

i386と返ってくれば、Intel CPUを使用していることと同じ挙動でターミナルを使用できます。

もし、Rosettaなんて見つからない！という場合には

これまで、Rossetaを使用する機会がなかった場合には、チェックボックスがないかもしれません（詳しくはわからないです）
そのような時には、以下のコードを実行してみてください。
（おまじないです）

softwareupdate --install-rosetta

Condaをインストール

今回は、「Mambaforge3」というものを導入してみます。
Rosetta2から起動したターミナルに以下のコードを入力します。
ワーキングディレクトリは、好きな場所で良いと思います。
今回、dropboxではなく、外付けSSDまたはSDXCを使用しています。
少しずつ変更がありますが、ご容赦ください。

# ワーキングディレクトリの変更: ここにダウンロードされます
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
# パスを通します
echo ". ~/mambaforge/etc/profile.d/conda.sh" >> ~/.bashrc && \
echo "conda activate base" >> ~/.bashrc
source ~/.bashrc

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

Chnannel設定

Biocondaが推奨するchannel設定をしておきます。
以下のコードを順番に実行します。

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

ツールのインストール

解析の前処理などに使用するツールをインストールします。
コンフリクトも起こるようで、わたしには解決ができません。
もし、不具合が発生するようでしたら、環境をリセットし、その時に必要なツールのみインストールすると良いと思います。

mamba install -c bioconda seqkit -y                    # 追記2/29
mamba install -c bioconda bowtie2 -y
mamba install -c bioconda sra-tools -y
mamba install -c bioconda samtools -y
mamba install -c bioconda fastp -y
mamba install -c bioconda fastqc -y
mamba install -c bioconda fastq-screen -y
mamba install -c bioconda hisat2 -y
mamba install -c bioconda stringtie -y
mamba install -c conda-forge pigz -y
mamba install libtiff pillow simplejson matplotlib -y  # <- multiqcで必要
mamba install -c bioconda multiqc -y                   # 実施しました2/9
mamba install -c bioconda igv -y                       # 起動確認2/9
mamba install -c bioconda igvtools -y                  # 起動確認2/9

multiqcを実施してみました。関連パッケージのインストールが必要になるとのことでした。（2024/2/9追記）
seqkitを実施してみました。（2024/2/29追記）

データの配置

データの配置は皆様それぞれと思いますが、コードを実行する前に、ワーキングディレクトリを確認すると良いと思います。
何も考えず、シンプルに動かす場合には、全てワーキングディレクトリに保存してしまいましょう。
また、必要なデータやファイルがワーキングディレクトリと別に保存されている場合には、その都度、指定するようにお願いします。

データ保存場所のイメージ

FASTAの簡易解析（2024/2/29追記）

seqkitを用いて、FASTAの簡易的な解析を行います。
以下のコードを実行します。

seqkit stats /Volumes/1TB/RNAseq/rawdata/*.fastq.gz

以下が表示されます。72データの解析に約2時間30分かかりました。
Fastpでトリムした後のデータも同じパスに含まれていたため、時間がかかってしまいました。
念の為、コピーしてテキストデータとして保存しておきます。

トリミング

Fastp

Fastpを用いてアダプター配列などのトリミングを行います。
深く理解して処理を進めているわけではないため、別にTrimmomaticでの処理も検討いたします（2024/2/29追記）。

まず、Fastpで処理を行います。
以下のコードを実行します。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata
# サンプルリストファイルのパス
SAMPLE_LIST_FILE="sample_list.txt"
## サンプルリストを読み込んで処理するループ
## NR > 1: 1行目(header)を無視する
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
    echo "Processing sample: $SAMPLE"
    
    # fastpの実行
    fastp --detect_adapter_for_pe \
        -i "${SAMPLE}_1.fastq.gz" \
        -I "${SAMPLE}_2.fastq.gz" \
        -o "${SAMPLE}_1.trim.fastq.gz" \
        -O "${SAMPLE}_2.trim.fastq.gz" \
        -h "${SAMPLE}_report.html" -w 1
    
    echo "Sample $SAMPLE processing complete."
done

以下のデータ・レポートが出力されます。

fastp report (.html) の例

Trimmomatic（2024/2/29追記）

Trimmomaticも実行してみました。
本法ではアダプター配列を指定して処理を進めることになります。
以下のコードを実行します。

# アダプター配列情報のダウンロード
cd /Users/user/Downloads
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip

# Trimmomaticのインストール
mamba install -c bioconda trimmomatic -y

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata

# サンプルリストファイルのパス
SAMPLE_LIST_FILE="sample_list.txt"

## サンプルリストを読み込んで処理するループ
## NR > 1: 1行目(header)を無視する
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
    echo "Processing sample: $SAMPLE"

    # trimmomaticの実行
    trimmomatic PE -threads 10 -trimlog ${SAMPLE}.txt \
         ${SAMPLE}_1.fastq.gz ${SAMPLE}_2.fastq.gz \
         paired_${SAMPLE}_1.trim.fastq.gz unpaired_${SAMPLE}_1.trim.fastq.gz \
         paired_${SAMPLE}_2.trim.fastq.gz unpaired_${SAMPLE}_2.trim.fastq.gz \
         ILLUMINACLIP:/Users/user/Downloads/Trimmomatic-0.39/adapters/TruSeq2-PE.fa:2:30:10

    echo "Sample $SAMPLE processing complete."
done

試しに1サンプルのみ処理してみましたが、約1時間40分かかりました。
以下のデータが出力されます。

この処理をループですべて実行するには、、、、約30時間（18データ）かかることになります意外にも約5時間で終わりました！
後述のマッピング処理では、コードの内容に「paired_」を加えて処理すれば良さそうですね。

（2024/2/29追記ここまで）

クオリティチェック（QC）

FastQC

FastQCを用いてデータの質（クオリティ）の確認をします。
以下のコードを実行します。

fastqc --nogroup -o . *.trim.fastq.gz

以下のデータ・レポートが出力されます。

fastqc report (.html) の例

この続きは実行中のため、後日追記いたします。
ようやく解析のスタートラインが見えてきました。

（2024/2/7 20:50）

MultiQC

MultitQCを用いてデータの質（クオリティ）の確認をします。
「/Volumes/1TB/RNAseq/rawdata」に格納されている「fastp report」と「fastqc report」を認識して、まとめて提示してくれます。
以下のコードを実行します。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata
# MultiQC実行
multiqc .

以下のフォルダ・レポートが出力されます。

multiqc report (.html) の例

マッピング

HISAT2

HISAT2を用いてリードデータのマッピングします。
まずリファレンスゲノム配列のインデックスを準備しておきます。
作成するディレクトリは、リードデータと同じ場所でも大丈夫です。
今回は「cd /Volumes/1TB/RNAseq/03」に格納しています。
以下のコードを実行します（約40分）。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/03
# リファレンスゲノム配列のダウンロード
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/10090/109/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.fna.gz
# 解凍
unpigz -d GCF_000001635.27_GRCm39_genomic.fna.gz
# index作成(約40分)
hisat2-build GCF_000001635.27_GRCm39_genomic.fna index_GRCm39_genomic

以下の8つのインデックスが出力されます。

作成したインデックスを使用して、リードデータをマッピングします。
以下のコードを実行します（時間がかかるためO/N：19:00〜）。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata
# サンプルリストファイルのパス
SAMPLE_LIST_FILE="sample_list.txt"
## サンプルリストを読み込んで処理するループ
## NR > 1: 1行目(header)を無視する
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
    echo "Processing sample: $SAMPLE"
    
    # hisat2 (mapping)の実行
    FILE1=${SAMPLE}_1.trim.fastq.gz
    FILE2=${SAMPLE}_2.trim.fastq.gz
    INDEX=/Volumes/1TB/RNAseq/03/index_GRCm39_genomic
    hisat2 -p 1 -x ${INDEX} -1 ${FILE1} -2 ${FILE2} \
        | samtools sort -@ 1 > ${SAMPLE}.sorted.bam
done

以下のデータ「.bam」が出力されます。

bamファイルのインデックス

bamファイルの情報を検索するための索引となります。
以下のコードを実行します。

# ワーキングディレクトリ
cd /Volumes/1TB/RNAseq/rawdata
# サンプルリストファイルのパス
SAMPLE_LIST_FILE="sample_list.txt"
## サンプルリストを読み込んで処理するループ
## NR > 1: 1行目(header)を無視する
awk 'NR > 1 {print $2}' "$SAMPLE_LIST_FILE" | while IFS= read -r SAMPLE; do
    echo "Processing sample: $SAMPLE"
    
    # SAMtools (index)の実行
    samtools index -@ 1 ${SAMPLE}.sorted.bam
done 

以下のデータ「.bai」が出力されます。

BAMとBAI

BAM（Binary Alignment Map）とBAI（Binary Alignment Index）は、DNAやRNAの配列データを保存するファイル形式
通常、BAMとBAIはペアで使用し、アラインメントデータの効率的な操作と解析に使用する

⚫︎BAM

• リード配列（DNA・RNA断片）のアラインメント情報のバイナリ

• リード配列、クオリティスコア、アラインメント位置、アラインメント方向などの情報を含む

⚫︎BAI

• BAMに含まれる情報を高速に検索するための索引（インデックス）

• BAMファイル内のリード位置情報を含む

• 特定のリードやリード範囲をBAMファイルから効率的に取得できる

マッピング結果の確認

IGV（Integrative Genomics Viewer）

IGVを用いて、リードデータの可視化を行います。
以下のコードを実行します。

ivg

IVGが起動しました。
試しにgenome情報をダウンロードしてみました。
Genomes > Select Hosted Genome… > 「Mouse mm39」

BAMとBAIを読み込ませると良いのだと思いますが、上手くいきませんでした。勉強不足です。いつの日か、更新したいと思います。

更新