使用NCBI-ePCR和Primer3进行引物批量化设计

NCBI已经不再维护并下架了ePCR软件,转而推荐使用其Primer-BLAST网页工具。这对于单个引物设计任务比较方便,但不利于基因组较大的非模式物种的特异性引物设计或批量化的引物设计。

本教程讲解使用NCBI-ePCR和Primer3进行引物批量化设计。

1.下载并安装NCBI-ePCR和Primer3软件

$ wget http://ftp.debian.org/debian/pool/main/e/epcr/epcr_2.3.12-1.orig.tar.gz -P ~/software
$ tar zxf ~/software/epcr_2.3.12-1.orig.tar.gz -C /opt/biosoft/
$ make -j 4
$ echo 'PATH=$PATH:/opt/biosoft/e-PCR-2.3.12/' >> ~/.bashrc
$ source ~/.bashrc
$ wget https://sourceforge.net/projects/primer3/files/primer3/2.4.0/primer3-2.4.0.tar.gz -P ~/software/
$ tar zxf ~/software/primer3-2.4.0.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/primer3-2.4.0/src/
$ make all
$ echo 'PATH=$PATH:/opt/biosoft/primer3-2.4.0/src/' >> ~/.bashrc
$ source ~/.bashrc

2. 使用ePCR进行引物验证

首先,使用famap命令和fahash命令分两步将基因组序列转换为哈希数据库。

$ famap -t N -b genome.famap genome.fasta
程序将FASTA格式的序列转换为famap数据库文件。

常用参数:
-t <string>
设置碱基类型。可以设置4种值:n,允许含有小写碱基atcgn,其它字符转换为n或N;nx,允许含有小写碱基和兼并碱基字符,其它字符转换为n或N;N,仅允许大写碱基,atcgn自动转换为大写,其它字符转换为N;NX,允许大写碱基和兼并碱基字符,其它字符转换为N。 -b <string>
设置输出的famap数据库文件路径。

$ fahash -b genome.hash -w 12 -f 3 genome.famap
程序进一步将famap文件转换为hash数据库文件。

常用参数:
-b <string>
输出hash数据库文件。
-w <int>
设置wordsize长度。
-f <int>
设置wordcnt长度。

然后,使用re-PCR将引物和数据库进行比对,寻找引物比对结果。

分别对多个引物进行比对,得到各个引物的匹配结果:
$ re-PCR -p genome.hash -n 1 -g 1 ACTATTGATGATGA AGGTAGATGTTTTT …

输入一对引物,并设置产物长度期望值,得到一对引物的匹配结果:
$ re-PCR -s genome.hash -n 1 -g 1 ACTATTGATGATGA AGGTAGATGTTTTT 50-1000
常用参数:
-p
输入hash数据库,在命令行中直接输入引物序列,将引物和数据库进行比对。
-s
输入hash数据库,在命令行中必须输入一对引物和产物长度期望范围,得到一对引物的匹配结果。
-n
设置允许的错配碱基数。
-g
设置允许的gap数。

3. Primer3软件的使用

Primer3是命令行形式的引物设计软件,能很好地用于引物的批量设计。 其主程序是primer3_core。 该命令的输入是Boulder-IO格式,适合于软件读入数据;输出文件默认下是适合人类阅读的格式,也可以是Boulder-IO格式,有助于下一步引物结果的批量操作。 Boulder-IO格式以文本形式记录着引物设计信息,且每个引物信息的结尾使用“等于换行符”分隔。每个记录由多个标签和对应的值构成,用于指定输入信息或输出结果。例如:

SEQUENCE_ID=example
SEQUENCE_TEMPLATE=GTAGTCAGTAGACNATGACNACTGACGATGCAGACNACACACACACACACAGCACACAGGTATTAGTGGGCCATTCGATCCCGACCCAAATCGATAGCTACGATGACG
SEQUENCE_TARGET=37,21
PRIMER_TASK=generic
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=1
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_OPT_SIZE=18
PRIMER_MIN_SIZE=15
PRIMER_MAX_SIZE=21
PRIMER_MAX_NS_ACCEPTED=1
PRIMER_PRODUCT_SIZE_RANGE=75-100
P3_FILE_FLAG=1
SEQUENCE_INTERNAL_EXCLUDED_REGION=37,21
PRIMER_EXPLAIN_FLAG=1

Primer3使用示例和参数:

$ primer3_core -p3_settings_file p3_settings_file -strict_tags -format_output input_file result.p3.out
程序的输入文件要求是Boulder-IO格式,如果没有输入文件,则会从标准输入读取数据。

常用参数:
-p3_settings_file
用于输入primer3的配置文件。这是因为primer3可设置的参数太多了,将大部分参数放入到该配置文件,有利于参数的输入。primer3的默认参数不好,特别是默认参数下没有开启热力学计算。推荐使用配置文件来根据自己的需求来设定相关参数。此外,输入文件中的参数设置能取代此文件中的设定值。
此文件格式:第1行固定为"Primer3 File - http://primer3.sourceforge.net";第2行固定为"P3_FILE_TYPE=settings";第3行是个空行;从第4行开始则是标准的Boulder-IO格式内容。
-format_output
让primer3_core产生人类易读的结果,否则产生机器易读的结果(Boulder-IO格式结果)。
-strict_tags
要求输入文件中的标签要全部正确。设置该参数后,如果有标签不能被程序识别,则报错并停止执行。不设置此参数,则软件忽略不识别的参数。
-p3_settings_file=file_path
指定 primer_core 的配置文件,该配置文件的设定会取代默认设置。当然,
-echo_settings_file
打印出p3_settings_file中的设置信息。如果没有指定设置文件,或含有-format_output,则该参数失效。
-output=file_path
指定输出文件路径,如果不指定,则输出到标准输出。
-error=file_path
指定错误信息输出路径,如果不指定,则输出到stderr中。

Primer3使用的难点是根据自身需求准备p3_settings_file文件。我的一个示例如下(使用时需要去掉#注释部分和空行,且必须第三行留空行):

Primer3 File - http://primer3.sourceforge.net
P3_FILE_TYPE=settings

P3_FILE_ID=P3 Settings from Lianfu Chen
PRIMER_FIRST_BASE_INDEX=1
PRIMER_TASK=generic
PRIMER_NUM_RETURN=5
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=0
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_PICK_ANYWAY=1
PRIMER_THERMODYNAMIC_PARAMETERS_PATH=/opt/biosoft/primer3-2.4.0/src/primer3_config/

### 引物 TM 值设定: ###
PRIMER_TM_FORMULA=1                 # TM 的计算方法, 1 表示使用 the SantaLucia parameters (Proc Natl Acad Sci 95:1460-65)
PRIMER_MIN_TM=55.0  
PRIMER_OPT_TM=60.0
PRIMER_MAX_TM=65.0
PRIMER_PAIR_MAX_DIFF_TM=5.0         # 两个引物之间的 TM 值最多相差 5 摄氏度
PRIMER_WT_TM_LT=0
PRIMER_WT_TM_GT=0
PRIMER_PAIR_WT_DIFF_TM=0.0

### 引物长度设定: ###
PRIMER_MIN_SIZE=18
PRIMER_OPT_SIZE=20
PRIMER_MAX_SIZE=22
PRIMER_WT_SIZE_LT=0
PRIMER_WT_SIZE_GT=0

### 引物 GC 含量设定:###
PRIMER_MIN_GC=30.0
PRIMER_MAX_GC=70.0
PRIMER_WT_GC_PERCENT_LT=0.0
PRIMER_WT_GC_PERCENT_GT=0.0

### 引物的热力学计算: ###
# 开启热力学计算,开启后,根据热力学的值 TH 值来计算罚分。 TH 的罚分方法为:罚分系数 * (1 / (引物TM - 4 - TH值))。
# 该方法好处是,TH 值越大,罚分力度越重。
PRIMER_THERMODYNAMIC_OLIGO_ALIGNMENT=1

# 引物自身进行反向互补,罚分系数计算为 9 / ( 1/(60-4-45) - 1/(60-4-0) ) = 123.2
# 这样计算罚分系数的原则是,若 TM 为最佳 60 摄氏度的时候,所允许的最大罚分值减去最小罚分值为 9,和 3' 端碱基的稳定性的罚分额度一致。
PRIMER_MAX_SELF_ANY=8.00
PRIMER_WT_SELF_ANY=0.0
PRIMER_MAX_SELF_ANY_TH=45.00
PRIMER_WT_SELF_ANY_TH=123.2

# 引物自身进行 3' 端反向互补形成引物二聚体,罚分系数计算为 9 / ( 1/(60-4-35) - 1/(60-4-0) ) = 302.4
PRIMER_MAX_SELF_END=3.00
PRIMER_WT_SELF_END=0.0
PRIMER_MAX_SELF_END_TH=35.00
PRIMER_WT_SELF_END_TH=302.4

# left primer 和 right primer 序列的反向互补
PRIMER_PAIR_MAX_COMPL_ANY=8.00
PRIMER_PAIR_WT_COMPL_ANY=0.0
PRIMER_PAIR_MAX_COMPL_ANY_TH=45.00
PRIMER_PAIR_WT_COMPL_ANY_TH=123.2

# left primer 和 right primer 进行 3' 端反向互补形成引物二聚体
PRIMER_PAIR_MAX_COMPL_END=3.00
PRIMER_PAIR_WT_COMPL_END=0.0
PRIMER_PAIR_MAX_COMPL_END_TH=35.00
PRIMER_PAIR_WT_COMPL_END_TH=302.4

# 发夹结构,罚分系数计算为 9 / ( 1/(60-4-24) - 1/(60-4-0) ) = 672
PRIMER_MAX_HAIRPIN_TH=24.00
PRIMER_WT_HAIRPIN_TH=672

# 3' 端碱基的稳定性
PRIMER_MAX_END_STABILITY=9.0
PRIMER_WT_END_STABILITY=1

### 碱基序列设定: ###
PRIMER_LOWERCASE_MASKING=0              # 模板序列中包含小写字符不影响引物设计
PRIMER_MAX_POLY_X=4                     # 引物序列中不能包含单核苷酸连续长度超过 4 bp
PRIMER_MAX_NS_ACCEPTED=0                # 引物中允许的 N 的数目
PRIMER_WT_NUM_NS=0.0                    # 每个 N 的罚分
PRIMER_MAX_END_GC=5                     # 引物中 3' 端 5bp 碱基中允许的最大 Gs 或 Cs 的数目
PRIMER_GC_CLAMP=0                       # 引物中 3' 端碱基中不能出现连续的 Gs 和 Cs 序列
PRIMER_LIBERAL_BASE=1                   # 是否允许有诸如 N A R Y 等类型的碱基。必须在设定PRIMER_MAX_NS_ACCEPTED 不为 0 后方有效。
PRIMER_LIB_AMBIGUITY_CODES_CONSENSUS=0  # 如果设置为 1,则 C 能与 S 完美匹配,任意碱基能和 N 完美匹配。

### 碱基质量设定: ###
PRIMER_MIN_QUALITY=0                    # primer 序列允许最小的碱基质量
PRIMER_MIN_END_QUALITY=0
PRIMER_QUALITY_RANGE_MIN=0
PRIMER_QUALITY_RANGE_MAX=100
PRIMER_WT_SEQ_QUAL=0.0
PRIMER_WT_END_QUAL=0.0

## 引物位置设置: ###
PRIMER_SEQUENCING_LEAD=50               # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 表明引物的 3' 端距目标区域有 50bp。
PRIMER_SEQUENCING_SPACING=500           # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了在同一条链上的两个引物的距离.
PRIMER_SEQUENCING_INTERVAL=250          # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了在不同链上的两个引物的距离。
PRIMER_SEQUENCING_ACCURACY=20           # 该参数仅在 PRIMER_TASK=pick_sequencing_primers 时有效. 该值决定了引物设计的区间.
PRIMER_OUTSIDE_PENALTY=0
PRIMER_INSIDE_PENALTY=-1.0
PRIMER_WT_POS_PENALTY=0.0

### PCR 反应体系设定: ###
PRIMER_SALT_MONOVALENT=50.0                 # 单价盐离子浓度(mM)
PRIMER_SALT_CORRECTIONS=1                   # PRIMER_SALT_CORRECTIONS=1 means use the salt correction in SantaLucia et al 1998
PRIMER_SALT_DIVALENT=1.5                    # 二价镁离子浓度(mM)
PRIMER_DNTP_CONC=0.6                        # 总dNTPs浓度(mM)
PRIMER_DNA_CONC=50.0                        # DNA产物浓度(mM)

### PCR 产物的设定: ###
PRIMER_PRODUCT_MIN_TM=-1000000.0
PRIMER_PRODUCT_OPT_TM=0.0
PRIMER_PRODUCT_MAX_TM=1000000.0
PRIMER_PAIR_WT_PRODUCT_TM_LT=0.0
PRIMER_PAIR_WT_PRODUCT_TM_GT=0.0
PRIMER_PRODUCT_OPT_SIZE=0
PRIMER_PAIR_WT_PRODUCT_SIZE_LT=0.0
PRIMER_PAIR_WT_PRODUCT_SIZE_GT=0.0

## 罚分因子: ###
PRIMER_PAIR_WT_PR_PENALTY=1.0                # left primer和right primer的罚分之和。此和乘以此系数,再加其它罚分作为最终罚分。
PRIMER_PAIR_WT_IO_PENALTY=0.0                # 将 internal oligo 的罚分乘以此系数,加入到引物的最终罚分中。

### internal oligo 的设定:###
# TM 值
PRIMER_INTERNAL_MIN_TM=57.0
PRIMER_INTERNAL_OPT_TM=60.0
PRIMER_INTERNAL_MAX_TM=63.0
PRIMER_INTERNAL_WT_TM_LT=1.0
PRIMER_INTERNAL_WT_TM_GT=1.0

# 长度
PRIMER_INTERNAL_MIN_SIZE=18
PRIMER_INTERNAL_OPT_SIZE=20
PRIMER_INTERNAL_MAX_SIZE=27
PRIMER_INTERNAL_WT_SIZE_LT=1.0
PRIMER_INTERNAL_WT_SIZE_GT=1.0

# GC 含量
PRIMER_INTERNAL_MIN_GC=20.0
PRIMER_INTERNAL_MAX_GC=80.0
PRIMER_INTERNAL_OPT_GC_PERCENT=50.0
PRIMER_INTERNAL_WT_GC_PERCENT_LT=0.0
PRIMER_INTERNAL_WT_GC_PERCENT_GT=0.0

# 热力学
PRIMER_INTERNAL_MAX_SELF_ANY=12.00
PRIMER_INTERNAL_WT_SELF_ANY=0.0
PRIMER_INTERNAL_MAX_SELF_END=12.00
PRIMER_INTERNAL_WT_SELF_END=0.0

# 碱基序列
PRIMER_INTERNAL_MAX_POLY_X=5
PRIMER_INTERNAL_MAX_NS_ACCEPTED=0

# 碱基质量
PRIMER_INTERNAL_MIN_QUALITY=0
PRIMER_INTERNAL_WT_END_QUAL=0.0
PRIMER_INTERNAL_WT_SEQ_QUAL=0.0

# 非引物数据库
#PRIMER_INTERNAL_MAX_LIBRARY_MISHYB=12.00
#PRIMER_INTERNAL_WT_LIBRARY_MISHYB=0.0

# PCR 反应体系
PRIMER_INTERNAL_SALT_MONOVALENT=50.0
PRIMER_INTERNAL_SALT_DIVALENT=1.5
PRIMER_INTERNAL_DNA_CONC=50.0
PRIMER_INTERNAL_DNTP_CONC=0.0
=

4. 编写程序调用primer3进行引物设计再调用e-PCR软件进行特异性验证

编写程序primer3_with_ePCR_validation.pl,输入模板序列的FASTA文件,即可对所有的模板序列批量化进行引物设计。若同时输入全基因组的序列,进一步可以对设计出的引物进行特异性筛选。

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] seq.fasta > primer3.out

    程序对seq.fasta文件中的序列调用primer3进行引物批量化设计。并将结果输出到标准输出。

    --min_product_length <INT>  default: 80
    设置引物得到的最小产物长度。

    --max_product_length <INT>  default: 150
    设置引物得到的最大产物长度。

    --primer-num <int>    default: 5
    设置使用primer3软件设计的引物对数量。

    --p3_setting_file <STRING>
    程序使用Primer3进行引物批量设计,该参数设置所使用的Primer3配置文件路径。该参数是必须的,以利于primer3设置正确有效的引物。

    --CPU <int>    default: 1
    设置并行化运行primer3的任务数。程序调用ParaFly命令进行并行化运算,需要能直接在terminal中直接运行ParaFly命令。

    --db <string>
    若使用该参数输入FASTA格式的数据库序列,程序能额外对primer3得到的引物进行e-PCR验证。

    --max-mismatch <int>    default: 1
    --max-gap <int>    default: 1
    设置e-PCR过程中允许引物和数据库序列最大的错配和gap数量。设置更高能得到更好的特异性。

    --max-hit <int>    default: 1
    设置允许e-PCR对数据库搜索得到的最大结果数量。若模板序列存在数据库中,推荐设置为1,若模板序列不存在数据库中,推荐设置为0。

    --no-primer-seq <string>
    若设置该参数的值,则生成相应的文件,输出没有引物结果的序列。

    --tmp-dir <string>    default: primer3_with_ePCR_validation.tmp
    设置临时文件夹,程序的中间文件放入该文件夹中。

USAGE
if (@ARGV==0){die $usage}

...

使用AnimalTFDB3.0对动物基因组的转录因子进行注释

1. 动物转录因子数据AnimalTFDB3.0简介

动物转录因子数据库AnimalTFDB3.0对97个动物基因组的转录因子(Transcription Factor)和转录辅助因子(Transcription cofactor)进行了归纳整理。基于DNA结合结构域,将动物转录因子分成了73个基因家族,将转录辅助因子分成了83个基因家族。此外,动物转录因子分为六大类(Basic Domain Group、Zinc-Coordinating Group、Beta-Scaffold Factors、Helix-turn-helix、Other Alpha-Helix Group和Unclassified Structure),动物转录辅助因子也分为六大类(Co-activator/repressors、Chromatin Remodeling Factors、General Cofactors、Histone-modifying Enzymes、Cell Cycle和Other Cofactors)。

动物转录因子数据库AnimalTFDB3.0提供了网页工具进行转录因子分析。该网页工具一次仅支持上传不超过1000条蛋白序列。该网页工具采用HMM算法进行计算,分析速度比较快。数据库作者对大部分转录因子家族采用了PFam数据库的HMM模型,并自己构建了一些转录因子家族的HMM模型。但作者没有提供HMM数据库的下载,仅提供了97个动物转录因子和转录辅助因子蛋白序列的FASTA格式文件下载。

使用AnimalTFDB3.0数据库的网页工具仅能进行转录因子分析,且由于仅能支持不超过1000条序列进行分析而比较麻烦,不利于动物全基因组的转录因子分析。以下讲解下载AnimalTFDB3.0数据库FASTA文件,并自行编写程序进行转录因子和转录辅助因子注释。

2. 下载AnimalTFDB3.0数据库蛋白序列

首先,进入AnimalTFDB3.0数据库下载页面,采用复制粘贴方式获取97个物种的物种名,保存到文件species.txt中。

Ailuropoda melanoleuca;Anas platyrhynchos;Anolis carolinensis;Aotus nancymaae;Astyanax mexicanus;Bos taurus;Caenorhabditis elegans;Callithrix jacchus;Canis familiaris;Capra hircus;Carlito syrichta;Cavia aperea;Cavia porcellus;Cebus capucinus;Cercocebus atys;Chinchilla lanigera;Chlorocebus sabaeus;Choloepus hoffmanni;Ciona intestinalis;Ciona savignyi;Colobus angolensis palliatus;Cricetulus griseus chok1gshd;Cricetulus griseus crigri;Danio rerio;Dasypus novemcinctus;Dipodomys ordii;Drosophila melanogaster;Echinops telfairi;Equus caballus;Erinaceus europaeus;Felis catus;Ficedula albicollis;Fukomys damarensis;Gadus morhua;Gallus gallus;Gasterosteus aculeatus;Gorilla gorilla;Heterocephalus glaber female;Heterocephalus glaber male;Homo sapiens;Ictidomys tridecemlineatus;Jaculus jaculus;Latimeria chalumnae;Lepisosteus oculatus;Loxodonta africana;Macaca fascicularis;Macaca mulatta;Macaca nemestrina;Mandrillus leucophaeus;Meleagris gallopavo;Mesocricetus auratus;Microcebus murinus;Microtus ochrogaster;Monodelphis domestica;Mus caroli;Mus musculus;Mus pahari;Mus spretus;Mustela putorius furo;Myotis lucifugus;Nannospalax galili;Nomascus leucogenys;Notamacropus eugenii;Ochotona princeps;Octodon degus;Oreochromis niloticus;Ornithorhynchus anatinus;Oryctolagus cuniculus;Oryzias latipes;Otolemur garnettii;Ovis aries;Pan paniscus;Pan troglodytes;Papio anubis;Pelodiscus sinensis;Peromyscus maniculatus bairdii;Petromyzon marinus;Poecilia formosa;Pongo abelii;Procavia capensis;Propithecus coquereli;Pteropus vampyrus;Rattus norvegicus;Rhinopithecus bieti;Rhinopithecus roxellana;Saimiri boliviensis boliviensis;Sarcophilus harrisii;Sorex araneus;Sus scrofa;Taeniopygia guttata;Takifugu rubripes;Tetraodon nigroviridis;Tupaia belangeri;Tursiops truncatus;Vicugna pacos;Xenopus tropicalis;Xiphophorus maculatus

然后,根据物种名,编写下载数据库文件的脚本(每个物种下载4个数据文件)并批量下载数据。

$ mkdir /opt/biosoft/AnimalTFDB3.0
$ cd /opt/biosoft/AnimalTFDB3.0

生成文件species.list,包含97个物种名,每行一个物种名。根据http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/download网>页中的信息,使用vim编辑和perl单行命令抓取全部的物种名信息。
$ perl -p -i -e 's/;/\n/g;' species.txt

$ mkdir data_of_97_species
$ cd data_of_97_species
$ perl -e 'while (<>) { chomp; s/\s+/_/g; print "wget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_Cof_protein.fa\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF\nwget http://bioinfo.life.hust.edu.cn/static/AnimalTFDB3/download/$_\_TF_cofactors\n"; }' ../species.list > download.list
$ ParaFly -c download.list -CPU 2

Caenorhabditis_elegans的FASTA文件有问题,存在序列名一致的多条序列,推荐进行修改。
$ perl -p -i -e 's/>(\S+)\t(\S+)/>$2\t$1/' data_of_97_species/Caenorhabditis_elegans_*_protein.fa
$ cd ..

分别合并转录因子数据库和转录辅助因子数据库的FASTA文件。

得到转录因子序列数据库,FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*TF_protein.fa > AnimalTFDB3.0_TF_protein.fasta
$ perl -p -i -e 's/^(>\S+)\t.*?\t.*?\t/$1\t/' AnimalTFDB3.0_TF_protein.fasta


得到转录辅助因子数据库,FASTA文件头部仅包含序列名和基因家族信息。
$ cat data_of_97_species/*_Cof_protein.fa > AnimalTFDB3.0_Cof_protein.fasta
由于Cof_protein.fa文件中不包含转录辅助因子基因家族信息,需要编写程序读取TF_cofactors文件信息并将基因家族信息写入FASTA文件
$ perl -e 'while (<>) { @_ = split /\t/; $gf{$_[2]} = $_[3]; } open IN, "AnimalTFDB3.0_Cof_protein.fasta"; while (<IN>) { if (m/^>(\S+)\t(\S+)/) { if (exists $gf{$1}) {print ">$1\t$gf{$1}\n";} elsif (exists $gf{$2}) {print ">$2\t$gf{$2}\n";} else {print STDERR "No type: $_\n";} } else { print; } }' data_of_97_species/*_TF_cofactors > aa; mv aa AnimalTFDB3.0_Cof_protein.fasta

3. 使用DIAMOND将数据库所有蛋白和数据库自身进行比对,计算每个基因家族的BLAST阈值

使用ALL_VS_ALL的方式分别对转录因子数据库和转录辅助因子数据库进行DIAMOND比对(evalue 1e-5、identity 30%且coverage 30%)。根据DIAMOND比对结果,得到每个转录因子家族中每个基因的比对信息(evalue、hitNum),编写程序进而计算每个基因家族的BLAST阈值。

$ diamond --in AnimalTFDB3.0_TF_protein.fasta --db AnimalTFDB3.0_TF_protein
$ diamond blastp --db AnimalTFDB3.0_TF_protein --query AnimalTFDB3.0_TF_protein.fasta --out TF.tab --outfmt 6 --sensitive --max-target-seqs 50000 --evalue 1e-5 --id 30 --subject-cover 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1

编写程序geneFamily_BLAST_threshold_of_sensitivity.pl,输入DIAMOND结果文件和数据库FASTA文件,得到各个sensitivity级别上每个基因家族的evalue和hitNum阈值。

4. 将全基因组蛋白序列和数据库进行比对,再根据阈值文件进行准确注释。

使用FALCON对三代测序数据进行基因组组装

1. FALCON软件简介

FALCON是PacBio公司开发的一款用于三代基因组De novo组装软件。相比于HGAP4软件,FALCON软件的基因组组装原理基本一致。但FALCON使用命令行运行,更适合于大基因组的组装,且能分析双倍体序列,并在基因组组装结果中给出包含变异位点信息的等位基因序列(alternative contigs / a-contigs)和主要的基因组序列(primary contig / p-contig)。每一条a-contig都有其对应的p-contig序列。因此,FALCON软件适合双倍体物种的基因组组装,能给出单倍的基因序列。其基因组组装结果中的p-contigs序列总长度要小于其它基因组组装软件(例如Canu和HGAP)的基因组序列。

FALCON-Unzip则是真正的单倍型组装软件,它能在FALCON或HGAP4软件的基因组组装结果基础上,利用较长的PacBio reads进行单倍型分析,对p-contigs序列向单倍型进行转换,同时输出单倍型序列(haplotig)区块。

2. FALCON软件下载与安装

FALCON软件在2018年在算法上有了较大的更改,将FALCON在内的多个软件整合到了pb-assembly软件(https://github.com/PacificBiosciences/pb-assembly)中。同时提供快捷的bioconda安装方法:

$ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ sh Miniconda3-latest-Linux-x86_64.sh
    进入交互式界面:首先,输入yes,按Enter键,表示同意软件的license;然后,输入conda程序安装路径/opt/biosoft/miniconda3_for_pb-assembly,按enter键; 再设置是否添加Minicoda3到PATH变量,直接按Enter键,表示选择no,以免和系统自带的软件冲突。若需要切换到miniconda3环境,输入命令export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH即可。
$ export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH
$ conda config --add channels defaults
$ conda config --add channels bioconda
$ conda config --add channels conda-forge
    配置bioconda channel(https://bioconda.github.io/index.html)
$ conda install pb-assembly
    安装pb-assembly软件
$ source activate /opt/biosoft/miniconda3_for_pb-assembly
    激活FALCON软件的运行环境

3. FALCON软件的使用

FALCON软件运行示例:

首先,载入FALCON软件的运行环境。

然后,准备三种输入文件:测序数据FASTA格式文件、包含测序数据文件路径的FOFN列表文件、程序运行的参数配置文件。下载基因组大小为4.6Mb的E. coli基因组PacBio测序数据。
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139\
_c100713652400000001823152404301535_s1_p0.subreads.tar | tar xvf -
$ ls */ecoli.?.fasta > input.fofn
$ wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg
当使用最新版的FALCON进行组装时,需要对此配置文件进行一些修改。

最后,对E. coli基因组进行组装
$ fc_run.py fc_run_ecoli_local.cfg

生成的最终主要结果文件为 2-asm-falcon/p_ctg.fa 。

运行FALCON最关键的步骤是准备配置文件fc_run.cfg。一个配置文件示例及其参数讲解如下:

[General]
输入参数:
input_fofn = input.fofn
设置数据输入文件为input.fofn,该文本文件中每行是一个PacBio测序数据的FASTA文件路径。可以使用相对路径。
input_type = raw
设置输入数据类型为原始测序数据(raw),或者为修正后的数据(preads)。

对数据进行分块:
pa_DBsplit_option = -x500 –s200
设置预组装(pre-assembly)过程中对raw reads数据进行分块的参数。这些参数会传递给Dbsplit命令来执行数据分块。
-s<int>    default: 200
表示每份数据含有50Mb的数据量,该参数默认值为200。默认参数情况下,使用daligner进行比对需要消耗约16G内存。推荐设置该参数的值为预测的基因组大小。此外,软件作者推荐对小于10Mb的基因组设置该参数值为50,对于大基因组设置该参数值为200。
-x<int>
忽略长度小于指定阈值的reads。
ovlp_DBsplit_option=-x500 –s200
设置OLC组装过程中对preads数据进行分块的参数。

对重复序列进行屏蔽:
pa_HPCTANmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。该命令是正常daligner命令的一个变种,适合对串联重复序列进行比对。
pa_HPCREPmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对散在重复进行屏蔽的参数。其参数会传递给HPC.REPmask命令。该分析步骤会迭代3次。HPC.TANmask命令进行数据分析时有两个关键参数-g和-c,用于设置分析的数据块数量和覆盖度阈值。pa_REPmask_code设置3次迭代过程中这两个参数的值。这两个参数值的设置与数据分块方法和测序数据量有关。
pa_REPmask_code = 1,666;4,266;8,53
设置三次HPC.REPmask命令中的-g和-c参数的值。以上参数值则表示不对散在重复进行屏蔽。该参数的设置可以参考网址https://dazzlerblog.wordpress.com/2016/04/01/detecting-and-soft-masking-repeats/。在该文章中,作者给出了一个示例:对30Gb的基因组测序30x并将数据分割成3600个大小为250Mb的数据块,使用参数值“1,20;10,15;100,10”。每个数据块的测序深度为0.008x,对1个数据块的数据进行比对后,对reads中覆盖度超过20x的位点进行屏蔽,即对高度重复(20 / 0.008 = 2500x)的基因组序列位点进行屏蔽;然后第二次迭代中,对10个数据块的数据进行比对,对reads中覆盖度超过15x的位点进行屏蔽,即对中度重复(15 / 0.08 = 187.5x)的基因组序列位点进行屏蔽;最后第三次迭代中,对100个数据块的数据进行比对,对reads中覆盖度超过10x的位点进行屏蔽,即对低度重复(10 / 0.8 = 12.5x)的基因组序列位点进行屏蔽。示例中基因组大小是30Gb,超级大,一般含有大量重复序列,使用上述参数对数据中的重复位点进行屏蔽。
对于常规的基因组组装,设置pa_REPmask_code参数后,在三次迭代过程中,推荐逐步对重复次数为1000x、100x和10x的序列位点进行屏蔽。例如:对300Mb大小的基因组测序100x,将数据分割成150个大小为200Mb的数据块,可以考虑设置pa_REPmask_code值为“1,666;4,266;8,53”。
ovlp_HPCTANmask_option = -k20 -h480 -w8 -e.80
设置在OLC组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。

预组装(pre-assembly)参数:
length_cutoff = -1
设置对种子序列的长度筛选阈值。若该参数值设置为-1,则程序自动计算种子序列的长度筛选阈值,根据如下两个参数挑选最长的reads序列作为种子序列,直到其数据量达到基因组指定覆盖度为止。
genome_size = 4652500
seed_coverage = 30
设置基因组大小和种子序列覆盖度。推荐设置seed_coverage参数的值为20~40x。当然是用更多的种子序列有利于基因组的组装效果,但是会消耗更多运行时间。
pa_daligner_option = -k14 -w6 -h35 -e.70 -l1000
在新版本FALCON中该参数额外再分出了pa_HPCdaligner_option参数,专门用于对数据块的分配。pa_daligner_option参数专门用于调用daligner软件对raw reads进行比对,进行overlap分析。其常用参数:
-k<int>    default: 14
设置k-mer大小。推荐设置为14~18。设置较低的值,能提高sensitivity,得到更多的reads重叠结果,但增加磁盘、内存、计算和时间消耗,适合数据质量较差的情况;设置较高的值,能提高specificity,系统资源消耗更少,运行速度更快,仅适合数据质量较高的情况。
-w<int>    default: 6
daligner命令对两个reads进行比对,得到一个斜对角线的比对带。该带的宽度设置为2的w次方。
-h<int>    default: 35
在斜对角线带的比对区域上,需要至少有指定数目的碱基能被k-mers覆盖。-k、-w和-h是daligner命令进行比对的主要参数。默认设置下,寻找两两reads之间重叠,先将reads分割成长度为14bp的k-mers序列;对其中一条序列宽度为64bp的区域进行分析,要求至少有35个碱基位点能和被另外一条序列的k-mers覆盖,则找到重叠。
-k、-w和-h的默认参数适合于PacBio测序的raw reads。对于修正后的read,推荐使用参数“daligner -k20 -h60 -e.95”。
-e<float>    default: 0.70
设置identity阈值。推荐设置为0.70~0.80。设置为较高的值有利于单倍型序列的组装。
-l<int>    default: 1000
忽略长度小于低于设定值的reads。
-T<int>    default: 4
设置每个daligner使用的CPU线程数。
-M<int>
尝试仅使用指定大小GB的内存。程序通过忽略覆盖度过高的k-mers来降低内存使用。一般情况下,对
pa_HPCdaligner_option = -v -B4 -M16
程序调用HPC.daligner进行分析,相比于daligner命令,其增加的常用参数:
-v
将-v参数传递给LAsort和LAmerge命令。
-B<int>    default: 4
设置每个daligner任务对指定数目的数据块进行分析。在最新版本的FALCON中,该参数不能再设置到pa_daligner_option参数中,否则程序运行出错。
falcon_sense_option = --output-multi --min-idt 0.70 --min-cov 4 --max-n-read 200
FALCON使用fc_consensus命令根据daligner进行Overlapping分析的结果进行一致性分析,从而对种子序列进行校正。常用参数:
--output-multi
输出每个种子序列多个修正后的区域。
--min-cov    default: 6
设置当种子序列上某位点覆盖度低于指定阈值时,则在该位点对序列打断或截短。
--min-cov-aln    default: 10
设置当种子序列的平均覆盖度低于指定阈值时,过滤掉该序列。
--min-n-read    default: 10
--max-n-read    default: 500
对覆盖度在--min_n_read和--max_n_read两个参数设定范围内的位点进行一致性分析,从而对种子序列进行校正。
--min-idt    default: 0.7
设置比对结果中reads重叠部分的最小identity阈值。
--n-core    default: 24
设置运行的线程数,默认值是24。
falcon_sense_greedy = False
falcon_sense_skip_contained = False

OLC算法进行组装的参数:
length_cutoff_pr = 12000
对预组装结果中长度超过此阈值的reads使用OLC算法进行下一步组装。
ovlp_HPCdaligner_option = -v -B4 -M16
ovlp_daligner_option = -k20 -w6 -h60 -e.96 -l1000
调用daligner对校正后的preads进行overlapping分析时,参数相对要设置更加严格。
-k<int>    default: 14
对preads进行比对,推荐该参数的值为18~24。
-e<float>    default: 0.70
对preads进行比对,推荐该参数的值为0.93~0.96。
overlap_filtering_setting = --max-diff 100 --max-cov 100 --min-cov 2 --bestn 10
过滤overlaps。若reads首尾两端的覆盖度比平均覆盖度大很多,则表明reads首尾是重复序列;若reads首尾两端的覆盖度比平均覆盖度相差较小很多,则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。
--bestn    default: 10
设置报告reads此数目的最优overlaps。
--min-cov和—max-cov
所允许的reads首尾的覆盖度范围。
--max-diff
设置所允许的首尾覆盖度值的最大差异。
fc_ovlp_to_graph_option = --min_len 4000 --min_idt 0.96

任务投递与计算资源消耗参数:
[job.defaults]
job_type = local
设置任务提交方式。local表示使用本地主机运行FALCON,运行相对稳定。也可以设置为集群任务提交方式,可以设置的值有:sge、lsf、pbs、torque和slurm。该参数设置不同的值,后面的submit也要设置对应的信息。
pwatcher_type=blocking
submit = bash -C ${CMD} >| ${STDOUT_FILE} 2>| ${STDERR_FILE}
MB=32768
NPROC=6
njobs=32
FALCON软件流程对数据进行了切割,再进行并行化运算,加快程序运行速度。MB设置任务实用的内存数;NPROC设置单个任务的CPU线程数;njobs设置并行任务数。后续分别对每个流程中的计算资源进行分配。da表示预组装过程中运行的daligner命令;la表示预组装过程中运行的Lasort和LAmerge命令;cns表示预组装过程中运行的fc_consensus命令;pda表示运行OLC组装过程中运行的的daligner命令;pla表示OLC组装过程中运行的Lasort和LAmerge命令;asm表示最终组装过程。
[job.step.da]
NPROC=4
MB=32768
njobs=240
[job.step.la]
NPROC=4
MB=32768
njobs=240
[job.step.cns]
NPROC=8
MB=65536
njobs=240
[job.step.pda]
NPROC=4
MB=32768
njobs=240
[job.step.pla]
NPROC=4
MB=32768
njobs=240
[job.step.asm]
NPROC=24
MB=196608
njobs=1

FALCON的结果文件:

0-rawreads/
该目录存放对raw subreads进行overlpping分析与校正的结果;
0-rawreads/cns-runs/cns_*/*/*.fasta存放校正后的序列信息。

1-preads_ovl/
该目录存放对校正后reads进行overlapping的结果;

2-asm-falcon/
该目录是最终结果目录,主要的结果文件是p_ctg.fa和a_ctg.fa。

使用SignalP对蛋白序列进行信号肽预测

1. 信号肽简介

信号肽是蛋白质N-末端一段编码长度为5-30的疏水性氨基酸序列,用于引导新合成蛋白质向通路转移的短肽链。信号肽存在于分泌蛋白、跨膜蛋白和真核生物细胞器内的蛋白中。

信号肽指引蛋白质转移的方式有两种:(1)常规的分泌(Sec/secretory)通路;(2)双精氨酸转移(Tat/twin-arginine)通路。前者存在于原核生物蛋白质转移到质膜过程中,以及真核生物蛋白质转移到内质网膜的过程中。后者存在于细菌、古菌、叶绿体和线粒体中,信号肽序列较长、疏水性较弱且尾部区含有两个连续精氨酸。相比于前者转运非折叠蛋白质,后者能转运折叠蛋白质跨越双层脂质膜 。

信号肽指引蛋白质转运后,将由信号肽酶进行切除。信号肽酶有三种:(1)一型信号肽酶(SPaseI);(2)二型信号肽酶(SPaseII);(3)三型信号肽酶(SPaseIII)。大部分信号肽由SPaseI进行移除,SPaseI存在古菌、细菌和真核生物中,且在真核生物的内质网膜上仅存在一型信号肽酶。细菌和古菌脂蛋白的信号肽C端含有一段称为 lipobox 的保守区域,由SPaseII切除其信号肽,且lipobox紧邻切除位点(CS/Cleavage Site)的氨基酸是半胱氨酸,这和锚定到膜的功能是相关的。细菌的四型菌毛蛋白信号肽由SPaseIII进行切除。此外:分泌通路(Sec)相关信号肽能由SPaseI、SPaseII和SPaseIII切除,但是双精氨酸转移(Tat)通路相关信号肽仅由 SPaseI和SPaseII切除。

使用SignalP 5.0能对原核生物的信号肽Sec/SPI、Sec/SPII和Tat/SPI,和对真核生物仅含有 Sec/SPI信号肽进行预测 。 SignalP 5.0目前不能对Tat/SPII进行预测。此外,由于没有足够大的数据进行训练,SignalP 5.0 也不能对Sec/SPIII进行分析。

SignalP是最为常用的信号肽分析软件,常用于分泌蛋白预测的第一步。到目前2019年4月,SignalP版本到了第5版。第一版本基于人工神经网络(Artificial Neural Network);第二版本引入隐马尔科夫模型(Hidden Markov Models);第三版本改进切除位点(Cleavage site)预测方法;第四版本改进对信号肽和跨膜螺旋的区分能力。这四个版本仅能对Sec/SPI类型的信号肽进行预测。而第五版本结合了深度神经网络(deep neural network)、条件随机分类(Conditional random field classification)和迁移学习(transfer learning)方法,能对信号肽进行更准确的预测。

可以使用SignalP网页工具进行分析。但一次仅支持最多5000条序列。以下讲解本地化运行SignalP软件。

2. SignalP软件下载和安装

# 需要填写edu邮箱和相关信息来获取下载地址
$ wget http://www.cbs.dtu.dk/download/6B91F6BC-5A05-11E9-8172-2ED6B9CD16B5/signalp-5.0.Linux.tar.gz -P ~/software/
$ tar zxf /home/train/software/signalp-5.0.Linux.tar.gz -C /opt/biosoft/
$ echo 'PATH=$PATH:/opt/biosoft/signalp-5.0/bin' >> ~/.bashrc
$ source ~/.bashrc

3. SignalP软件使用

对真核生物的全基因组蛋白序列进行信号肽预测:
$ signalp -batch 30000 -org euk -fasta proteins.fasta -gff3 -mature
signalp命令的参数:

-batch <int> default: 10000
程序并行化运行的序列条数。程序能多线程运行,速度很快。推荐设置该参数值大于FASTA文件的序列总条数,虽然增加内存消耗,但能加快程序运行。
-org <string> default: euk
设置分析的物种类型。该参数值有4种:arch,古菌;gram+,革兰氏阳性细菌,gram-,革兰氏阴性细菌;euk,真核生物。
-fasta <string>
输入FASTA格式的蛋白序列文件。
-prefix <string>
设置输出文件前缀。默认在程序运行目录下输出结果文件,和输入文件名的前缀相同,后缀为_summary.signalp5。
-gff3
添加该参数,输出GFF3格式的信号肽预测结果。
-mature
添加该参数,对含有信号肽的蛋白序列,切除信号肽后输出其成熟蛋白序列。可以用于下游的跨膜区分析。
-tmp <string> default: /tmp
设置临时文件夹路径。
-format <string> default: short
设置输出格式。该参数值有两个:short,仅输出一个信号肽预测的整合文本结果;long,额外输出每条序列的各位点预测文本结果和图片结果;
-plot <string> default: png
设置输出图片结果的类型。当--format参数为long时,该参数生效。该参数值有三个:png;eps;none表示仅得到表格文件。
-version
打印程序版本并退出。

网页服务器禁止IP访问

我的wordpress站点总是有很多垃圾评论。于是需要对相应的IP地址禁止访问。

我的方法是:在/etc/httpd/conf.d目录下生成后缀为多个.conf的文件,每个文件分别对一个目标文件夹进行保护。比如生成文件/etc/httpd/conf.d/deny_chenlianfu.conf,其内容为:

<Directory "/home/chenlianfu/wordpress">
    AllowOverride None
    Options MultiViews FollowSymLinks ExecCGI
    Order allow,deny
    Allow from all
    Deny from 100.42.17.90
    Deny from 101.4.136.34
    Deny from 101.91.215.188
    Deny from 101.98.247.14
    Deny from 95.85.80.227
    Deny from 95.85.80.82
    Deny from 95.85.80.86
    Deny from 98.174.90.36
</Directory>

然后重启httpd服务,使配置文件生效,从而禁止配置文件中的IP访问我的网站。

/etc/init.d/httpd restart

为了对多个文件夹进行保护,则要生成多个配置文件。我编写名为wordpress_deny_ips.pl的程序,输入含有准备禁止访问IP地址的文本文件,自动生成多个配置文件:

#!/usr/bin/perl
use strict;

my $usage = <<USAGE;
USAGE:
    perl $0 web.txt

USAGE
if (@ARGV==0){die $usage}

open IN, "/etc/httpd/conf.d/deny_chenlianfu.conf";
my %ip;
while (<IN>) {
    $ip{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
}
close IN;

open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
my %ip_new;
while (<IN>) {
    $ip_new{$1} = 1 if m/(\d+\.\d+\.\d+\.\d+)/;
}
close IN;

my $number = 0;
foreach (keys %ip_new) {
    if (exists $ip{$_}) {
        print STDERR "Duplicate $_\n";
    }
    else {
        $ip{$_} = 1;
        $number ++;
    }
}

my @num = keys %ip;
my $num = @num;
print STDERR "$number ips were add\n$num ips in total\n";

my $out = '
    AllowOverride None
    Options MultiViews FollowSymLinks ExecCGI
    Order allow,deny
    Allow from all
';
foreach (sort keys %ip) {
    $out .= "    Deny from $_\n";
}
$out .= "</Directory>\n";


open OUT, ">", "/etc/httpd/conf.d/deny_chenlianfu.conf" or die "Can not create file /etc/httpd/conf.d/deny_chenlianfu.conf, $!\n";
print OUT '<Directory "/home/chenlianfu/wordpress">' . $out;
close OUT;

open OUT, ">", "/etc/httpd/conf.d/deny_zhengyue.conf" or die "Can not create file /etc/httpd/conf.d/deny_zhengyue.conf, $!\n";
print OUT '<Directory "/home/zhengyue/wordpress">' . $out;
close OUT;

open OUT, ">", "/etc/httpd/conf.d/deny_wuchangsong.conf" or die "Can not create file /etc/httpd/conf.d/deny_wuchangsong.conf, $!\n";
print OUT '<Directory "/home/wuchangsong/wordpress">' . $out;
close OUT;

很多垃圾评论是全英文或含有日文。可以在主体对应的functinons.php文件中添加一些设置来禁止全英文或包含日文的评论。我在文件./wp-content/themes/twentyeleven/functions.php的尾部添加以下代码:

function refused_english_comments( $incoming_comment ) {
    $pattern = '/[一-龥]/u';
    //禁止全英文评论
    if(!preg_match($pattern, $incoming_comment['comment_content'])) {
        wp_die( "您的评论中必须包含汉字!" );
    }

    $pattern = '/[あ-んア-ン]/u';
    //禁止日文评论
    if(preg_match($pattern, $incoming_comment['comment_content'])) {
        wp_die( "评论禁止包含日文" );
    }
    return( $incoming_comment );
}
add_filter('preprocess_comment', 'refused_english_comments');

然后重启网页服务,即可生效。

检测全基因组序列中的端粒序列

使用三代测序数据能获得较好的、甚至完整的基因组序列。通过检测基因组序列两端的属于端粒(Telomere)的特定的重复序列,可以知道基因组组装是否得到完整的染色体水平的序列。若在序列中间检测到端粒序列,可以知道基因组组装过程中对Contigs有错误的连接。

在染色体序列首尾存在端粒序列。在人类中,端粒序列由重复单元TTAGGG,串联重复约2500次组成。 不同的物种,端粒的重复单元可能不一样,可以在端粒数据库中查询。 我对一种子囊菌使用PacBio测序数据进行了基因组组装,得到12条序列,发现大部分序列首尾均出现了TTAGGG/CCCTAA的串联重复序列。我对一种昆虫物种PacBio组装基因组序列进行分析,发现序列首尾出现一些长度较长微卫星重复序列,而没有固定的重复单元,这和端粒数据库中的结果一致。

编写Perl程序对全基因组序列的端粒序列进行搜索,查看基因组序列的完整情况:

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 genome.fasta > telomere_info.txt

    大部分物种端粒序列的重复单元是TTAGGG/CCCTAA。本程序能在基因组中寻找端粒重复单元的串联重复序列,并给出位点信息。

    --split-length <int>    default: 100000
    --overlap-length <int>    default: 10000
    程序会将每条序列打断后进行重复单元搜索。这两个参数设置打断的序列长度和相邻两序列之间的重叠长度。

    --repeat-unit <string>    default: TTAGGG
    设置重复单元碱基序列,该重复单元的反向互补也将作为重复单元进行搜索。可以在端粒数据库(http://telomerase.asu.edu/sequences_telomere.html)中寻找目标端粒重复单元。
    vertebrate sp.      TTAGGG
    plants sp.          TTTAGGG
    Pezizomycotina      TTAGGG

    --min-repeat-num <int>    default: 4
    设置重复单元最小重复次数。


USAGE
if (@ARGV==0){die $usage}

my ($splitLength, $overlapLength, $repeatunit, $minRepeatNum);
GetOptions(
    "split-length:i" => \$splitLength,
    "overlap-length:i" => \$overlapLength,
    "repeat-unit:s" => \$repeatunit,
    "min-repeat-num:s" => \$minRepeatNum,
);
$splitLength ||= 100000;
$overlapLength ||= 10000;
$repeatunit ||= "TTAGGG";
$repeatunit = uc($repeatunit);
my $repeatunit_reverse = reverse $repeatunit;
$repeatunit_reverse =~ tr/ATCG/TAGC/;
$minRepeatNum ||= 4;

# 读取基因组序列
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
my (%seq, $seq_id);
while (<IN>) {
    chomp;
    if (m/^>(\S+)/) {
        $seq_id = $1;
    }
    else {
        $_ = uc($_);
        $seq{$seq_id} .= $_;
    }
}
close IN;

# 将基因组序列打断
my (%seq_split, %seq_length);
foreach my $id (keys %seq) {
    my $seq = $seq{$id};
    my $length = length($seq);
    $seq_length{$id} = $length;
    my $pos = 0;
    while ($pos < $length) {
        $seq_split{$id}{$pos} = substr($seq, $pos, $splitLength + $overlapLength);
        $pos += $splitLength;
    }
}

print "SeqID\tSeqLength\tStart\tEnd\tLength\tType\n";
foreach my $id (sort keys %seq_split) {
    foreach my $pos (sort {$a <=> $b} keys %{$seq_split{$id}}) {
        my $seq = $seq_split{$id}{$pos};
        while ($seq =~ m/(($repeatunit){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit\n";
        }
        while ($seq =~ m/(($repeatunit_reverse){$minRepeatNum,})/g) {
            my $length = length($1);
            my $end = pos($seq);
            $end = $end + $pos;
            my $start = $end - $length + 1;
            print "$id\t$seq_length{$id}\t$start\t$end\t$length\t$repeatunit_reverse\n";
        }
    }
}

使用DIAMOND将全基因组蛋白序列比对到Nr数据库

1. DIAMOND简介

对全基因组的基因进行Nr注释是必不可少的一步。由于Nr数据库非常大,导致使用BLAST会消耗巨大的计算资源和时间。使用DIAMOND则能快500-20000倍,而获得和BLAST比较一致的结果。特别是对于长度为100-150bp,数量超过1M的核酸序列,DIAMOND的速度比BLASTX快20000倍;当e-value设置为1e-5时,DIAMOND能报告80~90%的BLASTX匹配结果;若设置sensitivie模式,DIAMOND速度是BLASTX的2500倍,DIAMOND能报告94%的BLASTX匹配结果。因此,对大批量的序列进行Nr注释或规模较大的蛋白数据库的比对,使用DIAMOND是优先选择的方式。

虽然DIAMOND软件的运行速度快,但是其结果中没有匹配的Subject序列的描述信息。若需要蛋白名称和物种的描述信息,需要自行编写软件从数据库的FASTA文件中提取并整合到DIAMOND结果中。

2. DIAMOND下载和安装

$ wget https://github.com/bbuchfink/diamond/releases/download/v0.9.24/diamond-linux64.tar.gz -P ~/software/
$ mkdir /opt/biosoft/diamond
$ tar zxf ~/software/diamond-linux64.tar.gz -C /opt/biosoft/diamond
$ echo 'PATH=$PATH:/opt/biosoft/diamond' >> ~/.bashrc
$ source ~/.bashrc

3. DIAMOND使用

DIAMOND软件的主命令是diamond,它的使用包含几个子命令。DIAMOND最常用的使用方法:

使用DIAMOND软件的子命令makedb将FASTA格式的蛋白序列创建成后缀为dmnd的数据库文件:
$ diamond makedb --in nr_eukaryon.fasta -d nr_eukaryon_20190405
    将Nr数据库(版本20190405)的真核子集创建成DIAMOND数据库,数据库包含37.4M条蛋白序列,服务器是24线程2.0GHz的CPU资源,耗时17min。

将Illumina测序数据比对到数据库:
$ diamond blastx --db nr_eukaryon_20190405 --query reads.fq.gz --out reads.tab

将全基因组蛋白序列比对到数据库:
$ diamond blastp --db nr_eukaryon_20190405 --query proteins.fasta --out nr.tab --outfmt 6 --sensitive --max-target-seqs 20 --evalue 1e-5 --id 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1
    对18700个蛋白序列进行Nr注释,服务器160线程2.4GHz的CPU资源,耗时91min(不加--sensitive参数只需要16min)。

此外,DIAMOND还包含其它几个子命令:

blastx
    将核酸序列比对到蛋白序列数据库
view
    从DAA文件生成格式化的输出结果
getseq
    从DIAMOND数据库文件获取FASTA序列
version
    打印DIAMOND软件的版本号
dbinfo
    打印数据库的信息
help
    打印DIAMOND软件的详细帮助信息

DIAMOND的参数:

数据库参数:
--in <string>    default: STDIN
    输入FASTA格式的蛋白序列数据库文件。
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--taxonmap <string>
--taxonnodes <string>
    这两个参数输入NCBI的物种分类信息。前者输入Nr数据库中蛋白编号和物种编号的对应关系信息( ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz);后者输入物种编号的层次信息(ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdmp.zip)。输入这两个文件有利于DIAMOND和Nr数据库进行比对时,可以选择和数据库中指定的物种大类的数据进行比对。

输入参数:
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--query | -q <string>    default: STDIN
    输入需要注释的FASTA或FASTQ格式的序列文件。可以是带.gz后缀的压缩文件。
--taxonlist <string>
    输入NCBI分类编号,仅对数据库中的目标子集进行比对。可以输入多个使用逗号分隔的编号ID。例如,古菌(2157)、细菌(2)、病毒(10239)、真菌(4751)、动物(33208)和植物(3193)。使用该参数,必须要使用--taxonmap和--taxonnodes参数构建数据库。
--query-gencode <int>    default:1
    设置在进行BLASTX模式进行比对时所使用的遗传密码。
--strand <string>    default: both
    设置对query序列的那条链进行比对。可以设置的值为:boh、plus和minus。
--min-orf | -l <int>
    在使用BLASTX模式进行比对时,若序列的某个ORF低于此值,则忽略该ORF。默认设置下:若核酸序列长度低于30,则值为1;若核酸序列长度低于100,则值为20;若核酸序列长度不低于100,则值为40。若该值设为1,则表示在BLASTX模式时,识别所有长度的ORF并进行比对。

比对参数:
--sensitive
    添加该参数,则能得到更多比对结果。该模式适合比对较长的序列。默认模式主要适用于比对short reads序列(Illumina reads),搜寻比对长度为30~40aa且bit得分大于50的匹配结果。
--more-sensitive
    相比于sensitive,能得到更全的比对结果。
--frameshifit | -F <int>
    在进行BLASTX模式时,设置对移码匹配的罚分,且一般推荐设置该参数值为15左右。设置该参数后,有利于比对较长的有INDEL错误的reads序列。默认情况是是禁用了此功能。
--gapopen <int>    default: BLOSUM62矩阵的默认值为1
    设置打开Gap罚分。
--gapextend <int>    default: BLOSUM62矩阵的默认值为1
    设置延长Gap罚分。
--matrix <string>    default: BLOSUM62
    设置计分矩阵。可以设置的值有:BLOSUM45,BLOSUM50,BLOSUM62,BLOSUM80,BLOSUM90,PAM250,PAM70和PAM30。
--comp-based-stats <bool>    default: 1
    设置是否对比对得分进行校正。
--algo <bool>
    设置种子序列搜索的算法。0表示使用double-indexed算法,1表示query-indexed算法。第一种算法是软件的主要算法,但不适合query数据较少的情况,这时需要使用第二种算法。默认设置下,软件会根据query和数据库的数据量自动选择相应的算法。两种算法的结果基本一致。
--freq-sd <int>
    设置种子序列最大的匹配频率。若种子序列在query或数据库种出现的频率过高,则忽略之。在一般模式下,该参数值为50;在sensitive模式下,该参数值为100,在more-sensitive模式下,该参数值为200。


输出参数:
--out | -o <string>    default:STDOUT
    设置输出文件。
--outfmt | -f <int>    default: 6
    设置输出格式。支持的格式有:0,两两比对格式;5,XML格式;6,BLAST表格格式;100,DIAMOND匹配存档(DDA)格式,该格式可以使用diamond view命令转换成其它格式;101,SAM格式;102,分类鉴定结果,该模式下结果文件分三列,QueryID、NCBI物种分类ID、最佳匹配evaule;103,PAF格式。
--salltitles
    在DDA格式结果中包含all subject titles。
--sallseqid
    在DDA格式结果中包含all subject ids。由于id比title要短,相比于--salltitles,本参数输出的结果文件更小。
--compress <bool>    default: 0
    设置是否对暑促文件进行压缩。0表示不压缩,1表示进行gzip压缩。
--max-target-seqs | -k <int>    default: 25
    设置每个query比对结果的最大匹配数目。
--top <int>
    若设置该参数的值,则报告一定数目的top匹配结果。例如设置该参数值为10,则报告匹配得分为top10的结果,即报告所有得分和最高得分相差10%以内的匹配结果。该参数会取代--max-target-seqs参数设置的值。
--range-culling
    添加该参数能剔除匹配到query序列相同位置的hit结果。该参数适合于较长的query DNA序列。当DNA序列较长时,它可能跨越多个基因,而一般情况下仅报告25个最佳匹配结果,若这25个hits结果可能evaule值很好但得分较低,且一起不能覆盖query序列的50%,会导致阵势的比对结果较差。因此,添加该参数能剔除匹配到query序列相同位置的hit结果,以使最终的比对结果对query序列的覆盖率超过50%。
--range-cover <float>    default: 50.0
    和参数--range-culling配合使用,表示剔除hits以使所有的比对结果能对query序列的覆盖率超过此参数设定的百分比。
--long-reads
    对于long read DNA测序数据,推荐添加该参数,等同参数--range-culling --top 10 -F 15。
--evalue | -e <float>    default: 0.001
    设置比对的evalue阈值。
--min-score <int>
    设置比对的bit score阈值。设置该参数会替代--evalue参数。
--id <int>
    设置identity阈值。
--query-cover <int>
    设置对query序列的覆盖度阈值。该阈值是HSP的阈值。
--subject-cover <int>
    设置对subject序列的覆盖度阈值。该阈值是HSP的阈值。
--max-hsps <int>
    设置报告的HSPs数量。默认是报告所有的HSPs,若设置成1,则仅报告最优HSP。
--unal <bool>
    设置是否报告没有匹配的比对结果。默认设置下,--outfmt为1、5和101时,该参数值为1,其它输出方式该参数值为0。
--no-self-hits
    添加该参数,则不报告序列和自己的比对结果。

性能参数:
--threads | -p <int>    default: Max
    设置程序运行所使用的CPU线程数。默认是服务器可用的最大CPU线程数。
--block-size | -b <float>    default: 2.0
    设置每次处理多少G的序列字符数。该参数控制程序消耗的内存量,一般内存消耗该值的6倍。设置更大的值,会消耗更多的内存和临时磁盘空间,但能提高性能。默认设置下,程序一次处理2G个序列字符,消耗内存12G。
--tmpdir | -t <string>    default: directory of --out
    设置临时文件夹路径。推荐该文件所在剩余磁盘空间有100G及以上。若将该参数设置为/dev/shm,则会将临时文件存放在内存,会增加内存消耗和计算性能。
--index-chunks | -c <int>    default: 4
    将seed index分成指定的份数。推荐将该参数值设置成1,能增加计算性能和内存使用量。

DIAMOND软件运行过程中,若使用参数 –block-size 20.0 –tmpdir /dev/shm –index-chunks 1 能使用较大内存,不会对数据进行分块,程序运行较快。程序运行过程中可能完全不会生成过程文件。但在程序运行过程中给出一些进度信息。当shape到15的时候,则表示程序快运行完毕。通过查看shape值来确定程序运行的进度。

CPU threads: 160
Scoring parameters: (Matrix=BLOSUM62 Lambda=0.267 K=0.041 Penalties=11/1)
Temporary directory: /dev/shm
Opening the database…
Target sequences to report alignments for: 20
Opening the input file…
Opening the output file…
Loading query sequences…
Masking queries…
Building query seed set…
Algorithm: Double-indexed
Building query histograms…
Allocating buffers…
Loading reference sequences…
Building reference histograms…
Allocating buffers…
Initializing temporary storage…
Processing query chunk 0, reference chunk 0, shape 0, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Processing query chunk 0, reference chunk 0, shape 1, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
...
Processing query chunk 0, reference chunk 0, shape 15, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Deallocating buffers…
Computing alignments…
Deallocating reference…
Loading reference sequences…
Deallocating buffers…
Deallocating queries…
Loading query sequences…
Closing the input file…
Closing the output file…
Closing the database file…
Deallocating taxonomy…
Total time = 659.515s
Reported 33888205 pairwise alignments, 33890340 HSPs.
122518 queries aligned.

DIAMOND默认设置下输出表格格式的结果。结果分12列,其结果信息和BLAST默认设置-outfmt 6输出的格式完全一致。

1. query序列ID
2. subject序列ID
3. Identity百分比
4. 匹配长度
5. 错配长度
6. 打开Gap的次数
7. query序列起始位点
8. query序列结束位点
9. subject序列起始位点
10. subject序列结束位点
11. E-vaule值
12. bitscore得分

CentOS7系统配置firewalld使用NAT方法将网络共享给内部局域网机器

学校实验室有台服务器申请了固定的公网IP,能连接外部网络,同时该机器和其它几台内部服务器连接在一个路由器上。需要将该服务器的网络共享给其它内网服务器。进行如下设置即可。

首先,外网服务器有两根网线连接,一根和外网相连,一根和内网路由器相连。在外网服务器进行NAT转发设置:

第一种方法:
开启NAT转发
# firewall-cmd --permanent --zone=public --add-masquerade
开放DNS使用的53端口,否则可能导致内网服务器虽然设置正确的DNS,但是依然无法进行域名解析。
# firewall-cmd --zone=public --add-port=53/tcp --permanent
重启防火墙
# systemctl restart firewalld.service
# 检查是否允许NAT转发
firewall-cmd --query-masquerade
# 关闭NAT转发
firewall-cmd --remove-masquerade

第二种方法:
开启ip_forward转发
在/etc/sysctl.conf配置文件尾部添加
net.ipv4.ip_forward=1
然后让其生效
# sysctl -p
执行firewalld命令进行转发:
firewall-cmd --permanent --direct --passthrough ipv4 -t nat -I POSTROUTING -o enoxxxxxx -j MASQUERADE -s 192.168.1.0/24
注意enoxxxxxx对应外网网口名称。
重启防火墙
# systemctl restart firewalld.service

设置外网服务器的IP地址为192.168.1.1,然后,在内网服务器中对配置文件/etc/sysconfig/network-scripts/ifcfg-enxxxxxx进行修改,设置网卡的网关为外网服务器的IP地址:

ONBOOT=yes
IPADDR=192.168.1.2
PREFIX=24
GATEWAY=192.168.1.1
DNS1=119.29.29.29

重启内网服务器的网络服务,即可访问外网了。

# systemctl restart network.service

创建NCBI的Nr数据库的子集数据库

1. 为什么要做Nr子集数据库

NCBI官网仅提供Nr全数据库。该数据库太大,将物种的蛋白序列使用Blastp比对到Nr数据库非常消耗计算和时间。对1个蛋白序列可能需要1个CPU计算半个小时。若对全基因组2万个基因分析,普通台式机8个CPU要计算2000/(2*8*24)=52天。这主要是由于Nr数据库太大导致的。为了能尽快得到Nr注释结果,可以按物种分类将Nr数据库分割成子集数据库,能得到更快的比对速度。以下是创建Nr子集数据库的步骤。

2. 创建Nr子集数据库的步骤

从NCBI下载Nr数据库和分类数据库文件

cd /opt/biosoft/Nr_database
# 下载Nr数据库(FASTA文件)
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /blast/db/FASTA/nr.gz ./

# 下载NCBI的分类数据库文件
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/taxdump.tar.gz ./
ascp -T -l 200M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --host=ftp.ncbi.nih.gov --user=anonftp --mode=recv /pub/taxonomy/accession2taxid/prot.accession2taxid.gz ./

# 解压缩两个NCBI的分类数据库文件
gzip -dc prot.accession2taxid.gz > prot.accession2taxid
mkdir ~/.taxonkit
tar zxf taxdump.tar.gz -C ~/.taxonkit
# 其主要有效文件有两个:
# names.dmp 记录物种名及其分类编号
# nodes.dmp 记录分类编号的节点信息
# 查看~/.taxonkit/names.dmp文件,使用关键词检索得到目标类的分类编号,例如:
# fungi 4751             # grep -P "\|\s+[fF]ungi\w*\s*\|" ~/.taxonkit/names.dmp
# plants 3193            # grep -P "\|\s+[pP]lant\w*\s*\|" ~/.taxonkit/names.dmp
# animals 33208          # grep -P "\|\s+[aA]nimal\w*\s*\|" ~/.taxonkit/names.dmp

下载并安装NCBI分类数据库解析软件TaxonKitTaxonKit,并解析nodes.dmp文件的物种节点信息,得到指定类的所有物种列表信息。

# 下载并安装NCBI分类数据库解析软件TaxonKit
wget https://github.com/shenwei356/taxonkit/releases/download/v0.2.4/taxonkit_linux_amd64.tar.gz
tar zxvf taxonkit_linux_amd64.tar.gz

# 提取古菌(2157)、细菌(2)和病毒(10239)这几个大类下的所有物种编号。
./taxonkit list -j 8 --ids 2,2157,10239 > sub.meta.list

# 再编写程序extract_sub_data_from_Nr.pl获得列表中物种在Nr数据库中的序列信息。
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.meta.list --acc2taxid prot.accession2taxid - > nr_meta.fasta

提取fungi/plants/animals子集

./taxonkit list -j 8 --ids 4751 > sub.fungi.list
./taxonkit list -j 8 --ids 3193 > sub.plants.list
./taxonkit list -j 8 --ids 33208 > sub.animals.list
./taxonkit list -j 8 --ids 10239 > sub.virus.list
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.fungi.list --acc2taxid prot.accession2taxid - > nr_fungi.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.plants.list --acc2taxid prot.accession2taxid - > nr_plants.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.animals.list --acc2taxid prot.accession2taxid - > nr_animals.fasta
gzip -dc nr.gz | perl extract_sub_data_from_Nr.pl --sub_taxon sub.virus.list --acc2taxid prot.accession2taxid - > nr_virus.fasta

使用makeblastdb创建blast本地数据库

makeblastdb -in nr_fungi.fasta -dbtype prot -title nr_fungi -parse_seqids -out nr_fungi_`date +%Y%m%d` -logfile nr_fungi_`date +%Y%m%d`.log
makeblastdb -in nr_plants.fasta -dbtype prot -title nr_plants -parse_seqids -out nr_plants_`date +%Y%m%d` -logfile nr_plants_`date +%Y%m%d`.log
makeblastdb -in nr_animals.fasta -dbtype prot -title nr_animals -parse_seqids -out nr_animals_`date +%Y%m%d` -logfile nr_animals_`date +%Y%m%d`.log
makeblastdb -in nr_virus.fasta -dbtype prot -title nr_virus -parse_seqids -out nr_virus_`date +%Y%m%d` -logfile nr_virus_`date +%Y%m%d`.log

CentOS7的服务开启方法和防火墙管理方法

CentOS6升级到CentOS7后,最明细的改变就是服务的启动方式和防火墙管理方式改变了。

1. CentOS的服务开启方法

在CentOS6中的服务开启方式(以httpd服务为例):

# /etc/init.d/httpd start      启动服务
# /etc/init.d/httpd stop       停止服务
# /etc/init.d/httpd restart    重启服务

# chkconfig httpd on           设置服务开机启动
# chkconfig httpd off          设置服务开机不启动
# chkconfig --list httpd       检测服务是否开机启动

在CentOS7中的服务器开启方式(以httpd服务为例):

# systemctl start httpd.service      启动服务
# systemctl stop httpd.service       停止服务
# systemctl restart httpd.service    重启服务

# systemctl enable httpd.service     设置服务开机启动
# systemctl disable httpd.service    设置服务开机不启动
# systemctl is-enabled httpd.service 检测服务是否开机启动

我常用的服务有:

防火墙服务: iptables firewalld.service
ssh远程联机:sshd     sshd.service    22端口
Apache网站: httpd    httpd.service   80端口
Mysql数据库:mysqld   mariadb.service 3306端口
网络服务:   network  network.service

CentOS7的服务都采用systemctl命令来进行操作。仅有network服务依旧可以和CentOS6一样进行操作。

2. CentOS的防火墙端口开放方法

在CentOS6中仅需要修改/etc/sysconfig/iptables配置文件内容,即可开放相应的端口:

-A INPUT -m state --state NEW -m tcp -p tcp --dport 22 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 80 -j ACCEPT
-A INPUT -m state --state NEW -m tcp -p tcp --dport 3306 -j ACCEPT

在CentOS7中需要使用firewall-cmd命令对端口进行管理:

查看已经开放的端口
# firewall-cmd --list-ports
添加端口
# firewall-cmd --add-port=80/tcp --permanent
加入--permanent参数,使永久生效。
重启防火墙服务
# systemctl restart firewalld.service
重启后,再查看端口,则生效了。
设置服务开机启动
# systemctl enable firewalld.service