annotate merge_and_filter.r @ 17:b95fa7e426c3 draft

Uploaded
author davidvanzessen
date Thu, 24 Nov 2016 10:24:19 -0500
parents 4e596473c25c
children 949a30f04d9b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
1 args <- commandArgs(trailingOnly = TRUE)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
2
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
3
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
4 summaryfile = args[1]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
5 sequencesfile = args[2]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
6 mutationanalysisfile = args[3]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
7 mutationstatsfile = args[4]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
8 hotspotsfile = args[5]
14
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
9 aafile = args[6]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
10 gene_identification_file= args[7]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
11 output = args[8]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
12 before.unique.file = args[9]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
13 unmatchedfile = args[10]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
14 method=args[11]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
15 functionality=args[12]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
16 unique.type=args[13]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
17 filter.unique=args[14]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
18 class.filter=args[15]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
19 empty.region.filter=args[16]
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
20
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
21 summ = read.table(summaryfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
22 sequences = read.table(sequencesfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
23 mutationanalysis = read.table(mutationanalysisfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
24 mutationstats = read.table(mutationstatsfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
25 hotspots = read.table(hotspotsfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
14
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
26 AAs = read.table(aafile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
27 gene_identification = read.table(gene_identification_file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
28
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
29 if(method == "blastn"){
10
4b695ca65213 Uploaded
davidvanzessen
parents: 9
diff changeset
30 #"qseqid\tsseqid\tpident\tlength\tmismatch\tgapopen\tqstart\tqend\tsstart\tsend\tevalue\tbitscore"
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
31 gene_identification = gene_identification[!duplicated(gene_identification$qseqid),]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
32 ref_length = data.frame(sseqid=c("ca1", "ca2", "cg1", "cg2", "cg3", "cg4", "cm"), ref.length=c(81,81,141,141,141,141,52))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
33 gene_identification = merge(gene_identification, ref_length, by="sseqid", all.x=T)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
34 gene_identification$chunk_hit_percentage = (gene_identification$length / gene_identification$ref.length) * 100
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
35 gene_identification = gene_identification[,c("qseqid", "chunk_hit_percentage", "pident", "qstart", "sseqid")]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
36 colnames(gene_identification) = c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
37 }
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
38
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
39 input.sequence.count = nrow(summ)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
40 print(paste("Number of sequences in summary file:", input.sequence.count))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
41
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
42 filtering.steps = data.frame(character(0), numeric(0))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
43
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
44 filtering.steps = rbind(filtering.steps, c("Input", input.sequence.count))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
45
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
46 filtering.steps[,1] = as.character(filtering.steps[,1])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
47 filtering.steps[,2] = as.character(filtering.steps[,2])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
48 #filtering.steps[,3] = as.numeric(filtering.steps[,3])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
49
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
50 summ = merge(summ, gene_identification, by="Sequence.ID")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
51
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
52 summ = summ[summ$Functionality != "No results",]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
53
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
54 print(paste("Number of sequences after 'No results' filter:", nrow(summ)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
55
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
56 filtering.steps = rbind(filtering.steps, c("After 'No results' filter", nrow(summ)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
57
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
58 if(functionality == "productive"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
59 summ = summ[summ$Functionality == "productive (see comment)" | summ$Functionality == "productive",]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
60 } else if (functionality == "unproductive"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
61 summ = summ[summ$Functionality == "unproductive (see comment)" | summ$Functionality == "unproductive",]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
62 } else if (functionality == "remove_unknown"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
63 summ = summ[summ$Functionality != "No results" & summ$Functionality != "unknown (see comment)" & summ$Functionality != "unknown",]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
64 }
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
65
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
66 print(paste("Number of sequences after functionality filter:", nrow(summ)))
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
67
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
68 filtering.steps = rbind(filtering.steps, c("After functionality filter", nrow(summ)))
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
69
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
70 result = merge(summ, mutationanalysis[,!(names(mutationanalysis) %in% names(summ)[-1])], by="Sequence.ID")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
71
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
72 print(paste("Number of sequences after merging with mutation analysis file:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
73
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
74 result = merge(result, mutationstats[,!(names(mutationstats) %in% names(result)[-1])], by="Sequence.ID")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
75
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
76 print(paste("Number of sequences after merging with mutation stats file:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
77
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
78 result = merge(result, hotspots[,!(names(hotspots) %in% names(result)[-1])], by="Sequence.ID")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
79
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
80 print(paste("Number of sequences after merging with hotspots file:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
81
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
82 sequences = sequences[,c("Sequence.ID", "FR1.IMGT", "CDR1.IMGT", "FR2.IMGT", "CDR2.IMGT", "FR3.IMGT", "CDR3.IMGT")]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
83 names(sequences) = c("Sequence.ID", "FR1.IMGT.seq", "CDR1.IMGT.seq", "FR2.IMGT.seq", "CDR2.IMGT.seq", "FR3.IMGT.seq", "CDR3.IMGT.seq")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
84 result = merge(result, sequences, by="Sequence.ID", all.x=T)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
85
14
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
86 AAs = AAs[,c("Sequence.ID", "CDR3.IMGT")]
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
87 names(AAs) = c("Sequence.ID", "CDR3.IMGT.AA")
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
88 result = merge(result, AAs, by="Sequence.ID", all.x=T)
59765d2c8890 Uploaded
davidvanzessen
parents: 13
diff changeset
89
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
90 print(paste("Number of sequences in result after merging with sequences:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
91
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
92 result$VGene = gsub("^Homsap ", "", result$V.GENE.and.allele)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
93 result$VGene = gsub("[*].*", "", result$VGene)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
94 result$DGene = gsub("^Homsap ", "", result$D.GENE.and.allele)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
95 result$DGene = gsub("[*].*", "", result$DGene)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
96 result$JGene = gsub("^Homsap ", "", result$J.GENE.and.allele)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
97 result$JGene = gsub("[*].*", "", result$JGene)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
98
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
99 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
100
12
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
101 splt = strsplit(class.filter, "_")[[1]]
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
102 chunk_hit_threshold = as.numeric(splt[1])
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
103 nt_hit_threshold = as.numeric(splt[2])
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
104
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
105 higher_than=(result$chunk_hit_percentage >= chunk_hit_threshold & result$nt_hit_percentage >= nt_hit_threshold)
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
106
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
107 if(!all(higher_than, na.rm=T)){ #check for no unmatched
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
108 result[!higher_than,"best_match"] = paste("unmatched,", result[!higher_than,"best_match"])
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
109 }
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
110
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
111 if(class.filter == "101_101"){
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
112 result$best_match = "all"
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
113 }
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
114
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
115 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
116
12
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
117 write.table(x=result, file=gsub("merged.txt$", "before_filters.txt", output), sep="\t",quote=F,row.names=F,col.names=T)
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
118
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
119 print(paste("Number of empty CDR1 sequences:", sum(result$CDR1.IMGT.seq == "")))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
120 print(paste("Number of empty FR2 sequences:", sum(result$FR2.IMGT.seq == "")))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
121 print(paste("Number of empty CDR2 sequences:", sum(result$CDR2.IMGT.seq == "")))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
122 print(paste("Number of empty FR3 sequences:", sum(result$FR3.IMGT.seq == "")))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
123
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
124 if(empty.region.filter == "leader"){
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
125 result = result[result$FR1.IMGT.seq != "" & result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
126 print(paste("Number of sequences after empty FR1, CDR1, FR2, CDR2 and FR3 column filter:", nrow(result)))
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
127 filtering.steps = rbind(filtering.steps, c("After empty FR1, CDR1, FR2, CDR2, FR3 filter", nrow(result)))
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
128 } else if(empty.region.filter == "FR1"){
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
129 result = result[result$CDR1.IMGT.seq != "" & result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
130 print(paste("Number of sequences after empty CDR1, FR2, CDR2 and FR3 column filter:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
131 filtering.steps = rbind(filtering.steps, c("After empty CDR1, FR2, CDR2, FR3 filter", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
132 } else if(empty.region.filter == "CDR1"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
133 result = result[result$FR2.IMGT.seq != "" & result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
134 print(paste("Number of sequences after empty FR2, CDR2 and FR3 column filter:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
135 filtering.steps = rbind(filtering.steps, c("After empty FR2, CDR2, FR3 filter", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
136 } else if(empty.region.filter == "FR2"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
137 result = result[result$CDR2.IMGT.seq != "" & result$FR3.IMGT.seq != "", ]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
138 print(paste("Number of sequences after empty CDR2 and FR3 column filter:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
139 filtering.steps = rbind(filtering.steps, c("After empty CDR2, FR3 filter", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
140 }
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
141
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
142 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
143
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
144 if(empty.region.filter == "leader"){
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
145 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)),]
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
146 } else if(empty.region.filter == "FR1"){
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
147 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)),]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
148 } else if(empty.region.filter == "CDR1"){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
149 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)),]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
150 } else if(empty.region.filter == "FR2"){
2
e85fec274cde Uploaded
davidvanzessen
parents: 1
diff changeset
151 result = result[!(grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq) | grepl("n|N", result$CDR3.IMGT.seq)),]
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
152 }
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
153
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
154 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
155 print(result[result$Sequence.ID == "JY8QFUQ01BNS72",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
156
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
157 print(paste("Number of sequences in result after n filtering:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
158 filtering.steps = rbind(filtering.steps, c("After N filter", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
159
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
160 cleanup_columns = c("FR1.IMGT.Nb.of.mutations",
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
161 "CDR1.IMGT.Nb.of.mutations",
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
162 "FR2.IMGT.Nb.of.mutations",
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
163 "CDR2.IMGT.Nb.of.mutations",
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
164 "FR3.IMGT.Nb.of.mutations")
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
165
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
166 for(col in cleanup_columns){
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
167 result[,col] = gsub("\\(.*\\)", "", result[,col])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
168 result[,col] = as.numeric(result[,col])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
169 result[is.na(result[,col]),] = 0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
170 }
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
171
5
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
172 write.table(result, before.unique.file, sep="\t", quote=F,row.names=F,col.names=T)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
173
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
174 if(filter.unique != "no"){
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
175 clmns = names(result)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
176
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
177 if(empty.region.filter == "leader"){
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
178 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)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
179 } else if(empty.region.filter == "FR1"){
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
180 result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
181 } else if(empty.region.filter == "CDR1"){
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
182 result$unique.def = paste(result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
183 } else if(empty.region.filter == "FR2"){
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
184 result$unique.def = paste(result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
185 }
12
6b66c1c57f22 Uploaded
davidvanzessen
parents: 10
diff changeset
186
15
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
187 if(filter.unique == "remove"){
5
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
188 result = result[duplicated(result$unique.def) | duplicated(result$unique.def, fromLast=T),]
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
189 }
15
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
190 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
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
191
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
192 result = result[!duplicated(result$unique.def),]
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
193 }
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
194
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
195 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
196 print(result[result$Sequence.ID == "JY8QFUQ01BNS72",c("Sequence.ID","best_match")])
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
197
5
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
198 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)
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
199
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
200 filtering.steps = rbind(filtering.steps, c("After filter unique sequences", nrow(result)))
012a738edf5a Uploaded
davidvanzessen
parents: 3
diff changeset
201
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
202 if(nrow(summ) == 0){
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
203 stop("No data remaining after filter")
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
204 }
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
205
13
933fb21568ce Uploaded
davidvanzessen
parents: 12
diff changeset
206 result$best_match_class = gsub(",.*", "", result$best_match) #gsub so the unmatched don't have a class after it
933fb21568ce Uploaded
davidvanzessen
parents: 12
diff changeset
207
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
208 result$past = do.call(paste, c(result[unlist(strsplit(unique.type, ","))], sep = ":"))
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
209
15
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
210
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
211
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
212
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
213 result.matched = result[!grepl("unmatched", result$best_match),]
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
214 result.unmatched = result[grepl("unmatched", result$best_match),]
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
215
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
216 result = rbind(result.matched, result.unmatched)
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
217
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
218 result = result[!(duplicated(result$past)), ]
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
219
15
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
220
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
221
16
4e596473c25c Uploaded
davidvanzessen
parents: 15
diff changeset
222 print(result[result$Sequence.ID == "JY8QFUQ01C310D",c("Sequence.ID","best_match")])
15
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
223
61d0a6318711 Uploaded
davidvanzessen
parents: 14
diff changeset
224
13
933fb21568ce Uploaded
davidvanzessen
parents: 12
diff changeset
225 result = result[,!(names(result) %in% c("past", "best_match_class"))]
1
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
226
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
227 print(paste("Number of sequences in result after", unique.type, "filtering:", nrow(result)))
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
228
faae21ba5c63 Uploaded
davidvanzessen
parents: 0
diff changeset
229 filtering.steps = rbind(filtering.steps, c("After remove duplicates based on filter", nrow(result)))
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
230
10
4b695ca65213 Uploaded
davidvanzessen
parents: 9
diff changeset
231 unmatched = result[grepl("^unmatched", result$best_match),c("Sequence.ID", "chunk_hit_percentage", "nt_hit_percentage", "start_locations", "best_match")]
4b695ca65213 Uploaded
davidvanzessen
parents: 9
diff changeset
232
0
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
233 print(paste("Number of rows in result:", nrow(result)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
234 print(paste("Number of rows in unmatched:", nrow(unmatched)))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
235
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
236 matched.sequences = result[!grepl("^unmatched", result$best_match),]
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
237
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
238 write.table(x=matched.sequences, file=gsub("merged.txt$", "filtered.txt", output), sep="\t",quote=F,row.names=F,col.names=T)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
239
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
240 matched.sequences.count = nrow(matched.sequences)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
241 unmatched.sequences.count = sum(grepl("^unmatched", result$best_match))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
242
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
243 filtering.steps = rbind(filtering.steps, c("Number of matched sequences", matched.sequences.count))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
244 filtering.steps = rbind(filtering.steps, c("Number of unmatched sequences", unmatched.sequences.count))
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
245 filtering.steps[,2] = as.numeric(filtering.steps[,2])
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
246 filtering.steps$perc = round(filtering.steps[,2] / input.sequence.count * 100, 2)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
247
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
248 write.table(x=filtering.steps, file=gsub("unmatched", "filtering_steps", unmatchedfile), sep="\t",quote=F,row.names=F,col.names=F)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
249
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
250 write.table(x=result, file=output, sep="\t",quote=F,row.names=F,col.names=T)
c33d93683a09 Uploaded
davidvanzessen
parents:
diff changeset
251 write.table(x=unmatched, file=unmatchedfile, sep="\t",quote=F,row.names=F,col.names=T)