利用primer3批量设计引物

1、Primer3的下载与安装

primer3官网上下载相应版本的程序,解压后安装,然后检测是否有安装错误。例如安装2.6.1版本:

wget https://github.com/primer3-org/primer3/archive/refs/tags/v2.6.1.tar.gz
tar xzvf primer3-2.6.1.tar.gz
cd primer3-2.6.1/src
make
make test

安装好之后只需要调用src文件夹中的primer3_core程序就可以了,可以拷贝到其它的地方使用。

primer3的各参数说明可以参考官方文档,中文的可以参考这篇博客

2、模板的准备

将模板序列单独放在一个文件,如果是fasta格式的可以命名为fasta格式的代码读取文件,如果是tsv格式的(第一列为序列id,第二列为模板序列),使用tsv格式对应的代码读取文件。然后用下面代码进行批量设计。如果对引物条件不敏感,可以在参数中加入这个参数PRIMER_PICK_ANYWAY=1“,如果想获得多对引物可以修改这个参数”PRIMER_NUM_RETURN=1“。其它参数也可以根据需要修改或者添加。

如果只是检测设计的引物是否满足条件,可以引入以下函数修改对应的参数。

#!/usr/bin/perl
use strict;
use warnings;
my $primer3_core = "./primer3_core";
my $Seq_file = "./input_file";
my $primer3out_file = "./primer3out.txt";
############Change the parameter if need #######
my $primer3input = <<SET;
PRIMER_TASK=generic
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=0
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_PRODUCT_SIZE_RANGE=75-100
PRIMER_EXPLAIN_FLAG=1
PRIMER_FIRST_BASE_INDEX=1
PRIMER_LIBERAL_BASE=1
SEQUENCE_FORCE_LEFT_END=41
PRIMER_NUM_RETURN=1
PRIMER_OPT_SIZE=20
PRIMER_MIN_SIZE=18
PRIMER_MAX_SIZE=25
PRIMER_MIN_GC=30.0
PRIMER_MAX_GC=70.0
PRIMER_MIN_TM=58.0
PRIMER_OPT_TM=60.0
PRIMER_MAX_TM=62.0
=
SET
open(F, $Seq_file) or die "Can not open the flankingSeq file $!\n";
open(POUT, ">$primer3out_file");
my $print_out_head = join("\t", ("primer_id", "Foward", "Reverse"));
print POUT $print_out_head."\n";
my %id2seq;
my $id="";
my $seq="";
############if the input file is the fasta format##################
while(<F>){
	if(/^>(\S+)/){
		if($id){
			$id2seq{$id}=$seq;
			$seq="";
		}
		$id=$1;
	}else{
		s/\s+$//;
		$seq.=$_;
	}
	if(eof){
		$id2seq{$id}=$seq;
	}
}
############if the input file is the tsv format, use this code################
#while(<F>){
#	s/\s+$//;
#	($id, $seq)=split/\t/;
#	$id2seq{$id}=$seq;
#}		
		
foreach my $primer_id(keys %id2seq){
	my (@primer_left, @primer_right)=&primerPick($primer_id, $id2seq{$primer_id});
	my $i=0;
	while($i<=$#primer_left){
		print POUT $primer_id."\t".$primer_left[$i]."\t". $primer_right[$i]."\n" if($primer_left[$i] and $primer_right[$i]);
	}
}
	
close(POUT);
system("rm ./primer3input.txt");
sub primerPick{
	my ($snp_id, $template)=@_;
	my $attachmentSet = <<ASET;
SEQUENCE_ID=$snp_id
SEQUENCE_TEMPLATE=$template
ASET
	my $primer3set=$attachmentSet.$primer3input;
	open(OUT, ">./primer3input.txt");
	print OUT $primer3set;
	close(OUT);	
	
	my $primer3output = `$primer3_core --default_version=2 --io_version=4 --strict_tags  <./primer3input.txt `;
	my @primer_left="";
	if($primer3output=~/PRIMER_LEFT_(\d+)_SEQUENCE=(\w+)/g){
		$primer_left[$1]=$2;
	}
	my @primer_right="";
	if($primer3output=~/PRIMER_RIGHT_(\d+)_SEQUENCE=(\w+)/g){
		$primer_right[$1]=$2;
	}
	
	return @primer_left,@primer_right;
}
sub checkPrimer{
	my ($primer_id, $primer_left, $primer_right, $template)=@_;
	my $primer3set=$primer3input;
	$primer3set=~s/PRIMER_TASK=generic/PRIMER_TASK=check_primers/;
	$primer3set=~s/PRIMER_PICK_LEFT_PRIMER=1/PRIMER_PICK_LEFT_PRIMER=0/;
	$primer3set=~s/PRIMER_PICK_RIGHT_PRIMER=1/PRIMER_PICK_RIGHT_PRIMER=0/;
	my $attachmentSet = <<ASET;
SEQUENCE_ID=$primer_id
SEQUENCE_TEMPLATE=$template
SEQUENCE_PRIMER=$primer_left
SEQUENCE_PRIMER_REVCOMP=$primer_right
ASET
	$primer3set = $attachmentSet.$primer3set;
	open(OUT, ">./primer3input.txt");
	print OUT $primer3set;
	close(OUT);
	my $primer3out = `$primer3_core --default_version=2 --io_version=4 --strict_tags  <./primer3input.txt `;
	if($primer3out=~/PRIMER_PAIR_EXPLAIN=considered 1, ok 1/){
		return 1;
	}else{
		return 0;
	}
}

perl删除字符串中的重复字符

在设计KASP引物时,需要知道双等位基因型对应的碱基,这就需要从GT中却除重复的基因型,可以参考下面的代码:

#!/usr/bin/perl
use strict;
use warnings;
my $string="A/C	A/C	A/C	A/C	A/C	A/A	A/C	A/A	A/C	A/C	A/C	A/C	A/A	A/C	A/C	A/A";
$string=~s/\s+//g;
while($string =~ s/((.).*)\2+/$1/g) {};
print $string."\n";

perl循环读取文件时判断是否到最后一行

当需要去掉序列的分行符或者是建立序列的索引时,当序列文件较大时,不好将文件读入内存,需要一行一行的去读取,但怎样判断到了文件最后一行呢,在perl中一般用eof函数,但 eof 和带空圆括弧对 () 的 eof()表示的意思是不一样的,特别容易混淆。带圆括弧的 eof() 只是检测一组文件中的最后一个文件的文件结束,而 eof(没有圆括弧)在 while (<>) 循环里检查每个文件的文件结束。下面是建立fasta文件索引的代码。

#!/usr/bin/perl
use strict;
use warnings;
my $genome_file = “genomic.fna";
open(F, $genome_file);
#open(OUT, ">outfile");
my %chr2seq;
my $id="";
my $seq="";
while(<F>){
	if(/^>(\S+)/){
		if($id){
			$chr2seq{$id}=$seq;
			$seq="";
		}
		$id=$1;
	}else{
		s/\s+$//;
		$seq.=$_;
	}
	if(eof){
		$chr2seq{$id}=$seq;
	}
}

perl批量导入手机号码和邮箱

当单位或者公司发布通讯录时,会增加很多新的联系人,一个一个的导入实在是太麻烦了,所以打算自己写个代码生成个联系人文件,然后导入手机中就可以了,费话不多说,开干。

首先 将联系人信息用EXCEL整理成三列 ,第一列为姓名,第二列为手机号码,第三列为邮箱,其它也是可以加一些分组信息和地址之类的信息,但考虑信息越详细,信息泄露的后果就越严重,因此只选择添加手机号码和邮箱。整理好后复制到一个”input.txt”的文本文件,可以用记事本创建。运行下面perl程序。

将生成的“contacts.vcf”文件传到手机,然后用“联系人”或“电话本”打开这个文件,就会提示导入联系人,成功导入后,选择手机中的整理联系人,合并名字重复的联系人就可以了。

#!/usr/bin/perl
##用perl处理Excel通讯录(从Excel复制到input.txt),最后形成vCard的
###############################################################
#input.txt文件格式如下:
#第一列姓名,第二列手机号码,第三email,中间制表符分开
###############################################################
#比如输入文件内容如下:
#张三 手机号码 email
#李四 手机号码 email
###############################################################
#输出内容如下:
#BEGIN:VCARD
#VERSION:2.1
#N;CHARSET=UTF-8;ENCODING=QUOTED-PRINTABLE:;=E6=9D=8E=E5=9B=9B;;;
#FN;CHARSET=UTF-8;ENCODING=QUOTED-PRINTABLE:=E6=9D=8E=E5=9B=9B
#TEL;CELL:1234567
#EMAIL;X-internet:yuyin110110110@163.com
#END:VCARD
###############################################################
##perl convert2vcf.pl 
##Eidit by XiaowenH 2020-05-20
###############################################################
use strict;
use warnings;

#先打开读取文件内容,保存到数组里面
open(FH, "input.txt");#读取通讯录
my @data=<FH>;
close(FH);
#输出结果
open(OUT,">contacts.vcf"); #输出通讯录
foreach my $line (@data){
    $line=~s/\s+$//;
    my @l=split/\t/,$line;
    my $name = $l[0];
    $name=~s/^\s+//g;  ###去掉名字开头和结尾空格
    $name=~s/\s+$//g;

    my $phone="";
    if(exists $l[1]){
        $l[1]=~s/^\s+//g; ###去掉电话号码开头和结尾空格
        $l[1]=~s/\s+$//g;
	$phone=$l[1];
    }

	
    my $email = "";
    if(exists $l[2]){
	$email = $l[2];
    }
    print OUT "BEGIN:VCARD\n";
    print OUT "VERSION:2.1\n";
    print OUT "N;CHARSET=UTF-8;ENCODING=QUOTED-PRINTABLE:;$name;;;\n";
    print OUT "FN;CHARSET=UTF-8;ENCODING=QUOTED-PRINTABLE:$name\n";  #姓名
    print OUT "COUNTRYISO:CN\n";
    print OUT "TEL;CELL:$phone\n";#手机号
    print OUT "EMAIL;X-internet:$email\n";#邮箱
    print OUT "END:VCARD\n";
}
close(OUT);

Bioperl的简单快速安装

Bioperl 由于涉及到的包太多,需要超过12000的单独的检测,安装起来比较麻烦,而且很多检测都通不过,对于新手过于繁杂,现在bioperl官网介绍了一种简单快速安装的方法————用第三方工具进行安装,经亲自测试,效果非常好,现简要介绍而下。
1 安装perlbrew(https://perlbrew.pl/)
在终端输入

curl -L https://install.perlbrew.pl | bash

如果系统没有curl命令
对于linux系统

wget -O - https://install.perlbrew.pl | bash

FreeBSD系统

fetch -o- https://install.perlbrew.pl | sh

也可以用cpan安装,在终端输入(注意不需要先打开cpan)

sudo cpan App::perlbrew
perlbrew init

如果是用cpan安装,perlbrew程序安装在/usr/bin/perlbrew 或者 /usr/local/bin/perlbrew,默认的perlbrew的根目录是~/perl5/perlbrew 。

2 安装cpanm

perlbrew install-cpanm

3 安装Bio::perl

cpanm Bio::Perl

如果找不到cpanm,就输入cpanm的绝对路径

/root/perl5/perlbrew/bin/cpanm Bio::Perl

如果其中报错缺少哪个模块或者哪个模块版本太老,就先安装报错的模块,我的显示是 Test::More 模块版本太老,我就用cpan “install Test::More”就好了。如果正确安装会显示

–> Working on Bio::Perl
Fetching http://www.cpan.org/authors/id/C/CJ/CJFIELDS/BioPerl-1.007001.tar.gz … OK
Configuring BioPerl-1.007001 … OK
Building and testing BioPerl-1.007001 … OK
Successfully installed BioPerl-1.007001
1 distribution installed

一般在第一步“Working on Bio::Perl”会等一段时间,因为是要下载安装包,只要没有报错,等着就可以了。

基因注释

很多人觉得基因注释很高深,特别对于那些不太懂编程的来说,感到无从下手,去年我也折腾了差不多半年,最后发现也就那么回事。一开始我是找各种基因注释的工具,例如研究geneontology自带的perl程序,里面涉及了众多模板,而且是模块里面又有模块,运行时不是少了这个模块就是那个模块有问题,最后折腾了两三个月也没有得到理想结果。后来用了blast2go这个程序用起来是没有问题,可实在是太慢,几个小时下去,泡都不冒一个。最后默默在躺在床上反思,我的目的是什么,注释的原理是什么。按照这个思路,慢慢的去找这些问题的答案,最后终于开窍了,想明白了真的好简单,现将思路简单归纳如下:
1、基因注释方法
1)根据已注释的信息,利用序列相似性原则,去注释未知序列。
2)查找文献资料注释基因功能。
基因批量注释的话都是采用第一种。
2、所需资料
1)将需要注释的序列翻译成氨基酸序列;
2)下载现有的已注释的蛋白序列(uniprot_sprot.fasta.gz,文件大小为78MB);
3)下载联系蛋白序号和注释号(GO的accession)的文件(idmapping.td.gz,文件大小为2.17G);
4)下载最新的GO数据库文件assocdb-data.gz
3、注释方法
1)利用formatdb格式化已注释的蛋白序列为blast库;
2)将需要注释的序列与blast库做blast;
3)找到与需注释序列相似度高的蛋白序列号;
4)找到这些蛋白序列号对应的注释号(GO号);
5)用GO号在GO数据库文件中的term表中找到注释信息或者在geneontology官网用GO号就可以查到详细的注释信息;

如果用NCBI数据库,方法也差不多,先下载已注释的蛋白序列,格式化数据成blast库,将需注释的序列与这个blast库比对,找到相似度高的已注释的蛋白序列号,用此序列号用NCBI中gene_info.gz和gene2access找到蛋白序列的gene号,然后将此号用gene2go数据找到go号,最后用go号找到其注释信息就行了。

如果按照以上方法操作的话,其实基因注释非常简单,几乎只要做个blast,然后将信息对应过来就可以了,远没有那么复杂。

Perl输出函数传递的二维数组

最近用DBI模块写个了mysql接入函数,但得到的return一直不能输出结果,后来经过思考,原来是二维数组传递的问题,终于解决了,特此为记。

下面这个是mysql接入函数:

sub mysql_connect(){
        use DBI;
	my ($sql,$db)=@_;
	my $usr="selectonly";
	my $passwd="yuyin_110";
	my $dbh=DBI->connect("DBI:mysql:$db",$usr,$passwd);

	if(!$dbh){
		print "connect lost\n";
	}else{
		my $sth=$dbh->prepare($sql)or die "Can’t prepare $sql: $dbh->errstr\n";
		$sth->execute()or die "Can’t execute the query: $sth->errstr";
		my $data=$sth->fetchall_arrayref();
		$sth->finish;
		return $data;
	}
}

现在的问题是输出$data的内容,而$data返回的是一个二维数组的引用地址,可以用以下方法分层输出这个数组:

	my $sql="select * from test limit 1000";
	my $b=mysql_connect($sql,"test");
	my @r=@$b;    #####因为$b获得的只是一个引用地址,可以用@变成数组,这个数组这储存了一系列一维数组地址
	foreach(@r){
		my @each=@$_;            ######道理与获得二维地址一样
		foreach(@each){
			print  $_."\t";
		}
		print  "\n";

	}

Bioperl的简单安装

按照Bioperl上介绍的方法在linux下安装Bioperl老是安不上,或者是安装上了,但不能用,上面介绍的几种方法都试了,全不行,后面自己想了个办法,就是利用cpan只对要用到的模块进行单独安装,简单适应,如果你和我碰到了同样的问题不妨试试。
0、用root用户登录,不能由于权限问题,安装不上去。
1、确定cpan能用。
>perl -MCPAN -e shell
cpan>install Bundle::CPAN
cpan>q
2、升级cpan,保证安装的模块是最新的。
>cpan
cpan>install Module::Build
cpan>o conf prefer_installer MB
cpan>o conf commit
cpan>q
3、安装Bioperl最重要的模块SeqIO(该模块可以实现文件格式转换,计算序列长度,blast信息提取等),中间会有些选项要求选择,一路回车采用默认的就行了。
cpan>install Bio::SeqIO
4、安装SeqFeature模块(序列特征信息的获取或解析)。
cpan>install Bio::SeqFeature
5、安装GenBank模块
cpan>install Bio::GenBank
6、安装AlignIO和AlignI模块(数据格式格式转换)。
cpan>install Bio::AlignIO
cpan>install Bio::AlignI
7、安装DNAstatistics模块(序列统计分析,进化距离计算)。
cpan>install Bio::DNAstatistics

上面是一引起常用的模块,至于其它的一些模块如果要用到就按这种方法安装就行了。
利用cpan安装的模块默认路径是当前用户的perl5/lib/perl5目录下面,如/root/perl5/lib/perl5/,如果运行perl -e “use Bio::SeqIO”还是找不到该模块的话,可以将 /当前用户/perl5/lib/perl5/下的文件拷贝到@INC目录中,如/usr/local/lib64/perl5/中,最好也将/当前用户/perl5/bin/中的文件拷贝到/usr/local/bin/perl5中,以后可在那里调用。