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

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

应用场景

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

输入文件

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

输入文件长酱紫

写脚本

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

#!/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” 即为输出文件。

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
本文目录
本文目录