Mercurial > repos > davidvanzessen > argalaxy_tools
changeset 26:28fbbdfd7a87 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Mon, 13 Feb 2017 09:08:46 -0500 |
| parents | 94765af0db1f |
| children | b539aeb75980 |
| files | experimental_design/experimental_design.py imgt_loader/imgt_loader.py report_clonality/RScript.r report_clonality/RScript.r~ report_clonality/circos/fonts.conf~ report_clonality/circos/housekeeping.conf~ report_clonality/r_wrapper.sh report_clonality/r_wrapper.sh~ |
| diffstat | 7 files changed, 998 insertions(+), 349 deletions(-) [+] |
line wrap: on
line diff
--- a/experimental_design/experimental_design.py Thu Feb 09 07:20:09 2017 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,50 +0,0 @@ -import sys -import pandas as pd - -def main(): - patients = {} - files = [] - sample_id = sys.argv[1] - imgt_files = 0 - blast_files = 0 - #organize files - for arg in sys.argv[2:-2]: - if arg.find("/") is -1: - patients[sample_id] = files - files = [] - sample_id = arg - else: - df = pd.read_csv(arg, sep="\t", dtype=object, error_bad_lines=False) - if "Functionality" in list(df.columns.values): - df["VDJ Frame"][df["Functionality"] != "productive"] = "In-frame with stop codon" - imgt_files += 1 - else: - blast_files += 1 - files.append(df) - patients[sample_id] = files - columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', - u'CDR3 Seq', u'CDR3 Length', u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'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', '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', u'Sample', u'Replicate'] - if "N-REGION-nt nb" in files[0].columns: - columns.insert(30, "N-REGION-nt nb") - if blast_files is not 0: - print "Has a parsed blastn file, using limited columns." - columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', u'CDR3 Seq', u'CDR3 Length', u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'Sample', u'Replicate'] - - result = None - for patient_id, samples in patients.iteritems(): - count = 1 - for sample in samples: - sample['Sample'] = patient_id - sample['Replicate'] = str(count) - count += 1 - if result is None: - result = sample[columns] - else: - result = result.append(sample[columns]) - result.to_csv(sys.argv[-1], sep="\t", index=False, index_label="index") - -if __name__ == "__main__": - main()
--- a/imgt_loader/imgt_loader.py Thu Feb 09 07:20:09 2017 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,147 +0,0 @@ -import pandas as pd -try: - pd.options.mode.chained_assignment = None # default='warn' -except: - pass -import re -import argparse -import os - -def stop_err( msg, ret=1 ): - sys.stderr.write( msg ) - sys.exit( ret ) - -#docs.python.org/dev/library/argparse.html -parser = argparse.ArgumentParser() -parser.add_argument("--summ", help="The 1_Summary file from the imgt output") -parser.add_argument("--aa", help="The 5_AA-Sequence file from the imgt output") -parser.add_argument("--junction", help="The 6_Junction file from the imgt output") -parser.add_argument("--output", help="Output file") - -args = parser.parse_args() - -old_summary_columns = [u'Sequence ID', u'JUNCTION frame', u'V-GENE and allele', u'D-GENE and allele', u'J-GENE and allele', u'CDR1-IMGT length', u'CDR2-IMGT length', u'CDR3-IMGT length', u'Orientation'] -old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT'] -old_junction_columns = [u'JUNCTION'] - -added_summary_columns = [u'Functionality', u'V-REGION identity %', u'V-REGION identity nt', u'D-REGION reading frame', u'AA JUNCTION', u'Functionality comment', u'Sequence'] -added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT'] -added_junction_columns = [u"P3'V-nt nb", u'N-REGION-nt nb', u'N1-REGION-nt nb', u"P5'D-nt nb", u"P3'D-nt nb", u'N2-REGION-nt nb', u"P5'J-nt nb", u"3'V-REGION trimmed-nt nb", - u"5'D-REGION trimmed-nt nb", u"3'D-REGION trimmed-nt nb", u"5'J-REGION trimmed-nt nb", u"N-REGION", u"N1-REGION", u"N2-REGION"] - -outFile = args.output - -#fSummary = pd.read_csv(triplets[0][0], sep="\t", low_memory=False) -fSummary = pd.read_csv(args.summ, sep="\t", dtype=object) -#fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False) -fSequence = pd.read_csv(args.aa, sep="\t", dtype=object) -#fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False) -fJunction = pd.read_csv(args.junction, sep="\t", dtype=object) -tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]] - -tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"] -tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"] - -tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"] -tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"] - -tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"] -tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"] - -tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"] -tmp["CDR3 Length DNA"] = '1' -tmp["Strand"] = fSummary["Orientation"] -tmp["CDR3 Found How"] = 'a' - -for col in added_summary_columns: - tmp[col] = fSummary[col] - -for col in added_sequence_columns: - tmp[col] = fSequence[col] - -for col in added_junction_columns: - tmp[col] = fJunction[col] - -outFrame = tmp - -outFrame.columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', u'CDR3 Seq', u'CDR3 Length', - u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'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"] - -""" -IGHV[0-9]-[0-9ab]+-?[0-9]?D? -TRBV[0-9]{1,2}-?[0-9]?-?[123]? -IGKV[0-3]D?-[0-9]{1,2} -IGLV[0-9]-[0-9]{1,2} -TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])? -TRGV[234589] -TRDV[1-3] - -IGHD[0-9]-[0-9ab]+ -TRBD[12] -TRDD[1-3] - -IGHJ[1-6] -TRBJ[12]-[1-7] -IGKJ[1-5] -IGLJ[12367] -TRAJ[0-9]{1,2} -TRGJP?[12] -TRDJ[1-4] -""" - -vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)", - r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)", - r"(IGKV[0-3]D?-[0-9]{1,2})", - r"(IGLV[0-9]-[0-9]{1,2})", - r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)", - r"(TRGV[234589])", - r"(TRDV[1-3])", - r"(IGHV[0-9]S[0-9]+)"] - -dPattern = [r"(IGHD[0-9]-[0-9ab]+)", - r"(TRBD[12])", - r"(TRDD[1-3])"] - -jPattern = [r"(IGHJ[1-6])", - r"(TRBJ[12]-[1-7])", - r"(IGKJ[1-5])", - r"(IGLJ[12367])", - r"(TRAJ[0-9]{1,2})", - r"(TRGJP?[12])", - r"(TRDJ[1-4])"] - -vPattern = re.compile(r"|".join(vPattern)) - -dPattern = re.compile(r"|".join(dPattern)) - -jPattern = re.compile(r"|".join(jPattern)) - - -def filterGenes(s, pattern): - if type(s) is not str: - return "NA" - res = pattern.search(s) - if res: - return res.group(0) - return "NA" - - - -outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern)) -outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern)) -outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern)) - - -tmp = outFrame["VDJ Frame"] -tmp = tmp.replace("in-frame", "In-frame") -tmp = tmp.replace("null", "Out-of-frame") -tmp = tmp.replace("out-of-frame", "Out-of-frame") -outFrame["VDJ Frame"] = tmp -outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len) -safeLength = lambda x: len(x) if type(x) == str else 0 -#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? -#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? -outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")
--- a/report_clonality/RScript.r Thu Feb 09 07:20:09 2017 -0500 +++ b/report_clonality/RScript.r Mon Feb 13 09:08:46 2017 -0500 @@ -614,86 +614,86 @@ colnames(coincidence.table) = c("Coincidence Type", "Raw Coincidence Freq") write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) } - } else if(clonality_method == "old") { - clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")]) - - #write files for every coincidence group of >1 - samples = unique(clonalFreq$Sample) - for(sample in samples){ - clonalFreqSample = clonalFreq[clonalFreq$Sample == sample,] - if(max(clonalFreqSample$Type) > 1){ - for(i in 2:max(clonalFreqSample$Type)){ - clonalFreqSampleType = clonalFreqSample[clonalFreqSample$Type == i,] - clonalityFrame.sub = clonalityFrame[clonalityFrame$clonaltype %in% clonalFreqSampleType$clonaltype,] - clonalityFrame.sub = clonalityFrame.sub[order(clonalityFrame.sub$clonaltype),] - write.table(clonalityFrame.sub, file=paste("coincidences_", sample, "_", i, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) - } - } - } - - clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")]) - clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count - clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")]) - clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample") - - ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15') - tcct = textConnection(ct) - CT = read.table(tcct, sep="\t", header=TRUE) - close(tcct) - clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T) - clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight - - ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")]) - ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) - clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads) - ReplicateReads$Reads = as.numeric(ReplicateReads$Reads) - ReplicateReads$squared = as.numeric(ReplicateReads$Reads * ReplicateReads$Reads) - - ReplicatePrint <- function(dat){ - write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) - } - - ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) - lapply(ReplicateSplit, FUN=ReplicatePrint) - - ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")]) - clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T) - - ReplicateSumPrint <- function(dat){ - write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) - } - - ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) - lapply(ReplicateSumSplit, FUN=ReplicateSumPrint) - - clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")]) - clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T) - clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow - clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2) - clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1) - - ClonalityScorePrint <- function(dat){ - write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) - } - - clonalityScore = clonalFreqCount[c("Sample", "Result")] - clonalityScore = unique(clonalityScore) - - clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"]) - lapply(clonalityScoreSplit, FUN=ClonalityScorePrint) - - clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")] - - - - ClonalityOverviewPrint <- function(dat){ - dat = dat[order(dat[,2]),] - write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) - } - - clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample) - lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint) + } + clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")]) + + #write files for every coincidence group of >1 + samples = unique(clonalFreq$Sample) + for(sample in samples){ + clonalFreqSample = clonalFreq[clonalFreq$Sample == sample,] + if(max(clonalFreqSample$Type) > 1){ + for(i in 2:max(clonalFreqSample$Type)){ + clonalFreqSampleType = clonalFreqSample[clonalFreqSample$Type == i,] + clonalityFrame.sub = clonalityFrame[clonalityFrame$clonaltype %in% clonalFreqSampleType$clonaltype,] + clonalityFrame.sub = clonalityFrame.sub[order(clonalityFrame.sub$clonaltype),] + write.table(clonalityFrame.sub, file=paste("coincidences_", sample, "_", i, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T) + } + } + } + + clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")]) + clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count + clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample") + + ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15') + tcct = textConnection(ct) + CT = read.table(tcct, sep="\t", header=TRUE) + close(tcct) + clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T) + clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight + + ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")]) + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) + clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads) + ReplicateReads$Reads = as.numeric(ReplicateReads$Reads) + ReplicateReads$squared = as.numeric(ReplicateReads$Reads * ReplicateReads$Reads) + + ReplicatePrint <- function(dat){ + write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F) } + + ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) + lapply(ReplicateSplit, FUN=ReplicatePrint) + + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T) + + ReplicateSumPrint <- function(dat){ + write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F) + } + + ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) + lapply(ReplicateSumSplit, FUN=ReplicateSumPrint) + + clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T) + clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow + clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2) + clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1) + + ClonalityScorePrint <- function(dat){ + write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F) + } + + clonalityScore = clonalFreqCount[c("Sample", "Result")] + clonalityScore = unique(clonalityScore) + + clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"]) + lapply(clonalityScoreSplit, FUN=ClonalityScorePrint) + + clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")] + + + + ClonalityOverviewPrint <- function(dat){ + dat = dat[order(dat[,2]),] + write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F) + } + + clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample) + lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint) + } bak = PRODF @@ -724,16 +724,14 @@ PRODF.with.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) > 2,] PRODF.no.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) < 4,] + write.table(PRODF.no.D, "productive_no_D.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T) UNPROD.with.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) > 2,] UNPROD.no.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) < 4,] + write.table(UNPROD.no.D, "unproductive_no_D.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T) num_median = function(x, na.rm=T) { as.numeric(median(x, na.rm=na.rm)) } - print("---- table prod.with.d cdr3.length ----") - print(table(PRODF.with.D$CDR3.Length, useNA="ifany")) - print(median(PRODF.with.D$CDR3.Length, na.rm=T)) - newData = data.frame(data.table(PRODF.with.D)[,list(unique=.N, VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=mean(.SD$P3V.nt.nb, na.rm=T), @@ -772,10 +770,6 @@ newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1) write.table(newData, "junctionAnalysisProd_median_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F) - print("---- table unprod.with.d cdr3.length ----") - print(table(UNPROD.with.D$CDR3.Length, useNA="ifany")) - print(median(UNPROD.with.D$CDR3.Length, na.rm=T)) - newData = data.frame(data.table(UNPROD.with.D)[,list(unique=.N, VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=mean(.SD$P3V.nt.nb, na.rm=T), @@ -819,11 +813,11 @@ newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=mean(.SD$P3V.nt.nb, na.rm=T), - N1=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + N1=mean(.SD$N.REGION.nt.nb, na.rm=T), P2=mean(.SD$P5J.nt.nb, na.rm=T), DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)), - Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + Total.N=mean(.SD$N.REGION.nt.nb, na.rm=T), Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)), Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))), by=c("Sample")]) @@ -833,39 +827,43 @@ newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=num_median(.SD$P3V.nt.nb, na.rm=T), - N1=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + N1=median(.SD$N.REGION.nt.nb, na.rm=T), P2=num_median(.SD$P5J.nt.nb, na.rm=T), DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)), - Total.N=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + Total.N=median(.SD$N.REGION.nt.nb, na.rm=T), Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)), Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))), by=c("Sample")]) newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1) write.table(newData, "junctionAnalysisProd_median_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F) + print(paste("mean N:", mean(UNPROD.no.D$N.REGION.nt.nb, na.rm=T))) + newData = data.frame(data.table(UNPROD.no.D)[,list(unique=.N, VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=mean(.SD$P3V.nt.nb, na.rm=T), - N1=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + N1=mean(.SD$N.REGION.nt.nb, na.rm=T), P2=mean(.SD$P5J.nt.nb, na.rm=T), DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)), - Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + Total.N=mean(.SD$N.REGION.nt.nb, na.rm=T), Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)), Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))), by=c("Sample")]) newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1) write.table(newData, "junctionAnalysisUnProd_mean_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F) + print(paste("median N:", num_median(UNPROD.no.D$N.REGION.nt.nb, na.rm=T))) + newData = data.frame(data.table(UNPROD.no.D)[,list(unique=.N, VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), P1=num_median(.SD$P3V.nt.nb, na.rm=T), - N1=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + N1=median(.SD$N.REGION.nt.nb, na.rm=T), P2=num_median(.SD$P5J.nt.nb, na.rm=T), DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)), - Total.N=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)), + Total.N=median(.SD$N.REGION.nt.nb, na.rm=T), Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)), Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))), by=c("Sample")]) @@ -960,20 +958,13 @@ clonaltype.in.replicates$clonaltype = do.call(paste, c(clonaltype.in.replicates[clonaltype], sep = ":")) clonaltype.in.replicates = clonaltype.in.replicates[,c("clonaltype","Replicate", "ID", "Sequence", "Sample")] -print(head(clonaltype.in.replicates)) - clonaltype.counts = data.frame(table(clonaltype.in.replicates$clonaltype)) names(clonaltype.counts) = c("clonaltype", "coincidence") -print(head(clonaltype.counts)) - clonaltype.counts = clonaltype.counts[clonaltype.counts$coincidence > 1,] -head(clonaltype.counts) - clonaltype.in.replicates = clonaltype.in.replicates[clonaltype.in.replicates$clonaltype %in% clonaltype.counts$clonaltype,] clonaltype.in.replicates = merge(clonaltype.in.replicates, clonaltype.counts, by="clonaltype") -print(head(clonaltype.in.replicates)) clonaltype.in.replicates = clonaltype.in.replicates[order(clonaltype.in.replicates$clonaltype),c("coincidence","clonaltype", "Sample", "Replicate", "ID", "Sequence")] write.table(clonaltype.in.replicates, "clonaltypes_replicates.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/report_clonality/RScript.r~ Mon Feb 13 09:08:46 2017 -0500 @@ -0,0 +1,658 @@ +# ---------------------- load/install packages ---------------------- + +if (!("gridExtra" %in% rownames(installed.packages()))) { + install.packages("gridExtra", repos="http://cran.xl-mirror.nl/") +} +library(gridExtra) +if (!("ggplot2" %in% rownames(installed.packages()))) { + install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") +} +library(ggplot2) +if (!("plyr" %in% rownames(installed.packages()))) { + install.packages("plyr", repos="http://cran.xl-mirror.nl/") +} +library(plyr) + +if (!("data.table" %in% rownames(installed.packages()))) { + install.packages("data.table", repos="http://cran.xl-mirror.nl/") +} +library(data.table) + +if (!("reshape2" %in% rownames(installed.packages()))) { + install.packages("reshape2", repos="http://cran.xl-mirror.nl/") +} +library(reshape2) + +if (!("lymphclon" %in% rownames(installed.packages()))) { + install.packages("lymphclon", repos="http://cran.xl-mirror.nl/") +} +library(lymphclon) + +# ---------------------- parameters ---------------------- + +args <- commandArgs(trailingOnly = TRUE) + +infile = args[1] #path to input file +outfile = args[2] #path to output file +outdir = args[3] #path to output folder (html/images/data) +clonaltype = args[4] #clonaltype definition, or 'none' for no unique filtering +ct = unlist(strsplit(clonaltype, ",")) +species = args[5] #human or mouse +locus = args[6] # IGH, IGK, IGL, TRB, TRA, TRG or TRD +filterproductive = ifelse(args[7] == "yes", T, F) #should unproductive sequences be filtered out? (yes/no) +clonality_method = args[8] + +# ---------------------- Data preperation ---------------------- + +inputdata = read.table(infile, sep="\t", header=TRUE, fill=T, comment.char="") + +setwd(outdir) + +# remove weird rows +inputdata = inputdata[inputdata$Sample != "",] + +#remove the allele from the V,D and J genes +inputdata$Top.V.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.V.Gene) +inputdata$Top.D.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.D.Gene) +inputdata$Top.J.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.J.Gene) + +inputdata$clonaltype = 1:nrow(inputdata) + +PRODF = inputdata +UNPROD = inputdata +if(filterproductive){ + if("Functionality" %in% colnames(inputdata)) { # "Functionality" is an IMGT column + PRODF = inputdata[inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)", ] + UNPROD = inputdata[!(inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)"), ] + } else { + PRODF = inputdata[inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" , ] + UNPROD = inputdata[!(inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" ), ] + } +} + +clonalityFrame = PRODF + +#remove duplicates based on the clonaltype +if(clonaltype != "none"){ + clonaltype = paste(clonaltype, ",Sample", sep="") #add sample column to clonaltype, unique within samples + PRODF$clonaltype = do.call(paste, c(PRODF[unlist(strsplit(clonaltype, ","))], sep = ":")) + PRODF = PRODF[!duplicated(PRODF$clonaltype), ] + + UNPROD$clonaltype = do.call(paste, c(UNPROD[unlist(strsplit(clonaltype, ","))], sep = ":")) + UNPROD = UNPROD[!duplicated(UNPROD$clonaltype), ] + + #again for clonalityFrame but with sample+replicate + clonalityFrame$clonaltype = do.call(paste, c(clonalityFrame[unlist(strsplit(clonaltype, ","))], sep = ":")) + clonalityFrame$clonality_clonaltype = do.call(paste, c(clonalityFrame[unlist(strsplit(paste(clonaltype, ",Replicate", sep=""), ","))], sep = ":")) + clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$clonality_clonaltype), ] +} + +PRODF$freq = 1 + +if(any(grepl(pattern="_", x=PRODF$ID))){ #the frequency can be stored in the ID with the pattern ".*_freq_.*" + PRODF$freq = gsub("^[0-9]+_", "", PRODF$ID) + PRODF$freq = gsub("_.*", "", PRODF$freq) + PRODF$freq = as.numeric(PRODF$freq) + if(any(is.na(PRODF$freq))){ #if there was an "_" in the ID, but not the frequency, go back to frequency of 1 for every sequence + PRODF$freq = 1 + } +} + + + +#write the complete dataset that is left over, will be the input if 'none' for clonaltype and 'no' for filterproductive +write.table(PRODF, "allUnique.txt", sep=",",quote=F,row.names=F,col.names=T) +write.table(PRODF, "allUnique.csv", sep="\t",quote=F,row.names=F,col.names=T) +write.table(UNPROD, "allUnproductive.csv", sep=",",quote=F,row.names=F,col.names=T) + +#write the samples to a file +sampleFile <- file("samples.txt") +un = unique(inputdata$Sample) +un = paste(un, sep="\n") +writeLines(un, sampleFile) +close(sampleFile) + +# ---------------------- Counting the productive/unproductive and unique sequences ---------------------- + +if(!("Functionality" %in% inputdata)){ #add a functionality column to the igblast data + inputdata$Functionality = "unproductive" + search = (inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND") + if(sum(search) > 0){ + inputdata[search,]$Functionality = "productive" + } +} + +inputdata.dt = data.table(inputdata) #for speed + +if(clonaltype == "none"){ + ct = c("clonaltype") +} + +inputdata.dt$samples_replicates = paste(inputdata.dt$Sample, inputdata.dt$Replicate, sep="_") +samples_replicates = c(unique(inputdata.dt$samples_replicates), unique(as.character(inputdata.dt$Sample))) +frequency_table = data.frame(ID = samples_replicates[order(samples_replicates)]) + + +sample_productive_count = inputdata.dt[, list(All=.N, + Productive = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",]), + perc_prod = 1, + Productive_unique = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",list(count=.N),by=ct]), + perc_prod_un = 1, + Unproductive= nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",]), + perc_unprod = 1, + Unproductive_unique =nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",list(count=.N),by=ct]), + perc_unprod_un = 1), + by=c("Sample")] + +sample_productive_count$perc_prod = round(sample_productive_count$Productive / sample_productive_count$All * 100) +sample_productive_count$perc_prod_un = round(sample_productive_count$Productive_unique / sample_productive_count$All * 100) + +sample_productive_count$perc_unprod = round(sample_productive_count$Unproductive / sample_productive_count$All * 100) +sample_productive_count$perc_unprod_un = round(sample_productive_count$Unproductive_unique / sample_productive_count$All * 100) + + +sample_replicate_productive_count = inputdata.dt[, list(All=.N, + Productive = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",]), + perc_prod = 1, + Productive_unique = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",list(count=.N),by=ct]), + perc_prod_un = 1, + Unproductive= nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",]), + perc_unprod = 1, + Unproductive_unique =nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",list(count=.N),by=ct]), + perc_unprod_un = 1), + by=c("samples_replicates")] + +sample_replicate_productive_count$perc_prod = round(sample_replicate_productive_count$Productive / sample_replicate_productive_count$All * 100) +sample_replicate_productive_count$perc_prod_un = round(sample_replicate_productive_count$Productive_unique / sample_replicate_productive_count$All * 100) + +sample_replicate_productive_count$perc_unprod = round(sample_replicate_productive_count$Unproductive / sample_replicate_productive_count$All * 100) +sample_replicate_productive_count$perc_unprod_un = round(sample_replicate_productive_count$Unproductive_unique / sample_replicate_productive_count$All * 100) + +setnames(sample_replicate_productive_count, colnames(sample_productive_count)) + +counts = rbind(sample_replicate_productive_count, sample_productive_count) +counts = counts[order(counts$Sample),] + +write.table(x=counts, file="productive_counting.txt", sep=",",quote=F,row.names=F,col.names=F) + +# ---------------------- Frequency calculation for V, D and J ---------------------- + +PRODFV = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.V.Gene")]) +Total = ddply(PRODFV, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFV = merge(PRODFV, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFV = ddply(PRODFV, c("Sample", "Top.V.Gene"), summarise, relFreq= (Length*100 / Total)) + +PRODFD = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.D.Gene")]) +Total = ddply(PRODFD, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFD = merge(PRODFD, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFD = ddply(PRODFD, c("Sample", "Top.D.Gene"), summarise, relFreq= (Length*100 / Total)) + +PRODFJ = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.J.Gene")]) +Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length))) +PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE) +PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total)) + +# ---------------------- Setting up the gene names for the different species/loci ---------------------- + +Vchain = "" +Dchain = "" +Jchain = "" + +if(species == "custom"){ + print("Custom genes: ") + splt = unlist(strsplit(locus, ";")) + print(paste("V:", splt[1])) + print(paste("D:", splt[2])) + print(paste("J:", splt[3])) + + Vchain = unlist(strsplit(splt[1], ",")) + Vchain = data.frame(v.name = Vchain, chr.orderV = 1:length(Vchain)) + + Dchain = unlist(strsplit(splt[2], ",")) + if(length(Dchain) > 0){ + Dchain = data.frame(v.name = Dchain, chr.orderD = 1:length(Dchain)) + } else { + Dchain = data.frame(v.name = character(0), chr.orderD = numeric(0)) + } + + Jchain = unlist(strsplit(splt[3], ",")) + Jchain = data.frame(v.name = Jchain, chr.orderJ = 1:length(Jchain)) + +} else { + genes = read.table("genes.txt", sep="\t", header=TRUE, fill=T, comment.char="") + + Vchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "V",c("IMGT.GENE.DB", "chr.order")] + colnames(Vchain) = c("v.name", "chr.orderV") + Dchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "D",c("IMGT.GENE.DB", "chr.order")] + colnames(Dchain) = c("v.name", "chr.orderD") + Jchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "J",c("IMGT.GENE.DB", "chr.order")] + colnames(Jchain) = c("v.name", "chr.orderJ") +} +useD = TRUE +if(nrow(Dchain) == 0){ + useD = FALSE + cat("No D Genes in this species/locus") +} +print(paste("useD:", useD)) + +# ---------------------- merge with the frequency count ---------------------- + +PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) + +PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE) + +PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) + +# ---------------------- Create the V, D and J frequency plots and write the data.frame for every plot to a file ---------------------- + +pV = ggplot(PRODFV) +pV = pV + geom_bar( aes( x=factor(reorder(Top.V.Gene, chr.orderV)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage") +write.table(x=PRODFV, file="VFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + +png("VPlot.png",width = 1280, height = 720) +pV +dev.off(); + +if(useD){ + pD = ggplot(PRODFD) + pD = pD + geom_bar( aes( x=factor(reorder(Top.D.Gene, chr.orderD)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + pD = pD + xlab("Summary of D gene") + ylab("Frequency") + ggtitle("Relative frequency of D gene usage") + write.table(x=PRODFD, file="DFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + + png("DPlot.png",width = 800, height = 600) + print(pD) + dev.off(); +} + +pJ = ggplot(PRODFJ) +pJ = pJ + geom_bar( aes( x=factor(reorder(Top.J.Gene, chr.orderJ)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage") +write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + +png("JPlot.png",width = 800, height = 600) +pJ +dev.off(); + +pJ = ggplot(PRODFJ) +pJ = pJ + geom_bar( aes( x=factor(reorder(Top.J.Gene, chr.orderJ)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage") +write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + +png("JPlot.png",width = 800, height = 600) +pJ +dev.off(); + +# ---------------------- Now the frequency plots of the V, D and J families ---------------------- + +VGenes = PRODF[,c("Sample", "Top.V.Gene")] +VGenes$Top.V.Gene = gsub("-.*", "", VGenes$Top.V.Gene) +VGenes = data.frame(data.table(VGenes)[, list(Count=.N), by=c("Sample", "Top.V.Gene")]) +TotalPerSample = data.frame(data.table(VGenes)[, list(total=sum(.SD$Count)), by=Sample]) +VGenes = merge(VGenes, TotalPerSample, by="Sample") +VGenes$Frequency = VGenes$Count * 100 / VGenes$total +VPlot = ggplot(VGenes) +VPlot = VPlot + geom_bar(aes( x = Top.V.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + ggtitle("Distribution of V gene families") + + ylab("Percentage of sequences") +png("VFPlot.png") +VPlot +dev.off(); +write.table(x=VGenes, file="VFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + +if(useD){ + DGenes = PRODF[,c("Sample", "Top.D.Gene")] + DGenes$Top.D.Gene = gsub("-.*", "", DGenes$Top.D.Gene) + DGenes = data.frame(data.table(DGenes)[, list(Count=.N), by=c("Sample", "Top.D.Gene")]) + TotalPerSample = data.frame(data.table(DGenes)[, list(total=sum(.SD$Count)), by=Sample]) + DGenes = merge(DGenes, TotalPerSample, by="Sample") + DGenes$Frequency = DGenes$Count * 100 / DGenes$total + DPlot = ggplot(DGenes) + DPlot = DPlot + geom_bar(aes( x = Top.D.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + ggtitle("Distribution of D gene families") + + ylab("Percentage of sequences") + png("DFPlot.png") + print(DPlot) + dev.off(); + write.table(x=DGenes, file="DFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) +} + +JGenes = PRODF[,c("Sample", "Top.J.Gene")] +JGenes$Top.J.Gene = gsub("-.*", "", JGenes$Top.J.Gene) +JGenes = data.frame(data.table(JGenes)[, list(Count=.N), by=c("Sample", "Top.J.Gene")]) +TotalPerSample = data.frame(data.table(JGenes)[, list(total=sum(.SD$Count)), by=Sample]) +JGenes = merge(JGenes, TotalPerSample, by="Sample") +JGenes$Frequency = JGenes$Count * 100 / JGenes$total +JPlot = ggplot(JGenes) +JPlot = JPlot + geom_bar(aes( x = Top.J.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + ggtitle("Distribution of J gene families") + + ylab("Percentage of sequences") +png("JFPlot.png") +JPlot +dev.off(); +write.table(x=JGenes, file="JFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) + +# ---------------------- Plotting the cdr3 length ---------------------- + +CDR3Length = data.frame(data.table(PRODF)[, list(Count=.N), by=c("Sample", "CDR3.Length.DNA")]) +TotalPerSample = data.frame(data.table(CDR3Length)[, list(total=sum(.SD$Count)), by=Sample]) +CDR3Length = merge(CDR3Length, TotalPerSample, by="Sample") +CDR3Length$Frequency = CDR3Length$Count * 100 / CDR3Length$total +CDR3LengthPlot = ggplot(CDR3Length) +CDR3LengthPlot = CDR3LengthPlot + geom_bar(aes( x = CDR3.Length.DNA, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + ggtitle("Length distribution of CDR3") + + xlab("CDR3 Length") + + ylab("Percentage of sequences") +png("CDR3LengthPlot.png",width = 1280, height = 720) +CDR3LengthPlot +dev.off() +write.table(x=CDR3Length, file="CDR3LengthPlot.csv", sep=",",quote=F,row.names=F,col.names=T) + +# ---------------------- Plot the heatmaps ---------------------- + + +#get the reverse order for the V and D genes +revVchain = Vchain +revDchain = Dchain +revVchain$chr.orderV = rev(revVchain$chr.orderV) +revDchain$chr.orderD = rev(revDchain$chr.orderD) + +if(useD){ + plotVD <- function(dat){ + if(length(dat[,1]) == 0){ + return() + } + img = ggplot() + + geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="gold", high="blue", na.value="white") + + ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + + xlab("D genes") + + ylab("V Genes") + + png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name))) + print(img) + dev.off() + write.table(x=acast(dat, Top.V.Gene~Top.D.Gene, value.var="Length"), file=paste("HeatmapVD_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA) + } + + VandDCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.D.Gene", "Sample")]) + + VandDCount$l = log(VandDCount$Length) + maxVD = data.frame(data.table(VandDCount)[, list(max=max(l)), by=c("Sample")]) + VandDCount = merge(VandDCount, maxVD, by.x="Sample", by.y="Sample", all.x=T) + VandDCount$relLength = VandDCount$l / VandDCount$max + + cartegianProductVD = expand.grid(Top.V.Gene = Vchain$v.name, Top.D.Gene = Dchain$v.name, Sample = unique(inputdata$Sample)) + + completeVD = merge(VandDCount, cartegianProductVD, all.y=TRUE) + completeVD = merge(completeVD, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE) + completeVD = merge(completeVD, Dchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE) + VDList = split(completeVD, f=completeVD[,"Sample"]) + + lapply(VDList, FUN=plotVD) +} + +plotVJ <- function(dat){ + if(length(dat[,1]) == 0){ + return() + } + cat(paste(unique(dat[3])[1,1])) + img = ggplot() + + geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="gold", high="blue", na.value="white") + + ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + + xlab("J genes") + + ylab("V Genes") + + png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name))) + print(img) + dev.off() + write.table(x=acast(dat, Top.V.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapVJ_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA) +} + +VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")]) + +VandJCount$l = log(VandJCount$Length) +maxVJ = data.frame(data.table(VandJCount)[, list(max=max(l)), by=c("Sample")]) +VandJCount = merge(VandJCount, maxVJ, by.x="Sample", by.y="Sample", all.x=T) +VandJCount$relLength = VandJCount$l / VandJCount$max + +cartegianProductVJ = expand.grid(Top.V.Gene = Vchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(inputdata$Sample)) + +completeVJ = merge(VandJCount, cartegianProductVJ, all.y=TRUE) +completeVJ = merge(completeVJ, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE) +completeVJ = merge(completeVJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE) +VJList = split(completeVJ, f=completeVJ[,"Sample"]) +lapply(VJList, FUN=plotVJ) + +if(useD){ + plotDJ <- function(dat){ + if(length(dat[,1]) == 0){ + return() + } + img = ggplot() + + geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_gradient(low="gold", high="blue", na.value="white") + + ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + + xlab("J genes") + + ylab("D Genes") + + png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name))) + print(img) + dev.off() + write.table(x=acast(dat, Top.D.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapDJ_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA) + } + + + DandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.D.Gene", "Top.J.Gene", "Sample")]) + + DandJCount$l = log(DandJCount$Length) + maxDJ = data.frame(data.table(DandJCount)[, list(max=max(l)), by=c("Sample")]) + DandJCount = merge(DandJCount, maxDJ, by.x="Sample", by.y="Sample", all.x=T) + DandJCount$relLength = DandJCount$l / DandJCount$max + + cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(inputdata$Sample)) + + completeDJ = merge(DandJCount, cartegianProductDJ, all.y=TRUE) + completeDJ = merge(completeDJ, revDchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE) + completeDJ = merge(completeDJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE) + DJList = split(completeDJ, f=completeDJ[,"Sample"]) + lapply(DJList, FUN=plotDJ) +} + + +# ---------------------- calculating the clonality score ---------------------- + +if("Replicate" %in% colnames(inputdata)) #can only calculate clonality score when replicate information is available +{ + if(clonality_method == "boyd"){ + samples = split(clonalityFrame, clonalityFrame$Sample, drop=T) + + for (sample in samples){ + res = data.frame(paste=character(0)) + sample_id = unique(sample$Sample)[[1]] + for(replicate in unique(sample$Replicate)){ + tmp = sample[sample$Replicate == replicate,] + clone_table = data.frame(table(tmp$clonaltype)) + clone_col_name = paste("V", replicate, sep="") + colnames(clone_table) = c("paste", clone_col_name) + res = merge(res, clone_table, by="paste", all=T) + } + + res[is.na(res)] = 0 + infer.result = infer.clonality(as.matrix(res[,2:ncol(res)])) + + write.table(data.table(infer.result[[12]]), file=paste("lymphclon_clonality_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=F) + + res$type = rowSums(res[,2:ncol(res)]) + + coincidence.table = data.frame(table(res$type)) + colnames(coincidence.table) = c("Coincidence Type", "Raw Coincidence Freq") + write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T) + } + } else { + write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T) + + clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")]) + clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")]) + clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count + clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample") + + ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15') + tcct = textConnection(ct) + CT = read.table(tcct, sep="\t", header=TRUE) + close(tcct) + clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T) + clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight + + ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")]) + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) + clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads) + ReplicateReads$squared = ReplicateReads$Reads * ReplicateReads$Reads + + ReplicatePrint <- function(dat){ + write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) + } + + ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) + lapply(ReplicateSplit, FUN=ReplicatePrint) + + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T) + + ReplicateSumPrint <- function(dat){ + write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) + } + + ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) + lapply(ReplicateSumSplit, FUN=ReplicateSumPrint) + + clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")]) + clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T) + clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow + clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2) + clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1) + + ClonalityScorePrint <- function(dat){ + write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) + } + + clonalityScore = clonalFreqCount[c("Sample", "Result")] + clonalityScore = unique(clonalityScore) + + clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"]) + lapply(clonalityScoreSplit, FUN=ClonalityScorePrint) + + clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")] + + + + ClonalityOverviewPrint <- function(dat){ + write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) + } + + clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample) + lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint) + } +} + +imgtcolumns = c("X3V.REGION.trimmed.nt.nb","P3V.nt.nb", "N1.REGION.nt.nb", "P5D.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "P3D.nt.nb", "N2.REGION.nt.nb", "P5J.nt.nb", "X5J.REGION.trimmed.nt.nb", "X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb") +if(all(imgtcolumns %in% colnames(inputdata))) +{ + print("found IMGT columns, running junction analysis") + newData = data.frame(data.table(PRODF)[,list(unique=.N, + VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), + P1=mean(.SD$P3V.nt.nb, na.rm=T), + N1=mean(.SD$N1.REGION.nt.nb, na.rm=T), + P2=mean(.SD$P5D.nt.nb, na.rm=T), + DEL.DH=mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T), + DH.DEL=mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T), + P3=mean(.SD$P3D.nt.nb, na.rm=T), + N2=mean(.SD$N2.REGION.nt.nb, na.rm=T), + P4=mean(.SD$P5J.nt.nb, na.rm=T), + DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), + Total.Del=( mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T)), + + Total.N=( mean(.SD$N1.REGION.nt.nb, na.rm=T) + + mean(.SD$N2.REGION.nt.nb, na.rm=T)), + + Total.P=( mean(.SD$P3V.nt.nb, na.rm=T) + + mean(.SD$P5D.nt.nb, na.rm=T) + + mean(.SD$P3D.nt.nb, na.rm=T) + + mean(.SD$P5J.nt.nb, na.rm=T))), + by=c("Sample")]) + print(newData) + newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1) + write.table(newData, "junctionAnalysisProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) + + newData = data.frame(data.table(UNPROD)[,list(unique=.N, + VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T), + P1=mean(.SD$P3V.nt.nb, na.rm=T), + N1=mean(.SD$N1.REGION.nt.nb, na.rm=T), + P2=mean(.SD$P5D.nt.nb, na.rm=T), + DEL.DH=mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T), + DH.DEL=mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T), + P3=mean(.SD$P3D.nt.nb, na.rm=T), + N2=mean(.SD$N2.REGION.nt.nb, na.rm=T), + P4=mean(.SD$P5J.nt.nb, na.rm=T), + DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T), + Total.Del=(mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T) + + mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T)), + Total.N=( mean(.SD$N1.REGION.nt.nb, na.rm=T) + + mean(.SD$N2.REGION.nt.nb, na.rm=T)), + Total.P=( mean(.SD$P3V.nt.nb, na.rm=T) + + mean(.SD$P5D.nt.nb, na.rm=T) + + mean(.SD$P3D.nt.nb, na.rm=T) + + mean(.SD$P5J.nt.nb, na.rm=T))), + by=c("Sample")]) + newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1) + write.table(newData, "junctionAnalysisUnProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) +} + +# ---------------------- AA composition in CDR3 ---------------------- + +AACDR3 = PRODF[,c("Sample", "CDR3.Seq")] + +TotalPerSample = data.frame(data.table(AACDR3)[, list(total=sum(nchar(as.character(.SD$CDR3.Seq)))), by=Sample]) + +AAfreq = list() + +for(i in 1:nrow(TotalPerSample)){ + sample = TotalPerSample$Sample[i] + AAfreq[[i]] = data.frame(table(unlist(strsplit(as.character(AACDR3[AACDR3$Sample == sample,c("CDR3.Seq")]), "")))) + AAfreq[[i]]$Sample = sample +} + +AAfreq = ldply(AAfreq, data.frame) +AAfreq = merge(AAfreq, TotalPerSample, by="Sample", all.x = T) +AAfreq$freq_perc = as.numeric(AAfreq$Freq / AAfreq$total * 100) + + +AAorder = read.table(sep="\t", header=TRUE, text="order.aa\tAA\n1\tR\n2\tK\n3\tN\n4\tD\n5\tQ\n6\tE\n7\tH\n8\tP\n9\tY\n10\tW\n11\tS\n12\tT\n13\tG\n14\tA\n15\tM\n16\tC\n17\tF\n18\tL\n19\tV\n20\tI") +AAfreq = merge(AAfreq, AAorder, by.x='Var1', by.y='AA', all.x=TRUE) + +AAfreq = AAfreq[!is.na(AAfreq$order.aa),] + +AAfreqplot = ggplot(AAfreq) +AAfreqplot = AAfreqplot + geom_bar(aes( x=factor(reorder(Var1, order.aa)), y = freq_perc, fill = Sample), stat='identity', position='dodge' ) +AAfreqplot = AAfreqplot + annotate("rect", xmin = 0.5, xmax = 2.5, ymin = 0, ymax = Inf, fill = "red", alpha = 0.2) +AAfreqplot = AAfreqplot + annotate("rect", xmin = 3.5, xmax = 4.5, ymin = 0, ymax = Inf, fill = "blue", alpha = 0.2) +AAfreqplot = AAfreqplot + annotate("rect", xmin = 5.5, xmax = 6.5, ymin = 0, ymax = Inf, fill = "blue", alpha = 0.2) +AAfreqplot = AAfreqplot + annotate("rect", xmin = 6.5, xmax = 7.5, ymin = 0, ymax = Inf, fill = "red", alpha = 0.2) +AAfreqplot = AAfreqplot + ggtitle("Amino Acid Composition in the CDR3") + xlab("Amino Acid, from Hydrophilic (left) to Hydrophobic (right)") + ylab("Percentage") + +png("AAComposition.png",width = 1280, height = 720) +AAfreqplot +dev.off() +write.table(AAfreq, "AAComposition.csv" , sep=",",quote=F,na="-",row.names=F,col.names=T) + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/report_clonality/circos/fonts.conf~ Mon Feb 13 09:08:46 2017 -0500 @@ -0,0 +1,8 @@ + +default = LTe50046.ttf +normal = LTe50046.ttf +bold = LTe50048.ttf +condensed = LTe50050.ttf +condensedbold = LTe50054.ttf +mono = pragmata.ttf +glyph = wingding.ttf
--- a/report_clonality/r_wrapper.sh Thu Feb 09 07:20:09 2017 -0500 +++ b/report_clonality/r_wrapper.sh Mon Feb 13 09:08:46 2017 -0500 @@ -222,54 +222,40 @@ for sample in $samples; do echo "${clonality_method}" - if [[ "${clonality_method}" == "old" ]] ; then - echo "in old" - clonalityScore="$(cat $outputDir/ClonalityScore_$sample.csv)" - echo "<div class='tabbertab' title='$sample'><table class='pure-table pure-table-striped'>" >> $outputFile - - if [[ "${clonality_method}" == "boyd" ]] ; then - echo "<thead><tr><th colspan='2'>Clonality Score: $clonalityScore</th></tr></thead>" >> $outputFile - fi - - #replicate,reads,squared - echo "<tr><td>Replicate ID</td><td>Number of Sequences</td></tr>" >> $outputFile - while IFS=, read replicate reads squared - do - echo "<tr><td>$replicate</td><td>$reads</td></tr>" >> $outputFile - done < $outputDir/ReplicateReads_$sample.csv - - #sum of reads and reads squared - while IFS=, read readsSum squaredSum - do - echo "<tr><td>Sum</td><td>$readsSum</td></tr>" >> $outputFile - done < $outputDir/ReplicateSumReads_$sample.csv - - echo "<tr><td></td><td></td></tr>" >> $outputFile - - #overview - echo "<tr><td>Number of replicates containing the coincidence</td><td>Number of sequences shared between replicates</td></tr>" >> $outputFile - while IFS=, read type count weight weightedCount - do - if [[ "$type" -eq "1" ]]; then - echo "<tr><td>$type</td><td>$count</td></tr>" >> $outputFile - else - echo "<tr><td><a href='coincidences_${sample}_${type}.txt'>$type</a></td><td>$count</td></tr>" >> $outputFile - fi - - done < $outputDir/ClonalityOverView_$sample.csv - echo "</table></div>" >> $outputFile - else - echo "in new" + + echo "<div class='tabbertab' title='$sample'><table class='pure-table pure-table-striped'>" >> $outputFile + + if [[ "${clonality_method}" == "boyd" ]] ; then clonalityScore="$(cat $outputDir/lymphclon_clonality_${sample}.txt)" - echo "<div class='tabbertab' title='$sample'>" >> $outputFile - echo "Lymphclon clonality score: <br />$clonalityScore<br /><br />" >> $outputFile - echo "<table border = 1>" >> $outputFile - while read type count - do - echo "<tr><td>$type</td><td>$count</td></tr>" >> $outputFile - done < $outputDir/lymphclon_coincidences_$sample.txt - echo "</table></div>" >> $outputFile + echo "<tr><td>Clonality Score: </td><td>$clonalityScore</td></tr>" >> $outputFile fi + + #replicate,reads,squared + echo "<tr><td>Replicate ID</td><td>Number of Sequences</td></tr>" >> $outputFile + while read replicate reads squared + do + echo "<tr><td>$replicate</td><td>$reads</td></tr>" >> $outputFile + done < $outputDir/ReplicateReads_$sample.txt + + #sum of reads and reads squared + while read readsSum squaredSum + do + echo "<tr><td>Sum</td><td>$readsSum</td></tr>" >> $outputFile + done < $outputDir/ReplicateSumReads_$sample.txt + + echo "<tr><td></td><td></td></tr>" >> $outputFile + + #overview + echo "<tr><td>Number of replicates containing the coincidence</td><td>Number of sequences shared between replicates</td></tr>" >> $outputFile + while read type count weight weightedCount + do + if [[ "$type" -eq "1" ]]; then + echo "<tr><td>$type</td><td>$count</td></tr>" >> $outputFile + else + echo "<tr><td><a href='coincidences_${sample}_${type}.txt'>$type</a></td><td>$count</td></tr>" >> $outputFile + fi + done < $outputDir/ClonalityOverView_$sample.txt + echo "</table></div>" >> $outputFile done cat $dir/naive_clonality.htm >> $outputFile @@ -295,7 +281,7 @@ echo "<table class='pure-table pure-table-striped' id='junction_table'> <caption>Unproductive mean</caption><thead><tr><th>Donor</th><th>Number of sequences</th><th>V.DEL</th><th>P1</th><th>N1</th><th>P2</th><th>DEL.D</th><th>D.DEL</th><th>P3</th><th>N2</th><th>P4</th><th>DEL.J</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><th>Median.CDR3</th><thead></tr><tbody>" >> $outputFile while read Sample unique VDEL P1 N1 P2 DELD DDEL P3 N2 P4 DELJ TotalDel TotalN TotalP median do - echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELD</td><td>$DDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>$median</td></tr>" >> $outputFile + echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELD</td><td>$DDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>-</td></tr>" >> $outputFile done < $outputDir/junctionAnalysisUnProd_mean_wD.txt echo "</tbody></table>" >> $outputFile @@ -309,7 +295,7 @@ echo "<table class='pure-table pure-table-striped' id='junction_table'> <caption>Unproductive median</caption><thead><tr><th>Donor</th><th>Number of sequences</th><th>V.DEL</th><th>P1</th><th>N1</th><th>P2</th><th>DEL.D</th><th>D.DEL</th><th>P3</th><th>N2</th><th>P4</th><th>DEL.J</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><th>Median.CDR3</th><thead></tr><tbody>" >> $outputFile while read Sample unique VDEL P1 N1 P2 DELD DDEL P3 N2 P4 DELJ TotalDel TotalN TotalP median do - echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELD</td><td>$DDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>$median</td></tr>" >> $outputFile + echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELD</td><td>$DDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>-</td></tr>" >> $outputFile done < $outputDir/junctionAnalysisUnProd_median_wD.txt echo "</tbody></table>" >> $outputFile @@ -325,7 +311,7 @@ echo "<table class='pure-table pure-table-striped' id='junction_table'> <caption>Unproductive mean</caption><thead><tr><th>Donor</th><th>Number of sequences</th><th>V.DEL</th><th>P1</th><th>N</th><th>P2</th><th>DEL.J</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><th>Median.CDR3</th><thead></tr><tbody>" >> $outputFile while read Sample unique VDEL P1 N1 P2 DELJ TotalDel TotalN TotalP median do - echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>$median</td></tr>" >> $outputFile + echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>-</td></tr>" >> $outputFile done < $outputDir/junctionAnalysisUnProd_mean_nD.txt echo "</tbody></table>" >> $outputFile @@ -339,7 +325,7 @@ echo "<table class='pure-table pure-table-striped' id='junction_table'> <caption>Unproductive median</caption><thead><tr><th>Donor</th><th>Number of sequences</th><th>V.DEL</th><th>P1</th><th>N</th><th>P2</th><th>DEL.J</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><th>Median.CDR3</th><thead></tr><tbody>" >> $outputFile while read Sample unique VDEL P1 N1 P2 DELJ TotalDel TotalN TotalP median do - echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>$median</td></tr>" >> $outputFile + echo "<tr><td>$Sample</td><td>$unique</td><td>$VDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELJ</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td><td>-/td></tr>" >> $outputFile done < $outputDir/junctionAnalysisUnProd_median_nD.txt echo "</tbody></table>" >> $outputFile
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/report_clonality/r_wrapper.sh~ Mon Feb 13 09:08:46 2017 -0500 @@ -0,0 +1,203 @@ +#!/bin/bash + +inputFile=$1 +outputDir=$3 +outputFile=$3/index.html #$2 +clonalType=$4 +species=$5 +locus=$6 +filterproductive=$7 +clonality_method=$8 +dir="$(cd "$(dirname "$0")" && pwd)" +useD="false" +if grep -q "$species.*${locus}D" "$dir/genes.txt" ; then + echo "species D region in reference db" + useD="true" +fi +echo "$species" +if [[ "$species" == *"custom"* ]] ; then + loci=(${locus//;/ }) + useD="true" + echo "${loci[@]}" + if [[ "${#loci[@]}" -eq "2" ]] ; then + useD="false" + fi +fi +mkdir $3 +cp $dir/genes.txt $outputDir +Rscript --verbose $dir/RScript.r $inputFile $outputDir $outputDir $clonalType "$species" "$locus" $filterproductive ${clonality_method} 2>&1 +cp $dir/tabber.js $outputDir +cp $dir/style.css $outputDir +cp $dir/script.js $outputDir +cp $dir/jquery-1.11.0.min.js $outputDir +samples=`cat $outputDir/samples.txt` +echo "<html><center><h1><a href='index.html'>Click here for the results</a></h1>Tip: Open it in a new tab (middle mouse button or right mouse button -> 'open in new tab' on the link above)<br />" > $2 +echo "<table border = 1>" >> $2 +echo "<thead><tr><th>Sample/Replicate</th><th>All</th><th>Productive</th><th>Unique Productive</th><th>Unproductive</th><th>Unique Unproductive</th></tr></thead>" >> $2 +while IFS=, read sample all productive perc_prod productive_unique perc_prod_un unproductive perc_unprod unproductive_unique perc_unprod_un + do + echo "<tr><td>$sample</td>" >> $2 + echo "<td>$all</td>" >> $2 + echo "<td>$productive (${perc_prod}%)</td>" >> $2 + echo "<td>$productive_unique (${perc_prod_un}%)</td>" >> $2 + echo "<td>$unproductive (${perc_unprod}%)</td>" >> $2 + echo "<td>$unproductive_unique (${perc_unprod_un}%)</td></tr>" >> $2 +done < $outputDir/productive_counting.txt +echo "</table border></center></html>" >> $2 + +echo "productive_counting.txt" +echo "<html><head><title>Report on:" >> $outputFile +for sample in $samples; do + echo " $sample" >> $outputFile +done +echo "</title><script type='text/javascript' src='jquery-1.11.0.min.js'></script>" >> $outputFile +echo "<script type='text/javascript' src='tabber.js'></script>" >> $outputFile +echo "<script type='text/javascript' src='script.js'></script>" >> $outputFile +echo "<link rel='stylesheet' type='text/css' href='style.css'></head>" >> $outputFile +echo "<div class='tabber'><div class='tabbertab' title='Gene frequencies'>" >> $outputFile + +echo "<img src='CDR3LengthPlot.png'/><br />" >> $outputFile +echo "<img src='VFPlot.png'/>" >> $outputFile +if [[ "$useD" == "true" ]] ; then + echo "<img src='DFPlot.png'/>" >> $outputFile +fi +echo "<img src='JFPlot.png'/>" >> $outputFile +echo "<img src='VPlot.png'/>" >> $outputFile +if [[ "$useD" == "true" ]] ; then + echo "<img src='DPlot.png'/>" >> $outputFile +fi +echo "<img src='JPlot.png'/>" >> $outputFile +echo "<img src='AAComposition.png'/></div>" >> $outputFile + +count=1 +echo "<div class='tabbertab' title='Heatmaps'><div class='tabber'>" >> $outputFile +for sample in $samples; do + echo "<div class='tabbertab' title='$sample'><table border='1'><tr>" >> $outputFile + if [[ "$useD" == "true" ]] ; then + echo "<td><img src='HeatmapVD_$sample.png'/></td>" >> $outputFile + fi + echo "<td><img src='HeatmapVJ_$sample.png'/></td>" >> $outputFile + if [[ "$useD" == "true" ]] ; then + echo "<td><img src='HeatmapDJ_$sample.png'/></td>" >> $outputFile + fi + echo "</tr></table></div>" >> $outputFile + count=$((count+1)) +done +echo "</div></div>" >> $outputFile + +#echo "<div class='tabbertab' title='Interactive'><svg class='chart'></svg><script src='http://d3js.org/d3.v3.min.js'></script></div>" >> $outputFile + +hasReplicateColumn="$(if head -n 1 $inputFile | grep -q 'Replicate'; then echo 'Yes'; else echo 'No'; fi)" +echo "$hasReplicateColumn" +#if its a 'new' merged file with replicate info +if [[ "$hasReplicateColumn" == "Yes" && "$clonalType" != "none" ]] ; then + echo "<div class='tabbertab' title='Clonality'><div class='tabber'>" >> $outputFile + for sample in $samples; do + echo "${clonality_method}" + if [[ "${clonality_method}" == "old" ]] ; then + echo "in old" + clonalityScore="$(cat $outputDir/ClonalityScore_$sample.csv)" + echo "<div class='tabbertab' title='$sample'><table border='1'>" >> $outputFile + echo "<tr><td colspan='4'>Clonality Score: $clonalityScore</td></tr>" >> $outputFile + + #replicate,reads,squared + echo "<tr><td>Replicate ID</td><td>Number of Reads</td><td>Reads Squared</td><td></td></tr>" >> $outputFile + while IFS=, read replicate reads squared + do + + echo "<tr><td>$replicate</td><td>$reads</td><td>$squared</td><td></td></tr>" >> $outputFile + done < $outputDir/ReplicateReads_$sample.csv + + #sum of reads and reads squared + while IFS=, read readsSum squaredSum + do + echo "<tr><td>Sum</td><td>$readsSum</td><td>$squaredSum</td></tr>" >> $outputFile + done < $outputDir/ReplicateSumReads_$sample.csv + + #overview + echo "<tr><td>Coincidence Type</td><td>Raw Coincidence Freq</td><td>Coincidence Weight</td><td>Coincidences, Weighted</td></tr>" >> $outputFile + while IFS=, read type count weight weightedCount + do + echo "<tr><td>$type</td><td>$count</td><td>$weight</td><td>$weightedCount</td></tr>" >> $outputFile + done < $outputDir/ClonalityOverView_$sample.csv + echo "</table></div>" >> $outputFile + else + echo "in new" + clonalityScore="$(cat $outputDir/lymphclon_clonality_${sample}.csv)" + echo "<div class='tabbertab' title='$sample'>" >> $outputFile + echo "Lymphclon clonality score: <br />$clonalityScore<br /><br />" >> $outputFile + echo "<table border = 1>" >> $outputFile + while IFS=, read type count + do + echo "<tr><td>$type</td><td>$count</td></tr>" >> $outputFile + done < $outputDir/lymphclon_coincidences_$sample.csv + echo "</table></div>" >> $outputFile + fi + done + echo "</div></div>" >> $outputFile +fi + +hasJunctionData="$(if head -n 1 $inputFile | grep -q '3V-REGION trimmed-nt'; then echo 'Yes'; else echo 'No'; fi)" + +if [[ "$hasJunctionData" == "Yes" ]] ; then + echo "<div class='tabbertab' title='Junction Analysis'>" >> $outputFile + echo "<table border='1' id='junction_table'> <caption>Productive</caption><thead><tr><th>Sample</th><th>count</th><th>VH.DEL</th><th>P1</th><th>N1</th><th>P2</th><th>DEL.DH</th><th>DH.DEL</th><th>P3</th><th>N2</th><th>P4</th><th>DEL.JH</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><thead></tr><tbody>" >> $outputFile + while IFS=, read Sample unique VHDEL P1 N1 P2 DELDH DHDEL P3 N2 P4 DELJH TotalDel TotalN TotalP + do + echo "<tr><td>$Sample</td><td>$unique</td><td>$VHDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELDH</td><td>$DHDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJH</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td></tr>" >> $outputFile + done < $outputDir/junctionAnalysisProd.csv + echo "</tbody></table>" >> $outputFile + + echo "<table border='1' id='junction_table'> <caption>Unproductive</caption><thead><tr><th>Sample</th><th>count</th><th>VH.DEL</th><th>P1</th><th>N1</th><th>P2</th><th>DEL.DH</th><th>DH.DEL</th><th>P3</th><th>N2</th><th>P4</th><th>DEL.JH</th><th>Total.Del</th><th>Total.N</th><th>Total.P</th><thead></tr><tbody>" >> $outputFile + while IFS=, read Sample unique VHDEL P1 N1 P2 DELDH DHDEL P3 N2 P4 DELJH TotalDel TotalN TotalP + do + echo "<tr><td>$Sample</td><td>$unique</td><td>$VHDEL</td><td>$P1</td><td>$N1</td><td>$P2</td><td>$DELDH</td><td>$DHDEL</td><td>$P3</td><td>$N2</td><td>$P4</td><td>$DELJH</td><td>$TotalDel</td><td>$TotalN</td><td>$TotalP</td></tr>" >> $outputFile + done < $outputDir/junctionAnalysisUnProd.csv + echo "</tbody></table>" >> $outputFile + + echo "</div>" >> $outputFile +fi + +echo "<div class='tabbertab' title='Comparison'><table border='1'><tr><th>ID</th><th>Include</th></tr>" >> $outputFile +for sample in $samples; do + echo "<tr><td>$sample</td><td><input type='checkbox' onchange=\"javascript:compareAdd('$sample')\" id='compare_checkbox_$sample'/></td></tr>" >> $outputFile +done +echo "</table><div name='comparisonarea'>" >> $outputFile +echo "<table><tr id='comparison_table_vd'></tr></table>" >> $outputFile +echo "<table><tr id='comparison_table_vj'></tr></table>" >> $outputFile +echo "<table><tr id='comparison_table_dj'></tr></table>" >> $outputFile +echo "</div></div>" >> $outputFile + +echo "<div class='tabbertab' title='Downloads'>" >> $outputFile +echo "<table border='1'>" >> $outputFile +echo "<tr><th>Description</th><th>Link</th></tr>" >> $outputFile +echo "<tr><td>The dataset used to generate the frequency graphs and the heatmaps (Unique based on clonaltype, $clonalType)</td><td><a href='allUnique.csv'>Download</a></td></tr>" >> $outputFile +echo "<tr><td>The dataset used to calculate clonality score (Unique based on clonaltype, $clonalType)</td><td><a href='clonalityComplete.csv'>Download</a></td></tr>" >> $outputFile + +echo "<tr><td>The dataset used to generate the CDR3 length frequency graph</td><td><a href='CDR3LengthPlot.csv'>Download</a></td></tr>" >> $outputFile + +echo "<tr><td>The dataset used to generate the V gene family frequency graph</td><td><a href='VFFrequency.csv'>Download</a></td></tr>" >> $outputFile +if [[ "$useD" == "true" ]] ; then + echo "<tr><td>The dataset used to generate the D gene family frequency graph</td><td><a href='DFFrequency.csv'>Download</a></td></tr>" >> $outputFile +fi +echo "<tr><td>The dataset used to generate the J gene family frequency graph</td><td><a href='JFFrequency.csv'>Download</a></td></tr>" >> $outputFile + +echo "<tr><td>The dataset used to generate the V gene frequency graph</td><td><a href='VFrequency.csv'>Download</a></td></tr>" >> $outputFile +if [[ "$useD" == "true" ]] ; then + echo "<tr><td>The dataset used to generate the D gene frequency graph</td><td><a href='DFrequency.csv'>Download</a></td></tr>" >> $outputFile +fi +echo "<tr><td>The dataset used to generate the J gene frequency graph</td><td><a href='JFrequency.csv'>Download</a></td></tr>" >> $outputFile +echo "<tr><td>The dataset used to generate the AA composition graph</td><td><a href='AAComposition.csv'>Download</a></td></tr>" >> $outputFile + +for sample in $samples; do + if [[ "$useD" == "true" ]] ; then + echo "<tr><td>The data used to generate the VD heatmap for $sample.</td><td><a href='HeatmapVD_$sample.csv'>Download</a></td></tr>" >> $outputFile + fi + echo "<tr><td>The data used to generate the VJ heatmap for $sample.</td><td><a href='HeatmapVJ_$sample.csv'>Download</a></td></tr>" >> $outputFile + if [[ "$useD" == "true" ]] ; then + echo "<tr><td>The data used to generate the DJ heatmap for $sample.</td><td><a href='HeatmapDJ_$sample.csv'>Download</a></td></tr>" >> $outputFile + fi +done + +echo "</table>" >> $outputFile +echo "</div></html>" >> $outputFile