macでインフォマティクス

macでインフォマティクス

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

FASTA/FASTQ関係のツールまとめ

summary fasta/fastqの操作

随時更新

 

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

 

 

 

samtools/bcftools/htslibの導入はこちらを参考にしてください。綺麗にまとまってます。

Install samtools, bcftools and htslib on linux · GitHub

f:id:kazumaxneo:20181005144558j:plain

 

#Anaconda環境ならバージョン指定するだけで導入できます。（バージョンチェック）

conda install -c bioconda samtools==1.9
conda install -c bioconda bcftools==1.9

#htslibだけ入れたいなら
conda install -c bioconda htslib==1.9

 

高速化のTIPS

たくさんのデータに同じ処理をするなら、GNU parallelを使って並列化することで簡単に高速化できる。

例えばseqkit statsを全てのfastqに対して実行することを考える。ワイルドカードでも良いが、データが多いとかなり時間がかかってしまう。

seqkit stats -a -T *fq.gz > stats.txt

GNU parallelで並列化する（"sudo apt-get install parallel"）。

#listを作成
ls *fq.gz > list

#10並列（tailで1行目を捨ててエラー出力を追記）
cat list | parallel -j 10 'seqkit stats -a -T {} |tail -n 1 1>> stats.txt 2>&1'

＊HDDだとI/Oが足を引っ張る可能性があります。並列処理数を抑えて下さい。

(parallel --citation => will cite、で引用するのを忘れない事)

 

FASTA。

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

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

#or
awk '/>/ {$0 = ">rename_" ++n} 1' input.fa > rename.fa

またはseqkit renameを使う

#rename　duplicateしたIDをユニークな名前に修正。duplicateしていると判断されると"_2"などがつく。
seqkit rename -n input.fa > renamed.fa
  •  -n　check duplicated by full name instead of just id.

fastq => fasta

awk '(NR - 1) % 4 < 2' test.fq | sed 's/@/>/' > test.fa

#seqkit 
seqkit fq2fa in.fq > out.fa

#8並列（他のコマンドでももちろん利用可能、CPUとメモリがたくさん利用できる時だけ）。linuxのGNU paralel使用。
ls *.fq | parallel -j 8 'seqkit fq2fa {} > {.}.fa'

 fastq => fasta変換して、さらに連番へのリネームも同時実行

awk 'NR % 4 == 2 {print ">" ++n "\n" $0}' input.fq > output.fa

multi-lineのfastaからシングルラインのfastaに 修正(*1より）。

perl -pe '/^>/ ? print "\n" : chomp' input.fasta > output.fasta

指定領域を取り出す。

samtools faidx ref.fa chr2:1000-25000 > output.fa

(multi-)fastaの配列の改行を除き配列は１行で出力。（Biostarスレッドより）。

awk '/^>/ { if(NR>1) print ""; printf("%s\n",$0); next; } { printf("%s",$0);} END {printf("\n");}' input.fa > output.fa

multi-fastaの配列を合体させて強制的に１つに配列にする（Biostarスレッドより）。

cat multi_seq.fasta | sed -e '1!{/^>.*/d;}' | sed ':a;N;$!ba;s/\n//2g' > one_seq.fasta

#またはgrep ーVで>以外の行を抽出、seqretで修復
grep -V ">" multi_seq.fasta > out
seqret out one_seq.fasta

N（NNN...）を含む配列だけ除く。

#Nを1塩基以上含む配列は除く
cutadapt --quiet --max-n 1 input.fa > out.fa

リード長分布

samtools faidx input.fasta
cut -f2 input.fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

f:id:kazumaxneo:20190913000251p:plain

フォーマット変換

NGSなどのリードは対応していないが、小さめのFASTAデータセットであれば、オンラインで様々なフォーマットに変換できる。何千配列もあるような大きなデータは投げないよう注意する。

１、Sequence conversion

http://sequenceconversion.bugaco.com/converter/biology/sequences/

f:id:kazumaxneo:20191124145452p:plain

FASTAからNEXUS、Philip、clustal、EMBLなど様々なフォーマットに変換可能。

 

２、Sequence format converter

Sequence format converter

f:id:kazumaxneo:20191124145909p:plain

こちらも様々な出力に対応。

 

SMS

https://www.bioinformatics.org/sms2/

f:id:kazumaxneo:20191124150010p:plain

こちらも昔から保守されている。

 

-- 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

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

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

指定した配列を取り出す。opera_contig_216907~というヘッダの配列とopera_contig_162971~というヘッダの２配列のみ取り出す。

cat input.fasta |seqkit grep -p opera_contig_216907 -p opera_contig_162971

gzip圧縮ファイルならcat => zcatに変更。 

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

seqkit sort --quiet -r -l conitg.fa -L 10000 > sorted.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含量は別ファイルに出力される。

 

たくさんファイルがあるなら、並列処理することを考える。たとえばGNU parallelを使って、カレントのfastaファイルを対象にstatsコマンドを走らせる（別々に保存したい）。

ls *fasta | parallel --jobs 8 'seqkit stats {} >{}_stas.txt' 2>error_log 

 出力

f:id:kazumaxneo:20200330101750p:plain

 

 

FASTQ。

-- parallel-fastq-dump --

SRAのfastqを高速ダウンロード、例えばペアエンドデータSRRxxxxxxをダウンロードしてout_dir/にペアエンドfastq.gz出力。

parallel-fastq-dump --threads 8 --split-files --gzip --outdir out_dir --sra-id SRRxxxxxx


 

-- 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 --

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

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

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

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や２がついていないペアリードファイルにペアを識別する1 or 2を入れる。

reformat.sh in1=pair1.fq in2=pair2.fq out1=renamed1.fq out2=renamed2.fq addslash=t int=f

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以上同一の配列を１塩基のミスマッチまで許容して除去。

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全長シーケンスデータのクラスタリング１。全長読めてない配列を除去。

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全長シーケンスデータのクラスタリング２。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マスキング１。

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

 

指定リファレンスファイルに対して31-merで探索を行い、1塩基のミスマッチを許容して（hdist=1）マッチするリードがあれば、分離する。イルミナのスパイク（バクテリオファージ、5.5-kb）に対して検索をかけている。マッチしたペアリードとマッチしなかったペアリードをそれぞれペアの順番を同期したまま出力する（=４つのfastqが出力される）。

bbduk.sh in1=pair_1.fq in2=pair_2.fq \
 out1=unmatched_1.fq out2=unmatched_2.fq \
 outm1=matched_1.fq  outm2=matched_2.fq \
 ref=phix.fa k=31 hdist=1 stats=stats.txt
  • outm=<file> (outmatch) Write reads here that contain kmers matching the database.
  • ref=<file,file> Comma-delimited list of reference files.
  • hammingdistance=0 (hdist) Maximum Hamming distance for ref kmers (subs only). Memory use is proportional to (3*K)^hdist.

連結したリファレンス配列を指定すれば、バテクリアの汚染除去などにも使用できる。ただし扱いには注意する。

 

-- fixfastq -- (github)

#build
git clone https://github.com/davidvi/cleanFastq.git
cd cleanFastq/
g++ -std=c++11 fixFastq.cpp -o fixFastq
sudo cp fixFastq /usr/local/bin/

#run
fixFastq broken.fastq > fixed.fastq

--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

 

sff2fastq

Roche 454の.sff => fastq変換

git clone git://github.com/indraniel/sff2fastq.git;
cd sff2fastq
make
ln -s sff2fastq /usr/local/bin/

#実行
cat input.sff | sff2fastq > file.fastq

 

 

Ming Tangさんの素晴らしいワンライナー集（bioinformatics-one-liners）

 

追記 

--seqmajic--（リンク）

fastaフォーマットをphyフォーマットに変換（アライメント実行後のfastaを使う）

seqmagick convert input.fa output.phy

stockholm format（wiki）に変換（アライメント実行後のfastaを使う）

seqmagick convert input.fa output.sto

multi-fastaから先頭5配列を出力

seqmagick convert --head 5 input.fa firs_five_seq.fa

multi-fastaからラスト5配列を出力。ただし、５配列のうち、100bp以下は出力されない（--head 5 => --min-length 100 ）。

seqmagick convert --tail 5 --min-length 100 input.fa last_five_seq.fa

multi-fastaの100bp以上の配列を対象に、最後から5配列を出力（--min-length 100 => --head 5 ）

seqmagick convert --min-length 100 --tail 5 input.fa last_five_seq.fa

multi-fastaの100bp以上の配列を対象に、それぞれの配列の1bpから5bpまでのポジションを抽出

seqmagick convert --min-length 100 --cut 1:5 input.fa pos1-5.fa

 DNA <=> RNA変換

#DNA => RNA
seqmagick convert --transcribe dna2rna input_dna.fa output_rna.fa 

#RNA => DNA
seqmagick convert --transcribe rna2dna input_rna.fa output_dna.fa
  • --transcribe {dna2rna,rna2dna}
  • --translate {dna2protein,rna2protein,dna2proteinstop,rna2proteinstop} 

 duplicationを除く（配列をチェックしている。ヘッダーが違っても配列が100%同一だと除かれる）

seqmagick convert --deduplicate-sequences input.fa output.fa

 

追記

fastaからpseudo-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

+

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

 

 --BBSketch--（リンク）データに含まれる生物種を超高速推定（fastqも使えます）

sendsketch.sh in=contigs.fa

sketchファイルで結果は問い合わされ、 数秒で結果が返ってくる。

f:id:kazumaxneo:20191218180519p:plain

  • KID (or kmer identity) is simply the percent of the kmers that matched
  • WKID is the weighted kmer identity, which has been adjusted to compensate for genome size.
  • matches is the number of matching kmers between the query and reference sketch

 

 

 --fastp--（リンク）

fastq出力を指定しなければレポートファイル（htmlとjson）のみ出力される。

fastp -i pair1.fq -I pair2.fq
  • -I     read2 input file name (string [=])
  • -O   read2 output file name (string [=])
  • -3    enable per read cutting by quality in tail (3'), default is disabled (WARNING: this will interfere deduplication for SE data)
  • -w    worker thread number, default is 3 (int [=3])

 出力の一部。

f:id:kazumaxneo:20180521111137j:plain

 

 --fastqp--（リンク）

fastqの分析。inputファイルだけ指定すればランできる。

fastqp input.fastq.gz

出力ディレクトリはgz圧縮される。解凍して中を開く。複数のmetricsで分析が行われている。

  出力の一部。

f:id:kazumaxneo:20180622090028p:plain

 

 

各コマンドの詳細は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.

 

GNU Parallel 20150322

http://dx.doi.org/10.5281/zenodo.16303

 

*1

Question: Multiline Fasta To Single Line Fasta

https://www.biostars.org/p/9262/