macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

fastqから素早くインサートサイズを計算する

 

 bamファイルをすでに作っているなら、ペアエンドのインサートサイズはPicard-tools等ですぐ出せますが、raw fastqしかない時にいちいちbamにして求めるのは少し面倒です。ワンランナーで出すスクリプト書きました。好みにあわせて修正して使ってください。手順ですが、fastqをランダムサンプリングし、minimap2&&samtoolsでbamを作り、goleftに渡してインサートサイズを出しています。リファレンスがなければ、de novo assemblyを行い、得られたcontigを使ってください。

 

インストール

依存

  • samtools (1.7で動作確認。古いバージョンでなければO.K)。
  • seqkit(紹介
  • minimap2(紹介
  • goleft(紹介

> 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

f:id:kazumaxneo:20180611215507j:plain

インサートサイズ(outer insert size)は418、standard deviationは73だった。

 

 

感想

goleftはピカード(picard-tools)あたりに変え(紹介)、グラフィック出力にしてもいいんじゃないでしょうか。それから、ラージゲノムを使っている人は、小さな染色体だけ取り出して計算してください (ヒトならchr20とか)。時間がかかります。その時は、seqkitのサンプリングサイズを、10万リードから100万リードくらいに増やしてください。マッピングされたリードがあまりに少ないとエラーを吐きます。

 

お粗末様でした。

 

引用

https://www.biostars.org/p/50783/