読者です 読者をやめる 読者になる 読者になる

NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

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