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 "$_"; }
相关文章推荐
- poj蚂蚁问题
- android实现层级式导航
- Pool进程池创建大量子进程,进程间通信
- 页面锚定
- 大道至简第三篇阅读笔记
- 如何将开源项目部分代码作为private放在github上?
- 指针和句柄的区别
- 设计模式之单例模式
- javascript && 和 || 最清晰的描述
- 大道至简第二篇阅读笔记
- DWR如何获得返回对象
- Java NIO 01=====概述
- Python高级编程–正则表达式(习题)
- PHP自动加载
- 大道至简第一篇阅读笔记
- 微信开放平台JS SDK接入sha1算法
- linux的page cache策略
- android 之 Camera使用示例
- Maven实战——远程仓库的配置
- 第十二周项目4-点、圆的关系