ゲノム解析でエラーが出ているため解決策を教えていただきたいです。(bwa samtools)

Question

前提・実現したいこと

初心者でどうしたらいいのかわからず詰まっております。

マッピングしたく、以下を実行したところ

#!/bin/bash
set -euo pipefail
bwa=bwa-0.7.17/bwa
id=DRR002191
fq1=${id}_1.fastq.gz
fq2=${id}_2.fastq.gz
ref=RefHg38/hg38.fasta
rg="@RG\tID:${id}\tSM:${id}\tPL:illumina\tLB:${id}"

${bwa} mem
-R ${rg}
${ref}
${fq1} ${fq2}
| samtools view -@4 -b -1 - > ${id}.bam

以下のエラーメッセージが発生しました。

発生している問題・エラーメッセージ

[M::bwa_idx_load_from_disk] read 0 ALT contigs
[main_samview] fail to read the header from "-".

試したこと

samtoolsのインストールしなおしたりしました。

補足情報（FW/ツールのバージョンなど）

samtoolsの情報は以下です。
Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10.2)

Usage: samtools <command> [options]

Commands:
-- Indexing
dict create a sequence dictionary file
faidx index/extract FASTA
fqidx index/extract FASTQ
index index alignment

-- Editing
calmd recalculate MD/NM tags and '=' bases
fixmate fix mate information
reheader replace BAM header
targetcut cut fosmid regions (for fosmid pool only)
addreplacerg adds or replaces RG tags
markdup mark duplicates

-- File operations
collate shuffle and group alignments by name
cat concatenate BAMs
merge merge sorted alignments
mpileup multi-way pileup
sort sort alignment file
split splits a file by read group
quickcheck quickly check if SAM/BAM/CRAM file appears intact
fastq converts a BAM to a FASTQ
fasta converts a BAM to a FASTA

-- Statistics
bedcov read depth per BED region
coverage alignment depth and percent coverage
depth compute the depth
flagstat simple stats
idxstats BAM index stats
phase phase heterozygotes
stats generate stats (former bamcheck)

-- Viewing
flags explain BAM flags
tview text alignment viewer
view SAM<->BAM<->CRAM conversion
depad convert padded BAM to unpadded BAM

Comments

[otn]

Linux OSについての質問と言うより、bwa(?)というツールの使い方の質問のようなので、適切なタグを付けた方が、有識者の目にとまると思います。

[tare]

アドバイスありがとうございます！　bwaはタグ付けできなかったので、タイトルに入れました！

Answer 1

bwaは正しい出力を生成してるんでしょうか？
samtoolsに渡す前のbwaの出力をファイルに出力して内容をチェックしてみるのが先決だと思います。

Comments

[tare]

bwaは出力できているとおもうのですが、この先どうしたらいいかわからず、
教えていただけますでしょうか？

bwaは以下のようになっています。
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <[e-mail]>

Usage: bwa <command> [options]

Command: index index sequences in the FASTA format
mem BWA-MEM algorithm
fastmap identify super-maximal exact matches
pemerge merge overlapping paired ends (EXPERIMENTAL)
aln gapped/ungapped alignment
samse generate alignment (single ended)
sampe generate alignment (paired ended)
bwasw BWA-SW for long queries

shm manage indices in shared memory
fa2pac convert FASTA to PAC format
pac2bwt generate BWT from PAC
pac2bwtgen alternative algorithm for generating BWT
bwtupdate update .bwt to the new format
bwt2sa generate SA from BWT and Occ

[KojiDoi]

「できているとおもう」では話になりません。
${bwa} mem -R ${rg} ${ref} ${fq1} ${fq2} > out.sam
これでまともなout.samが出来ているかを確認するのが先です。

[tare]

全然理解しておらず、教えていただきありがとうございます。
${bwa} mem -R ${rg} ${ref} ${fq1} ${fq2} > out.sam
をいれたところ、{bwa}: コマンドが見つかりませんとでてしまいました。
out.samのファイルも空っぽなのができてしまいました。
教えていただきたいです。

[tare]

また再度、
${bwa} mem -R ${rg} ${ref} ${fq1} ${fq2} > out.sam
をいれたところ、
[M::bwa_idx_load_from_disk] read 0 ALT contigs
強制終了
になってしまいました。

[KojiDoi]

echo ${bwa} mem -R ${rg} ${ref} ${fq1} ${fq2}
このコマンドで何が表示されますか？

[tare]

$ echo ${bwa} mem -R ${rg} ${ref} ${fq1} ${fq2}
をいれたところ、
mem -R
表示されました。

[KojiDoi]

失礼、スクリプトの中で変数を定義してからbwaを呼び出しているのでしたね。であれば多分問題はそこではないですね。
実はちゃんとbwaを使ったことがないのであやふやな記憶ですが、事前にインデックスファイルを作っておく必要があったのではなかったでしょうか。そのファイルはちゃんとできていますか？

[tare]

../bwa-0.7.17/bwa index hg38.fasta
を実行して、参照ゲノム配列を収めたfastaファイルに対してのインデックスは作成でき、
出来上がったものを
ls -lh hg38.fasta
rw-r--r-- 1 ユーザ名 3.0G Apr 15 07:24 hg38.fasta
ls -lh hg38.fasta.*
-rw-r--r-- 1 ユーザ名 317K Apr 15 19:38 hg38.fasta.amb
-rw-r--r-- 1 ユーザ名 1.1K Apr 15 19:38 hg38.fasta.ann
-rw-r--r-- 1 ユーザ名 2.9G Apr 15 19:37 hg38.fasta.bwt
-rw-r--r-- 1 ユーザ名 738M Apr 15 19:38 hg38.fasta.pac
-rw-r--r-- 1 ユーザ名 1.5G Apr 15 20:19 hg38.fasta.sa

で確認してました。

[KojiDoi]

確認ですが、それらのファイルはサブディレクトリRefHg38の下に存在するのでしょうか？

[tare]

RefHg38の下に存在します。