Mercurial > repos > davidvanzessen > shm_csr

annotate merge_and_filter.r @ 96:385dea3c6cb5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 423a48569c69301fdbf893ac3a649128404dfff5
author rhpvorderman
date Fri, 05 Jan 2024 08:53:22 +0000
parents 8fcf31272f6e
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
1 args <- commandArgs(trailingOnly = TRUE)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
2
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
3
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
4 summaryfile = args[1]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
5 sequencesfile = args[2]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
6 mutationanalysisfile = args[3]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
7 mutationstatsfile = args[4]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
8 hotspotsfile = args[5]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
9 aafile = args[6]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
10 gene_identification_file= args[7]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
11 output = args[8]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
12 before.unique.file = args[9]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
13 unmatchedfile = args[10]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
14 method=args[11]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
15 functionality=args[12]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
16 unique.type=args[13]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
17 filter.unique=args[14]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
18 filter.unique.count=as.numeric(args[15])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
19 class.filter=args[16]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
20 empty.region.filter=args[17]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
21
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
22 print(paste("filter.unique.count:", filter.unique.count))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
23
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
24 summ = read.table(summaryfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
25 sequences = read.table(sequencesfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
26 mutationanalysis = read.table(mutationanalysisfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
27 mutationstats = read.table(mutationstatsfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
28 hotspots = read.table(hotspotsfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
29 AAs = read.table(aafile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
30 gene_identification = read.table(gene_identification_file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
31
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
32 fix_column_names = function(df){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
33 if("V.DOMAIN.Functionality" %in% names(df)){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
34 names(df)[names(df) == "V.DOMAIN.Functionality"] = "Functionality"
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
35 print("found V.DOMAIN.Functionality, changed")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
36 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
37 if("V.DOMAIN.Functionality.comment" %in% names(df)){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
38 names(df)[names(df) == "V.DOMAIN.Functionality.comment"] = "Functionality.comment"
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
39 print("found V.DOMAIN.Functionality.comment, changed")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
40 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
41 return(df)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
42 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
43
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
44 fix_non_unique_ids = function(df){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
45 df$Sequence.ID = paste(df$Sequence.ID, 1:nrow(df))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
46 return(df)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
47 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
48
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
49 summ = fix_column_names(summ)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
50 sequences = fix_column_names(sequences)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
51 mutationanalysis = fix_column_names(mutationanalysis)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
52 mutationstats = fix_column_names(mutationstats)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
53 hotspots = fix_column_names(hotspots)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
54 AAs = fix_column_names(AAs)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
55
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
56 if(method == "blastn"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
57 #"qseqid\tsseqid\tpident\tlength\tmismatch\tgapopen\tqstart\tqend\tsstart\tsend\tevalue\tbitscore"
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
58 gene_identification = gene_identification[!duplicated(gene_identification$qseqid),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
59 ref_length = data.frame(sseqid=c("ca1", "ca2", "cg1", "cg2", "cg3", "cg4", "cm"), ref.length=c(81,81,141,141,141,141,52))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
60 gene_identification = merge(gene_identification, ref_length, by="sseqid", all.x=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
61 gene_identification$chunk_hit_percentage = (gene_identification$length / gene_identification$ref.length) * 100
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
62 gene_identification = gene_identification[,c("qseqid", "chunk_hit_percentage", "pident", "qstart", "sseqid")]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
63 colnames(gene_identification) = c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
64 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
65
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
66 #print("Summary analysis files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
67 #print(names(summ))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
68
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
69
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
70
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
71 input.sequence.count = nrow(summ)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
72 print(paste("Number of sequences in summary file:", input.sequence.count))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
73
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
74 filtering.steps = data.frame(character(0), numeric(0))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
75
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
76 filtering.steps = rbind(filtering.steps, c("Input", input.sequence.count))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
77
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
78 filtering.steps[,1] = as.character(filtering.steps[,1])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
79 filtering.steps[,2] = as.character(filtering.steps[,2])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
80 #filtering.steps[,3] = as.numeric(filtering.steps[,3])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
81
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
82 #print("summary files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
83 #print(names(summ))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
84
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
85 summ = merge(summ, gene_identification, by="Sequence.ID")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
86
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
87 print(paste("Number of sequences after merging with gene identification:", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
88
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
89 summ = summ[summ$Functionality != "No results",]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
90
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
91 print(paste("Number of sequences after 'No results' filter:", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
92
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
93 filtering.steps = rbind(filtering.steps, c("After 'No results' filter", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
94
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
95 if(functionality == "productive"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
96 summ = summ[summ$Functionality == "productive (see comment)" | summ$Functionality == "productive",]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
97 } else if (functionality == "unproductive"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
98 summ = summ[summ$Functionality == "unproductive (see comment)" | summ$Functionality == "unproductive",]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
99 } else if (functionality == "remove_unknown"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
100 summ = summ[summ$Functionality != "No results" & summ$Functionality != "unknown (see comment)" & summ$Functionality != "unknown",]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
101 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
102
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
103 print(paste("Number of sequences after functionality filter:", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
104
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
105 filtering.steps = rbind(filtering.steps, c("After functionality filter", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
106
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
107 if(F){ #to speed up debugging
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
108 set.seed(1)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
109 summ = summ[sample(nrow(summ), floor(nrow(summ) * 0.03)),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
110 print(paste("Number of sequences after sampling 3%:", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
111
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
112 filtering.steps = rbind(filtering.steps, c("Number of sequences after sampling 3%", nrow(summ)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
113 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
114
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
115 print("mutation analysis files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
116 print(names(mutationanalysis[,!(names(mutationanalysis) %in% names(summ)[-1])]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
117
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
118 result = merge(summ, mutationanalysis[,!(names(mutationanalysis) %in% names(summ)[-1])], by="Sequence.ID")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
119
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
120 print(paste("Number of sequences after merging with mutation analysis file:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
121
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
122 #print("mutation stats files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
123 #print(names(mutationstats[,!(names(mutationstats) %in% names(result)[-1])]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
124
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
125 result = merge(result, mutationstats[,!(names(mutationstats) %in% names(result)[-1])], by="Sequence.ID")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
126
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
127 print(paste("Number of sequences after merging with mutation stats file:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
128
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
129 print("hotspots files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
130 print(names(hotspots[,!(names(hotspots) %in% names(result)[-1])]))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
131
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
132 result = merge(result, hotspots[,!(names(hotspots) %in% names(result)[-1])], by="Sequence.ID")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
133
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
134 print(paste("Number of sequences after merging with hotspots file:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
135
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
136 print("sequences files columns")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
137 print(c("FR1.IMGT", "CDR1.IMGT", "FR2.IMGT", "CDR2.IMGT", "FR3.IMGT", "CDR3.IMGT"))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
138
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
139 sequences = sequences[,c("Sequence.ID", "FR1.IMGT", "CDR1.IMGT", "FR2.IMGT", "CDR2.IMGT", "FR3.IMGT", "CDR3.IMGT")]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
140 names(sequences) = c("Sequence.ID", "FR1.IMGT.seq", "CDR1.IMGT.seq", "FR2.IMGT.seq", "CDR2.IMGT.seq", "FR3.IMGT.seq", "CDR3.IMGT.seq")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
141 result = merge(result, sequences, by="Sequence.ID", all.x=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
142
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
143 AAs = AAs[,c("Sequence.ID", "CDR3.IMGT")]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
144 names(AAs) = c("Sequence.ID", "CDR3.IMGT.AA")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
145 result = merge(result, AAs, by="Sequence.ID", all.x=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
146
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
147 print(paste("Number of sequences in result after merging with sequences:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
148
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
149 result$VGene = gsub("^Homsap ", "", result$V.GENE.and.allele)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
150 result$VGene = gsub("[*].*", "", result$VGene)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
151 result$DGene = gsub("^Homsap ", "", result$D.GENE.and.allele)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
152 result$DGene = gsub("[*].*", "", result$DGene)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
153 result$JGene = gsub("^Homsap ", "", result$J.GENE.and.allele)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
154 result$JGene = gsub("[*].*", "", result$JGene)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
155
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
156 splt = strsplit(class.filter, "_")[[1]]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
157 chunk_hit_threshold = as.numeric(splt[1])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
158 nt_hit_threshold = as.numeric(splt[2])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
159
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
160 higher_than=(result$chunk_hit_percentage >= chunk_hit_threshold & result$nt_hit_percentage >= nt_hit_threshold)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
161
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
162 if(!all(higher_than, na.rm=T)){ #check for no unmatched
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
163 result[!higher_than,"best_match"] = paste("unmatched,", result[!higher_than,"best_match"])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
164 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
165
93
8fcf31272f6e planemo upload commit a43893724cc769bed8a1f19a5b19ec1ba20cb63c
rhpvorderman
parents: 83
diff changeset
166 if(splt[1] == "101" & splt[2] == "101"){
8fcf31272f6e planemo upload commit a43893724cc769bed8a1f19a5b19ec1ba20cb63c
rhpvorderman
parents: 83
diff changeset
167 result$best_match = splt[3]
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
168 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
169
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
170 write.table(x=result, file=gsub("merged.txt$", "before_filters.txt", output), sep="\t",quote=F,row.names=F,col.names=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
171
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
172 print(paste("Number of empty CDR1 sequences:", sum(result$CDR1.IMGT.seq == "", na.rm=T)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
173 print(paste("Number of empty FR2 sequences:", sum(result$FR2.IMGT.seq == "", na.rm=T)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
174 print(paste("Number of empty CDR2 sequences:", sum(result$CDR2.IMGT.seq == "", na.rm=T)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
175 print(paste("Number of empty FR3 sequences:", sum(result$FR3.IMGT.seq == "", na.rm=T)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
176
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
177 if(empty.region.filter == "leader"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
178 result = result[result$FR1.IMGT.seq != "" & result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
179 } else if(empty.region.filter == "FR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
180 result = result[result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
181 } else if(empty.region.filter == "CDR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
182 result = result[result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
183 } else if(empty.region.filter == "FR2"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
184 result = result[result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
185 }
96
385dea3c6cb5 planemo upload commit 423a48569c69301fdbf893ac3a649128404dfff5
rhpvorderman
parents: 93
diff changeset
186 # If empty region filter is None, nothing happens.
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
187
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
188 print(paste("After removal sequences that are missing a gene region:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
189 filtering.steps = rbind(filtering.steps, c("After removal sequences that are missing a gene region", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
190
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
191 if(empty.region.filter == "leader"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
192 result = result[!(grepl("n|N", result$FR1.IMGT.seq) | grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
193 } else if(empty.region.filter == "FR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
194 result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
195 } else if(empty.region.filter == "CDR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
196 result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
197 } else if(empty.region.filter == "FR2"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
198 result = result[!(grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
199 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
200
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
201 print(paste("Number of sequences in result after n filtering:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
202 filtering.steps = rbind(filtering.steps, c("After N filter", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
203
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
204 cleanup_columns = c("FR1.IMGT.Nb.of.mutations",
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
205 "CDR1.IMGT.Nb.of.mutations",
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
206 "FR2.IMGT.Nb.of.mutations",
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
207 "CDR2.IMGT.Nb.of.mutations",
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
208 "FR3.IMGT.Nb.of.mutations")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
209
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
210 for(col in cleanup_columns){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
211 result[,col] = gsub("\\(.*\\)", "", result[,col])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
212 result[,col] = as.numeric(result[,col])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
213 result[is.na(result[,col]),] = 0
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
214 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
215
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
216 write.table(result, before.unique.file, sep="\t", quote=F,row.names=F,col.names=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
217
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
218
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
219 if(filter.unique != "no"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
220 clmns = names(result)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
221 if(filter.unique == "remove_vjaa"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
222 result$unique.def = paste(result$VGene, result$JGene, result$CDR3.IMGT.AA)
96
385dea3c6cb5 planemo upload commit 423a48569c69301fdbf893ac3a649128404dfff5
rhpvorderman
parents: 93
diff changeset
223 } else if(empty.region.filter == "leader" || empty.region.filter == "None"){
83
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
224 result$unique.def = paste(result$FR1.IMGT.seq, result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
225 } else if(empty.region.filter == "FR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
226 result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
227 } else if(empty.region.filter == "CDR1"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
228 result$unique.def = paste(result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
229 } else if(empty.region.filter == "FR2"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
230 result$unique.def = paste(result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
231 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
232
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
233 if(grepl("remove", filter.unique)){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
234 result = result[duplicated(result$unique.def) | duplicated(result$unique.def, fromLast=T),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
235 unique.defs = data.frame(table(result$unique.def))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
236 unique.defs = unique.defs[unique.defs$Freq >= filter.unique.count,]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
237 result = result[result$unique.def %in% unique.defs$Var1,]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
238 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
239
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
240 if(filter.unique != "remove_vjaa"){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
241 result$unique.def = paste(result$unique.def, gsub(",.*", "", result$best_match)) #keep the unique sequences that are in multiple classes, gsub so the unmatched don't have a class after it
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
242 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
243
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
244 result = result[!duplicated(result$unique.def),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
245 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
246
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
247 write.table(result, gsub("before_unique_filter.txt", "after_unique_filter.txt", before.unique.file), sep="\t", quote=F,row.names=F,col.names=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
248
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
249 filtering.steps = rbind(filtering.steps, c("After filter unique sequences", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
250
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
251 print(paste("Number of sequences in result after unique filtering:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
252
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
253 if(nrow(summ) == 0){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
254 stop("No data remaining after filter")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
255 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
256
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
257 result$best_match_class = gsub(",.*", "", result$best_match) #gsub so the unmatched don't have a class after it
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
258
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
259 #result$past = ""
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
260 #cls = unlist(strsplit(unique.type, ","))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
261 #for (i in 1:nrow(result)){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
262 # result[i,"past"] = paste(result[i,cls], collapse=":")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
263 #}
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
264
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
265
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
266
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
267 result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":"))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
268
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
269 result.matched = result[!grepl("unmatched", result$best_match),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
270 result.unmatched = result[grepl("unmatched", result$best_match),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
271
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
272 result = rbind(result.matched, result.unmatched)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
273
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
274 result = result[!(duplicated(result$past)), ]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
275
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
276 result = result[,!(names(result) %in% c("past", "best_match_class"))]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
277
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
278 print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
279
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
280 filtering.steps = rbind(filtering.steps, c("After remove duplicates based on filter", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
281
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
282 unmatched = result[grepl("^unmatched", result$best_match),c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
283
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
284 print(paste("Number of rows in result:", nrow(result)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
285 print(paste("Number of rows in unmatched:", nrow(unmatched)))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
286
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
287 matched.sequences = result[!grepl("^unmatched", result$best_match),]
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
288
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
289 write.table(x=matched.sequences, file=gsub("merged.txt$", "filtered.txt", output), sep="\t",quote=F,row.names=F,col.names=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
290
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
291 matched.sequences.count = nrow(matched.sequences)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
292 unmatched.sequences.count = sum(grepl("^unmatched", result$best_match))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
293 if(matched.sequences.count <= unmatched.sequences.count){
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
294 print("WARNING NO MATCHED (SUB)CLASS SEQUENCES!!")
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
295 }
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
296
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
297 filtering.steps = rbind(filtering.steps, c("Number of matched sequences", matched.sequences.count))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
298 filtering.steps = rbind(filtering.steps, c("Number of unmatched sequences", unmatched.sequences.count))
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
299 filtering.steps[,2] = as.numeric(filtering.steps[,2])
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
300 filtering.steps$perc = round(filtering.steps[,2] / input.sequence.count * 100, 2)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
301
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
302 write.table(x=filtering.steps, file=gsub("unmatched", "filtering_steps", unmatchedfile), sep="\t",quote=F,row.names=F,col.names=F)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
303
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
304 write.table(x=result, file=output, sep="\t",quote=F,row.names=F,col.names=T)
729738462297 "planemo upload commit c0ffc68aec5836d5b20b543106493056a87edf57"
rhpvorderman
parents: 81
diff changeset
305 write.table(x=unmatched, file=unmatchedfile, sep="\t",quote=F,row.names=F,col.names=T)