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;
   }
  }
 }
}

私の質問: 読み取りごとのグローバル マッピングの数を直接取得でき、すべてのトランスクリプトを調べる必要がない他の可能性はありますか? $sam -> getNumberOfMappings($read_name); のような Bio::DB::Sam でサブが見つかりませんでした。

私は 5,000 万回を超える読み取りがマッピングされた bam ファイルを使用しているため、ハッシュには巨大なメモリ リソース (約 40 GB の場合もあります) が必要になります。より少ないメモリでデータを保存する他の可能性はありますか?

どうもありがとう！

4

1 に答える 1

1

BAM ファイルは通常、リードの名前ではなく染色体の位置でソートされるため、リード マッピングはファイル内のどこにでも配置できます。最も簡単な方法は、SAM ファイルに移動して単純なシェル コマンドを実行することです。

 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

最初の列はマッピングの数です。2 番目は読み取り名です。

于 2012-01-15T18:21:58.867 に答える