perl novel可变剪接识别（1）

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

首先，拿到的结果是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";
}

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

perl novel可变剪接识别（1）相关推荐

正则表达式matlab,正则表达式中一个word的匹配 @MATLAB - 优秀的Free OS(Linux)版 - 北大未名BBS...
我目前想做的就是判断一个str是否可以被认为是有效的MATLAB index. 最好的方法是直接运行,然后看运行结果或报错类型,但是我不打算在不知道是什么类型的东西之前运行它,所以可以预先parse ...
moses binarize-all问题
要想在EMS(Experiment Management System)中运用Word Lattice进行调参或解码,就不得不将配置文件中binarize-all开关打开,并且指定ttable-bin ...
Bioinformatics：吉林大学刘富组-深度学习从宏基因组序列中识别短病毒序列Virtifier...
Virtifier: 基于深度学习方法的宏基因组病毒序列分类方法 Virtifier: a deep learning-based identifier for viral sequences fro ...
AEM：众里寻他千百度，用了SIP不迷路-识别污水处理单元中砷氧化细菌的多样性及代谢潜力...
活性污泥中 As(III)-氧化细菌的多样性和代谢潜能 Diversity and metabolic potentials of As(III)-oxidizing bacteria in acti ...
NAR：浙大陈云/马忠华团队揭示病原真菌组蛋白H3K27甲基化识别新机制
2021年9月27日,Nucleic Acids Research在线发表了浙江大学农学院马忠华团队题为"Fusarium BP1 is a reader of H3K27 methylat ...
JGG：口腔微生物组助力妊娠糖尿病识别
2021年1月20日,JGG 在线发表了中国科学院北京生命科学研究院王金锋副研究员与温州市人民医院张红萍主任医师合作的题为"The Oral Microbiome of Pregnant W ...
BMC：幼年特发性关节炎患儿肠道菌群的特征、生物标记的识别及其在临床预测中的作用...
论文题目:Gut microbiota in children with juvenile idiopathic arthritis: characteristics, biomarker ident ...
PhiSpy：在细菌基因组中识别噬菌体
PhiSpy:在细菌基因组中识别噬菌体 PhiSpy: a novel algorithm for finding prophages in bacterial genomes that combin ...
BiLSTM+CRF的损失由发射矩阵和转移矩阵计算而得 BiLSTM+CRF命名实体识别：达观杯败走记（下篇
如果是训练,那么直接用发射矩阵和真实标签去计算Loss,用于更新梯度. 这需要用到CRF中的forward函数. 如果是预测,那么就用发射矩阵去进行维特比解码,得到最优路径(预测的标签). 这需要用到 ...
下载perl的mysql模块_安装用于操作MySQL的Perl模块
在我使用Webmin(version 1.480)管理FreeBSD主机上的MySQL数据库服务器时出现: "警告:您的系统未安装Perl 模块 DBI 和 DBD::mysql,Webmi ...

perl novel可变剪接识别（1）

perl novel可变剪接识别（1）相关推荐

最新文章

热门文章