写了一个perl脚本，调用了两个子程序，这两个子程序中分别都要遍历一个文件。
刚开始将子程序需要重复遍历的文件句柄没有放在子程序中重复打开，而是放在了脚本开始；如下：
open DMR,"$ARGV[0]";
open MATC,"$ARGV[1]";
open MATS,"$ARGV[2]";
while(<DMR>){
...
@a=&sub1(a,b,c);
@b=&sub2(a,b,c);
...
}
sub sub1{
...
while(<MATC>){...}
return ...
}
sub sub2{
...
while(<MATC>){...}
return ...
}
close DMR;
close MATC;
close MATS;
每次脚本运行的结果都显示子程序调用了一次，第二次调用的时候都无法返回想要的值；
折腾了一上午，各种测试，各种修改，各种百度，没有找到原因；
最后测试到子程序才恍然大悟，第一次调用子程序后,文件句柄MATC和MATS已经被遍历完了，无法再次遍历，导致
子程序调用了一次后就无法再次返回有效值了，哭死....，菜鸟级选手，发现的太晚了。
解决办法就是将文件句柄放在子程序中进行遍历；
但是作为新手，每次读取一行文件，子程序都要重复打开文件句柄然后关闭，然后下次继续打开关闭，太不效率了；
后来又发现了seek（）函数
seek()函数是通过文件句柄来移动文件读写指针的方式来读取或写入文件的，以字节为单位进行读取和写入：
seek FILEHANDLE, POSITION, WHENCE
参数说明:
FILEHANDLE：文件句柄，用于存放一个文件唯一标识符。
POSITION：表示文件句柄(读写位置指针)要移动的字节数。
WHENCE：表示文件句柄(读写位置指针)开始移动时的起始位置，可以取的值为0、1、2；分别表示文件开头、当前位置和文件尾。
在我们的子程序中遍历文件，我们可以用seek（FILEHANDLE, 0, 0）来代替不断的open 和close 文件句柄了。
这样在脚本开头open，脚本结束close就可以了。
原始文件（<DMR>）
$zxg head ./DMR/Stress0_VS_Control0.DMR.txt
Chr Start End Length N_C Stress0 Control0 Diff AreaStat Direction
chr1 3998120 3998220 100 8 0.8635 0.5784 0.2851 55.3665 Hyper
chr1 6916517 6917222 705 69 0.9187 0.5976 0.3211 695.3536 Hyper
chr1 9454505 9454633 128 8 0.5281 0.2002 0.3280 52.8468 Hyper
chr1 10970705 10970870 165 26 0.8999 0.5681 0.3317 307.8769 Hyper
chr1 15955916 15956232 316 28 0.9360 0.6386 0.2973 267.6975 Hyper
chr1 17793102 17793290 188 9 0.8104 0.4703 0.3402 78.3708 Hyper
chr1 19440293 19441256 963 42 0.9334 0.4415 0.4919 618.3510 Hyper
chr1 20705333 20705755 422 38 0.9415 0.2969 0.6446 886.2865 Hyper
chr1 20955313 20955493 180 15 0.3614 0.6616 -0.3002 -110.5831 Hypo
control组的matrix 9个样品
$zxg head matrix_control.txt
chr pos C03 C04 C05 C112 C13 C19 C210 C214 C21
chr1 100000060 1.000,13,0 0.909,10,1 0.909,10,1 0.889,8,1 1.000,17,0 0.333,1,2 1.000,13,0 0.929,13,1 1.000,8,0
chr1 100000061 0.947,18,1 0.786,11,3 0.917,11,1 0.923,12,1 0.944,17,1 0.812,13,3 0.905,19,2 0.947,18,1 1.000,17,0
chr1 100000094 0.750,9,3 0.833,10,2 0.917,11,1 0.875,7,1 0.933,14,1 0.800,4,1 0.900,9,1 0.857,12,2 1.000,10,0
chr1 100000095 0.952,20,1 0.765,13,4 0.786,11,3 0.765,13,4 1.000,16,0 0.882,15,2 0.929,26,2 0.900,18,2 0.765,13,4
chr1 100000185 0.957,22,1 1.000,16,0 1.000,17,0 1.000,13,0 0.846,11,2 1.000,10,0 1.000,14,0 1.000,15,0 1.000,12,0
chr1 100000186 1.000,24,0 1.000,23,0 1.000,21,0 1.000,21,0 1.000,24,0 0.944,17,1 0.920,23,2 1.000,12,0 0.938,15,1
chr1 100000190 1.000,23,0 0.938,15,1 0.882,15,2 1.000,13,0 0.917,11,1 0.917,11,1 0.929,13,1 1.000,12,0 1.000,12,0
chr1 100000191 1.000,24,0 1.000,21,0 0.952,20,1 1.000,21,0 1.000,25,0 0.941,16,1 1.000,24,0 0.917,11,1 1.000,17,0
chr1 100000477 1.000,16,0 1.000,10,0 0.900,9,1 1.000,10,0 0.864,19,3 1.000,18,0 0.933,14,1 1.000,10,0 1.000,14,0
stressl组的matrix 9个样品
$zxg head matrix_stress.txt
chr pos H202 H203 H209 H2110 H2111 H217 H224 H228 H229
chr1 100000060 1.000,15,0 0.909,10,1 1.000,9,0 0.917,11,1 1.000,5,0 0.833,5,1 0.875,14,2 1.000,5,0 0.727,8,3
chr1 100000061 0.818,18,4 0.826,19,4 0.931,27,2 0.700,7,3 0.846,11,2 0.882,15,2 0.917,11,1 1.000,15,0 1.000,12,0
chr1 100000094 0.778,14,4 1.000,11,0 0.800,8,2 0.909,10,1 1.000,5,0 1.000,7,0 0.941,16,1 1.000,8,0 0.900,9,1
chr1 100000095 0.875,21,3 0.857,18,3 0.939,31,2 0.727,8,3 1.000,10,0 0.842,16,3 0.929,13,1 0.875,14,2 1.000,10,0
chr1 100000185 0.955,21,1 0.909,10,1 1.000,12,0 0.889,8,1 1.000,9,0 1.000,7,0 0.867,13,2 1.000,9,0 0.875,7,1
chr1 100000186 1.000,29,0 0.960,24,1 1.000,33,0 0.818,9,2 1.000,10,0 1.000,12,0 0.938,15,1 1.000,17,0 0.909,10,1
chr1 100000190 1.000,21,0 0.917,11,1 1.000,12,0 0.900,9,1 1.000,9,0 1.000,8,0 1.000,14,0 1.000,9,0 0.875,7,1
chr1 100000191 1.000,29,0 1.000,25,0 0.970,32,1 0.917,11,1 1.000,10,0 1.000,12,0 1.000,16,0 0.941,16,1 1.000,10,0
chr1 100000477 1.000,12,0 1.000,19,0 1.000,19,0 1.000,18,0 1.000,10,0 0.941,16,1 1.000,15,0 1.000,12,0 0.909,10,1
程序前边的变量声明
#! /usr/bin/perl -w
use strict;
my $head;
my %hash_c;
my %hash_s;
open DMR,"$ARGV[0]";
$head=<DMR>;
chomp $head;
主程序
while (<DMR>){
chomp; my @aa=split/\t/,$_;
my @cc=&methyl_level_c($aa[0],$aa[1],$aa[2]);#向子程序1传递三个参数chr，start，end
my @dd=&methyl_level_s($aa[0],$aa[1],$aa[2]);#向子程序2传递三个参数chr，start，end
$hash_c{$aa[0]}{"$aa[1]"."-"."$aa[2]"}=[$aa[3],$aa[4],$aa[5],$aa[6],$aa[7],$aa[8],$aa[9],$cc[0],$cc[1],$cc[2]];
$hash_s{$aa[0]}{"$aa[1]"."-"."$aa[2]"}=[$dd[0],$dd[1],$dd[2]];
}
这个脚本的目的是给子程序传递染色体号，DMR起始和中止位置，然后获得这个DMR区间内每个样品的甲基化水平。
每个样品的甲基化水平存放在两个matrix中，分别为control和stress.
子程序1
sub methyl_level_c{
my @a;
my @b;
my $f0;
my $f1;
my $f2;
my $f0_cg=0;
my $f1_cg=0;
my $f2_cg=0;
my $f0_level=0;
my $f1_level=0;
my $f2_level=0;
open MATC,"$ARGV[1]";#打开control_matrix.txt
<MATC>;#去表头
while(<MATC>){
chomp; @a=split/\t/,$_;
if($_[0]=~/$a[0]/ && $a[1]>=$_[1] && $a[1]<=$_[2]){ #匹配染色体号，界定DMR起始和终止位置的CG位点
foreach my $pos(2..10){
@b=split/,/,$a[$pos];#依次分割每个样品，取第一个值为该样品的该CG位点甲基化水平
if($pos<=4){if($b[0] ne "-"){$f0_cg+=1;$f0_level+=$b[0];}} #第2 3 4 列为第一组3个实验重复
elsif($pos>=5 and $pos<=7){if($b[0] ne "-"){$f1_cg+=1;$f1_level+=$b[0];}}#第5 6 7 列为第二组3个实验重复
elsif($pos>=8 and $pos<=10){if($b[0] ne "-"){$f2_cg+=1;$f2_level+=$b[0];}}#第57 8 9 列为第三组3个实验重复
}
}
#遍历一个DMR区间内的每个样品三个实验重复的每个CG位点甲基化水平并相加存储到$(group)_level中，同时统计参与计算的CG位点数$(group)_cg。
#下方为计算这个DMR的平均甲基化水平，测序深度为0的不参与计算，用“-”表示，如下
if($f0_cg!=0){$f0=$f0_level/$f0_cg;}else{$f0="-";};
if($f1_cg!=0){$f1=$f1_level/$f1_cg;}else{$f1="-";};
if($f2_cg!=0){$f2=$f2_level/$f2_cg;}else{$f2="-";};
}
return $f0,$f1,$f2; #返回每个组的每个DMR的甲基化水平
close MATC;
}
子程序2 与子程序1完全相同，只不过参与遍历的文件句柄由MATC变为了MATS
sub methyl_level_s{
my @a;
my @b;
my $f0;
my $f1;
my $f2;
my $f0_cg=0;
my $f1_cg=0;
my $f2_cg=0;
my $f0_level=0;
my $f1_level=0;
my $f2_level=0;
open MATS,"$ARGV[2]";
<MATS>;
while(<MATS>){
chomp;
@a=split/\t/,$_;
if($_[0]=~/$a[0]/ && $a[1]>=$_[1] && $a[1]<=$_[2]){
foreach my $pos(2..10){
@b=split/,/,$a[$pos];
if($pos<=4){if($b[0] ne "-"){$f0_cg+=1;$f0_level+=$b[0];}}
elsif($pos>=5 and $pos<=7){if($b[0] ne "-"){$f1_cg+=1;$f1_level+=$b[0];}}
elsif($pos>=8 and $pos<=10){if($b[0] ne "-"){$f2_cg+=1;$f2_level+=$b[0];}}
}
}
if($f0_cg!=0){$f0=$f0_level/$f0_cg;}else{$f0="-";};
if($f1_cg!=0){$f1=$f1_level/$f1_cg;}else{$f1="-";};
if($f2_cg!=0){$f2=$f2_level/$f2_cg;}else{$f2="-";};
}
return $f0,$f1,$f2;
close MATS;}
网友评论