构建样本vs基因矩阵
发表于:2021-09-29 | 分类: 生物信息
字数统计: 655 | 阅读时长: 2分钟 | 阅读量:

承接上一篇文章Swissprot数据库的本地化与序列比对

应用场景

分别预测了多个样本/基因组中某些基因的存在与否即数量,需要将这些样本/基因组中的基因数量情况合并在一起构建矩阵,此时,手动是非常困难和无趣的。又该请出Perl神了。

输入文件

在上一篇文章中,我们将多个宏基因组蛋白序列与SwissProt数据库做了比对,并根据比对到的ID与其他数据库做了mapping,得到了多个输出文件,保存为“sample.mapped”。其中“sample”可以是样本名,也可以是基因组名,它将出现在最后构建的矩阵中。这些文件既可作为本例的输入文件,其内容大概是下面酱紫的,各列以制表符分隔。

输入文件长酱紫

写脚本

上图是其中一个文件的内容的一部分,接下来我们将提取第19列的GO number来构建矩阵。将以下代码保存到文件中,命名为“get_matrix.pl”。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#!/usr/bin/perl
use strict;
use warnings;
# Author: Liu hualin
# Date: Sep 29, 2021

my %hash;
my %ids;
my %samples;

my @files = glob("*.mapped");
foreach (@files) {
$_=~/(.+).mapped/;
my $sample = $1;
$samples{$1}++;
open IN, "$_" || die;
<IN>;# 忽略第一行,如果第一行不是标题行,请将该行注释掉
while (<IN>) {
chomp;
my @lines = split /\t/;
if ($lines[18]=~/.+\; /) {
my @terms = split /\; /, $lines[18];# 18代表文件的第19列,若想提取其他列,可以自行修改该数字为“列号-1”,因为第一列代号为0
foreach (@terms) {
$ids{$_}++;
$hash{$sample}{$_}++;
}
}elsif ($lines[18]=~/\S+/) {
$ids{$lines[18]}++;
$hash{$sample}{$lines[18]}++;
}
}
}
open OUT, ">Matrix.txt" || die;

my @samples = sort keys %samples;
my @ids = sort keys %ids;

print OUT "\t" . join("\t", @samples) . "\n";
for (my $i=0; $i<@ids ;$i++) {
print OUT $ids[$i];
for (my $j=0; $j<@samples ;$j++) {
if (exists $hash{$samples[$j]}{$ids[$i]}) {
print OUT "\t$hash{$samples[$j]}{$ids[$i]}";
}else {
print OUT "\t0";
}
}
print OUT "\n";
}
close OUT;

构建矩阵

将脚本与输入文件放在同一目录下,在终端或Windows命令行中运行如下命令,得到的“Matrix.txt”即为输出文件。

1
perl get_matrix.pl

输出文件长酱紫

脚本获取

关注公众号“生信之巅”,聊天窗口回复“412a”获取下载链接。

生信之巅微信公众号 生信之巅小程序码

敬告:使用文中脚本请引用本文网址,请尊重本人的劳动成果,谢谢!Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

上一篇:
CAZy碳水化合物活性酶预测
下一篇:
Swissprot数据库的本地化与序列比对并与其他数据库快速mapping