annotate rDiff/src/tools/convert_reads_to_region_indicators.m @ 3:29a698dc5c7e default tip

Merge multiple heads.
author Dave Bouvier <dave@bx.psu.edu>
date Mon, 27 Jan 2014 14:15:36 -0500
parents 0f80a5141704
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
1 function [NEW_READS,UNEXPLAINED_READS,UNEXPLAINED_INDEX] = convert_reads_to_region_indicators(READS, gene)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
2 % Convert the reads into counts of EIRS
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
3
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
4 %UNEXPLAINED_REGIONS
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
5 UNEXPLAINED_REGIONS = 1:(length(gene.splicingevents) - 1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
6 UNEXPLAINED_REGIONS(gene.sequence == 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
7
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
8 %Extend GRAPHNODES to also include introns
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
9 EXT_GRAPHNODES = zeros(size(gene.graphnodes,1), length(gene.splicingevents) - 1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
10 EXT_GRAPHNODES(:, gene.unique_new_exons) = gene.graphnodes;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
11 % This puts gene.graphnodes into the exonic positions of EXT_GRAPHNODES
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
12
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
13 %Mark for each read into which region it falls
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
14 TEMP_READS = sparse(zeros(size(READS,1), length(gene.splicingevents) - 1));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
15 for i = 1:(length(gene.splicingevents) - 1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
16 TEMP_READS(:,i) = sum(READS(:, gene.splicingevents(i):(gene.splicingevents(i + 1) - 1)), 2) > 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
17 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
18
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
19 %UNEXPLAINED_INDEX = (((sum(EXT_GRAPHNODES,1) == 0) * TEMP_READS1') > 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
20 UNEXPLAINED_INDEX = ((( gene.sequence== 0) * TEMP_READS') > 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
21 UNEXPLAINED_READS = TEMP_READS(UNEXPLAINED_INDEX,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
22
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
23 TEMP_READS = TEMP_READS(not(UNEXPLAINED_INDEX),:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
24
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
25 %find the row of EXT_GRAPHNODES minimizing the mismatch between the rows of EXT_GRAPHNODES
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
26 %and the reads
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
27 %[MAX_VAL,MAX_NODE] = max( (1+diag(1./sum(EXT_GRAPHNODES,2))) * (EXT_GRAPHNODES * TEMP_READS'),[],1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
28 [MAX_VAL,MAX_NODE] = max(EXT_GRAPHNODES*TEMP_READS'+diag(1./sum(EXT_GRAPHNODES,2))*EXT_GRAPHNODES*TEMP_READS',[],1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
29
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
30 %Create sparse read matrix
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
31 NEW_READS = spconvert([[1,1,1; size(EXT_GRAPHNODES,1),2,1]; [MAX_NODE',(3:(length(MAX_NODE)+2))', ones(size(MAX_VAL))']])';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
32 NEW_READS = NEW_READS(3:end,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
33
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
34
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
35 return
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
36