Mercurial > repos > davidvanzessen > argalaxy_tools

comparison imgt_loader/imgt_loader.r @ 40:106275b54470 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 30 May 2017 07:26:33 -0400
parents bad6a9a53ae7
children 4d2a8f98a502
comparison
equal deleted inserted replaced
39:bad6a9a53ae7 40:106275b54470
1 args <- commandArgs(trailingOnly = TRUE) 1 args <- commandArgs(trailingOnly = TRUE)
2 options(show.error.locations = TRUE)
2 3
3 summ.file = args[1] 4 summ.file = args[1]
4 sequences.file = args[2] 5 sequences.file = args[2]
5 aa.file = args[3] 6 aa.file = args[3]
6 junction.file = args[4] 7 junction.file = args[4]
7 gapped.aa.file = args[5] 8 gapped.aa.file = args[5]
8 out.file = args[6] 9 out.file = args[6]
9 10
11 print(summ.file)
12 print(out.file)
13
10 summ = read.table(summ.file, sep="\t", header=T, quote="", fill=T) 14 summ = read.table(summ.file, sep="\t", header=T, quote="", fill=T)
11 sequences = read.table(sequences.file, sep="\t", header=T, quote="", fill=T) 15 sequences = read.table(sequences.file, sep="\t", header=T, quote="", fill=T)
12 aa = read.table(aa.file, sep="\t", header=T, quote="", fill=T) 16 aa = read.table(aa.file, sep="\t", header=T, quote="", fill=T)
13 gapped.aa = read.table(gapped.aa.file, sep="\t", header=T, quote="", fill=T) 17 gapped.aa = read.table(gapped.aa.file, sep="\t", header=T, quote="", fill=T)
14 junction = read.table(junction.file, sep="\t", header=T, quote="", fill=T) 18 junction = read.table(junction.file, sep="\t", header=T, quote="", fill=T)
19
20 print(paste("nrow(summ)", nrow(summ)))
21
22 write.table(summ, out.file, sep="\t", quote=F, row.names=F, col.names=T)
15 23
16 fix_column_names = function(df){ 24 fix_column_names = function(df){
17 if("V.DOMAIN.Functionality" %in% names(df)){ 25 if("V.DOMAIN.Functionality" %in% names(df)){
18 names(df)[names(df) == "V.DOMAIN.Functionality"] = "Functionality" 26 names(df)[names(df) == "V.DOMAIN.Functionality"] = "Functionality"
19 print("found V.DOMAIN.Functionality, changed") 27 print("found V.DOMAIN.Functionality, changed")
29 sequences = fix_column_names(sequences) 37 sequences = fix_column_names(sequences)
30 aa = fix_column_names(aa) 38 aa = fix_column_names(aa)
31 gapped.aa = fix_column_names(gapped.aa) 39 gapped.aa = fix_column_names(gapped.aa)
32 junction = fix_column_names(junction) 40 junction = fix_column_names(junction)
33 41
42 print(paste("nrow(summ)", nrow(summ)))
43
34 old_summary_columns=c('Sequence.ID','JUNCTION.frame','V.GENE.and.allele','D.GENE.and.allele','J.GENE.and.allele','CDR1.IMGT.length','CDR2.IMGT.length','CDR3.IMGT.length','Orientation') 44 old_summary_columns=c('Sequence.ID','JUNCTION.frame','V.GENE.and.allele','D.GENE.and.allele','J.GENE.and.allele','CDR1.IMGT.length','CDR2.IMGT.length','CDR3.IMGT.length','Orientation')
35 old_sequence_columns=c('CDR1.IMGT','CDR2.IMGT','CDR3.IMGT') 45 old_sequence_columns=c('CDR1.IMGT','CDR2.IMGT','CDR3.IMGT')
36 old_junction_columns=c('JUNCTION') 46 old_junction_columns=c('JUNCTION')
37 47
38 added_summary_columns=c('Functionality','V.REGION.identity..','V.REGION.identity.nt','D.REGION.reading.frame','AA.JUNCTION','Functionality.comment','Sequence') 48 added_summary_columns=c('Functionality','V.REGION.identity..','V.REGION.identity.nt','D.REGION.reading.frame','AA.JUNCTION','Functionality.comment','Sequence')
40 50
41 added_junction_columns=c('P3.V.nt.nb','N.REGION.nt.nb','N1.REGION.nt.nb','P5.D.nt.nb','P3.D.nt.nb','N2.REGION.nt.nb','P5.J.nt.nb','X3.V.REGION.trimmed.nt.nb','X5.D.REGION.trimmed.nt.nb','X3.D.REGION.trimmed.nt.nb','X5.J.REGION.trimmed.nt.nb','N.REGION','N1.REGION','N2.REGION') 51 added_junction_columns=c('P3.V.nt.nb','N.REGION.nt.nb','N1.REGION.nt.nb','P5.D.nt.nb','P3.D.nt.nb','N2.REGION.nt.nb','P5.J.nt.nb','X3.V.REGION.trimmed.nt.nb','X5.D.REGION.trimmed.nt.nb','X3.D.REGION.trimmed.nt.nb','X5.J.REGION.trimmed.nt.nb','N.REGION','N1.REGION','N2.REGION')
42 added_junction_columns=c(added_junction_columns, 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb') 52 added_junction_columns=c(added_junction_columns, 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb')
43 53
44 out=summ[,c("Sequence.ID","JUNCTION.frame","V.GENE.and.allele","D.GENE.and.allele","J.GENE.and.allele")] 54 out=summ[,c("Sequence.ID","JUNCTION.frame","V.GENE.and.allele","D.GENE.and.allele","J.GENE.and.allele")]
55
56 print(paste("nrow(summ)", nrow(summ)))
57 print(paste("nrow(aa)", nrow(aa)))
45 58
46 out[,"CDR1.Seq"] = aa[,"CDR1.IMGT"] 59 out[,"CDR1.Seq"] = aa[,"CDR1.IMGT"]
47 out[,"CDR1.Length"] = summ[,"CDR1.IMGT.length"] 60 out[,"CDR1.Length"] = summ[,"CDR1.IMGT.length"]
48 61
49 out[,"CDR2.Seq"] = aa[,"CDR2.IMGT"] 62 out[,"CDR2.Seq"] = aa[,"CDR2.IMGT"]
76 out = out[,c('Sequence.ID','JUNCTION.frame','Top V Gene','Top D Gene','Top J Gene','CDR1.Seq','CDR1.Length','CDR2.Seq','CDR2.Length','CDR3.Seq','CDR3.Length','CDR3.Seq.DNA','CDR3.Length.DNA','Strand','CDR3.Found.How','Functionality','V.REGION.identity..','V.REGION.identity.nt','D.REGION.reading.frame','AA.JUNCTION','Functionality.comment','Sequence','FR1.IMGT','FR2.IMGT','FR3.IMGT','CDR3.IMGT','JUNCTION','J.REGION','FR4.IMGT','P3.V.nt.nb','N.REGION.nt.nb','N1.REGION.nt.nb','P5.D.nt.nb','P3.D.nt.nb','N2.REGION.nt.nb','P5.J.nt.nb','X3.V.REGION.trimmed.nt.nb','X5.D.REGION.trimmed.nt.nb','X3.D.REGION.trimmed.nt.nb','X5.J.REGION.trimmed.nt.nb','N.REGION','N1.REGION','N2.REGION', 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb')] 89 out = out[,c('Sequence.ID','JUNCTION.frame','Top V Gene','Top D Gene','Top J Gene','CDR1.Seq','CDR1.Length','CDR2.Seq','CDR2.Length','CDR3.Seq','CDR3.Length','CDR3.Seq.DNA','CDR3.Length.DNA','Strand','CDR3.Found.How','Functionality','V.REGION.identity..','V.REGION.identity.nt','D.REGION.reading.frame','AA.JUNCTION','Functionality.comment','Sequence','FR1.IMGT','FR2.IMGT','FR3.IMGT','CDR3.IMGT','JUNCTION','J.REGION','FR4.IMGT','P3.V.nt.nb','N.REGION.nt.nb','N1.REGION.nt.nb','P5.D.nt.nb','P3.D.nt.nb','N2.REGION.nt.nb','P5.J.nt.nb','X3.V.REGION.trimmed.nt.nb','X5.D.REGION.trimmed.nt.nb','X3.D.REGION.trimmed.nt.nb','X5.J.REGION.trimmed.nt.nb','N.REGION','N1.REGION','N2.REGION', 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb')]
77 90
78 names(out) = c('ID','VDJ Frame','Top V Gene','Top D Gene','Top J Gene','CDR1 Seq','CDR1 Length','CDR2 Seq','CDR2 Length','CDR3 Seq','CDR3 Length','CDR3 Seq DNA','CDR3 Length DNA','Strand','CDR3 Found How','Functionality','V-REGION identity %','V-REGION identity nt','D-REGION reading frame','AA JUNCTION','Functionality comment','Sequence','FR1-IMGT','FR2-IMGT','FR3-IMGT','CDR3-IMGT','JUNCTION','J-REGION','FR4-IMGT','P3V-nt nb','N-REGION-nt nb','N1-REGION-nt nb','P5D-nt nb','P3D-nt nb','N2-REGION-nt nb','P5J-nt nb','3V-REGION trimmed-nt nb','5D-REGION trimmed-nt nb','3D-REGION trimmed-nt nb','5J-REGION trimmed-nt nb','N-REGION','N1-REGION','N2-REGION', 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb') 91 names(out) = c('ID','VDJ Frame','Top V Gene','Top D Gene','Top J Gene','CDR1 Seq','CDR1 Length','CDR2 Seq','CDR2 Length','CDR3 Seq','CDR3 Length','CDR3 Seq DNA','CDR3 Length DNA','Strand','CDR3 Found How','Functionality','V-REGION identity %','V-REGION identity nt','D-REGION reading frame','AA JUNCTION','Functionality comment','Sequence','FR1-IMGT','FR2-IMGT','FR3-IMGT','CDR3-IMGT','JUNCTION','J-REGION','FR4-IMGT','P3V-nt nb','N-REGION-nt nb','N1-REGION-nt nb','P5D-nt nb','P3D-nt nb','N2-REGION-nt nb','P5J-nt nb','3V-REGION trimmed-nt nb','5D-REGION trimmed-nt nb','3D-REGION trimmed-nt nb','5J-REGION trimmed-nt nb','N-REGION','N1-REGION','N2-REGION', 'P5.D1.nt.nb', 'P3.D1.nt.nb', 'N2.REGION.nt.nb', 'P5.D2.nt.nb', 'P3.D2.nt.nb', 'N3.REGION.nt.nb', 'P5.D3.nt.nb', 'P3.D2.nt.nb', 'N4.REGION.nt.nb', 'X5.D1.REGION.trimmed.nt.nb', 'X3.D1.REGION.trimmed.nt.nb', 'X5.D2.REGION.trimmed.nt.nb', 'X3.D2.REGION.trimmed.nt.nb', 'X5.D3.REGION.trimmed.nt.nb', 'X3.D3.REGION.trimmed.nt.nb', 'D.REGION.nt.nb', 'D1.REGION.nt.nb', 'D2.REGION.nt.nb', 'D3.REGION.nt.nb')
79 92
80 out[,"VDJ Frame"] = as.character(out[,"VDJ Frame"]) 93 out[,"VDJ Frame"] = as.character(out[,"VDJ Frame"])
81
82 fltr = out[,"VDJ Frame"] == "in-frame" 94 fltr = out[,"VDJ Frame"] == "in-frame"
83 if(any(fltr)){ 95 if(any(fltr, na.rm=T)){
84 out[fltr, "VDJ Frame"] = "In-frame" 96 out[fltr, "VDJ Frame"] = "In-frame"
85 } 97 }
86
87 fltr = out[,"VDJ Frame"] == "null" 98 fltr = out[,"VDJ Frame"] == "null"
88 if(any(fltr)){ 99 if(any(fltr, na.rm = T)){
89 out[fltr, "VDJ Frame"] = "Out-of-frame" 100 out[fltr, "VDJ Frame"] = "Out-of-frame"
90 } 101 }
91
92 fltr = out[,"VDJ Frame"] == "out-of-frame" 102 fltr = out[,"VDJ Frame"] == "out-of-frame"
93 if(any(fltr)){ 103 if(any(fltr, na.rm = T)){
94 out[fltr, "VDJ Frame"] = "Out-of-frame" 104 out[fltr, "VDJ Frame"] = "Out-of-frame"
95 } 105 }
96
97 fltr = out[,"VDJ Frame"] == "" 106 fltr = out[,"VDJ Frame"] == ""
98 if(any(fltr)){ 107 if(any(fltr, na.rm = T)){
99 out[fltr, "VDJ Frame"] = "Out-of-frame" 108 out[fltr, "VDJ Frame"] = "Out-of-frame"
100 } 109 }
101 110
102 for(col in c('Top V Gene','Top D Gene','Top J Gene')){ 111 for(col in c('Top V Gene','Top D Gene','Top J Gene')){
103 out[,col] = as.character(out[,col]) 112 out[,col] = as.character(out[,col])
104 fltr = out[,col] == "" 113 fltr = out[,col] == ""
105 fltr[is.na(fltr)] = T 114 fltr[is.na(fltr)] = T
106 if(any(fltr)){ 115 if(any(fltr, na.rm = T)){
107 out[fltr,col] = "NA" 116 out[fltr,col] = "NA"
108 } 117 }
109 } 118 }
110 119
111 write.table(out, out.file, sep="\t", quote=F, row.names=F, col.names=T) 120 write.table(out, out.file, sep="\t", quote=F, row.names=F, col.names=T)