NGSデータのマッピングとデータ処理のための汎用スクリプト

次世代シーケンシング（NGS）データの処理は、ゲノム解析の重要なステップです。ここでは、Minimap2とSAMtoolsを使用してNGSデータのマッピングとデータ処理を行うための汎用スクリプトを紹介します。

背景

NGS技術は、短時間で大量のデータを生成するため、データ処理の効率化が求められます。Minimap2は、高速かつ正確なリードマッピングツールとして広く使用されており、SAMtoolsはSAM/BAMファイルの操作に便利なツールです。本スクリプトでは、これらのツールを組み合わせて、効率的なデータ処理を実現します。

使用ツール

本スクリプトでは以下のツールを使用します：

Minimap2 - 高速かつ正確なリードマッピングツール
SAMtools - SAM/BAMファイルの操作ツール

これらのツールは、以下のコマンドでConda経由でインストールできます：

$ conda create -n mapping_env
$ conda activate mapping_env
$ conda install -c bioconda minimap2 samtools

スクリプト

以下に、Minimap2とSAMtoolsを使用した汎用スクリプトを示します。

#!/bin/bash

# ログファイルの設定
LOGFILE="path/to/mapping.log"

# Minimap2を使用したマッピングとSAMtoolsの処理をバックグラウンドで実行し、ログファイルに記録
{
    echo "Indexing reference genome with minimap2"
    minimap2 -d path/to/consensus.mmi path/to/consensus.fasta

    echo "Mapping reads with minimap2"
    minimap2 -ax sr path/to/consensus.mmi path/to/sample_R1.fastq.gz path/to/sample_R2.fastq.gz > path/to/mapped.sam

    echo "Converting SAM to BAM"
    samtools view -S -b path/to/mapped.sam > path/to/mapped.bam

    echo "Sorting and indexing BAM"
    samtools sort path/to/mapped.bam -o path/to/sorted.bam
    samtools index path/to/sorted.bam

    echo "Getting read counts for BAM"
    samtools idxstats path/to/sorted.bam > path/to/idxstats.txt

    echo "All tasks completed"
} &> $LOGFILE &

echo "Script is running in the background. Check the log file for details: $LOGFILE"

スクリプトの簡単な説明と使用方法

このスクリプトは、以下の手順でNGSデータのマッピングとデータ処理を行います：

参照ゲノムのインデックス作成（Minimap2）
リードデータの参照ゲノムへのマッピング（Minimap2）
SAMファイルをBAMファイルに変換（SAMtools）
BAMファイルのソートとインデックス作成（SAMtools）
BAMファイルのリードカウント取得（SAMtools）

スクリプトを実行するには、必要なファイルパスを適切に置き換えた後、以下のコマンドを実行します：

$ bash path/to/run_mapping.sh

スクリプトはバックグラウンドで実行され、進捗状況はログファイルに記録されます。

出力ファイルの例と見方

スクリプトの実行後、以下の出力ファイルが生成されます：

mapped.sam - マッピング結果のSAMファイル
mapped.bam - SAMファイルをBAM形式に変換したもの
sorted.bam - ソートおよびインデックス付けされたBAMファイル
idxstats.txt - BAMファイルのリードカウント情報

例えば、idxstats.txt の内容は以下のようになります：

seq1    50000   12345   6789
seq2    60000   23456   7890
...

各行は以下の情報を含みます：

seq1 - シーケンス名
50000 - シーケンスの長さ
12345 - マッピングされたリード数
6789 - マッピングされなかったリード数

このようにして、マッピングおよびデータ処理結果を詳細に確認することができます。

まとめ

本記事では、NGSデータのマッピングとデータ処理のための汎用スクリプトを紹介しました。このスクリプトを利用することで、Minimap2とSAMtoolsを効果的に活用し、効率的なデータ解析が可能となります。是非、あなたのプロジェクトにも活用してみてください。