1

我想计算转录表达式，因此我需要获取 bam 文件中所有读取的映射数量。我目前的程序是使用 Bio::DB::Sam 获取整体转录本并获取映射到其上的读数。结果存储在以 read_name 作为键（10 个字母）和 number_of_mappings 作为值（整数）的哈希中。

这是我正在使用的代码：

use Bio::DB:Sam;
use strict;

my %global_read_occurrences;


sub getGlobalReadOccurrences {

 my ($ids, $bam_file) = @_;

 $sam = Bio::DB::Sam -> new (-bam => $bam_file);

 foreach my $id (@{$ids}){
   my $alignments = $sam -> get_features_by_location(-seq_id => $transcript_id, -iterator => 1);


  while (my $alignment = $alignments -> next_seq){

   my $read_name = $alignment -> query -> name;

   if (exists($global_read_occurrences{$read_name})){
    $global_read_occurrences{$read_name}++;
   }
   else {
    $global_read_occurrences{$read_name} = 1;
   }
  }
 }
}

我的问题：是否有任何其他可能性我可以直接获得每次读取的全局映射数量并且我不必查看所有成绩单？我在 Bio::DB::Sam 中找不到任何像 $sam -> getNumberOfMappings($read_name);

我正在使用具有超过 5000 万个映射读取的 bam 文件，因此哈希将需要大量内存资源（有时大约 40 GB）这实际上可能还是来自其他地方？还有其他方法可以用更少的内存来存储数据吗？

非常感谢！

4

1 回答 1

1

BAM 文件通常按染色体位置排序，而不是按读取的名称，因此读取映射可以位于文件中的任何位置。最简单的方法是转到 SAM 文件并运行一个简单的 shell 命令：

 cut -f1,1 myfile.sam | sort | uniq -c

这将产生一个像这样的文件：

  2 HWI-EAS299_4_30M2BAAXX:2:99:965:826
  2 HWI-EAS299_4_30M2BAAXX:2:99:966:1932
  2 HWI-EAS299_4_30M2BAAXX:2:99:971:146
  2 HWI-EAS299_4_30M2BAAXX:2:9:997:1263
  2 HWI-EAS299_4_30M2BAAXX:2:99:972:281
  2 HWI-EAS299_4_30M2BAAXX:2:99:973:1904
  1 HWI-EAS299_4_30M2BAAXX:2:99:976:186
  2 HWI-EAS299_4_30M2BAAXX:2:99:986:687
  6 HWI-EAS299_4_30M2BAAXX:2:99:987:165
  2 HWI-EAS299_4_30M2BAAXX:2:99:99:1582
  2 HWI-EAS299_4_30M2BAAXX:2:99:99:160
  2 HWI-EAS299_4_30M2BAAXX:2:99:998:1139

第一列是映射的计数。第二个是读取名称。

于 2012-01-15T18:21:58.867 回答