macでインフォマティクス

macでインフォマティクス

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

FASTA/FASTQの操作ツール

 

情報が増えてきたので、これまで紹介してきたfasta、fastqの分析、変換(圧縮)、修復ツールをまとめておく。

 

 追記

 

 

FASTA

ヘッダの連番へのリネーム

perl -ane 'if(/\>/){$a++;print ">name_$a\n"}else{print;}' input.fa > rename.fa

-- seqkit --

カレントディレクトリの全fa.gzファイルを分析。

seqkit stats *.fa.gz

reverse complementary (相補鎖にしてさらにリバースにする)

seqkit seq -pr input.fa > output.fa

1000bp以上の配列だけ取り出す。

seqkit seq -m 1000 scaffolds.fa > output.fasta

最初の100-bpを取り出す。

seqkit subseq -r 1:100 contigs.fa > left.fa

 最後の100-bpを取り出す(マイナスで後ろからカウント)。 

seqkit subseq -r -100:-1 contig.fa > right.fa

 gtfで指定した領域の配列だけ取り出す。

seqkit subseq --gtf cDNA.gtf genome.fa > cdna.fa

長さでソートし、10000-bp以上の配列だけ出力。

seqkit sort --quiet -r -l conitg.fa -L 10000 > sorted.fa

100-bp windowでゲノム配列を切り出す(GC含量カウントに使ったりする)。

seqkit sliding -s 1 -W 100 genome.fa > 100bp_windows.fa

GC含量をカウント。

seqkit fx2tab genome.fa| head -n 1 | seqkit tab2fx | seqkit sliding -s 1 -W 100 | seqkit fx2tab -n -g > GC.txt

tab形式に変換し、GCなどの情報を表示。

seqkit fx2tab input.fa > output.fa

環状ゲノムのスタート位置を変える。

seqkit restart -i 100 genome.fa > genome_restart.fa

-- BBtools --

contig.fastaを分析する。

stats.sh in=contigs.fa gc=gc.txt gcformat=4

N50などの基礎的な情報が出力される。

user$ stats.sh in=contigs.fa gc=gc.txt gcformat=4

A C G T N IUPAC Other GC GC_stdev

0.2841 0.2276 0.2531 0.2352 0.0000 0.0000 0.0000 0.4807 0.0583

 

Main genome scaffold total:         42

Main genome contig total:           42

Main genome scaffold sequence total: 0.004 MB

Main genome contig sequence total:  0.004 MB  0.000% gap

Main genome scaffold N/L50:         9/131

Main genome contig N/L50:           9/131

Main genome scaffold N/L90:         42/41

Main genome contig N/L90:           42/41

Max scaffold length:                352

Max contig length:                  352

Number of scaffolds > 50 KB:        0

% main genome in scaffolds > 50 KB: 0.00%

 

 

Minimum Number        Number        Total         Total         Scaffold

Scaffold of            of            Scaffold      Contig        Contig  

Length  Scaffolds     Contigs       Length        Length        Coverage

-------- -------------- -------------- -------------- -------------- --------

    All             42             42         3,520         3,520 100.00%

     50             19             19         2,547         2,547 100.00%

    100             11             11         2,009         2,009 100.00%

    250              1             1           352           352 100.00%

gc含量は別ファイルに出力される。

 

 

 

 

 

 

 

FASTQ。

-- seqkit --

カレントディレクトリの全fqファイルを分析。

seqkit stats *.fq

フルネームが重複したリードを除去。

seqkit rmdup -n input.fq  

順番をバラバラにする。

seqkit shuffle input.fastq > shuffled.fastq

fastqをfastaに変換

seqkit fq2fa input.fastq > output.fa

duplicateしたIDをユニークな名前に修復。

seqkit rename input.fastq > renamed.fastq

 カレントディレクトリのfastqファイルから、配列が一致する配列を抽出。

seqkit common -s *fq > common_seq.fq

 リードを10%ランダムサンプリング。

seqkit sample -p 0.1 input.fastq > output.fastq

 fastqを均等に10ファイルに分けて出力。

seqkit split -p 10 input.fastq

-- seqtk --

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

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

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

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

クオリティを分析。

seqtk fqchk input.fq > output.txt

-- BBtools --

・様々な統計データを出力。

bbmap.sh in1=reads1.fq in2=reads2.fq bhist=bhist.txt qhist=qhist.txt aqhist=aqhist.txt lhist=lhist.txt ihist=ihist.txt ehist=ehist.txt qahist=qahist.txt indelhist=indelhist.txt mhist=mhist.txt gchist=gchist.txt idhist=idhist.txt scafstats=scafstats.txt

phread quality scoreを一般的なASCII+33からASCII +64に変換する。

reformat.sh in=reads.fq out=reads.fq qin=33 qout=64

fastq => fasta+qualへの変換。

reformat.sh in=reads.fq out=reads.fa qfout=reads.qual

インターレースのペアリードデータを2つのファイルに分ける。

reformat.sh in=interlace.fq out1=read1.fq out2=read2.fq

ペアリードデータをインターレースファイルにする。

reformat.sh in1=read1.fq in2=read2.fq out=reads.fq

ペアリードファイルがずれていないか検証。

reformat.sh in=reads.fq vint

ペアリードファイルの修復。

repair.sh -Xmx20g in1=broken1.fq in2=broken2 out1=fixed1.fq out2=fixed2.fq outs=singletons.fq repair

ペアリードのマージ。

bbmerge.sh in1=reads1.fq in2=reads2.fq out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq ihist=ihist.txt

オーバーラップするペアリードのエラー修復。

bbmerge.sh in=interlace.fq out=corrected.fq ecco mix

1や2がついていないペアリードファイルにペアを識別する1 or 2を入れる。

reformat.sh in=reads.fq out=renamed.fq addslash int

70文字ごとに改行する。(foldコマンドと同じ)

reformat.sh in=reads.fa out=wrapped.fa fastawrap=70

小文字の塩基を大文字にする。

reformat.sh in=reads.fq out=out.fq tuc

名前が全てユニークか確認する。

reformat.sh in=reads.fq out=out.fq uniquenames

クオリテイスコアを特定の範囲に限定する(処理でクオリティがおかしくなってしまった時に使う)

reformat.sh in=reads.fq out=out.fq mincalledquality=2 maxcalledquality=41

k-merのgraphを利用し、マージされなかったfastqを延長しマージを試みる。

bbmerge-auto.sh in=reads.fq out=merged.fq outu=unmerged.fq ihist=ihist.txt ecct extend2=20 iterations=5

リード長よりインサートサイズが短いペアリードからアダプター配列を見つける。

bbmerge.sh in1=reads1.fq in2=reads2.fq outa=adapters.fa

続けて見つかったアダプター(e.g., GATCGGAAGAGCACACGTCTGAACTCCAGTC )を除く。

bbmerge.sh in1=reads1.fq in2=read2.fq out=merged.fq adapter1=GATCGGAAGAGCACACGTCTGAACTCCAGTC 

k-merカバレッジ100でノーマライズ

bbnorm.sh in1=reads1.fq in2=reads2.fq out=normalized.fq target=100 min=5

duplication sequencesの除去。

dedupe.sh in=X.fa out=Z.fa outd=duplicates.fa

99%以上同一の配列を除去。

dedupe.sh in=X.fa out=Y.fa minidentity=99

200bp以上同一の配列を1塩基のミスマッチまで許容して除去。

dedupe.sh in=X.fq pattern=cluster%.fq ac=f am=f s=1 mo=200 c pc csf=stats.txt outbest=best.fq fo c mcs=3 cc dot=graph.dot

 PacBioの16S全長シーケンスデータのクラスタリング1。全長読めてない配列を除去。

reformat.sh in=reads_of_insert.fastq out=filtered.fq minlen=1420 maxlen=1640 maq=20 qin=33 

 1420-bp以下のリードと1640-bp以上のリード(アーティファクト)を除去する。また平均クオリティが20以下のリードも除去する。 

  PacBioの16S全長シーケンスデータのクラスタリング2。99%以上同一の配列をクラスタリング

dedupe.sh in=filtered.fq csf=stats_e26.txt outbest=best_e26.fq qin=33 usejni=t am=f ac=f fo c rnc=f mcs=3 k=27 mo=1420 ow cc pto nam=4 e=26 pattern=cluster_%.fq dot=graph.dot

 3'側のアダプタートリミング。

bbduk.sh in=reads.fq out=clean.fq ref=adapters.fa ktrim=r k=23 mink=11 hdist=1 tpe tbo

平均クオリティが10以下のリードをフィルタリング(除去)。

bbduk.sh in=reads.fq out=clean.fq maq=10

k-merフィルタリング。

bbduk.sh in=reads.fq out=unmatched.fq outm=matched.fq ref=phix.fa k=31 hdist=1 stats=stats.txt

指定ファイルに対して31-merで探索を行い、1塩基のミスマッチを許容して(hdist=1)マッチするリードがあれば、除去される。ここではイルミナのスパイク(バクテリオファージ、5.5-kb)に対して検索している。

k-merマスキング1。

bbduk.sh in=ecoli.fa out=ecoli_masked.fa ref=salmonella.fa k=25 ktrim=N

Quality recalibration

bbduk.sh in=reads.fq out=recalibrated.fq recalibrate sam=mapped.sam

samファイルを元にベースクオリティスコアを再校正したfastqを出力する。

Matching degenerate sequences such as primers

bbduk.sh in=reads.fq out=matching.fq literal=ACGTTNNNNNGTC copyundefined k=13 mm=f 

degenerateな配列を探索する。例えばメタ16Sの可変性領域を網羅的に探す

ユニークなリード数を示したヒストグラムをかく。

bbcountunique.sh in=reads out=histogram.txt

追記

--fqtools--

basetab 塩基をカウント

fqtools basetab input.fq 

qualtab クオリティスコアをカウント

fqtools qualtab input.fq 

count リード数をカウント

fqtools count input.fq

sequence 配列だけ表示(配列のみのcat)

fqtools sequence input.fq > output

header ヘッダーだけ表示(ヘッダーのみのcat)

fqtools header input.fq > output

quality クオリティスコアだけ表示(クオリティスコアのみのcat)

fqtools quality input.fq > output

validate 配列ATGCATGCを探す(grep、agrepに近い)。

fqtools find -s ATGCATGC input.fq

trim トリミング

#先頭10bp
fqtools trim -s 10 input.fq > trimmmed.fq

#100bpになるまでカット
fqtools trim -l 100 input.fq > trimmmed.fq

-- ngs crumbs--

ペアリードからインターレースのfastqを作る。

interleave_pairs R1.fq R2.fq > merged.fq

インターレースからペアリードに分離する。

interleave_pairs merged.fq -o R1.fq R2.fq

ペアじゃないリードを除く。

pair_matcher interelace.fq -o output.fq -p orphan.fq

 fastqのフォーマットを変更する(fasta、fastq、fastq-illumina)。

 convert_format input.fastq -f fasta -o output.fa

 --fastqutils--

 properpairs(リンクペアエンドでペアになっていないリードを除く。

fastqutils properpairs pair1.fq pair2.fq output_pair1.fq output_pair2.fq

 revcomp(リンク逆相補鎖(reverse complementary)に変換。

fastqutils revcomp input.fq > output.fq

 sort(リンクnameやポジションでソートする。

fastqutils sort input.fq > sort.fq

 split(リンク指定数に分割する。

#3つに分割する。
fastqutils split input.fq output 3

 convertqual(リンクイルミナのクオリティ値をsanger scaleに変換。

fastqutils convertqual input.fq > output.fq

 csencode(リンクcolor spaceをATGCに変換。

fastqutils csencode input.fq > output.fq
0 -> A
1 -> C
2 -> G
3 -> T
4,5,6,. -> N

 fromfasta(リンクqualとfastaからfastqを作成。

fastqutils fromfasta input.csfasta input.qual  > output.fq

 fromqseq(リンクillumina export.txtからfastqを作成。

fastqutils fromqseq export.txt

 stats(リンクstatistics

fastqutils stats input.fq > summary

 出力。全部位の情報も出る。

$ fastqutils stats forward.fq

Done! (0:58)                                                                                  

Space: basespace

Pairing: Fragmented

Quality scale: Illumina

Number of reads: 224924

 

Length distribution

Mean: 230.287452651

StdDev: 49.7016332922

Min: 20

25 percentile: 204

Median: 235

75 percentile: 262

Max: 301

Total: 224924

 

Quality distribution

pos mean stdev min 25pct 50pct 75pct max count

1 33.9262595366 0.436774185711 30 34 34 34 38 224924

2 33.8028489623 1.54960128452 10 34 34 34 38 224924

3 33.8104426384 1.54959912975 10 34 34 34 38 224924

4 33.8467082214 1.44051755323 10 34 34 34 38 224924

5 33.8715655066 1.56917725475 10 34 34 34 38 224924

6 37.6461204674 1.84245487446 10 38 38 38 38 224924

7 37.622272412 1.94525249168 10 38 38 38 38 224924

8 37.6478232647 1.86618349415 10 38 38 38 38 224924

9 37.6512333055 1.85726711398 10 38 38 38 38 224924

10 37.6626060358 1.79494590177 10 38 38 38 38 224924

11 37.6324714126 1.89571282258 9 38 38 38 38 224924

12 37.6014342622 1.99703097792 10 38 38 38 38 224924

13 37.6435818321 1.8427123548 10 38 38 38 38 224924

14 37.6336940478 1.88442047255 10 38 38 38 38 224924

15 37.6378332237 1.85404372502 9 38 38 38 38 224924

16 37.6625215628 1.78640763134 9 38 38 38 38 224924

17 37.6524514947 1.80913909303 10 38 38 38 38 224924

18 37.6461649268 1.81338343748 9 38 38 38 38 224924

 

追記

fastaからフェイクfastqを作成。

UCSCのkent utilitiesのfaToFastqを使う。

git clone git://github.com/ENCODE-DCC/kentUtils.git
cd kentUtils
make
rsyncでbinのパスを通す
sudo rsync -a -P ./bin/ /usr/local/bin/kentUtils/

 faToFastqだけならcondaでも入る(リンク)。

クオリティ"<"のfastqに見せかける。

faToFastq -qual=\< input.fa output.fq

 >head output.fq

$ head output.fq 

@DEAMERNANOPORE_20161206_FNFAB49164_MN16450_sequencing_run_MA_821_R9_4_NA12878_12_06_16_71094_ch374_read7088_strand.fast5

ATTCGTTTGGCTGGGTTCAGGATGCGTGTTGCTCGCTCCCGTTCCATTCATCACTGCTAAATTCTGTCCTCCCTTCATCACTCTTTTCGTTTCCATTCCCGTTTCCACTCCGATTCCCGTTCCTTCCTGCGTTCCATCACTCCGCTCCGTTGCTCCATTGCTCCGTTCCGTTAGGCTCCCGTTTCGTTCCACTCATCCGTTTCTGATCCATCACTCTGCTCCATTGCCGAGTCCACTCCCATTCAGATTCAGTTCTGCCCTCCGTTCCATTTCCGCTCATCATTCCAATCCGCCGCATGGGTTCCATTCCACTCAGGTTCCGTTCCGTTGCTCCGTTCCATTGCCTTCCGTTCGATCGGTTAGTTCCTCCCGTTTCCAGGTTTCCATCCACTCCGTTTCTGTTTCATCTGCTCCTTTCCTCCGTTCCACTTCATTCATTCCTGTCTCTCATCACTACATCCCAGTCACTCCTTTCAGTTCCTCCACTCCACTCCCGATTTCCTCTCCAGTCGGCTCCATCCGTTCCATTCCACTCATCCG

+

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

 

 

 

 

各コマンドの詳細はTopのツールのリンクで確認してください。

 

引用

BBMap (aligner for DNA/RNAseq) is now open-source and available for download.

 

SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation

Wei Shen, Shuai Le, Yan Li , Fuquan Hu 

PLOS ONE. 2016. doi:10.1371/journal.pone.0163962.