在處理大量fasta/fastq時,有時候會需要找一些速度快的程式,其中會發現在python, perl中有許多設計專門用來prase fasta/fastq相關的工具,但速度上總是不甚理想,此時,要是有一個跟awk一樣底層處理字串的工具,那絕對會很強大,這個bioawk便是為此而生的,他是Heng Li(寫samtools的大神)開發的一個小工具,速度會是現行其他工具的十倍速,所以很適合用來整合成自己pipeline的前線。
到哪下載和安裝呢?
可以直接到Heng Li的github下載,然後直接用make指令來編譯,比較舊版的linux系統,可能會需要安裝一些相關用來compile的工具。
那他有什麼基本指令呢?
其實使用起來就是awk,只是他在語法上把本來NR等常見的variable,制換成這些format(fasta/fastq/bed…)中的值,如下面的表格:
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 |
這邊有一些使用範例
Fasta
#fasta 算序列長度 bioawk -c fastx '{ print $name, length($seq) }' input.fasta 算每個序列的gc百分比 bioawk -c fastx '{ print $name, gc($seq) }' input.fasta 轉換成互補序列 bioawk -c fastx '{ print ">"$name;print revcomp($seq) }' input.fasta 只留下序列長度大於100的 bioawk -c fastx 'length($seq) > 100{ print ">"$name; print $seq }' input.fasta 轉換fasta成tabular的格式 bioawk -t -c fastx '{ print $name, $seq }' input.fasta 選取特定id的序列 bioawk -cfastx 'BEGIN{while((getline k <"IDs.txt")>0)i[k]=1}{if(i[$name])print ">"$name"\n"$seq}' input.fasta
Fastq
計算有幾個reads bioawk -t -c fastx 'END {print NR}' input.fastq 轉換fastq成為fasta bioawk -c fastx '{print ">"$name; print $seq}' input.fastq 計算phred平均值 bioawk -c fastx '{print ">"$name; print meanqual($qual)}' input.fastq 選取小於10的reads片段 bioawk -cfastx 'length($seq) > 10 {print "@"$name"\n"$seq"\n+\n"$qual}' input.fastq
閱讀參考
Bioawk basics
對「Bioawk:專門處理定序相關格式的awk(fasta,fastq,bed,sam,vcf,gff)」的一則回應