seqtk的安装与使用

seqtk的简介

seqtk是一款快速的、轻量级的FASTA或FASTAQ格式文件的处理工具。它来自于生信大神李恒之手，被称为序列处理的瑞士军刀。它处理FASTA/Q文件十分方便，可以大大提高序列分析的效率。

seqtk的安装

(py3env) yu@yu-virtual-machine:~$ sudo apt-get install seqtk

报错：

按照提示中执行

(py3env) yu@yu-virtual-machine:~$ sudo apt --fix-broken install

修复过程中，提示是否希望继续执行，输入“y”

继续执行安装操作，安装正常进行了

检查seqtk是否安装成功

(py3env) yu@yu-virtual-machine:~$ seqtk

展示了seqtk的版本号以及基本用法：

Usage:   seqtk <command> <arguments>
Version: 1.3-r106 #版本号

Command: seq       common transformation of FASTA/Q 
         comp      get the nucleotide composition of FASTA/Q 
         sample    subsample sequences 
         subseq    extract subsequences from FASTA/Q 
         fqchk     fastq QC (base/quality summary)
         mergepe   interleave two PE FASTA/Q files
         trimfq    trim FASTQ using the Phred algorithm

         hety      regional heterozygosity
         gc        identify high- or low-GC regions
         mutfa     point mutate FASTA at specified positions
         mergefa   merge two FASTA/Q files
         famask    apply a X-coded FASTA to a source FASTA
         dropse    drop unpaired from interleaved PE FASTA/Q
         rename    rename sequence names
         randbase  choose a random base from hets
         cutN      cut sequence at long N
         listhet   extract the position of each het

seqtk的使用

以fastq格式文件Akle_TTAGGC_L004_R1_001.fastq.gz为例

1.随机抽取序列

(base) yu@yu-virtual-machine:~/biodata$ gunzip -c Akle_TTAGGC_L004_R1_001.fastq.gz |seqtk sample -s 60 - 500 >test500.fq

代码解读：
gunzip -c 表示将解压后的文件输出到标准输出设备。
seqtk sample 表示随机抽样，其中，-s 60 表示随机数种子为60，- 500 表示抽取500条序列。
>test500.fq 表示将抽取结果重定向至文件test500.qc中

统计test500.fq文件的行数

(base) yu@yu-virtual-machine:~/biodata$ wc -l test500.fq 
2000 test500.fq  #按照fastq文件格式规定一个序列占四行，500个序列正好是2000行

2.fq文件和fa文件相互转换

(base) yu@yu-virtual-machine:~/biodata$ seqtk seq test500.fq>test500.fa

结果显示为：

(base) yu@yu-virtual-machine:~/biodata$ ll
-rw-rw-r--  1 yu yu    131264  9月 30 15:37 test500.fa
-rw-rw-r--  1 yu yu    131264  9月 27 17:21 test500.fq

3.对reads的两端进行适当的修剪

(base) yu@yu-virtual-machine:~/biodata$ seqtk trimfq -b 5 -e 10 test500.fq>cut.fq 
#切除reads前5bp和后10bp的碱基

部剪切前后情况比对（以前三条序列为例）

剪切前：

剪切后：

4.获得反向互补序列

(base) yu@yu-virtual-machine:~/biodata$ seqtk seq -r test500.fq >reverse500.fq

查看反向互补结果（还是以前三条序列为例）

原序列：

互补序列：

值得注意的是，反向互补序列文件依旧沿用原序列的序列标识符，两文件要做好区分，以免混淆

5.统计序列的碱基组成

(base) yu@yu-virtual-machine:~/biodata$ seqtk comp test500.fq >agct500.fq

查看结果（以前10个为例）

结果解读：
第1列：序列名称（序列标识符）
第2列：序列总长度
第3~6列：A、C、G、T碱基数目
（其余几列的含义我还没有搞清楚，搞清楚之后会补充）

6.合并双端测序得到R1和R2序列，实现两两配对

(base) yu@yu-virtual-machine:~/biodata$ seqtk mergefa Akle_TTAGGC_L004_R1_001.fastq.gz Akle_TTAGGC_L004_R2_001.fastq.gz >merge12.fa
[stk_mergefa] (same,diff,hom-het,het-hom,het-het)=(101827379,302054713,0,0,0)

结果(以前四条为例)

参考：
https://www.jianshu.com/p/309b79238553
https://www.jianshu.com/p/856a96ba565f
https://cloud.tencent.com/developer/article/1613420
https://www.cnblogs.com/xudongliang/p/6409534.html
https://www.jianshu.com/p/8d032a29d5a1