Mercurial > repos > vipints > deseq_hts
annotate deseq-hts_2.0/src/get_read_counts.m @ 10:2fe512c7bfdf draft
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 |
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; |