perl novel可变剪接识别（3）

好吧，该到最后的时刻了，识别novel可变剪接的类型，首先要在原理上判断类型的情况，把类型画出来比较好一些：

以上即为大致的分类，不过有两点要说明的：intron retained 1用tophat的结果无法实现，这个看覆盖度即可；另外还有两个分类：intergenic和other，顾名思义基因间区域内的和无法识别的，当然也可以按照其他分类，就不多说了，就以此为例吧：

#!/usr/bin/env perl
use warnings;
use strict;my %hash;
open GTF, $ARGV[0] or die $!;
while(<GTF>)
{chomp;my @tmp = split;my ($a, $b) = (split /-/, $tmp[2])[0,1];$hash{$tmp[0]}{$tmp[1]}{'pn'} = $tmp[3];$hash{$tmp[0]}{$tmp[1]}{'rs'} = $a;$hash{$tmp[0]}{$tmp[1]}{'re'} = $b;$hash{$tmp[0]}{$tmp[1]}{'type'} = $tmp[4];my @a = split /,/, $tmp[5];my @b = split /,/, $tmp[6];@{$hash{$tmp[0]}{$tmp[1]}{'start'}} = @a;@{$hash{$tmp[0]}{$tmp[1]}{'end'}} = @b;
}
close GTF;open NOVEL, $ARGV[1] or die $!;
print "#chr\tsplice_start\tsplice_end\tsplice_strand\tsplice_reads\tsplice_type\tgene_id\tgene_strand\tgene_type\n";
while(<NOVEL>)
{chomp;my @tmp = split;my $js = $tmp[1] + $tmp[5];my $je = $tmp[2] - $tmp[6] + 1;my ($chr, $read, $jpn) = @tmp[0,3,4];my $flag = 0;foreach my $t(keys %{$hash{$chr}}){if($je < $hash{$chr}{$t}{'rs'} or $js > $hash{$chr}{$t}{'re'}){next;}else{$flag = 1;if($js < $hash{$chr}{$t}{'rs'} and $je >= $hash{$chr}{$t}{'rs'}){if($hash{$chr}{$t}{'pn'} eq '+'){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t+\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t5UTR\t$t\t+\tCoding\n";}}else{if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t-\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t3UTR\t$t\t-\tCoding\n";}}last;}elsif($js <= $hash{$chr}{$t}{'re'} and $je > $hash{$chr}{$t}{'re'}){if($hash{$chr}{$t}{'pn'} eq '+'){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t+\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t3UTR\t$t\t+\tCoding\n";}}else{if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t-\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t5UTR\t$t\t-\tCoding\n";}}last;}else{my ($js_stat, $je_stat) = (0, 0);foreach my $i(@{$hash{$chr}{$t}{'end'}}){$js_stat = 1 if($js eq $i);}foreach my $j(@{$hash{$chr}{$t}{'start'}}){$je_stat = 1 if($je eq $j);}if($js_stat eq 1 and $je_stat eq 1){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tES\t$t\t$hash{$chr}{$t}{'pn'}\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\tES\t$t\t$hash{$chr}{$t}{'pn'}\tCoding\n";}}elsif($js_stat eq 1){if($hash{$chr}{$t}{'pn'} eq '+'){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\t3S\t$t\t+\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t3S\t$t\t+\tCoding\n";}}else{if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\t5S\t$t\t-\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t5S\t$t\t-\tCoding\n";}}}elsif($je_stat eq 1){if($hash{$chr}{$t}{'pn'} eq '+'){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\t5S\t$t\t+\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t5S\t$t\t+\tCoding\n";}}else{if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\t3S\t$t\t-\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\t3S\t$t\t-\tCoding\n";}}}else{my $f2 = 0;for(my $i = 0; $i < @{$hash{$chr}{$t}{'end'}}; $i ++){if($js >= $hash{$chr}{$t}{'start'}->[$i] and $je <= $hash{$chr}{$t}{'end'}->[$i]){$f2 = 1;if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tIR\t$t\t$hash{$chr}{$t}{'pn'}\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\tIR\t$t\t$hash{$chr}{$t}{'pn'}\tCoding\n";}last;}}if($f2 eq 0){if($hash{$chr}{$t}{'type'} eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t$hash{$chr}{$t}{'pn'}\tNonCoding\n";}else{print "$chr\t$js\t$je\t$jpn\t$read\tOther\t$t\t$hash{$chr}{$t}{'pn'}\tCoding\n";}}}last;}}}if($flag eq 0){print "$chr\t$js\t$je\t$jpn\t$read\tIntergenic\t.\t.\t.\n";}
}
close NOVEL;

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

正则表达式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可变剪接识别（3）

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

最新文章

热门文章