|
0
|
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
|
