seqtk抽取reads

2020-08-06 00:05:57 浏览数 (1)

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

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

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

https://github.com/lh3/seqtk

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

代码语言:javascript复制
  seqtk sample -s100 read1.fq 10000 > sub1.fq
  seqtk sample -s100 read2.fq 10000 > sub2.fq

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

代码语言:javascript复制
  seqtk sample -s100 read1.fq 0.5 > sub1.fq
  seqtk sample -s100 read2.fq 0.5 > sub2.fq

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

代码语言:javascript复制
  seqtk sample -s100 read1.fq.gz 10000 | gzip > sub1.fq.gz
  seqtk sample -s100 read2.fq.gz 10000 | gzip > sub2.fq.gz

0 人点赞