Thursday, July 04, 2013

Not mine - I found this csfasta to fastq conversion script that worked best...

#!/usr/bin/perl -w
use strict;
use Getopt::Long;

my $usage = qq(
This script is used for converting SOLiD/ABi FASTQ to Standard(Sanger) FASTQ
Usage :  perl $0 [Option]
Option:  -seq           seq file
         -qual          qual file
         -fastq         FASTQ file
         -o             output file
         -help          Help Information
Example: perl $0 -fastq solid.fastq -o std.fastq

Author : BENM <BinxiaoFeng\@gmail.com>
Version: 1.2
Date   : 2009-10-06
Update : 2012-03-27
\n);
my ($Seq,$Qual,$Fastq,$Output,$Header,$Help);
my %opts;
GetOptions
(
    \%opts,
    "seq:s"=>\$Seq,
    "qual:s"=>\$Qual,
    "fastq:s"=>\$Fastq,
    "o:s"=>\$Output,
    "header:s"=>\$Header,
    "help"=>\$Help
);

die($usage) if (((!defined $Seq)||(!defined $Qual))&&(!defined $Fastq)&&(!defined $Output)||($Help));

# Solexa->Sanger quality conversion table
my @conv_table;
for (-64..64) {
  $conv_table[$_+64] = chr(int(33 + 10*log(1+10**($_/10.0))/log(10)+.499));
}


# SOLiD color code
my @code = ([0,1,2,3],[1,0,3,2],[2,3,0,1],[3,2,1,0]);
my @bases = qw(A C G T);
my %decode = ();
foreach my $i(0..3)
{
    foreach my $j(0..3)
    {
        $decode{$code[$i]->[$j]} -> {$bases[$i]} = $bases[$j];
    }
}

open (OUT,">$Output") || die "Fail to create FASTA file:$Output for writing\n";
if (defined $Fastq)
{
    open (IN,$Fastq) || die "Fail to open FASTQ file:$Fastq for reading\n";
    while (<IN>)
    {
        if (/^@/)
        {
            print OUT $_;
            $_ = <IN>;
            s/\s+$//;
            my $seq = ($_=~/\d+/) ? col2base($_) : $_;
            print OUT "$seq\n";
            $_ = <IN>; print OUT $_; $_ = <IN>;
            s/\s+$//;
            my @t = split;
            my $qual = '';
            if ($_=~/^\d+\s*/)
            {
                my @t=split;
                $qual .= $conv_table[$_+64] for (@t); #intfq: SOLiD -> Standard
            }
            else
            {
                $qual .= $_;
            }
            substr($qual,0,1,'') if (length($qual)>length($seq));
            print OUT "$qual\n";
        }
    }
    close IN;
}
elsif ((defined $Seq)&&(defined $Qual))
{
    my ($Title,$Name,$seq,$qual)=("","","","");
    open (IN1,$Seq) || die "Fail to open Seq file:$Seq for reading\n";
    open (IN2,$Qual) || die "Fail to open Quality file:$Qual for reading\n";
    while(<IN2>)
    {
        if ($_=~/^\#\s+Title\:\s+(\S+)/)
        {
            $Title=$1;
            my $in1=<IN1>;
            while($in1 !~ /^\#\s+Title\:\s+(\S+)/)
            {
                $in1=<IN1>;
            }
        }
        else
        {
            if ($_=~/^#\s\w+/)
            {
                <IN2>;<IN1>;
            }
            elsif (/^[\@\>](\S+)/)
            {
                chomp $seq;
                chomp $qual;
                substr($qual,0,1,"") if (length($qual)>length($seq));
                warn "The length of quality is not match with seqence!\nSeq: $seq\nQual: $qual\n" if (length($qual)!= length($seq));
                print OUT "$Name\n$seq\n+\n$qual\n" if ($seq ne "")&&($qual ne "");
                ($seq,$qual)=("","");
                $Name="@".$Title.$1;
                $_ = <IN1>;
                if (/^[\@\>]/)
                {
                    $_ = <IN1>;
                }
                my $tmp_seq;
                while (($_ ne "")&&($_ !~ /^[\@\>]/)&&(!eof))
                {
                    s/\s+$//;
                    $tmp_seq.=$_;
                    $_=<IN1>;
                    chomp $_;
                }
                $tmp_seq .= $_ if ((eof)&&($_ !~ /^[\@\>]/));
                chomp $tmp_seq;
                $seq = (($tmp_seq=~/\d+/)) ? col2base($tmp_seq) : $tmp_seq;
                $_ = <IN2>;
                s/^\s+//;
                s/\s+$//;
                if (($_=~/\d+\s+\d+/)||($_=~/^\d+$/))
                {
                    $qual .= " " if ($qual=~/\d+$/);
                    my @t = split /\s+/,$_;
                    $qual .= $conv_table[$_+64] for (@t);
                }
                else
                {
                    $qual .= $_;
                }
            }
            else
            {
                s/^\s+//;
                s/\s+$//;
                if (($_=~/\d+\s+\d+/)||($_=~/^\d+$/))
                {
                    $qual .= " " if ($qual=~/\d+$/);
                    my @t = split /\s+/,$_;
                    $qual .= $conv_table[$_+64] for (@t);
                }
                else
                {
                    $qual .= $_;
                }
            }
        }
    }
   
    substr($qual,0,1,"") if (length($qual)>length($seq));
    print OUT "$Name\n$seq\n+\n$qual\n" if ($seq ne "")&&($qual ne "");
    close IN1;
    close IN2;
}
close OUT;

sub col2base
{
    my $col = shift;
    my @colors = split '',$col;
    my $string = $colors[0];
    if($string !~ /[acgt]/i){
        warn "$col has no header base\n";
        return 0;
    }
    my $last_base = $string;
    my $current_base = '';
    for(my $i=1;$i<@colors;$i++)
    {
        if (($last_base=~/N/i)&&(($colors[$i] eq ".")||($colors[$i]==5)))
        {
            $current_base = $bases[int(rand(@bases))];
        }
        else
        {
            $current_base = (exists $decode{$colors[$i]}->{$last_base}) ? $decode{$colors[$i]}->{$last_base} : "N";
        }
        $string .= $current_base;
        $last_base = $current_base;
    }
    substr($string,0,1,"");
    return $string;
}

No comments: