利用perl一键生成符合LEFse差异分析的Table表

全栈程序员站长

发布于 2022-09-07 16:53:48

3650

发布于 2022-09-07 16:53:48

文章被收录于专栏：全栈程序员必看

大家好，又见面了，我是你们的朋友全栈君。

利用perl一键生成符合在线LEFse差异分析的Table表

LEfSe分析的在线+本地运行的详细教程参考刘尧博客

基于Picrust2进行宏基因预测后，我们往往需要对数据进行可视化话，其中LEFse就是非常不错的选择,这里通过perl实现对表的格式化。 LEFse –Galaxy平台：http://huttenhower.sph.harvard.edu/galaxy

use strict;
use warnings;

my $mapFile=$ARGV[0];
my $tableFile=$ARGV[1];

my %mapHash=();
open(MAP,"$mapFile")or die $!;
while(my $line=<MAP>)
{ 
   
	chomp($line);
	my @arr=split(/\s+/,$line);
	$mapHash{ 
   $arr[0]}=$arr[3];
}
close(MAP);

my %hash=();
my @samples=();
my %sumHash=();
open(RF,"$tableFile") or die $!;
while(my $line=<RF>)
{ 
   
	next if($.==1);
	my @arr=split(/\t/,$line);
	if($.==2)
	{ 
   
		@samples=@arr;
		next;
	}
	my $taxonomy=$arr[$#arr];
	$taxonomy=~s/;/|/g;
	$taxonomy=~s/(k|p|c|o|f|g|s)__//g;
	$taxonomy=~s/\s+|\[|\]//g;
	$taxonomy=~s/\|*$//g;
	my @taxoArr=split(/\|/,$taxonomy);
	my $taxSetp="";
	for(my $tax=0;$tax<=$#taxoArr;$tax++)
	{ 
   
		if($tax==0)
		{ 
   
			$taxSetp=$taxoArr[$tax];
		}
		else
		{ 
   
			$taxSetp=$taxSetp . "|" . $taxoArr[$tax];
		}
		for(my $i=1;$i<$#samples;$i++)
		{ 
   
			${ 
   $hash{ 
   $taxSetp}}{ 
   $samples[$i]}+=$arr[$i];
		}
	}
	for(my $i=1;$i<$#samples;$i++)
	{ 
   
		$sumHash{ 
   $samples[$i]}+=$arr[$i];
	}
}
close(RF);

open(WF,">lefse_input.txt")or die $!;
#print WF "Taxonomy";
#for(my $i=1;$i<$#samples;$i++)
#{ 
   
# print WF "\t" .$samples[$i];
#}
#print WF"\n";
print WF "sampleType";
for(my $i=1;$i<$#samples;$i++)
{ 
   
	print WF "\t" . $mapHash{ 
   $samples[$i]};
}
print WF "\n";

foreach my $key(sort (keys %hash))
{ 
   
	next if($key=~/Unassigned/);
	print WF $key;
	for(my $i=1;$i<$#samples;$i++)
	{ 
   
		print WF "\t" . ${ 
   $hash{ 
   $key}}{ 
   $samples[$i]}/$sumHash{ 
   $samples[$i]}*100;
	}
	print WF "\n";
}
close(WF);

将以上代码保存为lefse.pl, 然后cd到相应目录，将输入文件lefse_input.txt（以picrust2输出的meta表为例），和map文件 sample-meta.txt,放到该pl同一目录下，然后运行