Perl学习日记#杂记

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。

©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容