started 12/05/2020
之前开始就经常需要自己处理一些简单的文本数据,现在把相关的具体做每个小项的东西记录下来,以后方便用
处理bs-seq对比后的数据
最近拿到了一组甲基化测序之后探测到的发生甲基化的位点的数据,大概是截取了基因组的所有CHH、CG、CHG区段,然后将对应区域的甲基化的C和未甲基化的C的数量记录下来,形式如下
1 1001 + 0 0 CHH CTA
1 1006 + 0 0 CHH CCC
1 1007 + 0 0 CHH CCT
1 1008 + 0 0 CHH CTA
1 1013 + 0 0 CHH CCC
1 1014 + 0 0 CHH CCT
1 1015 + 0 0 CHH CTA
文件没header,但是意思其实很好猜
现在的需求是将对应区段的所有行筛选出来,只需要用perl进行一个简单的分割以及判断即可
(vim写代码还是有点不习惯orz)代码如下:
#!/usr/bin/perl
@ARGV==1 or die;
$infile1 = $ARGV[0];
print STDOUT "Please input target chromosome:";
$chr_t = <STDIN>;
print STDOUT "Please input target region(F):";
$region_f = <STDIN>;
$region_f =~ s/\s+$//;
print STDOUT "Please input target region(B):";
$region_b = <STDIN>;
$region_b =~ s/\s+$//;
$infile1 =~ /.*\./;
$outfile = $&;
$outfile =~ s/\.//;
$sample = $outfile;
$outfile .= "_result_";
$outfile .= $region_f;
$outfile .= "_";
$outfile .= $region_b;
$outfile .= ".txt";
open ($bs1, "<", $infile1) or die;
open ($out, ">", $outfile) or die;
print $out "sample:$sample region: $region_f-$region_b\n";
foreach(<$bs1>)
{
$_=~ s/\s+$//;
$now = $_;
$now1 = $now;
@arr = split(/ /,$now);
$region = $arr[1];
$chr = $arr[0];
print STDOUT "region=$region chr=$chr\n";
if($chr < $chr_t || $region < $region_f)
{
next;
}
elsif($region >= $region_f && $region <= $region_b && $chr == $chr_t)
{
print $out "$now1\n";
}
else
{
last;
}
}
close $bs1;
close $out;
本以为就这么轻松的结束了,然后发现,3.8个g的txt,我的16gb内存的小macbookpro有点搞不定,文件半天打不开,一直开在open那里。。。这就有点绝望了。。。然后我发现我的电脑艰难的打开了,然后我发现程序顺利地跑完了??
ok,fine。