1. 背景

在前面小节我们使用了这些软件,因为混合使用比较让人混乱,这里总结理清楚一下.

2. seqtk

功能总览如下图所示.

2.1 seq

这个功能主要是对\(.fasta\)和\(.fastq\)格式的文件进行格式化.

\(-l\)
主要是让序列每行显示多少个碱基

#每行显示60个氨基酸
seqtk seq -l 60 atha.fasta | less -SN

-L
将短于\(x\)的序列删掉.

# 将短于100的氨基酸序列删除
seqtk seq -L 100 atha.fasta | less -SN

-U
将序列里的小写字母转化为大写,这个应该不用演示了.

2.2 Sample

当我们序列条数过多,全部用于实验会让代码速度减慢,所以有时我们会将序列抽样

从上图理解,格式是seqtk sample in.fa 分数|具体数字,也就是说我们可以抽取多少含量,也可以抽具体数字.

seqtk sample atha.fasta 2

如果我们不改变随机种子,那么每次随机都会是一样的结果.

-s
改变随机种子

2.3 subseq

用此指令提取序列.

可以观察到第一个参数是源文件,第二个参数是对应键名文件,我们根据\(name.list\)去提取文件.

seqtk subseq genome.fa name.list | less -N

我们可以改变name.list的文件内容,让\(subseq\)提取不同位置的碱基.代码保持不变,获得的碱基不同了.

这里有一个点需要科普.我们在文件里列的是1~10,但最后展示出来的是2~10.虽然我们起的名是\(name.list\)但是实际标准输入文件是bed格式的文件.在bed文件中,染色体序列的序号是从\(0\)开始,而输出文件是从\(1\)开始,这就造成了两者偏差.

3. seqkit

\(seqkit\)的很多功能与\(seqtk\)重叠.但是还是有一些功能用\(seqkit\)比较多.

3.1 faidx

此指令用于构建Index.创建 FASTA 索引文件并提取子序列.这个指令类似于samtools faidx,网上关于这个指令的资料没多少,但是samtools的资料蛮多.下面是参考\(samtools\)的资料.
samtools faidx(seqkit faidx也可以) 能够对fasta 序列建立一个后缀为.fai 的文件，根据这个.fai 文件和原始的fastsa文件，能够快速的提取任意区域的序列.该命令对输入的fasta序列有一定要求：对于每条序列，除了最后一行外，其他行的长度必须相同.生成的\(fai\)文件共5列,\t分隔.写着写着感觉我好像在上一节写过,算了不管了,接下来搞懂makeblastdb.