Mercurial > repos > vipints > rdiff
comparison rDiff/src/tools/convert_reads_to_region_indicators.m @ 0:0f80a5141704
version 0.3 uploaded
author | vipints |
---|---|
date | Thu, 14 Feb 2013 23:38:36 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:0f80a5141704 |
---|---|
1 function [NEW_READS,UNEXPLAINED_READS,UNEXPLAINED_INDEX] = convert_reads_to_region_indicators(READS, gene) | |
2 % Convert the reads into counts of EIRS | |
3 | |
4 %UNEXPLAINED_REGIONS | |
5 UNEXPLAINED_REGIONS = 1:(length(gene.splicingevents) - 1); | |
6 UNEXPLAINED_REGIONS(gene.sequence == 0); | |
7 | |
8 %Extend GRAPHNODES to also include introns | |
9 EXT_GRAPHNODES = zeros(size(gene.graphnodes,1), length(gene.splicingevents) - 1); | |
10 EXT_GRAPHNODES(:, gene.unique_new_exons) = gene.graphnodes; | |
11 % This puts gene.graphnodes into the exonic positions of EXT_GRAPHNODES | |
12 | |
13 %Mark for each read into which region it falls | |
14 TEMP_READS = sparse(zeros(size(READS,1), length(gene.splicingevents) - 1)); | |
15 for i = 1:(length(gene.splicingevents) - 1) | |
16 TEMP_READS(:,i) = sum(READS(:, gene.splicingevents(i):(gene.splicingevents(i + 1) - 1)), 2) > 0; | |
17 end | |
18 | |
19 %UNEXPLAINED_INDEX = (((sum(EXT_GRAPHNODES,1) == 0) * TEMP_READS1') > 0); | |
20 UNEXPLAINED_INDEX = ((( gene.sequence== 0) * TEMP_READS') > 0); | |
21 UNEXPLAINED_READS = TEMP_READS(UNEXPLAINED_INDEX,:); | |
22 | |
23 TEMP_READS = TEMP_READS(not(UNEXPLAINED_INDEX),:); | |
24 | |
25 %find the row of EXT_GRAPHNODES minimizing the mismatch between the rows of EXT_GRAPHNODES | |
26 %and the reads | |
27 %[MAX_VAL,MAX_NODE] = max( (1+diag(1./sum(EXT_GRAPHNODES,2))) * (EXT_GRAPHNODES * TEMP_READS'),[],1); | |
28 [MAX_VAL,MAX_NODE] = max(EXT_GRAPHNODES*TEMP_READS'+diag(1./sum(EXT_GRAPHNODES,2))*EXT_GRAPHNODES*TEMP_READS',[],1); | |
29 | |
30 %Create sparse read matrix | |
31 NEW_READS = spconvert([[1,1,1; size(EXT_GRAPHNODES,1),2,1]; [MAX_NODE',(3:(length(MAX_NODE)+2))', ones(size(MAX_VAL))']])'; | |
32 NEW_READS = NEW_READS(3:end,:); | |
33 | |
34 | |
35 return | |
36 |