本文举例介绍生物信息学多文件数据处理的一些编程思路，希望对大家有所帮助。
如果你是想从gff格式提取CDS，请看另一篇《Extract cds fastas from a gff annotation + reference sequence》。
师妹对某物种测序获得很多scaffold，用glean预测了基因结构（CDS、mRNA等），另外还预测了snp位点，她想统计各mRNA中的snp位点。
文件信息
1) 基因结构预测结果文件a.gff如下，第一列是序列id，第二列是注释信息来源，第三列是注释信息类型，四五列是开始与结束位置，第七列是序列方向/链(strand），第八列phase，第九列是描述。详见gff格式说明。
scaf002 glean mRNA 11153 14469 0.459 + . somestring scaf002 glean CDS 11153 11326 . + 0 somestring
备注： 此文件约10万行，文件大小7M左右。scaffold名称数量大约几百，存在多种格式。一条scaffold的所有记录是按照位置从小到大排序的。
2) snp文件b.snp，有三列，如下所示，第一列是scaffold名称，第二列是snp位点。
nscaf285 181 C T T T - C C C - - - - - T C C C - C Y C C C - C - - - - - - - - C Y T Y C Y T C SWSNP00000001 nscaf285 201 T - - - - T T - - - - T T T T T T - T T T T T - T - - - - - T T - Y T - T T Y T Y SWSNP00000002
备注：此文件大约200万行，文件大小约200M。scaffold名称数量大约几十，且都包含在a.gff文件中的scaffold名称中。
思路
背景：
a.gff中有多条scaffold,每条scaffold上有多个mRNA；b.snp中有多条scaffold，每条scaffold中有多个snp位点。
目标：
找到各scaffold上每条mRNA上的snp位点，结果按照snp位点数量从多到少对各mRNA排序输出。
大致思路：
嵌套循环snp数据（外循环）和mRNA信息数据（内循环），看每个snp位点是否存在与某个scaffold的某条snp上。
细化：
snp数据文件较大，且每行的记录只需要使用一次（和mRNA位点信息比较），如果读入内存再比较势必浪费大量内存。故只需要从a.gff一次性读入mRNA信息，然后在逐行读取snp文件b.snp过程中，逐一比较snp位点与mRNA位点。
数据存储结构：
前面说了要从a.gff一次性读入mRNA信息，数据结构的好坏与运行速度密切相关。
HUGOMORE42
1）mRNA位点信息。由于a.gff有多条scaffold，且b.snp中每行记录提供了scaffold名称，则自然想到用字典<1>（perl中为hash，python为map）存储mRNA为信息，通过scaffold名称查询对应的mRNA位点信息，从而缩小了对mRNA位点信息比较的范围，避免了将所有mRNA信息逐行存入数组中与snp位点信息进行一对一匹配。
此外，由于一条scaffold上有多处mRNA，我们关心的是各mRNA的位置，则可以将各mRNA的位置信息存入一个字典<2>，并将所有mRNA的信息存入一个列表（perl中为数组array，python中为列表list），匹配的时候遍历列表一一匹配即可。由于必须要遍历一条scaffold上所有mRNA，则不必考虑用字典存储所有mRNA来提高速度。
2） 存储各mRNA的snp位点信息。可以单独存储到一个数据结构中，但仍需要通过scaffold名称+mRNA信息（起始位点）来作为唯一标示，这些在前面讲的mRNA位点信息中已经存在，则只需在上面一段中的字典<2>中增加一对键值：”snp”=>位点列表的引用。在逐一比较snp位点与mRNA位点的过程中，发现一条snp存在某条mRNA上，则将该snp位点增加到该mRNA的snp位点列表中去。
最终存储所有mRNA的变量的数据结构如下（perl的Data::Dumper打印，非常直观）,
$VAR1 = { # $mrnas={...}表示其为指向字典的引用 'nscaf3092' => [ # 'nscaf3092' => [...]表示其为指向列表的引用 { # 一个mRNA的所有信息存在字典中 'strand' => '-', 'end' => '4128', 'start' => '1536', 'snp' => [1, 2, 3], # 又是[] 'snp_amount' => '3' } ], 'scaffold838' => [ { 'strand' => '-', 'end' => '17649', 'start' => '16538', 'snp' => [1, 2, 3, 4], 'snp_amount' => '4' }, { 'strand' => '-', 'end' => '15974', 'start' => '6356', 'snp' => [1, 2, 3, 4, 5], 'snp_amount' => '5' } ] }
其它：
结果按照snp位点数量从多到少排序输出，可以使用施瓦茨变换排序。
b.snp行数较多，可以实时显示进度，这在同系列的另一篇文章中有例子。该文章中的文件更大，记录数更多。
如果为了减少内存使用，当记录非常非常多的时候，记录单个mRNA信息的字典的key可以更精简一些。
代码实现（Perl）
HUGOMORE42
#!/usr/bin/perl
use strict;
# 因为需要在终端显示一些进度信息，所以结果写入结果文件中
die "Usage:$0 <*.gff> <*.snp> <result file>\n"
unless @ARGV == 3;
my $gff_file = shift @ARGV;
my $snp_file = shift @ARGV;
my $result_file = shift @ARGV;
print "read gff file ..."; # 当运算量大的时候，进度信息很有用
# 读取mRNA信息
my $mrnas = read_data_from_gff_file( $gff_file, "mRNA" );
print "Done.\n";
print "parse snp file\n";
my ( $seq, $loc ); # 只记录序列名称和snp位点
my $mrna; # 一个临时变量，存放单个mRNA的信息
open IN, "<", $snp_file
or die "failed to open file: $snp_file\n";
my $count = 0; # 用于进度显示
my $| = 1; # 取消输出缓存设置，用于快速实时显示进度
while (<IN>) {
$count++;
# 显示进度，\r是回车，即删除整行的文字，并显示新的文字，实现动态显示进度
print "\r$count";
next unless /(.+?)\t(\d+)\t/; # 格式约束
( $seq, $loc ) = ( $1, $2 ); # 赋值给变量
next unless exists $$mrnas{$seq}; # 筛选出当前snp记录的序列对应的所有mRNA信息
# $$mrnas{$seq} 等价于 $mrnas->{$seq}，后面的类似
for $mrna ( @{ $$mrnas{$seq} } ) { # 遍历该序列所有的mRNA位点
# 检查snp位点是否在mRNA中，如果不在则检查下一个mRNA
next unless $loc >= $$mrna{start} and $loc <= $$mrna{end};
$$mrna{snp} = [] unless ref $$mrna{snp} eq 'ARRAY'; # 初始化存储一个mRNA所有snp位点的指向列表的引用
push @{ $$mrna{snp} }, $loc; # 记录当前的snp位点
}
}
close IN;
# store and sort
# 先把有snp位点的mRNA信息存入一个列表
my @mrnas_with_snp;
for my $seq ( sort keys %$mrnas ) {
for $mrna ( @{ $$mrnas{$seq} } ) {
next unless ref $$mrna{snp} eq 'ARRAY'; # 只存入有snp的mRNA
$$mrna{snp_amount} = scalar @{ $$mrna{snp} }; # 新加入一个snp数量的信息
$$mrna{seq} = $seq; # 把序列名称也加进来
push @mrnas_with_snp, $mrna;
}
}
# 用施瓦茨变换来按照snp位点个数从大到小进行排序。
print "\nsort result ...";
@mrnas_with_snp = map { $_->[0] }
sort { $b->[1] <=> $a->[1] }
map { [ $_, $$_{snp_amount} ] } @mrnas_with_snp;
print "Down.\n";
# 把结果写入文件
$| = 0;
print "write result to $result_file ...";
open OUT, ">", $result_file
or die "failed to open file: $result_file\n";
for $mrna (@mrnas_with_snp) {
print OUT "$$mrna{seq}\t$$mrna{strand}\t$$mrna{start}\t$$mrna{end}\t";
print OUT $$mrna{snp_amount}, "\t";
print OUT join ",", @{ $$mrna{snp} };
print OUT "\n";
}
close OUT;
print "Done.\n";
# 读取gff文件的函数
sub read_data_from_gff_file ($$) { # $$限制函数参数为两个标量
my ( $file, $TYPE ) = @_;
my $data = {}; # 储存信息的变量，初始化为指向hash的引用
my ( @tmp, $seq, $type, $start, $end, $strand ); # 需要存储的信息的变量
open IN, "<", $file
or die "failed to open file: $file\n";
while (<IN>) {
@tmp = split /\t/, $_;
next unless @tmp == 9;
( $seq, $type, $start, $end, $strand ) =
( $tmp[0], $tmp[2], $tmp[3], $tmp[4], $tmp[6] );
next unless $type eq $TYPE; # only mRNA
$$data{$seq} = [] unless ref $$data{$seq} eq 'ARRAY'; # 初始化变量
push @{ $$data{$seq} },
{ start => $start, end => $end, strand => $strand };
}
close IN;
return $data;
}
总结
简单的数据处理程序，主要是程序结构和数据结构的选择。当数据量很大的时候，程序结构尤为重要；数据结构主要就是列表和字典的灵活应用。
类似的文章
-EOF-