Shiny从入门到入定——16-Escaping the graph

承接上一篇文章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!