Mercurial > repos > davidvanzessen > argalaxy_tools

view imgt_loader/imgt_loader.r @ 58:4c6df851e262 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 1a7731bb4fe2e80fb1fdc88121783068af3edb08"
author rhpvorderman
date Fri, 21 Apr 2023 10:04:39 +0000
parents 81b3eb11ed2c
children
line wrap: on
line source

args <- commandArgs(trailingOnly = TRUE)
options(show.error.locations = TRUE)

summ.file = args[1]
sequences.file = args[2]
aa.file = args[3]
junction.file = args[4]
gapped.aa.file = args[5]
out.file = args[6]

print(summ.file)
print(out.file)

summ = read.table(summ.file, sep="\t", header=T, quote="", fill=T)
sequences = read.table(sequences.file, sep="\t", header=T, quote="", fill=T)
aa = read.table(aa.file, sep="\t", header=T, quote="", fill=T)
gapped.aa = read.table(gapped.aa.file, sep="\t", header=T, quote="", fill=T)
junction = read.table(junction.file, sep="\t", header=T, quote="", fill=T)

print(paste("nrow(summ)", nrow(summ)))

write.table(summ, out.file, sep="\t", quote=F, row.names=F, col.names=T)

fix_column_names = function(df){
    if("V.DOMAIN.Functionality" %in% names(df)){
        names(df)[names(df) == "V.DOMAIN.Functionality"] = "Functionality"
        print("found V.DOMAIN.Functionality, changed")
    }
    if("V.DOMAIN.Functionality.comment" %in% names(df)){
        names(df)[names(df) == "V.DOMAIN.Functionality.comment"] = "Functionality.comment"
        print("found V.DOMAIN.Functionality.comment, changed")
    }
    return(df)
}

print_missing_columns = function(df, cols, label=""){
	cols_in_df = cols %in% names(df)
	if(sum(!cols_in_df) > 0){
		print("Columns are missing from summary file, don't have:")
		print(cols[!cols_in_df])
	} else {
		print(paste("No missing columns for", label))
	}
}

summ = fix_column_names(summ)
sequences = fix_column_names(sequences)
aa = fix_column_names(aa)
gapped.aa = fix_column_names(gapped.aa)
junction = fix_column_names(junction)

print(paste("nrow(summ)", nrow(summ)))

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')
old_sequence_columns=c('CDR1.IMGT','CDR2.IMGT','CDR3.IMGT')
old_junction_columns=c('JUNCTION')

added_summary_columns=c('Functionality','V.REGION.identity..','V.REGION.identity.nt','D.REGION.reading.frame','AA.JUNCTION','Functionality.comment','Sequence')
added_sequence_columns=c('FR1.IMGT','FR2.IMGT','FR3.IMGT','CDR3.IMGT','JUNCTION','J.REGION','FR4.IMGT')

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')
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')

minimum_columns=c("Sequence.ID","JUNCTION.frame","V.GENE.and.allele","D.GENE.and.allele","J.GENE.and.allele")
print_missing_columns(summ, minimum_columns)
out=summ[,minimum_columns]

print(paste("nrow(summ)", nrow(summ)))
print(paste("nrow(aa)", nrow(aa)))

out[,"CDR1.Seq"] = aa[,"CDR1.IMGT"]
out[,"CDR1.Length"] = summ[,"CDR1.IMGT.length"]

out[,"CDR2.Seq"] = aa[,"CDR2.IMGT"]
out[,"CDR2.Length"] = summ[,"CDR2.IMGT.length"]

out[,"CDR3.Seq"] = gapped.aa[,"CDR3.IMGT"]
out[,"CDR3.Length"] = summ[,"CDR3.IMGT.length"]

out[,"CDR3.IMGT"] = out[,"CDR3.Seq"]

out[,"CDR3.Seq.DNA"] = sequences[,"CDR3.IMGT"]
out[,"CDR3.Length.DNA"] = nchar(as.character(out[,"CDR3.Seq.DNA"]))
out[,"Strand"] = summ[,"Orientation"]
out[,"CDR3.Found.How"] = "a"

out[,added_summary_columns] = summ[,added_summary_columns]

out[,added_sequence_columns] = aa[,added_sequence_columns]

out[,added_junction_columns] = junction[,added_junction_columns]

out[,"Top V Gene"] = gsub(".* ", "", gsub("\\*.*", "", summ[,"V.GENE.and.allele"]))
out[,"Top D Gene"] = gsub(".* ", "", gsub("\\*.*", "", summ[,"D.GENE.and.allele"]))
out[,"Top J Gene"] = gsub(".* ", "", gsub("\\*.*", "", summ[,"J.GENE.and.allele"]))

out = out[!grepl("Less than", summ[,"V.GENE.and.allele"]),]
out = out[!grepl("Less than", summ[,"D.GENE.and.allele"]),]
out = out[!grepl("Less than", summ[,"J.GENE.and.allele"]),]

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')]

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')

out[,"VDJ Frame"] = as.character(out[,"VDJ Frame"])
fltr = out[,"VDJ Frame"] == "in-frame"
if(any(fltr, na.rm=T)){
	out[fltr, "VDJ Frame"] = "In-frame"
}
fltr = out[,"VDJ Frame"] == "null"
if(any(fltr, na.rm = T)){
	out[fltr, "VDJ Frame"] = "Out-of-frame"
}
fltr = out[,"VDJ Frame"] == "out-of-frame"
if(any(fltr, na.rm = T)){
	out[fltr, "VDJ Frame"] = "Out-of-frame"
}
fltr = out[,"VDJ Frame"] == ""
if(any(fltr, na.rm = T)){
	out[fltr, "VDJ Frame"] = "Out-of-frame"
}

for(col in c('Top V Gene','Top D Gene','Top J Gene')){
	out[,col] = as.character(out[,col])
	fltr = out[,col] == ""
	fltr[is.na(fltr)] = T
	if(any(fltr, na.rm = T)){
		out[fltr,col] = "NA"
	}
}

write.table(out, out.file, sep="\t", quote=F, row.names=F, col.names=T)