seqtk抽取reads

做测序数据分析经常要从原始的raw reads里面抽取部分做分析。

比如说不同样本之间的比较，不同平台之间的比较，以及不同的产品之间的比较等等。只有相同的起始reads数进行后续的分析，这样的比较才是一个合理且公正的比较。那么怎么随机抽取一定的数目的reads呢？

今天给大家安利一个小工具，叫seqtk

https://github.com/lh3/seqtk

比如说我们要从pair end的原始fastq文件中抽取10000条reads，可以用下面的命令。其中-s是seed，控制随机抽取，但是要注意在抽R1和R2的时候，一定要用相同的seed，这样才能保证抽出来的R1和R2仍然是配对的，否则有可能会错位。后面10000表示抽取的reads数目。

seqtk sample -s100 read1.fq 10000 > sub1.fq
seqtk sample -s100 read2.fq 10000 > sub2.fq

除了可以指定抽取的reads条数以外，还可以指定抽取的百分比，比如下面的命令就是抽取原始reads的一半。

seqtk sample -s100 read1.fq 0.5 > sub1.fq
seqtk sample -s100 read2.fq 0.5 > sub2.fq

这里还有一个小技巧，如果原始文件是压缩文件，也可以直接使用seqtk进行抽取，不需要先解压。不过抽出来的reads需要使用管道，进行压缩。这样才能保证抽完还是压缩文件。

seqtk sample -s100 read1.fq.gz 10000 | gzip > sub1.fq.gz
seqtk sample -s100 read2.fq.gz 10000 | gzip > sub2.fq.gz