统计2.2.1cufflinks的矩阵和差异表达结果--大量数据结构

统计最新版的cufflinks的矩阵结果和差异表达结果(针对基因),涉及的数据结构比较复杂,有二维数组,hash的数组,数组的hash等等:

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

die "perl $0 <file: norm.fpkm/diff.fpkm> <option: interclass/interblock/de> <outprefix>\n" if(@ARGV ne 3);

open EXP, $ARGV[0] or die $!;
my $string = <EXP>;
chomp($string);
my @labels = split /\t/, $string;
shift @labels;
@labels = map {$_ =~ s/_\d+$// ? $_ : $_} @labels;

my ($sum, @gene, %hash, @rc, @final);

if($ARGV[1] eq "interblock")
{
	open OUT, ">> $ARGV[2].txt" or die $!;
	my ($x, $flag) = (0, 0);
	for(my $i = 0; $i < @labels; $i ++)
	{
		if(exists $hash{$labels[$i]})
		{
			push @{$rc[$x]}, $i;
			$flag = 1;
		}else{
			$x ++ if($flag == 1);
			$hash{$labels[$i]} = 0;
			push @final, $labels[$i];
			push @{$rc[$x]}, $i;
		}
	}
=cut
	for (my $i = 0; $i < @rc; $i ++)
	{
		for (my $j = 0; $j < @{$rc[$i]}; $j ++)
		{
			print "$i\t$rc[$i][$j]\n";
		}
	}
=cut
	while(<EXP>)
	{
		chomp;
		my @tmp = split;
		for(my $i = 0; $i < @rc; $i ++)
		{
			my $fpkm;
			for(my $j = 0; $j < @{$rc[$i]}; $j ++)
			{
				$fpkm += $tmp[$rc[$i][$j] + 1];
			}
			if($fpkm > 0 and $tmp[0] !~ /^XLOC/)
			{
				$gene[$i]{ref} ++;
			}elsif($fpkm > 0){
				$gene[$i]{new} ++;
			}
		}
	}
	print OUT "\nTYPE: interBlock\n";
	for(my $i = 0; $i < @final; $i ++)
	{
		print OUT "$final[$i]\t";
		foreach my $key(sort keys %{$gene[$i]})
		{
			print OUT "$key: $gene[$i]{$key}\t";
		}
		print OUT "\n";
	}

}elsif($ARGV[1] eq "interclass"){
	open OUT, "> $ARGV[2].txt" or die $!;
	my $all_sum;
	while(<EXP>)
	{
		chomp;
		my @tmp = split;
		my $all;
		for(my $i = 1; $i < @tmp; $i ++)
		{
			$all += $tmp[$i];
			if($tmp[$i] > 0 and $tmp[0] !~ /^XLOC/)
			{
				$gene[$i - 1]{ref} ++;
			}elsif($tmp[$i] > 0){
				$gene[$i - 1]{new} ++;
			}
		}
		$all_sum ++ if($all > 0);
		$sum ++;
	}
	print OUT "All genes: $sum\n";
	print OUT "All samples cover: $all_sum\n\n";
	print OUT "TYPE: interClass\n";
	for(my $i = 0; $i < @labels; $i ++)
	{
		print OUT "$labels[$i]\t";
		foreach my $key(sort keys %{$gene[$i]})
		{
			print OUT "$key: $gene[$i]{$key}\t";
		}
		print OUT "\n";
	}

}elsif($ARGV[1] eq "de"){
	open OUT, "> $ARGV[2].txt" or die $!;
	open DE, "> $ARGV[2].diffexpression.txt" or die $!;
	print DE "$string\n";
	my %gene;
	while(<EXP>)
	{
		chomp;
		my @tmp = split;
		if(abs($tmp[9]) > 1 and $tmp[13] eq "yes")
		{
			print DE "$_\n";
			$hash{$tmp[4]."-".$tmp[5]} ++;
			if($tmp[9] > 1)
			{
				$gene{$tmp[4]."-".$tmp[5]}[0] ++;
			}elsif($tmp[9] < -1){
				$gene{$tmp[4]."-".$tmp[5]}[1] ++;
			}
		}
	}
	foreach my $key(keys %hash)
	{
		$gene{$key}[0] = 0 if(!defined $gene{$key}[0]);
		$gene{$key}[1] = 0 if(!defined $gene{$key}[1]);
		print OUT "$key: $hash{$key}\tUP: $gene{$key}[0]\tDOWN: $gene{$key}[1]\n";
	}
	
}else{
	die "check correct option\n";
}


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值