bamファイルをすでに作っているなら、ペアエンドのインサートサイズはPicard-tools等ですぐ出せますが、raw fastqしかない時にいちいちbamにして求めるのは少し面倒です。ワンランナーで出すスクリプト書きました。好みにあわせて修正して使ってください。手順ですが、fastqをランダムサンプリングし、minimap2&&samtoolsでbamを作り、goleftに渡してインサートサイズを出しています。リファレンスがなければ、de novo assemblyを行い、得られたcontigを使ってください。
インストール
依存
> cat insert_size.sh
#!/bin/bash
#samtools,seqkit、minimap2、goleftを使う。パスが通っていること。
#version1 2018/06/11
CMDNAME=`basename $0`
if [ $# -ne 3 ]; then
echo "Usage: $CMDNAME ref.fa pair1.fq pair2.fq" 1>&2
exit 1
fi
#10万リードサンプリング (goleftはリードが少なすぎるとエラーになる)
seqkit sample -n 100000 $2 > sampling1.fq 2>/dev/null
seqkit sample -n 100000 $3 > sampling2.fq 2>/dev/null
#mapping
minimap2 -t 8 -ax sr $1 sampling1.fq sampling2.fq 2>/dev/null | samtools sort -@ 8 -O BAM -o sampling.bam - 2>/dev/null
#index
samtools index sampling.bam
#出力
echo
echo ********results********
goleft covstats sampling.bam |column -t
echo
#samplingで出たfastqとbamを消す。
rm sampling1.fq sampling2.fq sampling.bam*
exit
上記をコピーしてファイル名insert_size.shで保存。画面出力は2>/dev/nullで捨てている。進捗が欲しければ消してください。
実行権をつける。
chmod u+x cat insert_size.sh
> insert_size.sh
$ insert_size.sh
Usage: insert_size.sh ref.fa pair1.fq pair2.fq
引数は3つ必要( reference.fasta、pair1.fastq、pair2.fastq)。
パスの通ったディクトリに移動しておく。
ラン
ペアエンドシーケンスデータのインサートサイズを計算。引数はフルパスで与えてもOKです。
insert_size.sh ref.fa pair1.fq pair2.fq
インサートサイズ(outer insert size)は418、standard deviationは73だった。
感想
goleftはピカード(picard-tools)あたりに変え(紹介)、グラフィック出力にしてもいいんじゃないでしょうか。それから、ラージゲノムを使っている人は、小さな染色体だけ取り出して計算してください (ヒトならchr20とか)。時間がかかります。その時は、seqkitのサンプリングサイズを、10万リードから100万リードくらいに増やしてください。マッピングされたリードがあまりに少ないとエラーを吐きます。
お粗末様でした。
引用