基因家族分析中,在计算kaks时需要将每组串联重复基因序列提取出,然后分别进行多序列比对。
我们给大家提供的脚本只能一对一对的提取序列,那数量少还可以,可如果串联重复基因太多就有些麻烦了。这里提供一个批量提取串联重复基因序列的脚本,可以将每组串联重复基因序列提取到单独的文件中。
首先我们需要准备3个文件:
1.串联重复基因对文件。每一行是一对串联重复基因,中间用Tab键分割。如下图:
2.所有基因cds序列。
3.一个空的目录。
脚本运行命令如下:
perl get_fa_by_id_kaks.pl WRKY.tandem ../data/Sspon.v20180123.cds.fa test/
WRKY.tandem:串联重复基因对文件;
../data/Sspon.v20180123.cds.fa:所有基因cds序列文件;
test/:空目录。
get_fa_by_id_kaks.pl就是脚本了,其代码如下:
#email: huangls@biomics.com.cn
die "perl $0 <idlist> <fa> <OUT>" unless ( @ARGV == 3 );
use Math::BigFloat;
use Bio::SeqIO;
use Bio::Seq;
$in = Bio::SeqIO->new(
-file => "$ARGV[1]",
-format => 'Fasta'
);
my %gene;
while ( my $seq = $in->next_seq() ) {
my ( $id, $sequence, $desc ) = ( $seq->id, $seq->seq, $seq->desc );
$gene{$id} = $seq;
}
open IN, "$ARGV[0]" or die "$!";
my $n = 1;
while (<IN>) {
chomp;
next if /^#/;
my @a = split /\s+/;
my $out = Bio::SeqIO->new(
-file => ">$ARGV[2]/dup_gene_paired$n.fa",
-format => 'Fasta'
);
if ( exists $gene{$a[0]} ) {
my ( $id, $sequence, $desc ) = ( $gene{$a[0]}->id, $gene{$a[0]}->seq, $gene{$a[0]}->desc );
my $newSeqobj = Bio::Seq->new(-seq => $sequence,
-id =>$id,
);
$out->write_seq($newSeqobj);
}
if ( exists $gene{$a[1]} ) {
my ( $id, $sequence, $desc ) = ( $gene{$a[1]}->id, $gene{$a[1]}->seq, $gene{$a[1]}->desc );
my $newSeqobj = Bio::Seq->new(-seq => $sequence,
-desc => $desc,
-id =>$id,
);
$out->write_seq($newSeqobj);
}
$n++;
$out->close();
}
close(IN);
$in->close();
此外,我们在网易云课堂上有各种教学视频,有兴趣可以了解一下:
1. 文章越来越难发?是你没发现新思路,基因家族分析发2-4分文章简单快速,学习链接:基因家族分析实操课程
2. 转录组数据理解不深入?图表看不懂?点击链接学习深入解读数据结果文件,学习链接:转录组(有参)结果解读;转录组(无参)结果解读
3. 转录组数据深入挖掘技能-WGCNA,提升你的文章档次,学习链接:WGCNA-加权基因共表达网络分析
4. 转录组数据怎么挖掘?学习链接:转录组标准分析后的数据挖掘
6. 更多学习内容:linux、perl、R语言画图,更多免费课程请点击以下链接:
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!