Bio :: DB :: Sam - bamファイルのすべての読み込みのマッピング数を取得

Question

私は転写表現を計算したいので、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）のようなサブを見つけることができませんでした。

マップされた読み込み数が5,000万を超えるバムファイルを使用しているため、ハッシュに膨大なメモリリソース（時には約40 GB）が必要になることがあります。これは実際に可能なのでしょうか？そして、より少ないmemでデータを格納する他の可能性はありますか？

ありがとうございます！

Answer 1

BAMファイルは、通常、読み込みの名前ではなく、染色体の場所によってソートされるため、読み込みのマッピングはファイルのどこにでも配置できます。最初の列はのカウントです

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 

：これは、このようなファイルを生成します

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

：あなたがするための最も簡単な方法は、SAMファイルにアクセスし、簡単なシェルコマンドを実行することですマッピング。 2番目は読み込み名です。