ログイン 概要 よくある質問

先日からsamtoolsを使い始めて、わからない点が2つあったので教えて下さい。

1、リファレンスゲノム用のindexが作成できない samtools faidx ref.fastaとやると [fai_build_core] different line length in sequence '(null)'. Segmentation fault となってしまい作成できない。

2、samtools viewのオプションの使い方がわからない。 mappingされたものとmappingされてないものを Samtools view -f p -F P で抽出できるようですが、pやPの意味がわからない。また他にもよく使うオプションがあったら教えて下さい。

ずっとwetな研究をしてきたので、わかりやすく教えて頂けると幸いです。

質問日 May 31 '11 at 19:45

pan's gravatar image

pan
1157

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

(Jun 01 '11 at 14:26) nob_fj ♦ nob_fj's gravatar image

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

(Jun 01 '11 at 14:42) pan pan's gravatar image

予想するに、"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)'.
セグメンテーション違反です

回答日 May 31 '11 at 20:58

nob_fj's gravatar image

nob_fj ♦
50761328

edited May 31 '11 at 21:02

1は、ref.fastaの各行の長さが違ったりしてませんか?(配列の最後の行以外で)

(May 31 '11 at 20:20) t_a_b_e t_a_b_e's gravatar image

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

(Jun 01 '11 at 11:45) pan pan's gravatar image
1

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.
(Jun 01 '11 at 12:42) nob_fj ♦ nob_fj's gravatar image
あなたの回答
プレビューをトグルする

この質問をフォローする

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

回答

回答とコメント

タグ:

×47
×10
×4

質問日: May 31 '11 at 19:45

閲覧数: 13,883 回

最終更新日: Jun 05 '12 at 23:02

powered by OSQA