annotate igblast/igblast.r @ 2:4f8a1d121384 draft

Uploaded
author davidvanzessen
date Mon, 05 Sep 2016 11:15:22 -0400
parents afe85eb6572e
children 124b7fd92a3e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
1 args <- commandArgs(trailingOnly = TRUE)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
2
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
3 infile=args[1]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
4 outfile=args[2]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
5
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
6 blasted = read.table(infile, header=T, sep="\t", fill=T, stringsAsFactors=F, comment.char="")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
7
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
8 blasted$ID = 1:nrow(blasted)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
9 blasted$VDJ.Frame = "Out-of-frame"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
10
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
11 search = blasted$inFrame == "true" & blasted$noStop == "false"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
12 if(sum(search) > 0){
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
13 blasted[search ,]$VDJ.Frame = "In-frame with stop codon"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
14 }
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
15
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
16 search = blasted$inFrame == "true" & blasted$noStop == "true"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
17 if(sum(search) > 0){
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
18 blasted[search ,]$VDJ.Frame = "In-frame"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
19 }
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
20
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
21 blasted$Top.V.Gene = blasted$vSegment
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
22 blasted$Top.D.Gene = blasted$dSegment
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
23 blasted$Top.J.Gene = blasted$jSegment
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
24 blasted$CDR1.Seq = blasted$cdr1aa
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
25 blasted$CDR1.Length = nchar(blasted$CDR1.Seq)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
26 blasted$CDR2.Seq = blasted$cdr2aa
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
27 blasted$CDR2.Length = nchar(blasted$CDR2.Seq)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
28 blasted$CDR3.Seq = blasted$cdr3aa
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
29 blasted$CDR3.Length = nchar(blasted$CDR3.Seq)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
30 blasted$CDR3.Seq.DNA = blasted$cdr3nt
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
31 blasted$CDR3.Length.DNA = nchar(blasted$CDR3.Seq.DNA)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
32 blasted$Strand = "+/-"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
33 blasted$CDR3.Found.How = "found"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
34
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
35 search = blasted$cdr3nt == ""
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
36 if(sum(search) > 0){
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
37 blasted[search,]$CDR3.Found.How = "NOT_FOUND"
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
38 }
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
39
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
40 blasted$AA.JUNCTION = blasted$CDR3.Seq
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
41
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
42 n = c("X.reads_count", "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", "AA.JUNCTION")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
43
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
44 n[!(n %in% names(blasted))]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
45
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
46 blasted = blasted[,c("X.reads_count", "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", "AA.JUNCTION")]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
47
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
48 names(blasted) = c("frequency.count", "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", "AA JUNCTION")
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
49
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
50 #duplicate rows based on frequency.count
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
51 blasted = blasted[rep(seq_len(nrow(blasted)), blasted$frequency.count),]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
52 blasted$ID = 1:nrow(blasted)
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
53
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
54 blasted = blasted[,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", "AA JUNCTION")]
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
55
afe85eb6572e Uploaded
davidvanzessen
parents:
diff changeset
56 write.table(blasted, outfile, quote=F, sep="\t", row.names=F, col.names=T)