view cf_scripts/split_for_velvet.pl @ 3:032b13f5b087 draft

Uploaded

#!/usr/bin/perl

use strict;
my %read1;
my %read2;
my %qual1;
my %qual2;

while (my $line=<STDIN>) {
	chomp();
	my ($name, $flag, $chr1, $pos1, $mq, $cigar, $chr2, $pos2, $dist, $seq, $quals, $stuff)=split(/\s+/, $line);
	if($flag & 16) {	#reverse-complement it
		$seq=reverse($seq);
		$seq=~tr/ACGT/TGCA/;
		$quals=reverse($quals);
	}
	if($flag & 64) {	#first in pair
		$read1{$name}=$seq;
		$qual1{$name}=$quals;
	}
	else {
		$read2{$name}=$seq;
		$qual2{$name}=$quals;
	}
}

open(FQ1, ">$ARGV[0]") or die "Can't open $ARGV[0]";
open(FQ2, ">$ARGV[1]") or die "Can't open $ARGV[1]";
open(FQ3, ">$ARGV[2]") or die "Can't open $ARGV[2]";

foreach my $key(keys %read1) {
    
    	if(exists $read2{$key}) {
		print FQ1 "\@$key\n$read1{$key}\n+\n$qual1{$key}\n";
        	print FQ2 "\@$key\n$read2{$key}\n+\n$qual2{$key}\n";
	}
	else {
		print FQ3 "\@$key\n$read1{$key}\n+\n$qual1{$key}\n";	
	}
}
foreach my $key(keys %read2) {

	if(!(exists $read1{$key})) {
            	print FQ3 "\@$key\n$read2{$key}\n+\n$qual2{$key}\n";
	}
}
close(FQ1);
close(FQ2);
close(FQ3);