https://github.com/qqwang-berkeley/JUM
Tip revision: afeaab37960ac19b78511e8909b5e630ec355751 authored by Qingqing Wang on 10 April 2023, 14:50:26 UTC
Add files via upload
Add files via upload
Tip revision: afeaab3
count_reads_for_intron_retention_sam.pl
#!/usr/bin/perl
use strict;
use warnings;
#data_file_1=output_long_intron.gff; data_file_2=XXXXXXintersect_long_intron.sam; readlength= #;
my $len=scalar(@ARGV);
if ($len < 1) {
die "Usage:count_reads_for_intron_retention_sam.pl <data_file_1> <data_file_2> <readlength>\n";
}
my ($data_file_1, $data_file_2, $readlength) = @ARGV;
my %hash; #$hash{chr}=[start1, start2, start3...];
open(IN1, "<$data_file_1") or die "can't open input1 file: $!";
open(IN2, "<$data_file_2") or die "can't open input2 file: $!";
#while(<IN1>) {
# chomp;
# my @long=split(/\s+/,$_);
# my @longid=split(/=/,$long[8]);
# $hash{$long[0]}{$longid[1]}[0]=$long[3];
# $hash{$long[0]}{$longid[1]}[1]=$long[4]-$readlength;
# $hash{$long[0]}{$longid[1]}[2]=0;
#}
#close IN1;
while(<IN2>) {
chomp;
my @array=split(/\s+/,$_);
if (exists $hash{$array[0]}) {
push @{ $hash{$array[0]} }, $array[1];
}
else {
$hash{$array[0]}[0]=$array[1];
}
}
close IN2;
while(<IN1>) {
chomp;
my @long=split(/\s+/,$_);
my @longid=split(/=/,$long[8]);
if(exists $hash{$long[0]}) {
# my @sorted = sort @ { $hash{$long[0]} };
my @count = grep {($_ >= $long[3]) && ($_ <= $long[4]-$readlength)} @{ $hash{$long[0]} };
my $num = scalar @count;
print $longid[1]; print "\t"; print $num; print "\n";
}
else {
print $longid[1]; print "\t"; print "0"; print "\n";
}
}
close IN1;