# HG changeset patch
# User davidvanzessen
# Date 1496143593 14400
# Node ID 106275b54470a5250d4d9758e0dd0f214e6a3f58
# Parent bad6a9a53ae754673eb3123781d09326b8e602f5
Uploaded
diff -r bad6a9a53ae7 -r 106275b54470 experimental_design/experimental_design.r
--- a/experimental_design/experimental_design.r Mon Apr 24 09:39:59 2017 -0400
+++ b/experimental_design/experimental_design.r Tue May 30 07:26:33 2017 -0400
@@ -15,6 +15,11 @@
print(paste("Adding file", counter, "to", current.id))
dat = read.table(current, sep="\t", header=T, quote="", fill=T)
+ if(nrow(dat) == 0){
+ print(paste(counter, "of", current.id, "has no sequences, skipping"))
+ next
+ }
+
#IMGT check
dat$Sample = current.id
diff -r bad6a9a53ae7 -r 106275b54470 imgt_loader/imgt_loader.r
--- a/imgt_loader/imgt_loader.r Mon Apr 24 09:39:59 2017 -0400
+++ b/imgt_loader/imgt_loader.r Tue May 30 07:26:33 2017 -0400
@@ -1,4 +1,5 @@
args <- commandArgs(trailingOnly = TRUE)
+options(show.error.locations = TRUE)
summ.file = args[1]
sequences.file = args[2]
@@ -7,12 +8,19 @@
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"
@@ -31,6 +39,8 @@
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')
@@ -43,6 +53,9 @@
out=summ[,c("Sequence.ID","JUNCTION.frame","V.GENE.and.allele","D.GENE.and.allele","J.GENE.and.allele")]
+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"]
@@ -78,24 +91,20 @@
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)){
+if(any(fltr, na.rm=T)){
out[fltr, "VDJ Frame"] = "In-frame"
}
-
fltr = out[,"VDJ Frame"] == "null"
-if(any(fltr)){
+if(any(fltr, na.rm = T)){
out[fltr, "VDJ Frame"] = "Out-of-frame"
}
-
fltr = out[,"VDJ Frame"] == "out-of-frame"
-if(any(fltr)){
+if(any(fltr, na.rm = T)){
out[fltr, "VDJ Frame"] = "Out-of-frame"
}
-
fltr = out[,"VDJ Frame"] == ""
-if(any(fltr)){
+if(any(fltr, na.rm = T)){
out[fltr, "VDJ Frame"] = "Out-of-frame"
}
@@ -103,7 +112,7 @@
out[,col] = as.character(out[,col])
fltr = out[,col] == ""
fltr[is.na(fltr)] = T
- if(any(fltr)){
+ if(any(fltr, na.rm = T)){
out[fltr,col] = "NA"
}
}
diff -r bad6a9a53ae7 -r 106275b54470 imgt_loader/imgt_loader.sh
--- a/imgt_loader/imgt_loader.sh Mon Apr 24 09:39:59 2017 -0400
+++ b/imgt_loader/imgt_loader.sh Tue May 30 07:26:33 2017 -0400
@@ -66,6 +66,9 @@
find $PWD/$name/files -iname "5_*" -exec cat {} + > $PWD/$name/aa.txt
find $PWD/$name/files -iname "6_*" -exec cat {} + > $PWD/$name/junction.txt
+echo "summ.txt `cat $PWD/$name/summ.txt | wc -l`"
+echo "aa.txt `cat $PWD/$name/aa.txt | wc -l`"
+
#python $dir/imgt_loader.py --summ $PWD/$name/summ.txt --aa $PWD/$name/aa.txt --junction $PWD/$name/junction.txt --output $output
Rscript --verbose $dir/imgt_loader.r $PWD/$name/summ.txt $PWD/$name/sequences.txt $PWD/$name/aa.txt $PWD/$name/junction.txt $PWD/$name/gapped_aa.txt $output 2>&1
diff -r bad6a9a53ae7 -r 106275b54470 report_clonality/RScript.r
--- a/report_clonality/RScript.r Mon Apr 24 09:39:59 2017 -0400
+++ b/report_clonality/RScript.r Tue May 30 07:26:33 2017 -0400
@@ -122,12 +122,17 @@
clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$clonality_clonaltype), ]
}
+if(nrow(PRODF) == 0){
+ stop("No sequences left after filtering")
+}
+
prod.unique.sample.count = data.frame(data.table(PRODF)[, list(Productive_unique=.N), by=c("Sample")])
prod.unique.rep.count = data.frame(data.table(PRODF)[, list(Productive_unique=.N), by=c("Sample", "Replicate")])
unprod.unique.sample.count = data.frame(data.table(UNPROD)[, list(Unproductive_unique=.N), by=c("Sample")])
unprod.unique.rep.count = data.frame(data.table(UNPROD)[, list(Unproductive_unique=.N), by=c("Sample", "Replicate")])
+
PRODF$freq = 1
if(any(grepl(pattern="_", x=PRODF$ID))){ #the frequency can be stored in the ID with the pattern ".*_freq_.*"
@@ -176,6 +181,7 @@
#create the table on the overview page with the productive/unique counts per sample/replicate
#first for sample
+
sample.count = merge(input.sample.count, prod.sample.count, by="Sample", all.x=T)
sample.count$perc_prod = round(sample.count$Productive / sample.count$All * 100)
sample.count = merge(sample.count, prod.unique.sample.count, by="Sample", all.x=T)
@@ -188,6 +194,16 @@
#then sample/replicate
rep.count = merge(input.rep.count, prod.rep.count, by=c("Sample", "Replicate"), all.x=T)
+
+print(rep.count)
+
+fltr = is.na(rep.count$Productive)
+if(any(fltr)){
+ rep.count[fltr,"Productive"] = 0
+}
+
+print(rep.count)
+
rep.count$perc_prod = round(rep.count$Productive / rep.count$All * 100)
rep.count = merge(rep.count, prod.unique.rep.count, by=c("Sample", "Replicate"), all.x=T)
rep.count$perc_prod_un = round(rep.count$Productive_unique / rep.count$All * 100)
diff -r bad6a9a53ae7 -r 106275b54470 report_clonality/r_wrapper.sh
--- a/report_clonality/r_wrapper.sh Mon Apr 24 09:39:59 2017 -0400
+++ b/report_clonality/r_wrapper.sh Tue May 30 07:26:33 2017 -0400
@@ -42,10 +42,14 @@
do
echo "| $sample | " >> $2
echo "$all | " >> $2
- echo "$productive (${perc_prod}%) | " >> $2
- echo "$productive_unique (${perc_prod_un}%) | " >> $2
- echo "$unproductive (${perc_unprod}%) | " >> $2
- echo "$unproductive_unique (${perc_unprod_un}%) |
" >> $2
+ if [[ "$productive" != "0" ]] ; then
+ echo "$productive (${perc_prod}%) | " >> $2
+ echo "$productive_unique (${perc_prod_un}%) | " >> $2
+ echo "$unproductive (${perc_unprod}%) | " >> $2
+ echo "$unproductive_unique (${perc_unprod_un}%) | " >> $2
+ else
+ echo "No productive sequences! | " >> $2
+ fi
done < $outputDir/productive_counting.txt
echo "
" >> $2
echo "Table showing the number and percentage of (unique) productive and unproductive sequences per donor and per replicate.
" >> $2