Archives

Categories

对扩增子序列进行引物匹配: usearch 和 primersearch

primer
经典问题场景:
数据:16S 扩增子序列,已完成双端合并,amplicon.fasta 文件, 一对引物序列:primers.txt

我们想做:

  1. 需要知道引物是否匹配,过滤掉一些噪音序列;
  2. 多对引物在同一个文库,需要根据引物把目标序列提取出来;

下面介绍两个比较简单的引物识别程序(双端引物序列模式)

一、USEARCH

强大的数据处理工具集合,但是: 只能免费获取32位版本,最大内存使用4G, 大部分的数据分析任务还是可以应付。

下载地址:32Bit , 选择需要下载的版本,地址会通过邮件进行发送。

如果想购买64bit 版本可以移步:buy64bit , 也可通过代理商购买,现在的销售模式是大版本免费升级,当前大版本为10.0

这里只涉及到引物匹配子命令:search_pcr

命令行接口:

-db, 引物序列文件, fasta格式;
-strand, 支持 plus 和 both 两种模式
-maxdiffs, 单端引物去最大差异
-minamp,匹配区域的最小长度;
-maxamp,匹配区域的最大长度;
-pcrout, 制表符分隔的文件,记录比对的详细信息;
-threads,  线程数;
-ampout,    匹配区域的序列文件;
-log, 输出log 信息;

-pcr_strip_primers 可以切除引物区域,如果引物中存在简并碱基,可能程序会报错,使用需要谨慎。

cat primers.txt 

>fw
GTGCCAGCMGCCGCGG
>rw
CCGTCAATTCMTTTRAGTTT

现在执行:

usearch10.0.240_i86linux32  -search_pcr amplicon.fasta   -db  primers.txt -strand both -threads 5 \
-maxdiffs 2  -minamp 390  -maxamp 430  --pcrout amplicon-pcrout.txt -log  amplicon.log            \
-ampout  amplicon-pcrout.fasta;

现在:amplicon-pcrout.fasta就是我们需要的文件,对于16S序列,可能存在正向和方向序列共存的情况,以及正确切除引物的需求, 这个问题留到 seqtk_primer_strip 介绍。

pcrout 输出文件格式见: pcrout, 后面对序列的匹配情况进行过滤需要使用这个文件,因为同一个序列可能在不同位置都匹配,这样都会report结果。

二、EMBOSS

EMBOSS 是老牌的数据处理工具箱,内容很丰度,小工具相当多, 达到300个之多,今天需要介绍的是: primersearch 做引物匹配:

primersearch 命令行接口如下:

$ primersearch -h
Search DNA sequences for matches with primer pairs
Version: EMBOSS:6.6.0.0

   Standard (Mandatory) qualifiers:
  [-seqall]            seqall     Nucleotide sequence(s) filename and optional
                                  format, or reference (input USA)
  [-infile]            infile     Primer pairs file
  [-mismatchpercent]   integer    [0] Allowed percent mismatch (Any integer
                                  value)
  [-outfile]           outfile    [*.primersearch] Whitehead primer3_core
                                  program output file

   Additional (Optional) qualifiers: (none)
   Advanced (Unprompted) qualifiers: (none)
   General qualifiers:
   -help               boolean    Report command line options and exit. More
                                  information on associated and general
                                  qualifiers can be found with -help -verbose

主要参数:

-seqall 扩增子序列文件;
-infile 引物文件;
-mismatchpercent 错配百分比; 引物长度为20的话,设置0.5,就是一个错配;
-outfile 输出文件;

引物文件序列格式,制表符分隔:

V4-V5   GTGCCAGCMGCCGCGG        CCGTCAATTCMTTTRAGTTT

程序执行:

primersearch  -infile  primers.txt   -seqall  amplicon.fasta   -outfile  amplicon.txt -mismatchpercent  5

结果有点小麻烦:

Primer name V4-V5
Amplimer 1
        Sequence: 2226

        GTGCCAGC[CAM]GCCGCGG hits forward strand at 1 with 0 mismatches
        CCGTCAATTC[CAM]TTT[GAR]AGTTT hits reverse strand at [1] with 0 mismatches
        Amplimer length: 412 bp

序列的标识符丢掉了,对后续的解析程序会造成很大麻烦,使用小工具primersearch_hack, 可以将序列标识符放在序列注释部分,这样就可以解决问题了:

primersearch_hack   amplicon.fasta  >amplicon_hack.fasta

再重新执行,并查看结果,是不是很愉悦了。

Primer name V4-V5
Amplimer 1
        Sequence: 2226
        HISEQ:593:HW27CBCXY:2:1101:18543:2226
        GTGCCAGC[CAM]GCCGCGG hits forward strand at 1 with 0 mismatches
        CCGTCAATTC[CAM]TTT[GAR]AGTTT hits reverse strand at [1] with 0 mismatches
        Amplimer length: 412 bp

后续还是需要根据匹配的位置,错配等进行过滤,以后再讲。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有

Last Upate: 2017-09-29 4:17 PM

Comments are closed.