एक फाइल से क्षेत्रों है कि अन्य फ़ाइल में क्षेत्रों का हिस्सा हैं जाओ (लूप के बिना)

Question

मैं दो फ़ाइलों है:

regions.txt: पहले स्तंभ गुणसूत्र नाम, दूसरे और तीसरे आरंभ और समाप्ति स्थिति हो रहा है।

1 100 200 
1 400 600 
2 600 700 

कवरेज.txt: पहला कॉलम गुणसूत्र नाम है, फिर दूसरा और तीसरा प्रारंभ और अंत स्थिति है, और अंतिम कॉलम स्कोर है।

1 100 101 5 
1 101 102 7 
1 103 105 8 
2 600 601 10 
2 601 602 15 

यह फ़ाइल लगभग 300 जीबी लाइनों के साथ लगभग 15 जीबी है।

मैं मूल रूप से covers.txt में सभी क्षेत्रों का अर्थ प्राप्त करना चाहता हूं जो क्षेत्रों.txt में प्रत्येक क्षेत्र में हैं।

दूसरे शब्दों में, regions.txt में पहली पंक्ति से शुरू करें, यदि cover.txt में कोई पंक्ति है, जिसमें समान गुणसूत्र है, प्रारंभ-कवरेज> = प्रारंभ-क्षेत्र है, और अंत-कवरेज < = अंत क्षेत्र, फिर अपने स्कोर को एक नई सरणी में सहेजें। सभी coverages.txt में खोज खत्म करने के बाद क्षेत्र गुणसूत्र, प्रारंभ, अंत, और पाया गया सभी स्कोर का मतलब प्रिंट करें।

अपेक्षित उत्पादन:

1 100 200 14.6 which is (5+7+8)/3 
1 400 600 0  no match at coverages.txt 
2 600 700 12.5 which is (10+15)/2 

मैं निम्नलिखित MATLAB स्क्रिप्ट है जो बहुत लंबा समय लग के बाद से मैं coverage.txt कई समय के साथ पाश के लिए है का निर्माण किया। मुझे नहीं पता कि एक तेज अजीब समान स्क्रिप्ट कैसे बनाएं।

मेरे matlab स्क्रिप्ट

fc = fopen('coverage.txt', 'r'); 
ft = fopen('regions.txt', 'r'); 
fw = fopen('out.txt', 'w'); 

while feof(ft) == 0 

linet = fgetl(ft); 
scant = textscan(linet, '%d%d%d'); 
tchr = scant{1}; 
tx = scant{2}; 
ty = scant{3}; 
coverages = []; 

    frewind(fc); 
    while feof(fc) == 0 

    linec = fgetl(fc); 
    scanc = textscan(linec, '%d%d%d%d'); 
    cchr = scanc{1}; 
    cx = scanc{2}; 
    cy = scanc{3}; 
    cov = scanc{4}; 


     if (cchr == tchr) && (cx >= tx) && (cy <= ty) 

      coverages = cat(2, coverages, cov); 

     end 

    end 

    covmed = median(coverages); 
    fprintf(fw, '%d\t%d\t%d\t%d\n', tchr, tx, ty, covmed); 

end  

कोई सुझाव का उपयोग कर एक विकल्प बनाने की AWK, पर्ल, या, ... आदि मैं कृपा विस्फोट अगर कोई मुझे सिखाने कैसे में सभी छोरों से छुटकारा पाने के कर सकते हैं मेरे matlab लिपि।

धन्यवाद

Answer 1

यहां एक पर्ल समाधान है। मैं क्रोमोसोम के माध्यम से विभिन्न श्रेणियों तक पहुंचने के लिए हैश (उर्फ शब्दकोश) का उपयोग करता हूं, इस प्रकार लूप पुनरावृत्तियों की संख्या को कम करता है।

यह संभावित रूप से कुशल है, क्योंकि मैं प्रत्येक इनपुट लाइन पर regions.txt पर पूर्ण लूप नहीं करता हूं। मल्टीथ्रेडिंग का उपयोग होने पर क्षमता को और बढ़ाया जा सकता है।

#!/usr/bin/perl 

my ($rangefile) = @ARGV; 
open my $rFH, '<', $rangefile or die "Can't open $rangefile"; 

# construct the ranges. The chromosome is used as range key. 
my %ranges; 
while (<$rFH>) { 
    chomp; 
    my @field = split /\s+/; 
    push @{$ranges{$field[0]}}, [@field[1,2], 0, 0]; 
} 
close $rFH; 

# iterate over all the input 
while (my $line = <STDIN>) { 
    chomp $line; 
    my ($chrom, $lower, $upper, $value) = split /\s+/, $line; 
    # only loop over ranges with matching chromosome 
    foreach my $range (@{$ranges{$chrom}}) { 
     if ($$range[0] <= $lower and $upper <= $$range[1]) { 
      $$range[2]++; 
      $$range[3] += $value; 
      last; # break out of foreach early because ranges don't overlap 
     } 
    } 
} 

# create the report 
foreach my $chrom (sort {$a <=> $b} keys %ranges) { 
    foreach my $range (@{$ranges{$chrom}}) { 
     my $value = $$range[2] ? $$range[3]/$$range[2] : 0; 
     printf "%d %d %d %.1f\n", $chrom, @$range[0,1], $value; 
    } 
} 

उदाहरण मंगलाचरण:

1 100 200 6.7 
1 400 600 0.0 
2 600 700 12.5 

(क्योंकि (5 + 7 + 8)/3 = 6.66 ...)

Answer 2

आम तौर पर, मैं आर में फ़ाइलें लोड होगा और यह गणना, लेकिन यह देखते हुए कि उनमें से एक तो बहुत बड़ा है, यह एक समस्या बन जाएगा। यहां कुछ विचार दिए गए हैं जो इसे हल करने में आपकी सहायता कर सकते हैं।

1. क्रोमोसोम द्वारा coverage.txt को विभाजित करने पर विचार करें। इससे गणना कम मांग होगी।

2. इसके बजाय coverage.txt से अधिक पाशन के

   , आपको पहले स्मृति में regions.txt पूरा पढ़ा (मुझे लगता है यह ज्यादा छोटा होता है)। प्रत्येक क्षेत्र के लिए, आप स्कोर और एक संख्या रखते हैं।

3. प्रक्रिया coverage.txt लाइन द्वारा लाइन। प्रत्येक पंक्ति के लिए, आप गुणसूत्र और क्षेत्र निर्धारित करते हैं जो इस विशेष खिंचाव से संबंधित है। इसके लिए कुछ फुटवर्क की आवश्यकता होगी, लेकिन यदि regions.txt बहुत बड़ा नहीं है, तो यह अधिक कुशल हो सकता है। क्षेत्र के स्कोर और स्कोर संख्या में स्कोर जोड़ें।

एक वैकल्पिक, सबसे प्रभावी तरीका दोनों फ़ाइलों को क्रोमोसोम द्वारा पहले क्रमबद्ध करने की आवश्यकता होती है, फिर स्थिति के अनुसार।

1. regions.txt से एक लाइन लें।गुणसूत्र और पदों को रिकॉर्ड करें। यदि पिछले लूप से शेष रेखा है, तो 3 पर जाएं .; अन्यथा 2.

2. coverage.txt से एक लाइन लें।

3. जांचें कि यह वर्तमान क्षेत्र में है या नहीं।

   • हाँ: क्षेत्र में वृद्धि, वृद्धि संख्या जोड़ें। संख्या से विभाजित स्कोर, उत्पादन के लिए वर्तमान क्षेत्र लिखते हैं, करने के लिए 1.

जाना यह पिछले विधि कुछ ठीक ट्यूनिंग आवश्यकता है, लेकिन सबसे कारगर होगा - यह: 2.

नहीं पर ले जाएं प्रत्येक फ़ाइल को केवल एक बार जाने की आवश्यकता है और स्मृति में लगभग कुछ भी स्टोर करने की आवश्यकता नहीं है।

Answer 3

यहाँ क्षेत्रों में बिन अपने कवरेज के लिए एक सरल तरीका है MATLAB:

% extract the regions extents 
bins = regions(:,2:3)'; 
bins = bins(:); 

% extract the coverage - only the start is needed 
covs = coverage(:,2); 

% use histc to place the coverage start into proper regions 
% this line counts how many coverages there are in a region 
% and assigns them proper region ids. 
[h, i]= histc(covs(:), bins(:)); 

% sum the scores into correct regions (second output of histc gives this) 
total = accumarray(i, coverage(:,4), [numel(bins),1]); 

% average the score in regions (first output of histc is useful) 
avg = total./h; 

% remove every second entry - our regions are defined by start/end 
avg = avg(1:2:end); 

अब इस यह सोचते हैं कि क्षेत्रों गैर अतिव्यापी हैं काम करता है, लेकिन मुझे लगता है कि यह मामला है। साथ ही, coverage फ़ाइल में प्रत्येक प्रविष्टि को कुछ क्षेत्र में गिरना पड़ता है।

इसके अलावा, अगर आप पूरी फाइल में पढ़ने से बचना चाहते हैं, तो कवरेज पर इस दृष्टिकोण को 'ब्लॉक' करने के लिए तुच्छ है। आपको केवल bins की आवश्यकता है, आपकी क्षेत्र फ़ाइल, जो संभवतः छोटी है। आप ब्लॉक में कवरेज को संसाधित कर सकते हैं, total में वृद्धिशील रूप से जोड़ सकते हैं और अंत में औसत गणना कर सकते हैं।

Answer 4

यहाँ: उदाहरण के इनपुट पर

$ perl script.pl regions.txt <coverage.txt >output.txt 

आउटपुट एक तरीका join और awk का उपयोग कर। script.awk की

join regions.txt coverage.txt | awk -f script.awk - regions.txt 

सामग्री::

FNR==NR && $4>=$2 && $5<=$3 { 
    sum[$1 FS $2 FS $3]+=$6 
    cnt[$1 FS $2 FS $3]++ 
    next 
} 

{ 
    if ($1 FS $2 FS $3 in sum) { 
     printf "%s %.1f\n", $0, sum[$1 FS $2 FS $3]/cnt[$1 FS $2 FS $3] 
    } 
    else if (NF == 3) { 
     print $0 " 0" 
    } 
} 

परिणाम: की तरह चलाने के लिए

1 100 200 6.7 
1 400 600 0 
2 600 700 12.5 

वैकल्पिक रूप से, यहां एक लाइनर है:

join regions.txt coverage.txt | awk 'FNR==NR && $4>=$2 && $5<=$3 { sum[$1 FS $2 FS $3]+=$6; cnt[$1 FS $2 FS $3]++; next } { if ($1 FS $2 FS $3 in sum) printf "%s %.1f\n", $0, sum[$1 FS $2 FS $3]/cnt[$1 FS $2 FS $3]; else if (NF == 3) print $0 " 0" }' - regions.txt