MiSeqデータのMHC領域へのマッピング（5）：bwa, samtools indexまでを同時に行うスクリプト

2015-03-05

記事(3)-(4)の作業を一度に行うshell scriptを作りました。カレントディレクトリにあるすべての.fastqファイルに対して、ヒトChromosome 6へのマッピングを行い、.sort.bamファイル、index fileを作成します。

使い方：.fastqファイルおよびreferenceのindexがあるディレクトリで以下のスクリプト（mapping_mhc.sh）をを実行する。結果はmapping_resultsフォルダ内の各.fastq file名のフォルダ内に作られる。

（追記150412）スクリプトは以下のものをコピーしてテキストエディタ（CotEditorなど）に貼り付け、"mapping_mhc.sh"という名前で保存してください。

sh mapping_mhc.sh

#!/usr/bin/sh
# mapping_mhc.sh

# MHC領域のMiSeqデータ（DRA000908）をhuman genome referenceにmappingする
# bwa, samtools, bcftools, vcfutils.plを使用

# current directoryに存在するすべての.fastq fileについて行う
# paired end非対応

mkdir mapping_results	# 解析結果を入れるdirectory

for fname in *.fastq; do

# 拡張子 (.fastq)を除く
fname=`echo "$fname" | sed -E 's/\.fastq//g'`

mkdir ./mapping_results/$fname	# 個別のデータについての解析結果 directory

# mapping using bwa mem
bwa mem -t 4 Homo_sapiens.GRCh38.dna_rm.chromosome.6.fa $fname.fastq > ./mapping_results/$fname/$fname.sam

# samtools view .sam -> .bam
samtools view -bS ./mapping_results/$fname/$fname.sam > ./mapping_results/$fname/$fname.bam

# samtools sort
samtools sort ./mapping_results/$fname/$fname.bam ./mapping_results/$fname/$fname.sort

# samtools index
samtools index ./mapping_results/$fname/$fname.sort.bam

done