Mercurial > repos > vipints > rdiff

comparison rDiff/src/tools/read_utils/get_reads_for_gene.m @ 0:0f80a5141704

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 [reads1] = get_reads_for_gene(CFG,gene)
2
3
4 % Get the reads from the bam-file
5 if strcmp(gene.strand,'-')
6 [mask1, read_intron_list] = get_reads(CFG.curr_bamfile, [CFG.chr_prefix gene.chr],gene.start+1,gene.stop+1, '0');
7 else
8 [mask1, read_intron_list] = get_reads(CFG.curr_bamfile, [CFG.chr_prefix gene.chr],gene.start,gene.stop, '0');
9 end
10
11 % Bring the reads into a matrix form
12 if isempty(mask1)
13 reads1=zeros(0,gene.stop-gene.start+1);
14 return
15 else
16 reads1=sparse(mask1(1,:)',mask1(2,:)',ones(size(mask1,2),1),max(mask1(1,:)),gene.stop-gene.start+1);
17 end
18
19 % remove reads which are shorter than CFG.min_read_length
20 reads1=reads1(sum(reads1,2)>CFG.min_read_length,:);
21
22 if isempty(reads1)
23 return
24 end
25 %remove reads which could stem from other genes
26 [reads1,FLAG]=remove_reads_from_other_genes(reads1,gene);
27
28 %Subsample reads
29 if isempty(reads1)
30 return
31 end
32 if CFG.rDiff_subsample>0
33 if size(reads1,1)>CFG.rDiff_subsample
34 RP=randperm(size(reads1,1));
35 reads1=reads1(RP(1:CFG.rDiff_subsample),:);
36 end
37 end
38
39 %Clip reads
40 if isempty(reads1)
41 return
42 end
43 if CFG.bases_to_clip>0
44 CFG.bases_to_clip=3;
45 [reads1]=clip_reads(reads1,CFG.bases_to_clip);
46 end
47