Perl学习19之生信简单运用

Posted 2021-04-17 pythonic生物人

"pythonic生物人"的第33篇分享。

---

原创不易，点个“赞“或"在看"鼓励下呗！

摘要

Perl计算SAM文件中reads落在，每个1M bin区间中的reads数目，GC碱基数目，GC比。

正文开始啦

---

解决问题： Perl计算SAM文件中reads落在，每个1M bin区间中的reads数目，GC碱基数目，GC比。

• 输入文件1

file1，kemr区间文件，每个bin区间长1M。

chr1 7000001

chr1 8000001

chr1 9000001

chr1 10000001

chr1 11000001

chr1 12000001

chr1 14000001

chr1 15000001

chr1 16000001

chr1 18000001

chr1 19000001

chr1 20000001

• 输入文件2

file2，sam文件

• 问题解决代码：reads_inbin.pl

#!/usr/bin/perluse strict;use warnings;
#存储每个bin的区间{chr1=>{1000001=>0,2000001>0.......}}my (%bin,%GC,%ATGC);open BIN,"<",$ARGV[0];while(<BIN>){ chomp; my @line = split /	/; $bin{$line[0]}{$line[1]}=0;}close BIN;

open IN,"<",$ARGV[1];while(<IN>){ chomp; next if /^@/; my @line = split/	/,$_; next if $line[2] eq "*"; next if $line[4] < 25; #求每个bin区间例如，[1000001,1000001+1000000)上的reads数目 #keys %{$hash{line[2]}}取出所有1000001,二维哈希的循环写法 #foreach my $start (keys $bin{$line[2]}){#该写法报如下错误 #Type of argument to keys on reference must be unblessed hashref or arrayref foreach my $s (sort keys %{$bin{$line[2]}}){ if($line[3]>=$s && $line[3]<($s+1000000)){ $bin{$line[2]}{$s} ++; my $G = ($line[9] =~ s/G/G/g); my $C = ($line[9] =~ s/C/C/g); my $A = ($line[9] =~ s/A/A/g); my $T = ($line[9] =~ s/T/T/g); $GC{$line[2]}{$s} += ($G + $C); $ATGC{$line[2]}{$s} += ($G + $C + $A + $T);
 } }}close IN;
foreach (keys %bin){ my $chr=$_; foreach (keys $bin{$chr}){ printf "$chr	$_	$bin{$chr}{$_}	$GC{$chr}{$_}	$ATGC{$chr}{$_}	%.2f%%
",($GC{$chr}{$_}/$ATGC{$chr}{$_})*100; }}

perl reads_inbin.pl file1 test.sam

chr7 131000001 790 17650 39500 44.68%

chr7 79000001 706 12653 35300 35.84%

chr7 36000001 710 15303 35500 43.11%

chr7 24000001 725 14647 36250 40.41%

chr7 97000001 681 14593 34050 42.86%

chr7 47000001 796 18134 39800 45.56%

chr7 64000001 512 10840 25600 42.34%

chr7 33000001 749 14887 37450 39.75%

chr7 19000001 652 12071 32600 37.03%

。。。。。。。。。。。。。。。。。。。

同系列文章

---

原创不易"点赞"、"在看"鼓励下呗！