20251125
This commit is contained in:
@@ -0,0 +1,85 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
|
||||
use FindBin;
|
||||
use lib ("$FindBin::RealBin/../../../PerlLib");
|
||||
use Nuc_translator;
|
||||
|
||||
my $usage = "usage: $0 blat.psl.nameSorted.sam reads.tab.nameSorted\n\n";
|
||||
|
||||
|
||||
my $sam = $ARGV[0] or die $usage;
|
||||
my $seqs = $ARGV[1] or die $usage;
|
||||
|
||||
|
||||
main: {
|
||||
|
||||
open (my $sam_fh, "$sam") or die "Error, cannot open file $sam";
|
||||
open (my $seqs_fh, "$seqs") or die "Error, cannot open file $seqs";
|
||||
|
||||
my $seq_line = <$seqs_fh>;
|
||||
chomp $seq_line;
|
||||
my ($seq_acc, $seq, $qual) = split(/\t/, $seq_line);
|
||||
$seq_acc =~ s/\s//g; # rid any ws from acc name
|
||||
|
||||
unless ($qual) {
|
||||
$qual = 'B' x length($seq);
|
||||
}
|
||||
|
||||
while (my $sam_line = <$sam_fh>) {
|
||||
|
||||
my @x = split(/\t/, $sam_line);
|
||||
|
||||
my $acc = $x[0];
|
||||
my $flag = $x[1];
|
||||
|
||||
my $aligned_orient = ($flag & 0x0010) ? '-' : '+';
|
||||
|
||||
|
||||
while ($seq_acc lt $acc) {
|
||||
$seq_line = <$seqs_fh>;
|
||||
chomp $seq_line;
|
||||
($seq_acc, $seq, $qual) = split(/\t/, $seq_line);
|
||||
$seq_acc =~ s/\s//g; # no ws in acc name
|
||||
unless (defined $qual) {
|
||||
$qual = 'B' x length($seq); #$qual = "*";
|
||||
}
|
||||
}
|
||||
|
||||
if ($acc eq $seq_acc) {
|
||||
|
||||
if ($aligned_orient eq '-') {
|
||||
my $revseq = &reverse_complement($seq);
|
||||
my @q = split(//, $qual);
|
||||
my $revqual = join("", reverse(@q));
|
||||
$x[9] = $revseq;
|
||||
$x[10] = $revqual;
|
||||
}
|
||||
else {
|
||||
$x[9] = $seq;
|
||||
$x[10] = $qual;
|
||||
}
|
||||
}
|
||||
else {
|
||||
die "Error,\n[$acc]\nnot encountered in file: $seqs,\ncurrently cursor is at seq:\n[$seq_acc]\n";
|
||||
}
|
||||
|
||||
## set read quality score to a high value so that Scripture will use it.
|
||||
$x[4] = 255;
|
||||
$x[5] =~ s/H/S/gi; # convert hard to soft clips; samtools doesn't like H's
|
||||
|
||||
foreach my $val (@x) {
|
||||
unless (defined $val) {
|
||||
$val = "*";
|
||||
}
|
||||
}
|
||||
|
||||
print join("\t", @x);
|
||||
}
|
||||
|
||||
exit(0);
|
||||
}
|
||||
|
||||
|
||||
@@ -0,0 +1,125 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
|
||||
use FindBin;
|
||||
use Cwd;
|
||||
|
||||
use Getopt::Long qw(:config no_ignore_case bundling);
|
||||
|
||||
|
||||
$ENV{LC_ALL} = 'C';
|
||||
|
||||
my $usage = <<_EOUSAGE_;
|
||||
|
||||
################################################################################################
|
||||
#
|
||||
# Required:
|
||||
#
|
||||
# --genome genome in fasta format
|
||||
# --reads reads in fasta format
|
||||
#
|
||||
# Optional:
|
||||
#
|
||||
# --blat_params quote-delimited params to pass to blat, eg. "-q=rna -t=dna -maxIntron=1000" (default: "-q=rna -t=dna")
|
||||
#
|
||||
# --top number of top hits (default: 10)
|
||||
# --min_per_ID minimum percent identity (default: 95)
|
||||
#
|
||||
##############################################################################################
|
||||
|
||||
|
||||
_EOUSAGE_
|
||||
|
||||
;
|
||||
|
||||
|
||||
my ($genome_fa, $reads_fa);
|
||||
|
||||
my $blat_params = "-q=rna -t=dna";
|
||||
my $top_hits = 10;
|
||||
my $min_per_ID = 95;
|
||||
|
||||
|
||||
|
||||
&GetOptions( 'genome=s' => \$genome_fa,
|
||||
'reads=s' => \$reads_fa,
|
||||
|
||||
'blat_params=s' => \$blat_params,
|
||||
|
||||
'top=i' => \$top_hits,
|
||||
'min_per_ID=i' => \$min_per_ID,
|
||||
);
|
||||
|
||||
unless ($genome_fa && $reads_fa) {
|
||||
die $usage;
|
||||
}
|
||||
|
||||
|
||||
{
|
||||
my @required_progs = qw (blat psl2sam.pl);
|
||||
|
||||
foreach my $prog (@required_progs) {
|
||||
my $path = `sh -c "command -v $prog"`;
|
||||
unless ($path =~ /^\//) {
|
||||
die "Error, cannot locate required program: $prog";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
main: {
|
||||
|
||||
my $util_dir = "$FindBin::RealBin/../util";
|
||||
|
||||
my $cmd = "$util_dir/fasta_to_tab.pl $reads_fa > $reads_fa.tab";
|
||||
&process_cmd($cmd) unless (-s "$reads_fa.tab");
|
||||
|
||||
$cmd = "sort -T . -S 2G -k 1,1 $reads_fa.tab > $reads_fa.tab.sort";
|
||||
&process_cmd($cmd) unless (-s "$reads_fa.tab.sort");
|
||||
|
||||
# run blat
|
||||
$cmd = "blat $blat_params $genome_fa $reads_fa $reads_fa.psl";
|
||||
&process_cmd($cmd);
|
||||
|
||||
# convert to sam
|
||||
$cmd = "psl2sam.pl -q 0 -r 0 $reads_fa.psl | sort -T . -S 2G -k 1,1 > $reads_fa.psl.sam";
|
||||
&process_cmd($cmd);
|
||||
|
||||
|
||||
# add the reads
|
||||
$cmd = "$util_dir/blat_sam_add_reads2.pl $reads_fa.psl.sam $reads_fa.tab.sort > $reads_fa.psl.sam.wReads";
|
||||
&process_cmd($cmd);
|
||||
|
||||
## prune output to top matches:
|
||||
$cmd = "$util_dir/top_blat_sam_extractor.pl $reads_fa.psl.sam.wReads $top_hits $min_per_ID > $reads_fa.psl.sam.wReads.top";
|
||||
&process_cmd($cmd);
|
||||
|
||||
$cmd = "$FindBin::RealBin/cigar_tweaker $reads_fa.psl.sam.wReads.top $genome_fa > $reads_fa.psl.sam.wReads.top.tweaked";
|
||||
&process_cmd($cmd);
|
||||
|
||||
$cmd = "sort -T . -S 2G -k 3,3 -k 4,4n $reads_fa.psl.sam.wReads.top.tweaked > $reads_fa.psl.sam.wReads.top.tweaked.coordSorted.sam";
|
||||
&process_cmd($cmd);
|
||||
|
||||
|
||||
|
||||
exit(0);
|
||||
}
|
||||
|
||||
####
|
||||
sub process_cmd {
|
||||
my ($cmd) = @_;
|
||||
|
||||
print "CMD: $cmd\n";
|
||||
my $ret = system($cmd);
|
||||
|
||||
if ($ret) {
|
||||
die "Error, cmd: $cmd died with ret $ret";
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
@@ -0,0 +1,121 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
|
||||
my $SEE = 0;
|
||||
|
||||
###############
|
||||
# blat format: # Q=cDNA T=genomic
|
||||
################
|
||||
|
||||
# 0: match
|
||||
# 1: mis-match
|
||||
# 2: rep. match
|
||||
# 3: N's
|
||||
# 4: Q gap count
|
||||
# 5: Q gap bases
|
||||
# 6: T gap count
|
||||
# 7: T gap bases
|
||||
# 8: strand
|
||||
# 9: Q name
|
||||
# 10: Q size
|
||||
# 11: Q start
|
||||
# 12: Q end
|
||||
# 13: T name
|
||||
# 14: T size
|
||||
# 15: T start
|
||||
# 16: T end
|
||||
# 17: block count
|
||||
# 18: block Sizes
|
||||
# 19: Q starts
|
||||
# 20: T starts
|
||||
# 21: Q seqs (pslx format)
|
||||
# 22: T seqs (pslx format)
|
||||
|
||||
#############################################################
|
||||
## This script filters blat output and extracts only ##
|
||||
## the top scoring alignment chain for each accession. ##
|
||||
#############################################################
|
||||
|
||||
my $line_num = 0;
|
||||
my %data;
|
||||
|
||||
my $filename = $ARGV[0] or die "\n\nusage: $0 outputfile.psl [num_top_hits=1]\n\n";
|
||||
my $num_top_hits = $ARGV[1] || 1;
|
||||
|
||||
## First Pass, assign scores to entries associated with accessions.
|
||||
open (FILE, "$filename");
|
||||
my $line_num = 0;
|
||||
while (<FILE>) {
|
||||
$line_num++;
|
||||
my @x = split (/\t/);
|
||||
unless ($x[0] =~ /^\d/) {next;}
|
||||
my ($matches, $mismatches, $q_gap_num, $t_gap_num, $q_insert, $t_insert) = ($x[0], $x[1], $x[4], $x[6], $x[5], abs($x[7]));
|
||||
|
||||
my $score = &calculate_score($matches, $mismatches, $q_gap_num, $t_gap_num, $q_insert, $t_insert);
|
||||
|
||||
my $accession = $x[9];
|
||||
my $score_struct = {score=>$score,
|
||||
line_num=>$line_num};
|
||||
push (@{$data{$accession}}, $score_struct);
|
||||
}
|
||||
|
||||
close FILE;
|
||||
|
||||
## Identify each highest scoring match
|
||||
my %line_nums_to_print;
|
||||
|
||||
foreach my $accession (keys %data) {
|
||||
my @hits = @{$data{$accession}};
|
||||
@hits = reverse sort {$a->{score}<=>$b->{score}} @hits; #sort in reverse order of score.
|
||||
|
||||
if ($SEE) { #verify the top hit is chosen.
|
||||
foreach my $hit (@hits) {
|
||||
print $hit->{line_num} . ":" . $hit->{score} . " ";
|
||||
}
|
||||
|
||||
print "\n chose ";
|
||||
}
|
||||
|
||||
for (my $i = 0; $i < $num_top_hits && $i <= $#hits; $i++) {
|
||||
|
||||
my $top_hit = $hits[$i];
|
||||
print $top_hit->{line_num} . ":" . $top_hit->{score} . "\n" if $SEE;
|
||||
my $line_num = $top_hit->{line_num};
|
||||
$line_nums_to_print{$line_num} = 1;
|
||||
}
|
||||
}
|
||||
|
||||
open (FILE, "$filename");
|
||||
$line_num = 0;
|
||||
while (<FILE>) {
|
||||
$line_num++;
|
||||
if ($line_nums_to_print{$line_num}) {
|
||||
print;
|
||||
}
|
||||
}
|
||||
close FILE;
|
||||
|
||||
exit(0);
|
||||
|
||||
|
||||
####
|
||||
sub calculate_score {
|
||||
my ($match, $mismatch, $q_gap_num, $t_gap_num, $q_insert, $t_insert) = @_;
|
||||
|
||||
## JKent's code in pslFilter:
|
||||
# score = (psl->match + psl->repMatch)*reward - psl->misMatch*cost
|
||||
# - (psl->qNumInsert + psl->tNumInsert + 1) * gapOpenCost
|
||||
# - log(psl->qBaseInsert + psl->tBaseInsert + 1) * gapSizeLogMod;
|
||||
|
||||
#use jkent's default score parameters:
|
||||
my $reward = 1;
|
||||
my $cost = 1;
|
||||
my $gapOpenCost = 4;
|
||||
my $gapSizeLogMod = 1;
|
||||
|
||||
return ( ($match * $reward) - ($mismatch * $cost) -
|
||||
( ($q_gap_num + $t_gap_num) * $gapOpenCost) -
|
||||
( log ($q_insert + $t_insert + 1) * $gapSizeLogMod) );
|
||||
|
||||
}
|
||||
@@ -0,0 +1,277 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use FindBin;
|
||||
use lib ("$FindBin::RealBin/../../../PerlLib");
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
use threads;
|
||||
|
||||
use Fasta_reader;
|
||||
use Process_cmd;
|
||||
use Thread_helper;
|
||||
|
||||
use Getopt::Long qw (:config no_ignore_case bundling);
|
||||
use vars qw ($DEBUG $opt_h $opt_g $opt_t $opt_c $opt_o $opt_B $opt_I);
|
||||
|
||||
my $CPU = 1;
|
||||
my $num_top_hits = 1;
|
||||
my $KEEP_PSLX = 0;
|
||||
|
||||
&GetOptions( 'g=s' => \$opt_g,
|
||||
'd' => \$DEBUG,
|
||||
'h' => \$opt_h,
|
||||
'c=s' => \$opt_c,
|
||||
'o=s' => \$opt_o,
|
||||
't=s' => \$opt_t,
|
||||
'I=i' => \$opt_I,
|
||||
'N=i' => \$num_top_hits,
|
||||
'CPU=i' => \$CPU,
|
||||
'KEEP_PSLX' => \$KEEP_PSLX,
|
||||
);
|
||||
|
||||
our $SEE = 0;
|
||||
|
||||
$|++;
|
||||
|
||||
my $MAX_INTRON = $opt_I || 100000;
|
||||
|
||||
my $usage = <<_EOH_;
|
||||
|
||||
Script chunks EST alignments into more manageable data sets.
|
||||
|
||||
############################# Options ###############################
|
||||
#
|
||||
# -g <string> genomic_seq.db
|
||||
# -t <string> transcripts database
|
||||
# -I <int> maximum intron length (default: 100000)
|
||||
# -o <string> prefix for output file (default: 'blat')
|
||||
# -N <int> number of top hits (default: $num_top_hits)
|
||||
#
|
||||
# --CPU <int> number of threads (default: 1)
|
||||
#
|
||||
# -h this help menu
|
||||
# -d debug mode
|
||||
# --KEEP_PSLX retain the raw blat output files
|
||||
#
|
||||
###################### Process Args and Options #####################
|
||||
|
||||
|
||||
_EOH_
|
||||
|
||||
;
|
||||
|
||||
|
||||
my $genome_db = $opt_g;
|
||||
my $transcript_db = $opt_t;
|
||||
my $output_prefix = $opt_o || "blat";
|
||||
my $blat_path = "blat";
|
||||
my $util_dir = $FindBin::RealBin;
|
||||
|
||||
unless ($genome_db && $transcript_db) {
|
||||
die "$usage\n";
|
||||
}
|
||||
|
||||
my $ooc_cmd = "$blat_path $genome_db $transcript_db -q=rna -dots=100 -maxIntron=$MAX_INTRON -makeOoc=11.ooc nada";
|
||||
|
||||
my $blat_thr;
|
||||
|
||||
unless (-s "11.ooc") {
|
||||
$blat_thr = threads->create('process_cmd', $ooc_cmd);
|
||||
}
|
||||
|
||||
my $blat_partitions_dir = "blat_out_dir";
|
||||
unless (-d $blat_partitions_dir) {
|
||||
mkdir($blat_partitions_dir) or die "Error, cannot mkdir $blat_partitions_dir";
|
||||
}
|
||||
|
||||
my $num_seqs = `grep '>' $transcript_db | wc -l `;
|
||||
$num_seqs =~ s/\s//g;
|
||||
unless ($num_seqs && $num_seqs =~ /^\d+$/) {
|
||||
die "Error, cannot determine number of fasta entries in $transcript_db";
|
||||
}
|
||||
|
||||
my $seqs_per_partition = int($num_seqs/$CPU + 0.5);
|
||||
unless ($seqs_per_partition > 1) {
|
||||
$seqs_per_partition = 1;
|
||||
}
|
||||
|
||||
my @transcript_files = &partition_transcript_db($transcript_db, $seqs_per_partition, $blat_partitions_dir);
|
||||
|
||||
if ($blat_thr) {
|
||||
$blat_thr->join();
|
||||
if ($blat_thr->error()) {
|
||||
die "Error, $ooc_cmd died ...";
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
###############################
|
||||
## Run BLAT
|
||||
###############################
|
||||
my @pslx_files;
|
||||
|
||||
{
|
||||
|
||||
my $thread_helper = new Thread_helper($CPU);
|
||||
my %thread_to_checkpoint;
|
||||
|
||||
foreach my $transcript_file (@transcript_files) {
|
||||
|
||||
|
||||
## process blat search:
|
||||
my $cmd = "$blat_path $genome_db $transcript_file -q=rna -dots=100 "
|
||||
. " -maxIntron=$MAX_INTRON -out=pslx -ooc=11.ooc $transcript_file.pslx";
|
||||
|
||||
|
||||
my $checkpoint_file = "$transcript_file.pslx.completed";
|
||||
unless (-e $checkpoint_file) {
|
||||
|
||||
my $thread = threads->create('process_cmd', $cmd);
|
||||
$thread_helper->add_thread($thread);
|
||||
|
||||
my $thread_id = $thread->tid();
|
||||
$thread_to_checkpoint{$thread_id} = {thread => $thread,
|
||||
checkpoint => $checkpoint_file,
|
||||
};
|
||||
}
|
||||
push (@pslx_files, "$transcript_file.pslx");
|
||||
|
||||
}
|
||||
|
||||
$thread_helper->wait_for_all_threads_to_complete();
|
||||
|
||||
## write checkpoints for successful threads:
|
||||
foreach my $thread_id (keys %thread_to_checkpoint) {
|
||||
my $thread_info_href = $thread_to_checkpoint{$thread_id};
|
||||
my $thread = $thread_info_href->{thread};
|
||||
unless ($thread->error()) {
|
||||
my $checkpoint = $thread_info_href->{checkpoint};
|
||||
system("touch $checkpoint");
|
||||
}
|
||||
}
|
||||
|
||||
if (my @failures = $thread_helper->get_failed_threads()) {
|
||||
die "Error, ". scalar(@failures) . " blat searches failed. ";
|
||||
}
|
||||
}
|
||||
|
||||
######################
|
||||
## get top hits.
|
||||
######################
|
||||
|
||||
|
||||
my @top_hits_files;
|
||||
|
||||
|
||||
{
|
||||
|
||||
my %thread_to_checkpoint;
|
||||
my $thread_helper = new Thread_helper($CPU);
|
||||
|
||||
foreach my $pslx_file (@pslx_files) {
|
||||
|
||||
my $cmd = "$util_dir/blat_top_hit_extractor.pl $pslx_file $num_top_hits > $pslx_file.top_${num_top_hits}";
|
||||
|
||||
my $completed_checkpoint_file = "$pslx_file.top_${num_top_hits}.completed";
|
||||
unless (-e $completed_checkpoint_file) {
|
||||
|
||||
my $thread = threads->create('process_cmd', $cmd);
|
||||
$thread_helper->add_thread($thread);
|
||||
|
||||
my $thread_id = $thread->tid();
|
||||
$thread_to_checkpoint{$thread_id} = { thread => $thread,
|
||||
checkpoint => $completed_checkpoint_file,
|
||||
};
|
||||
|
||||
}
|
||||
push (@top_hits_files, "$pslx_file.top_${num_top_hits}");
|
||||
}
|
||||
|
||||
$thread_helper->wait_for_all_threads_to_complete();
|
||||
|
||||
## write checkpoints for successful threads:
|
||||
foreach my $thread_id (keys %thread_to_checkpoint) {
|
||||
my $thread_info_href = $thread_to_checkpoint{$thread_id};
|
||||
my $thread = $thread_info_href->{thread};
|
||||
unless ($thread->error()) {
|
||||
my $checkpoint = $thread_info_href->{checkpoint};
|
||||
system("touch $checkpoint");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (my @failures = $thread_helper->get_failed_threads()) {
|
||||
die "Error, ". scalar(@failures) . " blat top hit selectors failed. ";
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
if (-s "$output_prefix.gff3") {
|
||||
print STDERR "WARNING, REPLACING EXISTING FILE: $output_prefix.gff3. KILL THIS NOW TO PREVENT THIS. (you have 10 seconds)\n";
|
||||
sleep(10);
|
||||
print STDERR "OK, too late. replacing it now.\n";
|
||||
unlink("$output_prefix.gff3");
|
||||
}
|
||||
|
||||
foreach my $top_hits_file (@top_hits_files) {
|
||||
|
||||
# convert to gff3 format
|
||||
print STDERR "-converting $top_hits_file to gff3\n";
|
||||
|
||||
my $cmd = "$util_dir/pslx_to_gff3.pl < $top_hits_file >> $output_prefix.gff3";
|
||||
&process_cmd($cmd);
|
||||
}
|
||||
|
||||
## clean up the pslx files we no longer need.
|
||||
foreach my $pslx_file (@pslx_files) {
|
||||
unlink($pslx_file) unless $KEEP_PSLX; # these files can be huge. Once have top hits, no longer need all hits (hopefully).
|
||||
}
|
||||
|
||||
print STDERR "done.\n";
|
||||
|
||||
exit(0);
|
||||
|
||||
|
||||
####
|
||||
sub partition_transcript_db {
|
||||
my ($transcript_db, $seqs_per_partition, $blat_partitions_dir) = @_;
|
||||
|
||||
my @files;
|
||||
|
||||
my $checkpoint_file = "$blat_partitions_dir/partitions.completed";
|
||||
|
||||
if (-e $checkpoint_file) {
|
||||
my @files = <$blat_partitions_dir/partition.*.fa>;
|
||||
return(@files);
|
||||
}
|
||||
else {
|
||||
|
||||
|
||||
my $fasta_reader = new Fasta_reader($transcript_db);
|
||||
my $partition_counter = 0;
|
||||
|
||||
my $counter = 0;
|
||||
my $ofh;
|
||||
|
||||
while (my $seq_obj = $fasta_reader->next()) {
|
||||
my $fasta_entry = $seq_obj->get_FASTA_format();
|
||||
if ($counter % $seqs_per_partition == 0) {
|
||||
close $ofh if $ofh;
|
||||
$partition_counter++;
|
||||
my $outfile = "$blat_partitions_dir/partition.$counter.fa";
|
||||
open ($ofh, ">$outfile") or die "Error, cannot write to outfile: $outfile";
|
||||
push (@files, $outfile);
|
||||
}
|
||||
print $ofh $fasta_entry;
|
||||
$counter++;
|
||||
}
|
||||
|
||||
close $ofh if $ofh;
|
||||
|
||||
system("touch $checkpoint_file");
|
||||
|
||||
return(@files);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,242 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
|
||||
|
||||
################
|
||||
# blat format: # Q=cDNA T=genomic
|
||||
################
|
||||
|
||||
# 0: match
|
||||
# 1: mis-match
|
||||
# 2: rep. match
|
||||
# 3: N's
|
||||
# 4: Q gap count
|
||||
# 5: Q gap bases
|
||||
# 6: T gap count
|
||||
# 7: T gap bases
|
||||
# 8: strand
|
||||
# 9: Q name
|
||||
# 10: Q size
|
||||
# 11: Q start
|
||||
# 12: Q end
|
||||
# 13: T name
|
||||
# 14: T size
|
||||
# 15: T start
|
||||
# 16: T end
|
||||
# 17: block count
|
||||
# 18: block Sizes
|
||||
# 19: Q starts
|
||||
# 20: T starts
|
||||
# 21: Q seqs (pslx format)
|
||||
# 22: T seqs (pslx format)
|
||||
|
||||
## All sequences start at 0 here; array-based.
|
||||
|
||||
my $JOIN_GAP = 9; #join alignment segments if within this gap along the genomic sequence.
|
||||
|
||||
my $chain_number = 0;
|
||||
|
||||
while (<STDIN>) {
|
||||
unless (/\w/) { next; }
|
||||
chomp;
|
||||
my @x = split (/\t/);
|
||||
unless ($x[0] =~ /^\d/) {next;} #eliminate headers if present.
|
||||
my @alignment_segments;
|
||||
$chain_number++;
|
||||
my $strand = $x[8];
|
||||
my $cdna_name = $x[9];
|
||||
my $genomic_name = $x[13];
|
||||
my $genomic_length = $x[14];
|
||||
my $cdna_length = $x[10];
|
||||
my $cDNA_seqs = $x[21];
|
||||
my $genomic_seqs = $x[22];
|
||||
my $num_segs = $x[17];
|
||||
|
||||
my @cdna_coords = split (/,/, $x[19]);
|
||||
my @genomic_coords = split (/,/, $x[20]);
|
||||
my @lengths = split (/,/, $x[18]);
|
||||
my @cDNA_seqs = split (/,/, $x[21]);
|
||||
my @genomic_seqs = split (/,/, $x[22]);
|
||||
|
||||
## going to implement score as chain_score + segment score
|
||||
## Chain score = matches_num - mismatch_num
|
||||
my $chain_score = $x[0] - $x[1];
|
||||
|
||||
## report each segment match as a separate btab entry:
|
||||
my $segment_number = 0;
|
||||
for (my $i = 0; $i < $num_segs; $i++) {
|
||||
$segment_number++;
|
||||
my $length = $lengths[$i];
|
||||
my $cdna_coord = $cdna_coords[$i];
|
||||
my $genomic_coord = $genomic_coords[$i];
|
||||
my ($cdna_end5, $cdna_end3) = (++$cdna_coord, $cdna_coord + $length - 1);
|
||||
my ($genomic_end5, $genomic_end3) = (++$genomic_coord, $genomic_coord + $length -1);
|
||||
if ($strand eq "-") {
|
||||
($genomic_end5, $genomic_end3) = ($genomic_end3, $genomic_end5);
|
||||
($cdna_end5, $cdna_end3) = sort {$a<=>$b} ($cdna_length - $cdna_end5 + 1, $cdna_length - $cdna_end3 + 1);
|
||||
}
|
||||
my $segment_score = $chain_score + $length;
|
||||
my $per_id = -1;
|
||||
my ($gseq, $cseq);
|
||||
if ( ($gseq = $genomic_seqs[$i]) && ($cseq = $cDNA_seqs[$i])) {
|
||||
if ($gseq eq $cseq) {
|
||||
$per_id = 100;
|
||||
} else {
|
||||
## walk thru and determine num ids
|
||||
my $num_id = 0;
|
||||
my @gseq_array = split (//, $gseq);
|
||||
my @cseq_array = split (//, $cseq);
|
||||
for (my $j = 0; $j < $length; $j++) {
|
||||
if ($gseq_array[$j] eq $cseq_array[$j]) {
|
||||
$num_id++;
|
||||
}
|
||||
}
|
||||
$per_id = ($num_id/$length) * 100;
|
||||
}
|
||||
}
|
||||
|
||||
## Create btab line.
|
||||
my @btab;
|
||||
$btab[0] = $genomic_name;
|
||||
$btab[2] = $genomic_length;
|
||||
$btab[3] = "blat";
|
||||
$btab[5] = $cdna_name;
|
||||
$btab[6] = $genomic_end5;
|
||||
$btab[7] = $genomic_end3;
|
||||
$btab[8] = $cdna_end5;
|
||||
$btab[9] = $cdna_end3;
|
||||
|
||||
$btab[10] = $per_id;
|
||||
$btab[12] = $segment_score;
|
||||
|
||||
$btab[13] = $chain_number;
|
||||
$btab[14] = $segment_number;
|
||||
|
||||
$btab[18] = $length;
|
||||
push (@alignment_segments, [@btab]);
|
||||
|
||||
}
|
||||
|
||||
&process_alignment_chain(\@alignment_segments, $chain_number);
|
||||
}
|
||||
|
||||
exit(0);
|
||||
|
||||
|
||||
## Join alignment segments if within 5 bp along the genomic sequence
|
||||
sub process_alignment_chain {
|
||||
my $alignment_segments_aref = shift;
|
||||
my $chain_number = shift;
|
||||
my @segments = sort {$a->[8]<=>$b->[8]} @$alignment_segments_aref;
|
||||
if ($#segments > 0) {
|
||||
my @new_segments = ($segments[0]); # always holds the last segment analyzed.
|
||||
for (my $i=1; $i <= $#segments; $i++) {
|
||||
my $last_segment = $new_segments[$#new_segments];
|
||||
my $current_segment = $segments[$i];
|
||||
my $prev_end3 = $last_segment->[7];
|
||||
my $curr_end5 = $current_segment->[6];
|
||||
my $gap_length = abs ($prev_end3 - $curr_end5) - 1;
|
||||
|
||||
if ($gap_length <= $JOIN_GAP) {
|
||||
## Must join prev and current segment
|
||||
my $prev_seg_length = abs ($last_segment->[7] - $last_segment->[6]) + 1;
|
||||
my $curr_seg_length = abs ($current_segment->[7] - $current_segment->[6]) + 1;
|
||||
|
||||
## make prev end3 the new end3 for both genomic and cdna coordinates
|
||||
$last_segment->[7] = $current_segment->[7];
|
||||
$last_segment->[9] = $current_segment->[9];
|
||||
|
||||
## recalculate the percent ID
|
||||
my $prev_per_id = $last_segment->[10];
|
||||
my $curr_per_id = $current_segment->[10];
|
||||
|
||||
my $new_per_id = ($prev_seg_length * $prev_per_id + $curr_seg_length * $curr_per_id) /
|
||||
($prev_seg_length + $curr_seg_length + $gap_length);
|
||||
|
||||
$last_segment->[10] = $new_per_id;
|
||||
$last_segment->[18] = abs($last_segment->[9] - $last_segment->[8]) + 1;
|
||||
$last_segment->[12] += $curr_seg_length;
|
||||
|
||||
} else {
|
||||
|
||||
#make the current segment the last segment
|
||||
push (@new_segments, $current_segment);
|
||||
|
||||
}
|
||||
}
|
||||
@segments = @new_segments;
|
||||
|
||||
## renumber the segment numbers
|
||||
my $segnum = 0;
|
||||
foreach my $segment (@segments) {
|
||||
$segnum++;
|
||||
$segment->[14] = $segnum;
|
||||
}
|
||||
}
|
||||
|
||||
## write GFF format.
|
||||
my $orient;
|
||||
foreach my $segment (@segments) {
|
||||
|
||||
my $genomic_end5 = $segment->[6];
|
||||
my $genomic_end3 = $segment->[7];
|
||||
|
||||
unless ($orient) {
|
||||
if ($genomic_end5 < $genomic_end3) {
|
||||
$orient = '+';
|
||||
}
|
||||
elsif ($genomic_end3 < $genomic_end5) {
|
||||
$orient = '-';
|
||||
}
|
||||
}
|
||||
if ($orient) {
|
||||
last;
|
||||
}
|
||||
}
|
||||
|
||||
unless ($orient) {
|
||||
$orient = '+'; ## set a default
|
||||
}
|
||||
|
||||
|
||||
foreach my $segment (@segments) {
|
||||
|
||||
my $genomic_contig = $segment->[0];
|
||||
my $cdna_name = $segment->[5];
|
||||
|
||||
my $genomic_end5 = $segment->[6];
|
||||
my $genomic_end3 = $segment->[7];
|
||||
|
||||
my ($genomic_lend, $genomic_rend) = sort {$a<=>$b} ($genomic_end3, $genomic_end5);
|
||||
|
||||
my $cdna_end5 = $segment->[8];
|
||||
my $cdna_end3 = $segment->[9];
|
||||
|
||||
my $per_id = sprintf("%.2f", $segment->[10]);
|
||||
|
||||
my $chain_ID = "blat.proc$$.chain_" . $chain_number;
|
||||
|
||||
print join("\t", $genomic_contig, "BLAT", "cDNA_match",
|
||||
$genomic_lend, $genomic_rend, $per_id, $orient, ".",
|
||||
"ID=$chain_ID;Target=$cdna_name $cdna_end5 $cdna_end3 +") . "\n";
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,36 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
|
||||
my $usage = "usage: $0 genome_db transcript_db maxIntron outFile\n\n";
|
||||
|
||||
my $genome_db = $ARGV[0] or die $usage;
|
||||
my $transcript_db = $ARGV[1] or die $usage;
|
||||
my $max_intron = $ARGV[2] or die $usage;
|
||||
my $outFile = $ARGV[3] or die $usage;
|
||||
|
||||
main: {
|
||||
|
||||
my $ooc_cmd = "blat -t=dna -q=rna -maxIntron=$max_intron -makeOoc=11.ooc $genome_db $transcript_db $outFile";
|
||||
unless (-s "11.ooc") {
|
||||
&process_cmd($ooc_cmd);
|
||||
}
|
||||
|
||||
|
||||
my $blat_cmd = "blat -t=dna -q=rna -maxIntron=$max_intron -ooc=11.ooc $genome_db $transcript_db $outFile";
|
||||
&process_cmd($blat_cmd);
|
||||
|
||||
exit(0);
|
||||
}
|
||||
|
||||
sub process_cmd {
|
||||
my ($cmd) = @_;
|
||||
|
||||
my $ret = system($cmd);
|
||||
if ($ret) {
|
||||
die "Error, $cmd died with ret $ret";
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
@@ -0,0 +1,132 @@
|
||||
#!/usr/bin/env perl
|
||||
|
||||
use strict;
|
||||
use warnings;
|
||||
|
||||
use FindBin;
|
||||
use lib ("$FindBin::RealBin/../../../PerlLib");
|
||||
|
||||
use SAM_reader;
|
||||
use SAM_entry;
|
||||
|
||||
|
||||
my $usage = "usage: $0 blat.nameSorted.sam [num_top_hits=20] [min_per_ID=0]\n\n";
|
||||
|
||||
# note min perID is based on read length and blat2sam.pl default scoring with -q 0 -r 0
|
||||
|
||||
my $blat_sam = $ARGV[0] or die $usage;
|
||||
my $num_top_hits = $ARGV[1] || 20;
|
||||
my $min_per_ID = $ARGV[2] || 0;
|
||||
|
||||
main: {
|
||||
|
||||
my $sam_reader = new SAM_reader($blat_sam);
|
||||
|
||||
while ($sam_reader->has_next()) {
|
||||
|
||||
my @entries;
|
||||
|
||||
my $sam_entry = $sam_reader->get_next();
|
||||
push (@entries, $sam_entry);
|
||||
|
||||
while ($sam_reader->has_next()
|
||||
&&
|
||||
$sam_reader->preview_next()->get_read_name() eq $sam_entry->get_read_name()) {
|
||||
|
||||
push (@entries, $sam_reader->get_next());
|
||||
}
|
||||
|
||||
&report_top_hits(@entries);
|
||||
}
|
||||
|
||||
exit(0);
|
||||
|
||||
}
|
||||
|
||||
|
||||
####
|
||||
sub report_top_hits {
|
||||
my @entries = @_;
|
||||
|
||||
@entries = &get_top_hits(@entries);
|
||||
|
||||
|
||||
foreach my $entry (@entries) {
|
||||
print $entry->toString() . "\n";
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
|
||||
####
|
||||
sub get_top_hits {
|
||||
my @entries = @_;
|
||||
|
||||
my @structs;
|
||||
|
||||
foreach my $entry (@entries) {
|
||||
my @fields = $entry->get_fields();
|
||||
|
||||
my $seq_length = length($entry->get_sequence());
|
||||
|
||||
my $min_score = $seq_length - ( ( 1 - ($min_per_ID / 100)) * $seq_length * 3); # perfect_score - (num_mismatches * mismatch_penalty)
|
||||
|
||||
my $alignment_field = $fields[11];
|
||||
|
||||
$alignment_field =~ /AS:i:(\d+)/ or die "Error, no score reported for " . $entry->toString();
|
||||
|
||||
my $score = $1;
|
||||
unless (defined $score) {
|
||||
die "Error, no score for " . $entry->toString();
|
||||
}
|
||||
|
||||
#print "MIN_score: $min_score vs. score: $score\n";
|
||||
|
||||
if ($score >= $min_score) {
|
||||
|
||||
push (@structs, { entry => $entry,
|
||||
score => $score,
|
||||
}
|
||||
);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@structs = reverse sort {$a->{score}<=>$b->{score}} @structs;
|
||||
|
||||
|
||||
|
||||
if ($num_top_hits < 0 && scalar(@structs) > 1) {
|
||||
# multiply mapped reads
|
||||
# ignoring entry.
|
||||
return();
|
||||
}
|
||||
|
||||
elsif ($num_top_hits > 0 && scalar (@structs) > $num_top_hits) {
|
||||
@structs = @structs[0..$num_top_hits-1];
|
||||
}
|
||||
|
||||
## unwrap:
|
||||
my @ret;
|
||||
if (@structs) {
|
||||
my $top_struct = shift @structs;
|
||||
|
||||
my $top_score = $top_struct->{score};
|
||||
|
||||
push (@ret, $top_struct->{entry});
|
||||
|
||||
foreach my $struct (@structs) {
|
||||
if ($struct->{score} == $top_score) {
|
||||
push (@ret, $struct->{entry});
|
||||
}
|
||||
else {
|
||||
last;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return(@ret);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user