perl novel可變剪接識別(1)
阿新 • • 發佈:2019-01-28
想把之前做的可變剪接模型給大家說一下,看看有什麼遺漏的沒有,由於當時想法比較複雜,所以程式有點多,大致分三個部分來進行。
首先,拿到的結果是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_s\t$block_e\n"; 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_s\t$block_e\n"; } } close JUNC; &showtime("Junctions is done..."); sub showtime(){ my ($str) = @_; my $time = localtime; print STDERR "$str\t$time\n"; }
這樣就得到兩部分資料,已知與新的——新的將進行接下來的處理工作。