macでインフォマティクス

macでインフォマティクス

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

fastq / fastaの操作ツール seqtk

 

seqtkはfastqをfastaに変換したり、相補鎖に変換できるツール。ランダムサンプリング機能ももち、de novo transcriptome解析でアセンブルに有利なリードデプスに間引くツールとして用いられることもある(ペーパー)。動作が非常に高速のため使いやすい。似たツールにseqkitがある。導入して機能をチェックしてみる。

 

 インストール

mac os10.12でテストした。

brewで導入できる。

brew install seqtk

#Anaconda環境なら
conda install -c bioconda seqtk

 Github


 

 ラン

seq  - common transformation of FASTA/Q

fastqをfastaに変換。(.gzも扱える。)

seqtk seq -a input.fq > output.fa

 ILLUMINA 1.3+フォーマットをfastaに変換し、同時にquality20以下を小文字にする。

seqtk seq -aQ64 -q20 input.fq > output.fa

相補鎖(reverse complement)に変換。

seqtk seq -r input.fq > output.fq

様々なオプションがある。

  • -q mask bases with quality lower than INT [0]
  • -X mask bases with quality higher than INT [255]
  • -Q quality shift: ASCII-INT gives base quality [33]
  • -M  mask regions in BED or name list FILE [null]
  • -L drop sequences with length shorter than INT [0]
  • -c  mask complement region (effective with -M)
  • -r reverse complement -A force FASTA output (discard quality)
  • -C drop comments at the header lines
  • -N drop sequences containing ambiguous bases
  • -1 output the 2n-1 reads only
  • -2 output the 2n reads only
  • -V shift quality by '(-Q) - 33'
  • -U convert all bases to uppercases
  • -S strip of white spaces in sequences

 

 mergepe  - interleave two PE FASTA/Q files

ペアリードを1つのファイルにまとめインターレースにする。

seqtk mergepe R1.fq R2.fq > paired.fq

 インターレースファイルからペアじゃないリードを除去するならdropseを使う。

 

trimfq - trim FASTQ using the Phred algorithm

 低クオリティな塩基の除去。Phread quality scoreに基づく。

eqtk trimfq input.fa > output.fa
  • -q error rate threshold (disabled by -b/-e) [0.05]
  • -l maximally trim down to INT bp (disabled by -b/-e) [30]

デフォルトだと5%のエラー率を超えるとトリムされる(-q 0.05)。

  指定範囲の塩基のトリミング。

eqtk trimfq -b 5 -e 10 input.fa > output.fa
  • -b  trim INT bp from left (non-zero to disable -q/-l) [0]
  • -e trim INT bp from right (non-zero to disable -q/-l) [0]
  • -L retain at most INT bp from the 5'-end (non-zero to disable -q/-l) [0] 

 

subseq    - extract subsequences from FASTA/Q

fastqをランダムサンプリングする。-sで乱数発生を制御している。

seqtk sample -s100 read1.fq 10000 > sub1.fq 
seqtk sample -s100 read2.fq 10000 > sub2.fq
  • -s RNG seed [11]
  • -2 2-pass mode: twice as slow but with much reduced memory 

上記では10000リードずつ取り出される。ただし-sの値を変えるとペアじゃないリードが10000リード取り出されるので注意。ペアリードは同じパラメータで行う。

 

fqchk    - fastq QC (base/quality summary)

クオリティの分析。

seqtk fqchk input.fq > output.txt

$ head -5 output.txt

min_len: 101; max_len: 101; avg_len: 101.00; 38 distinct quality values

POS     #bases  %A      %C      %G      %T      %N      avgQ    errQ    %low    %high

ALL     252500  26.2    24.3    23.7    25.8    0.0     35.9    20.4    2.8     97.2

1       2500    13.0    44.2    25.4    17.4    0.0     32.6    30.7    0.4     99.6

2       2500    21.0    22.1    21.7    35.2    0.0     32.9    31.7    0.4     99.6

3       2500    25.4    25.2    22.1    27.2    0.0     33.0    31.8    0.2     99.8

 Rに持ち込んでグラフ化すれば見やすい。

 他にもいくつかの機能がある。

 

引用

GitHub - lh3/seqtk: Toolkit for processing sequences in FASTA/Q formats

 

Using seqtk to trim and process reads at an insanely high speed — ANGUS 2.0 documentation

 

FASTA/Q sequence processing toolkit -- seqtk