我有一个小程序,基本上可以处理爆炸命中列表,并通过迭代包含每个爆炸列表的哈希结果(作为哈希键)来检查爆炸结果之间是否存在重叠。
这涉及以相同方式将每个爆炸输入文件处理为$ ARGV。根据我要达到的目标,我可能想比较2、3或4个blast列表中的基因重叠。我想知道如何将基本处理块作为子例程编写,无论存在多少$ ARGV参数,都可以对其进行迭代。
例如,如果我输入2个爆炸列表,则下面的命令会正常工作:
#!/usr/bin/perl -w
use strict;
use File::Slurp;
use Data::Dumper;
$Data::Dumper::Sortkeys = 1;
if ($#ARGV != 1){
die "Usage: intersect.pl <de gene list 1><de gene list 2>\n"
}
my $input1 = $ARGV[0];
open my $blast1, '<', $input1 or die $!;
my $results1 = 0;
my (@blast1ID, @blast1_info, @split);
while (<$blast1>) {
chomp;
@split = split('\t');
push @blast1_info, $split[0];
push @blast1ID, $split[2];
$results1++;
}
print "$results1 blast hits in $input1\n";
my %blast1;
push @{$blast1{$blast1ID[$_]} }, [ $blast1_info[$_] ] for 0 .. $#blast1ID;
#print Dumper (\%blast1);
my $input2 = $ARGV[1];
open my $blast2, '<', $input2 or die $!;
my $results2 = 0;
my (@blast2ID, @blast2_info);
while (<$blast2>) {
chomp;
@split = split('\t');
push @blast2_info, $split[0];
push @blast2ID, $split[2];
$results2++;
}
my %blast2;
push @{$blast2{$blast2ID[$_]} }, [ $blast2_info[$_] ] for 0 .. $#blast2ID;
#print Dumper (\%blast2);
print "$results2 blast hits in $input2\n";
但我希望能够对其进行调整,以适应3或4个爆炸列表的输入。我想象一个子例程对此将是最好的,它会为每个输入调用,并且可能看起来像这样:
sub process {
my $input$i = $ARGV[$i-1];
open my $blast$i, '<', $input[$i] or die $!;
my $results$i = 0;
my (@blast$iID, @blast$i_info, @split);
while (<$blast$i>) {
chomp;
@split = split('\t');
push @blast$i_info, $split[0];
push @blast$iID, $split[2];
$results$i++;
}
print "$results$i blast hits in $input$i\n";
print Dumper (\@blast$i_info);
print Dumper (\@blast$iID);
# Call sub 'process for every ARGV...
&process for 0 .. $#ARGV;
更新:
我已删除了最后一个代码片段的哈希部分。
结果数据结构应为:
4次爆炸 <$input$i>
$VAR1 = [
'TCONS_00001332(XLOC_000827),_4.60257:9.53943,_Change:1.05146,_p:0.03605,_q:0.998852',
'TCONS_00001348(XLOC_000833),_0.569771:6.50403,_Change:3.51288,_p:0.0331,_q:0.998852',
'TCONS_00001355(XLOC_000837),_10.8634:24.3785,_Change:1.16613,_p:0.001,_q:0.998852',
'TCONS_00002204(XLOC_001374),_0.316322:5.32111,_Change:4.07226,_p:0.00485,_q:0.998852',
];
$VAR1 = [
'gi|50418055|gb|BC078036.1|_Xenopus_laevis_cDNA_clone_MGC:82763_IMAGE:5156829,_complete_cds',
'gi|283799550|emb|FN550108.1|_Xenopus_(Silurana)_tropicalis_mRNA_for_alpha-2,3-sialyltransferase_ST3Gal_V_(st3gal5_gene)',
'gi|147903202|ref|NM_001097651.1|_Xenopus_laevis_forkhead_box_I4,_gene_1_(foxi4.1),_mRNA',
'gi|2598062|emb|AJ001730.1|_Xenopus_laevis_mRNA_for_Xsox17-alpha_protein',
];
和输入:
TCONS_00001332(XLOC_000827),_4.60257:9.53943,_Change:1.05146,_p:0.03605,_q:0.998852 0.0 gi|50418055|gb|BC078036.1|_Xenopus_laevis_cDNA_clone_MGC:82763_IMAGE:5156829,_complete_cds
TCONS_00001348(XLOC_000833),_0.569771:6.50403,_Change:3.51288,_p:0.0331,_q:0.998852 0.0 gi|283799550|emb|FN550108.1|_Xenopus_(Silurana)_tropicalis_mRNA_for_alpha-2,3-sialyltransferase_ST3Gal_V_(st3gal5_gene)
TCONS_00001355(XLOC_000837),_10.8634:24.3785,_Change:1.16613,_p:0.001,_q:0.998852 0.0 gi|147903202|ref|NM_001097651.1|_Xenopus_laevis_forkhead_box_I4,_gene_1_(foxi4.1),_mRNA
TCONS_00002204(XLOC_001374),_0.316322:5.32111,_Change:4.07226,_p:0.00485,_q:0.998852 0.0 gi|2598062|emb|AJ001730.1|_Xenopus_laevis_mRNA_for_Xsox17-alpha_protein
您不能在变量名的中间插入变量值。(好吧,您可以,但是不应该。即使那样,您也不能在名称中间使用数组索引。)
这些名称无效:
@blast[$i]_info
@blast[$i]_ID
您需要将索引移到末尾:
@blast_info[$i]
@blast_ID[$i]
就是说,我会完全摆脱数组,而使用哈希代替。
您的第二个代码段未显示对子例程的调用。除非明确调用它,否则它将永远不会运行,并且您的程序将什么也不做。我将修改process
sub以采用单个参数,并为的每个元素调用它@ARGV
。例如
process($_) foreach @ARGV;
这是我编写程序的方式:
use strict;
use warnings;
use Data::Dumper;
my @blast;
push @blast, process($_) foreach @ARGV;
print Dumper(\@blast);
sub process {
my $file = shift;
open my $fh, '<', $file or die "Can't read file '$file' [$!]\n";
my %data;
while (<$fh>) {
chomp;
my ($id, undef, $info) = split '\t';
$data{$id} = $info;
}
return \%data;
}
尚不清楚最终的数据结构应该是什么样。(我做出了最大的猜测。)我建议阅读perlreftut以获得对引用的更好的基本了解,并使用它们在Perl中构建数据结构。
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句