您的位置:首页 > 其它

unmapbam to fastq和自己的annovar格式~~~

2015-05-31 19:01 302 查看
#!perl
use warnings;
use strict;

die "perl $0 <unmaped.bam> <outprefix>\n" if @ARGV != 2;

my %hash;
open BAM, "samtools view $ARGV[0] |" or die $!;
while(<BAM>)
{
chomp;
my @tmp = split;
push @{$hash{$tmp[0]}}, "$tmp[1]\t$tmp[9]\t$tmp[10]";
}

open OUT1, "| gzip > $ARGV[1]_1.fq.gz" or die $!;
open OUT2, "| gzip > $ARGV[1]_2.fq.gz" or die $!;
foreach(keys %hash)
{
next if(@{$hash{$_}} == 1);
foreach my $i(@{$hash{$_}})
{
my @str = split /\t/, $i;
if($str[0] == 69)
{
print OUT1 "\@$_/1\n$str[1]\n+\n$str[2]\n";
}elsif($str[0] == 133){
print OUT2 "\@$_/2\n$str[1]\n+\n$str[2]\n";
}else{
last;
}
}
}


#!perl
use warnings;
use strict;

die "perl $0 <VCF>
Note: VCF to new annovar format.\n" if(@ARGV != 1);

my (@lines, @names, @out);
open VCF, $ARGV[0] or die $!;
while(<VCF>)
{
chomp;
next if(/^##/);
@lines = split;
next if($lines[4] =~ /,/);
if(/^#{1}/)
{
@names = @lines[9..$#lines];
next;
}
my (@type, @num);
for(my $i = 0; $i < @names; $i ++)
{
my @hh = split /:/, $lines[$i + 9];
$hh[0] =~ s/\//\|/;
push @type, $hh[0];
if(defined $hh[1])
{
push @num, "$hh[1]";
}else{
push @num, ".,.";
}
}
my $t = join "\t", @type;
my $n = join "\t", @num;
if(length($lines[3]) > 1)
{
my $del = length($lines[3]) - 1;
my $start = $lines[1] + 1;
my $end = $lines[1] + $del;
my @str = split //, $lines[3];
shift @str;
my $change = join "", @str;
push @out, "$lines[0]\t$start\t$end\t$change\t-\tDeletion\t$lines[2]\t$lines[5]\t$t\t$n\n";
}elsif(length($lines[4]) > 1){
my @str = split //, $lines[4];
shift @str;
my $change = join "", @str;
push @out, "$lines[0]\t$lines[1]\t$lines[1]\t-\t$change\tInsertion\t$lines[2]\t$lines[5]\t$t\t$n\n";
}else{
push @out, "$lines[0]\t$lines[1]\t$lines[1]\t$lines[3]\t$lines[4]\tSNP\t$lines[2]\t$lines[5]\t$t\t$n\n";
}
}
my @ts = map {$_ . "_type"} @names;
my @rs = map {$_ . "_read"} @names;
my $t = join "\t", @ts;
my $r = join "\t", @rs;
print "#chr\tstart\tend\tref\tmut\tvariation_type\tdatabase\tquality\t$t\t$r\n";
foreach(@out)
{
print "$_";
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: