Mercurial > repos > vipints > deseq_hts

annotate deseq-hts_2.0/src/get_read_counts.m @ 10:2fe512c7bfdf draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DESeq2 version 1.0.19 added to the repo
author vipints <vipin@cbio.mskcc.org>
date Tue, 08 Oct 2013 08:15:34 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
10
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
1 function get_read_counts(anno_dir, outfile, varargin)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
2 %
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
3 % -- input --
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
4 % anno_dir: directory of genes
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
5 % outfile: output file
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
6 % varargin: list of BAM files (at least two)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
7
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
8 % DESeq paths
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
9 global DESEQ2_PATH DESEQ2_SRC_PATH
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
10
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
11 % interpreter paths
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
12 global INTERPRETER MATLAB_BIN_PATH OCTAVE_BIN_PATH
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
13
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
14 % SAMTools path
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
15 global SAMTOOLS_DIR
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
16
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
17 %%%% paths
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
18 addpath(sprintf('%s/tools', DESEQ2_PATH));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
19 addpath(sprintf('%s/mex', DESEQ2_PATH));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
20 addpath(sprintf('%s', DESEQ2_SRC_PATH));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
21
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
22 deseq_config;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
23
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
24 %%% read list of replicate groups from variable length argument list
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
25 rg_list = cell(1,size(varargin, 2));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
26 file_list = cell();
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
27 file_cond_ids = [];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
28 file_rep_ids = [];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
29 for idx = 1:size(varargin, 2)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
30 rg_list(idx) = varargin(idx);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
31 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
32 idx = strmatch('', rg_list, 'exact');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
33 rg_list(idx) = [];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
34 for idx = 1:length(rg_list),
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
35 items = separate(rg_list{idx}, ':');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
36 for idx2 = 1:length(items)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
37 if isempty(deblank(items{idx2})),
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
38 continue;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
39 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
40 file_list{end + 1} = items{idx2};
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
41 file_cond_ids(end + 1) = idx;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
42 file_rep_ids(end + 1) = idx2;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
43 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
44 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
45 clear idx idx2;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
46
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
47 %%%%% adapt to number of input arguments
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
48 file_num = length(file_list);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
49 RESULTS = cell(1, file_num);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
50
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
51 %%%% get annotation file
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
52 load(sprintf('%s', anno_dir));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
53
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
54 %%%%% mask overlapping gene regions -> later not counted
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
55 [genes] = mask_dubl(genes,0);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
56
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
57 %%%% remove genes with no annotated exons or where no
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
58 idx = find(arrayfun(@(x)(~isempty(x.exons)*~isempty(x.start)*~isempty(x.stop)), genes));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
59 fprintf('removed %i of %i genes, which had either no exons annotated or lacked a start or stop position\n', size(genes, 2) - size(idx, 2), size(genes, 2))
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
60 genes = genes(idx);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
61 clear idx;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
62
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
63 %%%% check if genes have field chr_num
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
64 if ~isfield(genes, 'chr_num')
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
65 chrms = unique({genes(:).chr});
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
66 for i = 1:length(genes)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
67 genes(i).chr_num = strmatch(genes(i).chr, chrms, 'exact');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
68 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
69 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
70
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
71 %%%% iterate over all given bam files
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
72 for f_idx = 1:file_num
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
73 expr1_bam = fullfile('', file_list{f_idx});
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
74 STAT = cell(size(genes, 2),1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
75 for i=1:size(genes,2)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
76 RESULT = cell(1,7);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
77 gene = genes(i);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
78 RESULT{4} = f_idx;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
79 RESULT{1} = gene.name;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
80 if isempty(gene.exons)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
81 RESULT{2} = inf;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
82 RESULT{3} = inf;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
83 RESULT{5} = [inf,inf];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
84 STAT{i} = RESULT;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
85 continue;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
86 elseif or(isempty(gene.start),isempty(gene.stop))
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
87 RESULT{2} = inf;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
88 RESULT{3} = inf;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
89 RESULT{5} = [inf,inf];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
90 STAT{i} = RESULT;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
91 continue;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
92 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
93 if ~isempty(gene.chr_num),
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
94 [mask1, read_intron_list] = get_reads(expr1_bam, gene.chr, gene.start, gene.stop, '0');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
95 clear read_intron_list;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
96 else
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
97 mask1 = [];
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
98 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
99
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
100 if isempty(mask1)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
101 reads1 = zeros(0,gene.stop-gene.start+1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
102 else
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
103 reads1 = sparse(mask1(1,:)',mask1(2,:)',ones(size(mask1,2),1),max(mask1(1,:)),gene.stop-gene.start+1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
104 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
105 if ~isempty(reads1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
106 [reads1,FLAG] = remove_reads_from_other_genes(reads1,gene);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
107 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
108 L = size(reads1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
109 RESULT{2}=[size(reads1,1)]; % number of all reads falling in that gene
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
110 EXON_IDX=zeros(1,gene.stop-gene.start+1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
111 for t=1:size(gene.transcripts,2)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
112 for e=1:size(gene.exons{t},1)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
113 EXON_IDX((gene.exons{t}(e,1)-gene.start+1):(gene.exons{t}(e,2)-gene.start+1))=1;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
114 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
115 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
116 reads1 = reads1(sum(reads1(:,find(EXON_IDX)),2)>0,:);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
117 L1 = sum(EXON_IDX);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
118 RESULT{3}=[size(reads1,1)]; % number of reads overlapping to exons
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
119 RESULT{5}=[L, L1]; % size of reads1, number of exonic positions
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
120 % old and weighted poisson new ,weighted regions reads and
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
121 % unexplained reads
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
122 clear reads1;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
123 STAT{i} = RESULT;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
124 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
125 RESULTS{f_idx} = STAT;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
126 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
127
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
128 S=size(genes,2);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
129 READCOUNTS_ALL=zeros(S, file_num);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
130 READCOUNTS_EXON=zeros(S, file_num);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
131 LENGTH_ALL=zeros(S,file_num);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
132 LEN_EXON=zeros(S, file_num);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
133
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
134 for j=1:file_num,
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
135 for i=1:S
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
136 T=RESULTS{j}{i};
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
137 if isempty(T)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
138 continue
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
139 else
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
140 READCOUNTS_ALL(i,j)=T{2};
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
141 READCOUNTS_EXON(i,j)=T{3};
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
142 LENGTH_ALL(i,j)=T{5}(1);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
143 LEN_EXON(i,j)=T{5}(2);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
144 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
145 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
146 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
147
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
148 %%%%% write results for all bam files
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
149 fid_conditions = fopen(sprintf('%s_CONDITIONS.tab', outfile), 'w');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
150 fid_counts = fopen(sprintf('%s_COUNTS.tab', outfile) ,'w');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
151 fprintf(fid_counts,'gene');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
152 fprintf(fid_conditions, 'file\tcondition\treplicate\n');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
153 for j = 1:length(file_list)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
154 fname = file_list{j} ;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
155 fname = separate(fname, '/');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
156 fname = fname{end};
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
157 fname = strrep(fname, '.bam', '') ;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
158 fprintf(fid_counts,'\t%s', fname);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
159 fprintf(fid_conditions, '%s\t%i\t%i\n', fname, file_cond_ids(j), file_rep_ids(j));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
160 end;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
161 fprintf(fid_counts,'\n') ;
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
162
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
163 for i = 1:size(genes,2)
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
164 fprintf(fid_counts,'%s',genes(i).name);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
165 for j = 1:length(file_list),
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
166 fprintf(fid_counts,'\t%i', READCOUNTS_EXON(i,j));
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
167 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
168 fprintf(fid_counts,'\n');
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
169 end
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
170 fclose(fid_counts);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
171 fclose(fid_conditions);
2fe512c7bfdf DESeq2 version 1.0.19 added to the repo
vipints <vipin@cbio.mskcc.org>
parents:
diff changeset
172 exit;