python与perl的矩阵转换及多样品的PCA

最新推荐文章于 2022-08-27 16:51:57 发布

浮生终有醒

最新推荐文章于 2022-08-27 16:51:57 发布

阅读量1.2k

点赞数

分类专栏：生物信息学

本文链接：https://blog.csdn.net/skenoy/article/details/8669157

版权

生物信息学专栏收录该内容

117 篇文章 9 订阅

订阅专栏

之前看过前辈用python转换矩阵，但python一直没系统学过，所以从网络中学到perl的矩阵转换，两个做了一下比较：

import sys
file = open(sys.argv[1], 'r')

arr = []
for line in file:
    info = line.rstrip().split()
    arr.append(info)
    
tarr = [[r[col] for r in arr] for col in range(len(arr[0]))]

for x in tarr:
    print '\t'.join(x)

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

my @matrix;
my $i=0;
open FA, $ARGV[0] || die $!;
while(<FA>)
{	
	foreach my $elem (split/\s+/, $_)
	{
		push @{$matrix[$i]}, $elem;
	}
	$i++;
}
close FA;

my @transition = map{my $x=$_;[map{$matrix[$_][$x]}0..$#matrix];}0..$#{$matrix[0]};

foreach my $row(@transition)
{
        foreach my $col (@{$row})
	{
                print $col, "\t";
        }
        print "\n";
}

关于perl的矩阵转换加入了一点点数据结构的知识了，一开始也没看懂数据结构，用的多了就熟练了，不过最常用的还是哈希的数组。
PCA，Principal Component Analysis，主成分分析，主要用于数据降维，但博主是生物信息领域，那么就看一下这方面的分析吧。在进行基因表达数据分析时，一个重要问题是确定每个实验数据是否是独立的，如果每次实验数据之间不是独立的，则会影响基因表达数据分析结果的准确性。对于利用二代测序所检测到的基因表达数据，如果用 PCA 方法进行分析，将基因作为变量通过分析确定一组“主要基因元素”，如果能够解释清楚基因的特征，那么就可以解释实验现象。拿之前的程序出来看看吧：

#!/usr/bin/env perl
use warnings;
use strict;
@ARGV == 0 and die "$0 <Matrix> <all> <need>\n";

open HEAD, shift or die $!;
open LST, shift or die $!;
open MATX, shift or die $!;

open RCMD , '>run_pca.R' or die $!;
open MX , '>mat.tab' or die $!;
open PF , '>pch_col.txt' or die $!;
open GT , '>genotypes.txt' or die $!;

my @ids = ();
while(<HEAD>){
	chomp;
	push @ids, $_;
}
close HEAD;

my %good_ids = ();
while(<LST>){
	chomp;
	$_ or next;
	my @F = split;
	$good_ids{$F[0]}  = $F[1];
}
close LST;

print RCMD 
'library("pls")
library("pcaMethods")
data<-read.table("mat.tab")
para <- read.table("pch_col.txt")
pdf("PCA.pdf")
res <-pca(data, method="ppca", nPcs=3, center=TRUE)
write.table(res@loadings, "loadings")
write.table(res@scores, "scores")
plotPcs(res, type = "scores", pch=as.character(para[,1]), col=as.integer(para[,2]))
';

for my $i (@ids){
	$_ = <MATX>;
	if(exists $good_ids{$i}){
		print MX $_;	
		if( $good_ids{$i} eq 'A'){
			print PF "@\t1\t$i\n";
		}elsif($good_ids{$i} eq 'B'){
			print PF "^\t2\t$i\n";
		}elsif($good_ids{$i} eq 'C'){
			print PF "&\t3\t$i\n";
		}elsif($good_ids{$i} eq 'D'){
			print PF "*\t4\t$i\n";
		}elsif($good_ids{$i} eq 'E'){
			print PF "-\t5\t$i\n";
		}elsif($good_ids{$i} eq 'F'){
			print PF "+\t6\t$i\n";
		}elsif($good_ids{$i} eq 'G'){
			print PF "#\t7\t$i\n";
		}elsif($good_ids{$i} eq 'H'){
			print PF ".\t8\t$i\n";
		}
		print GT "$i\t$good_ids{$i}\t$_";
	}
}

system ("R CMD BATCH run_pca.R");