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

标签：

#!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 "$_";
}

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

标签：

原文地址：http://blog.csdn.net/skenoy/article/details/46290737