-
Notifications
You must be signed in to change notification settings - Fork 1
/
generateProbeList.m
68 lines (59 loc) · 2.31 KB
/
generateProbeList.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
function probeList=generateProbeList(adapterList,probeHeader,probeSequence,params)
% params = struct('species','Mouse','verbose',1,'number',48,...
% 'specialTranscripts','C:\FISHerMan\Db\Mouse.STList.fas');
%% Remove transcripts without enough probes
if params(1).verbose
disp('generating the list of probes to order');
disp(' removing transcripts without enough probes');
end
pos = regexp(probeHeader, ':');
trimHeader = {};
for n = 1:length(probeHeader)
trimHeader{end+1,1} = probeHeader{n,1}(1:pos{n,1}(1)-1);
end
uniqueHeader = unique(trimHeader, 'stable');
indexTotal = zeros(length(trimHeader),1);
for n = 1:length(uniqueHeader)
index = ismember(trimHeader, uniqueHeader{n,1});
if sum(index) < params(1).number &&...
checkSpecialTranscripts(uniqueHeader{n,1},params) % check for Bin's special sequences
indexTotal = indexTotal+index;
if params(1).verbose
disp([' transcript ' uniqueHeader{n,1} ...
' has less than ' num2str(params(1).number) ' probes']);
end
end
end
indexTotal = logical(indexTotal);
probeHeader(indexTotal) = [];
probeSequence(indexTotal) = [];
trimHeader(indexTotal) = [];
%% Trim the number of probes of each transcript to be 48
uniqueHeader = unique(trimHeader, 'stable');
indexTotal = [];
for n = 1:length(uniqueHeader)
index = ismember(trimHeader, uniqueHeader{n,1});
if sum(index) > params(1).number
index=find(index>0);
index=index(params(1).number+1:end);
indexTotal=[indexTotal;index];
end
end
probeHeader(indexTotal) = [];
probeSequence(indexTotal) = [];
trimHeader(indexTotal) = [];
% disp(' randomizing and saving the list of probes');
% indexTotal = randperm(length(probeHeader))';
% probeHeader = probeHeader(indexTotal);
% probeSequence = probeSequence(indexTotal);
probeList = [params(1).species '.probes.nr.txt'];
if exist(probeList, 'file')
delete(probeList);
end
fastawrite(probeList,probeHeader,probeSequence);
%% Check how many transcripts are left after this step of screening
[geneNumLeft,geneNumDelete]=checkTranscriptsLeft(adapterList,probeHeader);
if params(1).verbose
disp([num2str(geneNumDelete) ' out of ' num2str(geneNumLeft+geneNumDelete)...
' FISH escaped FISHerMan''s net']);
end