搜索
查看: 124|回复: 2

[求助] perl小白,一个与gz有关的问题

[复制链接]

3

主题

4

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2020-4-24 10:07:20 | 显示全部楼层 |阅读模式
  为了得出两个vcf在捕获区域下的一致性写出了以下代码,但神奇的是在vcf非压缩状态时可以正常工作,在gz压缩时就产生了bug,具体为运算失效:以下统计的SNP/INDEL总数为零;
  各位有时间的话请看看,多谢!



#!/usr/bin/perl
use File::Basename;
if($ARGV[0]=~/gz$/){
        open A,"gzip -cd $ARGV[0] | "  or die "$!";}
else{
        open A,$ARGV[0] or die "$!";}
if($ARGV[1]=~/gz$/){
        open B,"gzip -cd $ARGV[1] | "  or die "$!";}
else{
        open B,$ARGV[1] or die "$!";}            #AB是vcf文件
if($ARGV[2]=~/gz$/){
        open C,"gzip -cd $ARGV[2] | "  or die "$!";}
else{
        open C,$ARGV[2] or die "$!";}           #C是捕获区域
open OUT,">$ARGV[3]" or die "$!";

my %hash_s,%hash_i;

while(<C>){
    chomp;
    my @s=split;
    while(<A>){
        chomp;
        next if(/#/);
        next if(/^(?!.*PASS)/);
        my @t=split;
        next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
        my $id=join "-",$t[0],$t[1],$t[3],$t[4];
        if(/Dels=0.00/){                                                                  #     运
            $hash_s{$id}=0;                                                             #     算
            $A_snp++;}                                                                   #      异
        elsif(/^(?!.*Dels)/){                                                             #      常
            $hash_i{$id}=0;                                                             #       区
            $A_indel++;}                                                                 #       域
}
seek A,0,0;
}
seek C,0,0;
while(<C>){
    chomp;
    my @s=split;
    while(<B>){
        chomp;
        next if(/#/);
        next if(/^(?!.*PASS)/);
        my @t=split;
        next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
        my $id=join "-",$t[0],$t[1],$t[3],$t[4];
        if(/Dels=0.00/){                                                             #    同
            $B_snp++;}                                                              #
        if(/^(?!.*Dels)/){                                                           #
            $B_indel++;}                                                            #     上
        if(exists $hash_s{$id}){                                                  
                $S_share++;
        }
        elsif(exists $hash_i{$id}){
                $I_share++;
        }
}
seek B,0,0;
}
$AS_novel=$A_snp-$S_share;
$AI_novel=$A_indel-$I_share;
$BS_novel=$B_snp-$S_share;
$BI_novel=$B_indel-$I_share;
if($A_snp==0or$A_indel==0or$B_snp==0or$B_indel==0){
print "SNP/INDEL = 0!\n" ;}
else{
$AS_share=sprintf("%.2f",$S_share/$A_snp*100);
$AI_share=sprintf("%.2f",$I_share/$A_indel*100);
$BS_share=sprintf("%.2f",$S_share/$B_snp*100);
$BI_share=sprintf("%.2f",$I_share/$B_indel*100);

print OUT "SNP\nS_share\tAS_share\tBS_share\tAS_novel\tBS_novel\n";
print OUT "$S_share\t$AS_share%\t$BS_share%\t$AS_novel\t$BS_novel\n";
print OUT "INDEL\nI_share\tAI_share\tBI_share\tAI_novel\tBI_novel\n";
print OUT "$I_share\t$AI_share%\t$BI_share%\t$AI_novel\t$BI_novel\n";
}
close A;
close B;
close C;
close OUT;




上一篇:次级图例要怎么加?
下一篇:遇到arguments must have same length
回复

使用道具 举报

3

主题

4

帖子

53

积分

注册会员

Rank: 2

积分
53
 楼主| 发表于 2020-4-24 11:16:27 | 显示全部楼层
我自己的猜测,有可能是gzip打开的文件无法进行多次循环,如果确实是这样的话应该如何解决呢?
回复 支持 反对

使用道具 举报

3

主题

4

帖子

53

积分

注册会员

Rank: 2

积分
53
 楼主| 发表于 2020-4-24 11:50:46 | 显示全部楼层

[求助]perl小白,一个关于gz的问题

为了得出两个vcf在捕获区域下的一致性写出了以下代码,但神奇的是在vcf非压缩状态时可以正常工作,在gz压缩时就产生了bug,具体为运算失效:以下统计的SNP/INDEL总数为零;
  各位有时间的话请看看,多谢!



#!/usr/bin/perl
use File::Basename;
if($ARGV[0]=~/gz$/){
        open A,"gzip -cd $ARGV[0] | "  or die "$!";}
else{
        open A,$ARGV[0] or die "$!";}
if($ARGV[1]=~/gz$/){
        open B,"gzip -cd $ARGV[1] | "  or die "$!";}
else{
        open B,$ARGV[1] or die "$!";}            #AB是vcf文件
if($ARGV[2]=~/gz$/){
        open C,"gzip -cd $ARGV[2] | "  or die "$!";}
else{
        open C,$ARGV[2] or die "$!";}           #C是捕获区域
open OUT,">$ARGV[3]" or die "$!";

my %hash_s,%hash_i;

while(<C>){
    chomp;
    my @s=split;
    while(<A>){
        chomp;
        next if(/#/);
        next if(/^(?!.*PASS)/);
        my @t=split;
        next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
        my $id=join "-",$t[0],$t[1],$t[3],$t[4];
        if(/Dels=0.00/){                                                                  #     运
            $hash_s{$id}=0;                                                             #     算
            $A_snp++;}                                                                   #      异
        elsif(/^(?!.*Dels)/){                                                             #      常
            $hash_i{$id}=0;                                                             #       区
            $A_indel++;}                                                                 #       域
}
seek A,0,0;
}
seek C,0,0;
while(<C>){
    chomp;
    my @s=split;
    while(<B>){
        chomp;
        next if(/#/);
        next if(/^(?!.*PASS)/);
        my @t=split;
        next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
        my $id=join "-",$t[0],$t[1],$t[3],$t[4];
        if(/Dels=0.00/){                                                             #    同
            $B_snp++;}                                                              #
        if(/^(?!.*Dels)/){                                                           #
            $B_indel++;}                                                            #     上
        if(exists $hash_s{$id}){                                                  
                $S_share++;
        }
        elsif(exists $hash_i{$id}){
                $I_share++;
        }
}
seek B,0,0;
}
$AS_novel=$A_snp-$S_share;
$AI_novel=$A_indel-$I_share;
$BS_novel=$B_snp-$S_share;
$BI_novel=$B_indel-$I_share;
if($A_snp==0or$A_indel==0or$B_snp==0or$B_indel==0){
print "SNP/INDEL = 0!\n" ;}
else{
$AS_share=sprintf("%.2f",$S_share/$A_snp*100);
$AI_share=sprintf("%.2f",$I_share/$A_indel*100);
$BS_share=sprintf("%.2f",$S_share/$B_snp*100);
$BI_share=sprintf("%.2f",$I_share/$B_indel*100);

print OUT "SNP\nS_share\tAS_share\tBS_share\tAS_novel\tBS_novel\n";
print OUT "$S_share\t$AS_share%\t$BS_share%\t$AS_novel\t$BS_novel\n";
print OUT "INDEL\nI_share\tAI_share\tBI_share\tAI_novel\tBI_novel\n";
print OUT "$I_share\t$AI_share%\t$BI_share%\t$AI_novel\t$BI_novel\n";
}
close A;
close B;
close C;
close OUT;

回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-6 09:53 , Processed in 0.027514 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.