登录
当前位置： 首页 > 工具软件 > Novel > 使用案例 >

perl novel可变剪接识别（2）

公西翼
2023-12-01

博主其实对未知的可变剪接分类有些困惑，但想了很久才使用了一项比较复杂的算法，接下来并不是分类，而是先转换数据库，因为ensembl|gencode数据库的格式并不能满足博主的需求，转换后更能方便地处理接下来的工作：

#!/usr/bin/env perl
use warnings;
use strict;

my (%gene);
open GTF, $ARGV[0] or die $!;
while(<GTF>)
{
	chomp;
	next if(/^#/);
	my @tmp = split;
	my ($gid) = $_ =~ /gene_id "([^;]+)";/; #匹配基因名称
	if($_ =~ /gene_type "protein_coding";/) #分coding和noncoding分别处理，因为noncoding没有UTR结构，不存在UTR的可变剪接
	{
		if($tmp[2] =~ /exon/)
		{
			push @{$gene{$tmp[0]}{$gid}{$tmp[6]}}, "$tmp[3],$tmp[4],1";
		}
	}else{
		if($tmp[2] =~ /exon/)
		{
			push @{$gene{$tmp[0]}{$gid}{$tmp[6]}}, "$tmp[3],$tmp[4],0";
		}
	}
}
close GTF;

open OUT, ">$ARGV[0].myformat" or die $!;
foreach my $c(keys %gene)
{
	foreach my $g(keys %{$gene{$c}})
	{
			foreach my $ps(keys %{$gene{$c}{$g}})
			{
				if($ps eq '+')
				{
					my ($s, $e, $type);
					if(@{$gene{$c}{$g}{$ps}} eq 1)
					{
						($s,$e, $type) = (split /,/, $gene{$c}{$g}{$ps}->[0])[0,1,2];
					}else{
						$s = (split /,/, $gene{$c}{$g}{$ps}->[0])[0];
						$e = (split /,/, $gene{$c}{$g}{$ps}->[-1])[1];
						$type = (split /,/, $gene{$c}{$g}{$ps}->[0])[2];
					}
					my (@start, @end);
					foreach my $loci(@{$gene{$c}{$g}{$ps}})
					{
						push @start, (split /,/, $loci)[0];
						push @end, (split /,/, $loci)[1];
					}
					my $ss = join ",", @start;
					my $ee = join ",", @end;
					print OUT join "\t", $c, $g, "$s-$e", $ps, $type, $ss, $ee, "\n";
				}else{
					my @ps_a = reverse @{$gene{$c}{$g}{$ps}};
					my ($s, $e, $type);
					if(@ps_a eq 1)
					{
						($s,$e,$type) = (split /,/, $ps_a[0])[0,1,2];
					}else{
						$s = (split /,/, $ps_a[0])[0];
						$e = (split /,/, $ps_a[-1])[1];
						$type = (split /,/, $ps_a[0])[2];
					}
					my (@start, @end);
					foreach my $loci(@ps_a)
					{
						push @start, (split /,/, $loci)[0];
						push @end, (split /,/, $loci)[1];
					}
					my $ss = join ",", @start;
					my $ee = join ",", @end;
					print OUT join "\t", $c, $g, "$s-$e", $ps, $type, $ss, $ee, "\n";
				}
			}
		
	}
}

这个格式类似ucsc的refseq格式，当然只是类似而已，而博主所需求的是gene类型和exon起始与终止而已。

说到数据库的转换，博主还想起来一件事情，就是refseq转gff格式，其实两者差的有些多。

有个妹子写了一个转换的，还不错，展示的内容非常的详细想要啥结果都可以，当然在这之中博主也贡献了一点功劳，嘿嘿！就拿来在这展示一下吧：

#!perl -w
use strict; 
die "Usage : perl $0 <in.refGene.lst> <out.gff>" unless (@ARGV == 2);
my ($in, $out) = @ARGV;

my($pre, $insert, @inserts, $utr, $utr_o, $i, $j, $cds, $nm, $chr, $direction, $start_exon, $end_exon, $start_cds, $end_cds, $cds_num, $start, $end, $tmp, $gene, @starts, @ends, $up, $down);

if ($in =~ /\.gz/){open IN, " gzip -dc $in | " || die $!;}
else{open IN, $in || die $!;}
if ($out =~ /\.gz/){open OUT, "| gzip > $out" || die $!;}
else {open OUT , "> $out" || die $!;}

while (<IN>){
		chomp;
		($nm, $chr, $direction, $start_exon, $end_exon, $start_cds, $end_cds, $cds_num, $start, $end, $tmp, $gene, $insert) = (split)[1..12,15];
		if ($nm =~ /NM/){
				$pre = 'mRNA';
		}else{
				$pre = 'ncRNA';
		}
		$start =~ s/,$//;
		$end =~ s/,$//;
		$insert =~ s/,$//;
		$insert =~ s/\-1/\./g;
		if ($cds_num > 1){
				@starts = split /,/, $start;
				@ends = split /,/, $end;
				@inserts = split /,/, $insert;
		}else{
				@starts = ($start);
				@ends = ($end);
				@inserts = ($insert);
		}
		print OUT "$chr\trefGene\t$pre\t$start_exon\t$end_exon\t.\t$direction\t.\tID=$nm; name=$gene;\n";
		if ($direction eq '+'){
				$utr = 5;
				$utr_o = 3;
				#print OUT "$chr\trefGene\t5-UTR\t$start_exon\t",$start_cds-1,"\t.\t$direction\t.\tParent=$nm;\n";
		} else {
				$utr = 3;
				$utr_o = 5;
				#print OUT "$chr\trefGene\t3-UTR\t$start_exon\t",$start_cds-1,"\t.\t$direction\t.\tParent=$nm;\n";
		}
		for ($i = 0; $i < @starts; $i ++){
				print OUT "$chr\trefGene\tintron\t",$ends[$i-1]+1,"\t",$starts[$i]-1,"\t.\t$direction\t.\tParent=$nm;\n" if ($i > 0);
				if ($pre eq 'ncRNA'){
						print OUT "$chr\trefGene\tCDS\t$starts[$i]\t$ends[$i]\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
						next;
				}
				if ($ends[$i] < $start_cds){
						print OUT "$chr\trefGene\t$utr-UTR\t$starts[$i]\t$ends[$i]\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
				}elsif($starts[$i] < $start_cds and $ends[$i] > $start_cds){
						print OUT "$chr\trefGene\t$utr-UTR\t$starts[$i]\t", $start_cds - 1, "\t.\t$direction\t.\tParent=$nm;\n";
						print OUT "$chr\trefGene\tCDS\t$start_cds\t$ends[$i]\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
				}elsif($starts[$i] >= $start_cds and $ends[$i] <= $end_cds){
						print OUT "$chr\trefGene\tCDS\t$starts[$i]\t$ends[$i]\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
				}elsif($starts[$i] < $end_cds and $ends[$i] > $end_cds){
						print OUT "$chr\trefGene\tCDS\t$starts[$i]\t$end_cds\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
						print OUT "$chr\trefGene\t$utr_o-UTR\t", $end_cds + 1, "\t$ends[$i]\t.\t$direction\t.\tParent=$nm;\n";
				}else{
						print OUT "$chr\trefGene\t$utr_o-UTR\t$starts[$i]\t$ends[$i]\t.\t$direction\t$inserts[$i]\tParent=$nm;\n";
				}
		}
}
close IN;
close OUT;

上面是正常的gff格式，下面将添加一些别的东西：

#!perl -w
use strict;
die "Usage : perl $0 <in.file> <out.file> " if (@ARGV > 2);
my ($in, $out) = @ARGV;

$in ||= 'refGene.sort.2.gz';
$out ||= 'refGene.ann.gz';

if($in =~ /\.gz/){
		open IN, "gzip -dc $in |" || die $!;
}else{
		open IN, $in || die $!;
}

if ($out =~ /\.gz/){
		open OUT, "| gzip > $out " || die $!;
}else{
		open OUT, "> $out" || die $!;
}


my (@tmps, $chr, $ref, $region, $start, $end, $num, $dot, $id, $count_cds, $count_intron, $gene, $nm);
my (%infos, %genes);
my @chrs = (1..22, 'X', 'Y');

$/ = "\n>";
chomp (my $line = <IN>);
@tmps = split /\n/, $line;

my ($chr_b, $end_b, $orientation, $id_b) = $tmps[0] =~ /(chr\w+)\t\w+\t\w+RNA\t\d+\t(\d+)\t\.\t([\+\-])\t\.\t(ID=\S+; name=\S+)$/;
$tmps[0] =~ s/>//;
print OUT "$tmps[0]\n";
@tmps[1..$#tmps] = &add_num ($orientation, @tmps[1..$#tmps]);
print OUT join "\n", @tmps[1..$#tmps];
print OUT "\n";

while (<IN>){
		chomp;
		@tmps = split /\n/, $_;
		($chr, $start, $end, $orientation, $id) = $tmps[0] =~ /(chr\w+)\t\w+\t\w+RNA\t(\d+)\t(\d+)\t\.\t([\+\-])\t\.\t(ID=\S+; name=\S+)$/;
		if ($chr eq $chr_b and $end_b < $start){
				print OUT "$chr\trefGene\tintergenic\t",$end_b + 1, "\t", $start - 1, "\t.\t.\t.\t$id_b|$id\n";
		}
		($chr_b, $end_b ,$id_b) = ($chr, $end, $id);
		$tmps[0] =~ s/>//;
		print OUT "$tmps[0]\n";
		@tmps[1..$#tmps] = &add_num ($orientation, @tmps[1..$#tmps]);
		print OUT join "\n", @tmps[1..$#tmps];
		print OUT "\n";
}
close IN;
close OUT;
#$/ = '\n';


sub add_num {
		my @tmps = @_;
		my ($count_cds, $count_intron) = (0, 0);
		if ($tmps[0] eq '+'){
				for my $tmp (@tmps[1..$#tmps]){
						if ($tmp =~ /CDS/){
								$count_cds ++;
								$tmp =~ s/\.\t\+/$count_cds\t\+/;
								#$count_cds ++;
						}elsif($tmp =~ /intron/){
								$count_intron ++;
								$tmp =~ s/\.\t\+/$count_intron\t\+/;
								#$count_intron ++;
						}
				}
				return @tmps[1..$#tmps];
		}
		#return @tmps[1..$#tmps];
		for my $tmp (@tmps[1..$#tmps]){
			if ($tmp =~ /CDS/){
					$count_cds ++;
			}elsif($tmp =~ /intron/){
					$count_intron ++;
			}
		}
		for my $tmp (@tmps[1..$#tmps]){
				if ($tmp =~ /CDS/){
						$tmp =~ s/\.\t\-/$count_cds\t\-/;
						$count_cds --;
				}elsif($tmp =~ /intron/){
						$tmp =~ s/\.\t\-/$count_intron\t\-/;
						$count_intron --;
				}
		}
		return @tmps[1..$#tmps];
}

脚本如下：

less refGene.txt.gz | sort -k3,3 -k5,5n |gzip > refGene.txt.sort.gz
perl format.pl refGene.txt.sort.gz refGene.gff.gz
less refGene.gff.gz  |perl -lane 'if($F[2] =~ /RNA/){$F[0] = ">$F[0]";}print join "\t",@F[0..7], "$F[8] $F[9]";'  |gzip > refGene.gff.2.gz
perl ann.sort.pl refGene.gff.2.gz refGene.ann.gz

最终的转换格式如下：

chr1	refGene	intron	763230	764381	1	+	.	Parent=NR_047525; 
chr1	refGene	CDS	764382	764484	2	+	.	Parent=NR_047525; 
chr1	refGene	intron	764485	787305	2	+	.	Parent=NR_047525; 
chr1	refGene	CDS	787306	787490	3	+	.	Parent=NR_047525; 
chr1	refGene	intron	787491	788049	3	+	.	Parent=NR_047525; 
chr1	refGene	CDS	788050	788146	4	+	.	Parent=NR_047525; 
chr1	refGene	intron	788147	788769	4	+	.	Parent=NR_047525; 
chr1	refGene	CDS	788770	794826	5	+	.	Parent=NR_047525; 
chr1	refGene	intergenic	794827	803449	.	.	.	ID=NR_047525; name=LOC643837;|ID=NR_027055; name=FAM41C;
chr1	refGene	ncRNA	803450	812182	.	-	.	ID=NR_027055; name=FAM41C;
chr1	refGene	CDS	803450	804055	3	-	.	Parent=NR_027055; 
chr1	refGene	intron	804056	809490	2	-	.	Parent=NR_027055; 
chr1	refGene	CDS	809491	810535	2	-	.	Parent=NR_027055; 
chr1	refGene	intron	810536	812124	1	-	.	Parent=NR_027055; 
chr1	refGene	CDS	812125	812182	1	-	.	Parent=NR_027055; 
chr1	refGene	intergenic	812183	852951	.	.	.	ID=NR_027055; name=FAM41C;|ID=NR_026874; name=LOC100130417;
chr1	refGene	ncRNA	852952	854817	.	-	.	ID=NR_026874; name=LOC100130417;
chr1	refGene	CDS	852952	853100	4	-	.	Parent=NR_026874; 
chr1	refGene	intron	853101	853400	3	-	.	Parent=NR_026874; 
chr1	refGene	CDS	853401	853555	3	-	.	Parent=NR_026874; 
chr1	refGene	intron	853556	854203	2	-	.	Parent=NR_026874; 
chr1	refGene	CDS	854204	854295	2	-	.	Parent=NR_026874; 
chr1	refGene	intron	854296	854713	1	-	.	Parent=NR_026874; 
chr1	refGene	CDS	854714	854817	1	-	.	Parent=NR_026874; 
chr1	refGene	intergenic	854818	861119	.	.	.	ID=NR_026874; name=LOC100130417;|ID=NM_152486; name=SAMD11;
chr1	refGene	mRNA	861120	879961	.	+	.	ID=NM_152486; name=SAMD11;
chr1	refGene	5-UTR	861120	861180	.	+	.	Parent=NM_152486; 
chr1	refGene	intron	861181	861300	1	+	.	Parent=NM_152486; 
chr1	refGene	5-UTR	861301	861320	.	+	.	Parent=NM_152486; 
chr1	refGene	CDS	861321	861393	1	+	0	Parent=NM_152486;

第6列为exon和intron的排列顺序，第8列为翻译偏移量，intergenic为添加两个gene或transcript之间的距离等，加强版的gff格式~~~


 

 类似资料：

相关阅读

SQLite中的Select语句可识别行号如何识别Java套接字中的EOF？龙卷风：识别/跟踪websocket的连接？python 实现性别识别在识别为9的地方无法识别09

相关文章

虎牙 语音识别算法一面面经Go语言可变参数TensorFlow搭建神经网络实现手写识别GCC自动识别的文件扩展名（一览表）快手视觉识别搜索算法一面凉经

相关问答

用java_home环境变量识别问题源日食不识别渐变依赖关系识别和解决javax.el.PropertyNotFoundException：目标不可达图QL突变中无法识别的参数无法识别局部变量类型推断

相关文档

Linkface 人脸识别技术 帮助文档变现猫接入文档 v1.4略知 知识付费与内容变现 帮助文档初识 PhalApi：探索接口服务开发的技艺初识 Spring Security
快捷导航： 新手教程 算法原理 架构设计 Java进阶 数据库进阶 大厂专栏 面试经验 编程笔记 编程问答 所有专题 文档资料 工具软件 电子书籍 小牛导航
在线工具： 房贷计算器 个税计算器 Linux命令查询 Json格式化 正则表达式 颜色转换 AES加解密 SHA1加密 MD5加密 毒鸡汤 字数统计 随机密码生成 进制转换 Base64编解码 励志句子
Copyright © 2019-2024 小牛知识库@xnip.cn. All Rights Reserved.