處理Fastq和Fasta檔案是進入次世代定序後常會遇到的情境，且當你有特殊需求時，可能就必須要從底層開始來設計你的序列分析方法。

之前寫過 Bioawk：專門處理定序相關格式的awk（fasta,fastq,bed,sam,vcf,gff） 和 awk 進階筆記：字串處理, awk裡好用的變數：FS, OFS, RS, ORS, NR, NF, FILENAME, FNR。 最近開始又要處理序列，所以來稍微複習一下。

推薦觀看美國猶他大學(University of Utah)的教授Aaron R. Quinlan的應用基因體學課程中分享AWK和Bioawk的影片：Applied Computational Genomics – awk and bioawk，裡面的講解簡單且易上手。

AWK是一個古老的Unix程式，在1970年代在當時資通訊最前沿的貝爾實驗室，由三位工程師Alfred Aho, Peter Weinberger, Brian Kernighan共同開發的，主要是一款字串處理程式。

AWK reads the input a line at a time. A line is scanned for each pattern in the program, and for each pattern that matches, the associated action is executed.

— Alfred V. Aho

使用bioawk和seqkit基本上，就可以解決百分之90以上，在fasta/fastq的前處理和基本分析需求，另外，fastp也是另一個蠻不錯的軟體，蠻適合用來修剪(trim)序列使用，速度也快。

bioawk基本上是建立在awk的語法上，但提供不錯原生的變數可使用，基本涵蓋次世代定序相關的資料結構，所以非常方便，減少許多重新發明輪子的問題，如下表：

BED	SAM	VCF	GFF	FASTX(FASTA/FASTQ)
chrom	qname	chrom	seqname	name
start	flaq	pos	source	seq
end	rname	id	feature	qual
name	pos	ref	start	comment
score	mapq	alt	end
strand	cigar	qual	score
thickstart	rnext	filter	filter
thickend	pnext	info	strand
rgb	tlen		group
blockcount	seq		attribute
blocksizes	qual
blockstarts

使用情境

計算序列檔案中有幾個reads？

bioawk -c fastx '{print $name}' input.fasta | wc -l

計算序列檔案中，每一個reads有多長，並且輸出成一個tab相隔的資料

bioawk -c fastx '{prin $name, length($seq)}' input.fasta > reads_length.txt

篩選序列檔案中，長度大於3000，小於4000的reads

bioawk -c fastx 'length($seq)>3000 && length($seq)<4000{print ">"$name;print $seq}' input.fasta > filtered_input.fasta

進階情境

篩選序列檔案中，含有在序列ID檔案中的reads

seqkit grep -f curated_id.txt input.fastq -o curated_input.fastq

篩選序列檔案中，用特定引子對去做對照，其可以夾出來的核酸區域

seqkit amplicon -F ATACTTCCTATCGCGTA -R TCTATCGCTCGATC input.fasta

閱讀參考：

生信小工具：awk的升级版bioawk

bioawk github

fastp gitbhub

seqkit 使用说明

2016. SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. PLos One

Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, https://doi.org/10.1093/bioinformatics/bty560

Bioawk basic