2019年12月11日水曜日

fastaファイルを取り扱う

大量にデータを含むfastaファイルを取り扱う必要があった。長さに応じて選別したり配列を取り出したりしたいのだが、そういえば何使えばよかったかなと思った。自分で最初から書いてもそこまで難しくはないだろうが、明らかに車輪の再発明っぽいしバグも紛れてしまうだろう。binの中をのぞいて、ucscのプログラムがたくさん入っている中から、faSizeなど使えそうなものを探した。

http://hgdownload.soe.ucsc.edu/admin/exe/macOSX.x86_64/README.txt

 また、fastaとfastqを扱う多機能ツールであるseqkitも以前インストールしていた。

https://bioinf.shenwei.me/seqkit/usage/

 どちらも多機能で、とりあえず使いたい機能を使うためにどのオプションを使おうかというところで悩む感じだが、思いついた時にパイプラインとして記録しておいたり、何がしかwrapper scriptを作っておいたらいいかもしれない。

色々やり方があるだろうし、もっと良いやり方もあると思うが、とりあえず思いついたものを記録していきたい。

(1)fastaやfastqを長い方から表示する

簡単な方法として、seqkitを使う。
 seqkit sort --by-length --reverse sequence.fa

ファイルに保存する場合は、
seqkit sort --by-length --reverse sequence.fa > result.fa

上だと全部が表示されるが、一定以上の長さのものだけ表示したい場合は、
seqkit sort --by-length --reverse sequence.fa | seqkit seq --min-len 500
seqkit sort --by-length --reverse sequence.fa | seqkit seq --min-len 500 > result.fa

これで、ある一定以上の長さの配列について、長いものから見て検討することができる。

(2)配列の名前と長さだけ表示する

faSize -detailed sequence.fa

この場合、自動的にseq IDだけになるようだ。

同じことをseqkitでやりたい場合
seqkit fx2tab --length --name --only-id sequence.fa

長さでソートされ、一定以上の長さの配列だけにして、配列長を表示したい場合
seqkit sort --by-length --reverse sequence.fa | seqkit seq --min-len 500 | seqkit fx2tab --length --name --only-id


思いついた都度、続く

(3) fastqファイルをfastaファイルに変換する

seqkitを使う場合、
seqkit fq2fa input.fastq

スレッドの数を増やす場合、
seqkit fq2fa --threads 8 input.fastq

ファイルを出力する場合
seqkit fq2fa input.fastq --out-file output.fasta

16スレッドを使い、fastqファイルをfastaファイルとして長さでソート、50Kb以上の配列を出力する場合、
seqkit fq2fa --threads 16 input.fastq | seqkit seq --min-len 50000 | seqkit sort --by-length --reverse
seqkit fq2fa --threads 16 input.fastq | seqkit seq --min-len 50000 | seqkit sort --by-length --reverse --out-file output.fa

 

(4) seqkitのグローバルフラッグ

Global Flags:
      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)
      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...
      --id-regexp string                regular expression for parsing ID (default "^(\\S+)\\s?")
  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)
  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")
      --quiet                           be quiet and do not show extra information
  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")
  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

 
思いついた都度、続く