Answers to: samtoolsの使い方

Comment by pan on pan's 質問

pan — Wed, 01 Jun 2011 14:42:01 +0900

かうぱーと（http://miyahn.blog80.fc2.com/）のサイトにはBAMからunmapped readを抽出方法が書いてあります。これが出来るなら逆の事も出来るだろうと思い調べ、あるPDF(http://bioinf.eva.mpg.de/ngsa_2010/data/samtools_talk.pdf)に行き着いたのですが、結局pやPの意味がわからなくて質問させて頂きました。

Comment by nob_fj on pan's 質問

nob_fj — Wed, 01 Jun 2011 14:26:33 +0900

"mappingされたものとmappingされてないものを Samtools view -f p -F P で抽出できる"の情報源もあれば、より親切な質問かもしれません。

Comment by nob_fj on nob_fj's 回答

nob_fj — Wed, 01 Jun 2011 12:42:46 +0900

2の質問の参考情報です。私自身は、あんまりSAMtoolsでフィルタリングしていない(スマートじゃなくawkを使っている)ので、精通していませんが、-fと-Fは共通の"FLAG"を使っているようですね。

下のメッセージの"Options:"の"-f"と"-F"の記述、および、最後の節の6番の記述が答えになるでしょうか。

解釈が正しいか自身は無いですが、 "-f p"が"required flag"に"p=0x1 (paired), "を指定、 "-F P"が"filtering flag"に"P=0x2 (properly paired), "を指定という意味になるでしょうか。

この指定が、"mappingされたものとmappingされてないものを"抽出することになるのかは、私は分かりません。実際に、小さいデータで細かい動作を確認すると良いのではないでしょうか。

$ /toolPath/samtools-0.1.12a/samtools view -F
view: option requires an argument -- F

Usage:   samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b       output BAM
         -h       print header for the SAM output
         -H       print header only (no alignments)
         -S       input is SAM
         -u       uncompressed BAM output (force -b)
         -x       output FLAG in HEX (samtools-C specific)
         -X       output FLAG in string (samtools-C specific)
         -c       print only the count of matching records
         -t FILE  list of reference names and lengths (force -S) [null]
         -T FILE  reference sequence file (force -S) [null]
         -o FILE  output file name [stdout]
         -R FILE  list of read groups to be outputted [null]
         -f INT   required flag, 0 for unset [0]
         -F INT   filtering flag, 0 for unset [0]
         -q INT   minimum mapping quality [0]
         -l STR   only output reads in library STR [null]
         -r STR   only output reads in read group STR [null]
         -?       longer help

Notes:

  1. By default, this command assumes the file on the command line is in
     the BAM format and it prints the alignments in SAM. If `-t' is
     applied, the input file is assumed to be in the SAM format. The
     file supplied with `-t' is SPACE/TAB delimited with the first two
     fields of each line consisting of the reference name and the
     corresponding sequence length. The `.fai' file generated by `faidx'
     can be used here. This file may be empty if reads are unaligned.

  2. SAM->BAM conversion: `samtools view -bT ref.fa in.sam.gz'.

  3. BAM->SAM conversion: `samtools view in.bam'.

  4. A region should be presented in one of the following formats:
     `chr1', `chr2:1,000' and `chr3:1000-2,000'. When a region is
     specified, the input alignment file must be an indexed BAM file.

  5. Option `-u' is preferred over `-b' when the output is piped to
     another samtools command.

  6. In a string FLAG, each character represents one bit with
     p=0x1 (paired), P=0x2 (properly paired), u=0x4 (unmapped),
     U=0x8 (mate unmapped), r=0x10 (reverse), R=0x20 (mate reverse)
     1=0x40 (first), 2=0x80 (second), s=0x100 (not primary),
     f=0x200 (failure) and d=0x400 (duplicate). Note that `-x' and
     `-X' are samtools-C specific. Picard and older samtools do not
     support HEX or string flags.

Comment by pan on nob_fj's 回答

pan — Wed, 01 Jun 2011 11:45:48 +0900

いろいろ試して例３であることがわかりました。配列名も配列の長さも問題なかったのですが、根本的に間違っていました。テキスト上に配列名書いて、配列コピペして、拡張子を.fastaにしてと外見上fastaファイルにしてたのですが、中身が違うんですね……今まで問題無かったのが不思議です。ありがとうございました。

Answer by nob_fj

nob_fj — Tue, 31 May 2011 20:58:50 +0900

予想するに、"ref.fasta"に">"で始まる配列名が無いのではないかと思います。

簡単なfastaファイルと、fastaでないファイルを入力として試しましたが、以下のメッセージが出ました。これの3と同じ(言語は違うようですが)メッセージが出ているので、配列名が無いのではないのでしょうか。

あまり分かりやすく無いかもしれませんが、こちらも参考に下さい。NGS Surfer's WikiのFASTAの説明ページ

例1正しいfastaファイル。

$ cat tmp1.fasta
>chr1
aaaa
aaaa
aaa
$ /toolpath/samtools-0.1.12a/samtools faidx tmp1.fasta
$ ls tmp1.fasta*
tmp1.fasta  tmp1.fasta.fai

例2 fastaとしては正しいが、samtoolsのfaidxが許容しない書式。

$ cat tmp2.fasta
>chr1
aaaa
aaaaa
aaa
$ /toolpath/samtools-0.1.12a/samtools faidx tmp2.fasta
[fai_build_core] different line length in sequence 'chr1'.
セグメンテーション違反です

例3 fastaではないファイル。

$ cat tmp3.fasta
aaaa
aaaa
aaa
$ /toolpath/samtools-0.1.12a/samtools faidx tmp3.fasta
[fai_build_core] different line length in sequence '(null)'.
セグメンテーション違反です

Comment by t_a_b_e on nob_fj's 回答

t_a_b_e — Tue, 31 May 2011 20:20:39 +0900

１は、ref.fastaの各行の長さが違ったりしてませんか？（配列の最後の行以外で）