perl novel可变剪接识别（1）

概述想把之前做的可变剪接模型给大家说一下，看看有什么遗漏的没有，由于当时想法比较复杂，所以程序有点多，大致分三个部分来进行。 首先，拿到的结果是tophat给出的junction的数据，其次博主使用的数据库是ensembl的数据库，gencode也可以，先得到已知的参考junction： #!/usr/bin/env perluse warnings;use strict;die "perl

想把之前做的可变剪接模型给大家说一下，看看有什么遗漏的没有，由于当时想法比较复杂，所以程序有点多，大致分三个部分来进行。

首先，拿到的结果是tophat给出的junction的数据，其次博主使用的数据库是ensembl的数据库，gencode也可以，先得到已知的参考junction：

#!/usr/bin/env perluse warnings;use strict;dIE "perl  <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";}

这样就得到两部分数据，已知与新的——新的将进行接下来的处理工作。

总结

以上是内存溢出为你收集整理的perl novel可变剪接识别（1）全部内容，希望文章能够帮你解决perl novel可变剪接识别（1）所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。