perl novel可变剪接识别(1)
发布时间:2020-12-16 00:15:11 所属栏目:大数据 来源:网络整理
导读:想把之前做的可变剪接模型给大家说一下,看看有什么遗漏的没有,由于当时想法比较复杂,所以程序有点多,大致分三个部分来进行。 首先,拿到的结果是tophat给出的junction的数据,其次博主使用的数据库是ensembl的数据库,gencode也可以,先得到已知的参考ju
想把之前做的可变剪接模型给大家说一下,看看有什么遗漏的没有,由于当时想法比较复杂,所以程序有点多,大致分三个部分来进行。 首先,拿到的结果是tophat给出的junction的数据,其次博主使用的数据库是ensembl的数据库,gencode也可以,先得到已知的参考junction: #!/usr/bin/env perl use warnings; use strict; die "perl $0 <junction.bed> <ensembl|gencode> <bp>n" unless @ARGV eq 3; &showtime("Start..."); my $in = shift; my %gene; open GTF,$in or die $!; while(<GTF>){ chomp; my @tmp = split; next if(/^#/); my ($gid) = $_ =~ /gene_name "([^;]+)";/; my ($tid) = $_ =~ /transcript_id "([^;]+)";/; if($tmp[2] =~ /exon/){ push @{$gene{$gid}{$tid}{$tmp[0]}},"$tmp[3]t$tmp[4]"; } } close GTF; my %hash; foreach my $g(keys %gene){ foreach my $t(keys %{$gene{$g}}){ foreach my $chr(keys %{$gene{$g}{$t}}){ my $flag = 0; my $end; foreach my $str(sort @{$gene{$g}{$t}{$chr}}){ if($flag == 0){ $end = (split /t/,$str)[1]; $flag = 1; }else{ my $f_end = (split /t/,$str)[0]; $hash{$chr}{$end}{$f_end} = 0; $end = (split /t/,$str)[1]; } } } } } &showtime("GTF is done..."); #上面两步将所有的EXON连接方式给存起来 my $junc = shift; my ($id) = $junc =~ /^(w+)_alignment/; open OUT1,">$id.novel" or die $!; open OUT2,">$id.ref" or die $!; my $bp = shift; #设置精确范围,可以在参考junction的范围内都识别成已知的剪接,博主设置为0为了精确 $bp ||= 0; open JUNC,$junc or die $!; while(<JUNC>){ chomp; next if($. == 1); my @tmp = split; my $flag = 0; my ($block_s,$block_e) = (split /,/,$tmp[10])[0,1]; my $intron_s = $tmp[1] + $block_s; my $intron_e = $tmp[2] - $block_e + 1; foreach my $s(keys %{$hash{$tmp[0]}}){ foreach my $e(keys %{$hash{$tmp[0]}{$s}}){ if(abs($intron_s - $s) <= $bp and abs($intron_e - $e) <= $bp){ $flag = 1; print OUT2 "$tmp[0]t$tmp[1]t$tmp[2]t$tmp[4]t$tmp[5]t$block_st$block_en"; delete $hash{$tmp[0]}{$s}{$e}; last; } } } if($flag eq 1) { next; }else{ print OUT1 "$tmp[0]t$tmp[1]t$tmp[2]t$tmp[4]t$tmp[5]t$block_st$block_en"; } } close JUNC; &showtime("Junctions is done..."); sub showtime(){ my ($str) = @_; my $time = localtime; print STDERR "$strt$timen"; } 这样就得到两部分数据,已知与新的——新的将进行接下来的处理工作。 (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |