Mercurial > repos > davidvanzessen > combined_immune_repertoire_imgt
changeset 1:d2b3bcabb478 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Mon, 09 Dec 2013 06:08:52 -0500 |
| parents | e71c59b72669 |
| children | 1c5927d0a4ce |
| files | RScript.r combined.sh combined_imgt.xml igblastmerge.py igparse.pl imgtconvert.py imgtconvert.sh immunerepertoirecombined_imgt/RScript.r immunerepertoirecombined_imgt/combined.sh immunerepertoirecombined_imgt/combined_imgt.xml immunerepertoirecombined_imgt/igblastmerge.py immunerepertoirecombined_imgt/igparse.pl immunerepertoirecombined_imgt/imgtconvert.py immunerepertoirecombined_imgt/imgtconvert.sh |
| diffstat | 14 files changed, 1939 insertions(+), 1908 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RScript.r Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,295 @@ +#options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) + +args <- commandArgs(trailingOnly = TRUE) + +inFile = args[1] +outFile = args[2] +outDir = args[3] +clonalType = args[4] + +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/") +} +require(ggplot2) +if (!("plyr" %in% rownames(installed.packages()))) { + install.packages("plyr", repos="http://cran.xl-mirror.nl/") +} +require(plyr) + +if (!("data.table" %in% rownames(installed.packages()))) { + install.packages("data.table", repos="http://cran.xl-mirror.nl/") +} +library(data.table) + + +test = read.table(inFile, sep="\t", header=TRUE, fill=T) + +test = test[test$Sample != "",] + +test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene) +test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene) +test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene) + +#test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":")) +test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":")) + +PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ] + +NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ] + +#PRODF = PROD[ -1] + +PRODF = PROD + +#PRODF = unique(PRODF) +PRODF = PRODF[!duplicated(PRODF$VDJCDR3), ] + +PRODFV = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Sample", "Top.V.Gene")]) +PRODFV$Length = as.numeric(PRODFV$Length) +Total = 0 +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=.N), by=c("Sample", "Top.D.Gene")]) +PRODFD$Length = as.numeric(PRODFD$Length) +Total = 0 +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=.N), by=c("Sample", "Top.J.Gene")]) +PRODFJ$Length = as.numeric(PRODFJ$Length) +Total = 0 +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)) + +V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53") +tcV = textConnection(V) +Vchain = read.table(tcV, sep="\t", header=TRUE) +PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) +close(tcV) + +D = c("v.name\tchr.orderD\nIGHD1-1\t1\nIGHD2-2\t2\nIGHD3-3\t3\nIGHD6-6\t4\nIGHD1-7\t5\nIGHD2-8\t6\nIGHD3-9\t7\nIGHD3-10\t8\nIGHD4-11\t9\nIGHD5-12\t10\nIGHD6-13\t11\nIGHD1-14\t12\nIGHD2-15\t13\nIGHD3-16\t14\nIGHD4-17\t15\nIGHD5-18\t16\nIGHD6-19\t17\nIGHD1-20\t18\nIGHD2-21\t19\nIGHD3-22\t20\nIGHD4-23\t21\nIGHD5-24\t22\nIGHD6-25\t23\nIGHD1-26\t24\nIGHD7-27\t25") +tcD = textConnection(D) +Dchain = read.table(tcD, sep="\t", header=TRUE) +PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE) +close(tcD) + + +J = c("v.name\tchr.orderJ\nIGHJ1\t1\nIGHJ2\t2\nIGHJ3\t3\nIGHJ4\t4\nIGHJ5\t5\nIGHJ6\t6") +tcJ = textConnection(J) +Jchain = read.table(tcJ, sep="\t", header=TRUE) +PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) +close(tcJ) + +setwd(outDir) + +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") + +png("VPlot.png",width = 1280, height = 720) +pV +dev.off(); + +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") + +png("DPlot.png",width = 800, height = 600) +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") + +png("JPlot.png",width = 800, height = 600) +pJ +dev.off(); + +revVchain = Vchain +revDchain = Dchain +revVchain$chr.orderV = rev(revVchain$chr.orderV) +revDchain$chr.orderD = rev(revDchain$chr.orderD) + +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() +} + +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(test$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() + } + 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() +} + +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(test$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) + +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() +} + +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(test$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) + + +sampleFile <- file("samples.txt") +un = unique(test$Sample) +un = paste(un, sep="\n") +writeLines(un, sampleFile) +close(sampleFile) + + +if("Replicate" %in% colnames(test)) +{ + clonalityFrame = PROD + clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":")) + clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] + clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")]) + 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", "VDJCDR3")]) + ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) + 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(Reads), ReadsSquaredSum=sum(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$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) +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/combined.sh Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,107 @@ +#!/bin/bash + +export IGDATA=/home/david/tmp/ncbi-igblast-1.0.0; + +clonalType=${@:(-3):1} +html=${@:(-2):1} +imageDir=${@:(-1):1} +dataCount=`expr $# - 3` +inputData=${@:(1):dataCount} +dir="$(cd "$(dirname "$0")" && pwd)" +array=("$@") +echo "<html><h3>Progress</h3><table><tr><td>info</td></tr>" > $html +echo "<tr><td>-----------------------------------</td></tr>" >> $html + +id=${inputData[0]} +forwardSlash="/" +mergerInput=() +count=0 +for current in $inputData; do + if [[ "$current" != *"$forwardSlash"* ]]; then + id=$current + count=0 + mergerInput+=($id) + continue + fi + fileName=$(basename $current) + convertedFileName="${fileName%.*}" + convertedFileName="$PWD/$convertedFileName.converted" + bash $dir/imgtconvert.sh $current $id $count $convertedFileName + mergerInput+=($convertedFileName) + count=$((count+1)) +done +wait + + + +echo "<tr><td>-----------------------------------</td></tr>" >> $html +echo "<tr><td>merging</td></tr>" >> $html + +python $dir/igblastmerge.py ${mergerInput[*]} --output $PWD/merged.txt + +echo "<tr><td>done</td></tr>" >> $html +echo "<tr><td>-----------------------------------</td></tr>" >> $html +echo "<tr><td>plotting</td></tr>" >> $html + + +inputFile=$PWD/merged.txt +outputFile=$html +outputDir=$imageDir +mkdir $outputDir +Rscript --verbose $dir/RScript.r $inputFile $outputDir $outputDir $clonalType 2>&1 +echo "<html>" > $outputFile +echo "<img src='VPlot.png'/>" >> $outputFile +echo "<img src='DPlot.png'/>" >> $outputFile +echo "<img src='JPlot.png'/>" >> $outputFile + +samples=`cat $outputDir/samples.txt` +count=1 +echo "<table border='1'><caption><a href='allUnique.tsv'><h3>$clonalType</h3></a></caption>" >> $outputFile +hasReplicateColumn="$(if head -n 1 $inputFile | grep -q 'Replicate'; then echo 'Yes'; else echo 'No'; fi)" +for sample in $samples; do + clonalityScore="$(cat $outputDir/ClonalityScore_$sample.csv)" + echo "<tr><td colspan='3' height='100'></td></tr>" >> $outputFile + echo "<tr><td colspan='3'><h1>$sample</h1></td></tr>" >> $outputFile + + echo "$hasReplicateColumn" + #if its a 'new' merged file with replicate info + if [[ "$hasReplicateColumn" == "Yes" ]] ; then + echo "<tr><td colspan='3'><a href='clonality_$sample.tsv'><h2>Clonality Score: $clonalityScore</h2></a></td></tr>" >> $outputFile + + #replicate,reads,squared + echo "<tr><td colspan='3'><table border='1'><tr><th>Replicate ID</th><th>Number of Reads</th><th>Reads Squared</th></tr>" >> $outputFile + while IFS=, read replicate reads squared + do + + echo "<tr><td><a href='clonality_${sample}_$replicate.tsv'>$replicate</a></td><td>$reads</td><td>$squared</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 + + echo "</table></td></tr>" >> $outputFile + + #overview + echo "<tr><td colspan='3'><table border='1'><tr><th>Coincidence Type</th><th>Raw Coincidence Freq</th><th>Coincidence Weight</th><th>Coincidences, Weighted</th></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></td></tr>" >> $outputFile + fi + + echo "<tr><td><h2>V-D Heatmap:</h2></td><td><h2>V-J Heatmap:</h2></td><td><h2>D-J Heatmap:</h2></td></tr><tr>" >> $outputFile + mv "$outputDir/HeatmapVD_$sample.png" "$outputDir/VD_$sample.png" + echo "<td><img src='VD_$sample.png'/></td>" >> $outputFile + mv "$outputDir/HeatmapVJ_$sample.png" "$outputDir/VJ_$sample.png" + echo "<td><img src='VJ_$sample.png'/></td>" >> $outputFile + mv "$outputDir/HeatmapDJ_$sample.png" "$outputDir/DJ_$sample.png" + echo "<td><img src='DJ_$sample.png'/></td></tr>" >> $outputFile + count=$((count+1)) +done +echo "</table>" >> $outputFile + +echo "</html>" >> $outputFile
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/combined_imgt.xml Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,38 @@ +<tool id="immunerepertoirecombined_imgt" name="Immunerepertoire combined" version="1.0"> + <description>pipeline for IMGT</description> + <command interpreter="bash"> + combined.sh + #for $i, $f in enumerate($patients) + ${f.id} + #for $j, $g in enumerate($f.samples) + ${g.sample} + #end for + #end for + "$clonaltype_select" $out_file $out_file.files_path + </command> + <inputs> + <repeat name="patients" title="Patients" min="1" default="1"> + <repeat name="samples" title="Sample" min="1" default="1"> + <param name="sample" type="data" label="Sample to Process" /> + </repeat> + <param name="id" type="text" label="ID" /> + </repeat> + <param name="clonaltype_select" type="select" label="Clonal Type Definition"> + <option value="Top.V.Gene,CDR3.Seq">Top.V.Gene, CDR3.Seq</option> + <option value="Top.V.Gene,CDR3.Seq.DNA">Top.V.Gene, CDR3.Seq.DNA</option> + <option value="Top.V.Gene,Top.J.Gene,CDR3.Seq">Top.V.Gene, Top.J.Gene, CDR3.Seq</option> + <option value="Top.V.Gene,Top.J.Gene,CDR3.Seq.DNA">Top.V.Gene, Top.J.Gene, CDR3.Seq.DNA</option> + <option value="Top.V.Gene,Top.D.Gene,Top.J.Gene,CDR3.Seq.DNA">Top.V.Gene, Top.D.Gene, Top.J.Gene, CDR3.Seq.DNA</option> + </param> + </inputs> + <outputs> + <data format="html" name="out_file" /> + </outputs> + <help> + The entire Immune Repertoire pipeline as a single tool, input several FASTA files, give them an ID and it will BLAST, parse, merge and plot them. + </help> + <requirements> + <requirement type="package" version="1.0.0">igBlastn</requirement> + </requirements> +</tool> +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/igblastmerge.py Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,51 @@ +import sys +# error +def stop_err( msg ): + sys.stderr.write( "%s\n" % msg ) + sys.exit() + +# main +def main(): + args = sys.argv[1:-2] + print args + try: + o = open(sys.argv[-1], 'w') + i = open(args[1], 'r') + separator = "\t" + newline = "\n" + line = i.readline() + #write the header + o.write(line[:line.rfind(newline)] + separator + "Sample" + separator + "Replicate" + newline) + i.close() + + current = 1 + sampleID = args[0] + count = 1 + + while True: + f = open(args[current], 'r') + line = f.readline() + line = f.readline() + while line: + o.write(line[:line.rfind(newline)] + separator + sampleID + separator + str(count) + newline) + line = f.readline() + f.close() + + if current >= (len(args) - 1): + break + if args[current + 1].find("/") is -1: + sampleID = args[current + 1] + current += 1 + count = 1 + else: + count += 1 + current += 1 + o.close() + + except Exception, ex: + stop_err('Error running new_column.py\n' + str(ex)) + sys.exit(0) + +if __name__ == "__main__": + print sys.argv + main()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/igparse.pl Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,1252 @@ +#!/usr/bin/perl +=head1 IGBLAST_simple.pl + +This version (1.4) has been heavily adapted since the original program was first created back in October 2012. +Bas Horsman (EMC, Rotterdam, The Netherlands) has contributed with minor - though important - code changes. + +From V 1.2 onwards a 'Change Log' is included at the end of the program + +=head2 Usage + +Requires no modules in general use; the Data::Dumper (supplied as part of the Perl Core module set) might be useful for debugging/adjustment +as it allows inspection of the data stores. + +The program takes a text file of the + + ./IGBLAST_simple.pl igBLASTOutput.txt <-optional: index of record to process-> + +Supply the text version of the igBLAST report in the format as in the example below. +The extra command line arugment is the record number (aka. BLAST report) to process. +If 0 or absent all are processed, if supplied that record (base 1) is processed and the program dies afterwards. + +=head2 Example Input + +A standard igBLAST record or set of them in a file; this being typical: + + BLASTN 2.2.27+ + + +Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A. +Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. +Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of +protein database search programs", Nucleic Acids Res. 25:3389-3402. + + + +Database: human_gl_V; human_gl_D; human_gl_J + 674 sequences; 179,480 total letters + + + +Query= HL67IUI01D26LR length=433 xy=1559_1437 region=1 +run=R_2012_04_10_11_57_56_ + +Length=433 + Score E +Sequences producing significant alignments: (Bits) Value + +lcl|IGHV3-30*04 330 2e-92 +lcl|IGHV3-30-3*01 330 2e-92 +lcl|IGHV3-30*01 327 2e-91 +lcl|IGHD3-16*01 14.4 11 +lcl|IGHD3-16*02 14.4 11 +lcl|IGHD1-14*01 12.4 43 +lcl|IGHJ4*02 78.3 1e-18 +lcl|IGHJ5*02 70.3 4e-16 +lcl|IGHJ4*01 68.3 2e-15 + + +Domain classification requested: imgt + + +V(D)J rearrangement summary for query sequence (Top V gene match, Top D gene match, Top J gene match, Chain type, V-J Frame, Strand): +IGHV3-30*04 IGHD3-16*01 IGHJ4*02 VH In-frame + + +V(D)J junction details (V end, V-D junction, D region, D-J junction, J start). Note that possible overlapping nucleotides at VDJ junction (i.e, nucleotides that could be assigned to either joining gene segment) are indicated in parentheses (i.e., (TACT)) but are not included under V, D, or J gene itself +AGAGA TATGAGCCCCATCATGACA ACGTTTG CCGGAA ACTAC + +Alignment summary between query and top germline V gene hit (from, to, length, matches, mismatches, gaps, percent identity) +FWR1 27 38 12 11 1 0 91.7 +CDR1 39 62 24 22 2 0 91.7 +FWR2 63 113 51 50 1 0 98 +CDR2 114 137 24 23 1 0 95.8 +FWR3 138 251 114 109 5 0 95.6 +CDR3 (V region only) 252 259 8 7 1 0 87.5 +Total N/A N/A 233 222 11 0 95.3 + + +Alignments + + <----FWR1--><----------CDR1--------><-----------------------FWR2------ + W A A S G F T F N T Y A V H W V R Q A P G K G + Query_1 27 TGGGCAGCCTCTGGATTCACCTTCAATACCTATGCTGTGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGC 96 +V 95.3% (222/233) IGHV3-30*04 64 ..T......................G..G.......A................................. 133 + C A A S G F T F S S Y A M H W V R Q A P G K G +V 95.7% (221/231) IGHV3-30-3*01 64 ..T......................G..G.......A................................. 133 +V 94.8% (221/233) IGHV3-30*01 64 ..T......................G..G.......A................................. 133 + + ----------------><----------CDR2--------><---------------------------- + L E W V A V I S Y D G S N K N Y A D S V K G R F + Query_1 97 TGGAGTGGGTGGCAGTTATATCATATGATGGAAGCAATAAAAACTACGCAGACTCCGTGAAGGGCCGATT 166 +V 95.3% (222/233) IGHV3-30*04 134 ..................................T......T............................ 203 + L E W V A V I S Y D G S N K Y Y A D S V K G R F +V 95.7% (221/231) IGHV3-30-3*01 134 .........................................T............................ 203 +V 94.8% (221/233) IGHV3-30*01 134 .A................................T......T............................ 203 + + ---------------------------FWR3--------------------------------------- + T I S R D N S K N T L Y L Q M N S L R V E D T + Query_1 167 CACCATCTCCAGAGACAATTCCAAGAACACGTTATATCTGCAAATGAACAGCCTGAGAGTTGAGGACACG 236 +V 95.3% (222/233) IGHV3-30*04 204 ...............................C.G.........................C.......... 273 + T I S R D N S K N T L Y L Q M N S L R A E D T +V 95.7% (221/231) IGHV3-30-3*01 204 ...............................C.G.........................C.......... 273 +V 94.8% (221/233) IGHV3-30*01 204 ...............................C.G.........................C.......... 273 + + --------------> + A V Y Y C T R D M S P I M T T F A G N Y W G Q + Query_1 237 GCTGTTTATTACTGTACGAGAGATATGAGCCCCATCATGACAACGTTTGCCGGAAACTACTGGGGCCAGG 306 +V 95.3% (222/233) IGHV3-30*04 274 .....G.........G.......----------------------------------------------- 296 + A V Y Y C A R +V 95.7% (221/231) IGHV3-30-3*01 274 .....G.........G.....------------------------------------------------- 294 +V 94.8% (221/233) IGHV3-30*01 274 .....G.........G.......----------------------------------------------- 296 +D 100.0% (7/7) IGHD3-16*01 12 ------------------------------------------.......--------------------- 18 +D 100.0% (7/7) IGHD3-16*02 12 ------------------------------------------.......--------------------- 18 +D 100.0% (6/6) IGHD1-14*01 8 -------------------------------------------------......--------------- 13 +J 100.0% (39/39) IGHJ4*02 10 -------------------------------------------------------............... 24 +J 100.0% (35/35) IGHJ5*02 17 -----------------------------------------------------------........... 27 +J 97.4% (38/39) IGHJ4*01 10 -------------------------------------------------------.............A. 24 + + + G T L V T V S S + Query_1 307 GAACCCTGGTCACCGTCTCCTCAG 330 +J 100.0% (39/39) IGHJ4*02 25 ........................ 48 +J 100.0% (35/35) IGHJ5*02 28 ........................ 51 +J 97.4% (38/39) IGHJ4*01 25 ........................ 48 + + +Lambda K H + 1.10 0.333 0.549 + +Gapped +Lambda K H + 1.08 0.280 0.540 + +Effective search space used: 64847385 + + +Query= HL67IUI01EQMLY length=609 xy=1826_1636 region=1 +run=R_2012_04_10_11_57_56_ + + +...etc... + +=head2 Example Output + + +Example output from the data above sent: + $ ./IGBLAST_simple.pl igBLASTOutput.txt 1 + D: Request to process just record '1' received + D: printOUTPUTData: Running + D: printOUTPUTData: HEADER Printout requested '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 Found How' + OUTPUT: # 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 Found How + D: ID is: 'HL67IUI01D26LR' + D: Minimum base marked-up (27) - aka. $AlignmentStart; maximum: (259) + D: Starting Search for CDR3 + D: markUpCDR3: Passed Parameters '251, 27, TGGGG....GG., WG.G' (& AA & DNA sequence) + D: markUpCDR3: returning: 223, 282, MOTIF_FOUND_IN_BOTH, (3) [NB: offset of :'+ 27' + D: CDR3 was found by pattern matching: 'MOTIF_FOUND_IN_BOTH' (250, 309) + D: Top Hits (raw)= 'IGHV3-30*04 IGHD3-16*01 IGHJ4*02 VH In-frame +' + D: Top Hits (parsed)= 'IGHV3-30*04, IGHD3-16*01, IGHJ4*02, VH, In-frame, +' + D: printOUTPUTData: Running + OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH + +=head4 Usage notes: + +Designed to be easy to "grep -v D:" or "grep OUTPUT:" for to select the parts you need: + + ./IGBLAST_simple.pl igBLASTOutput.txt 1 | grep OUTPUT: + + OUTPUT: # 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 Found How + OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01EQMLY In-frame IGHV4-39*01 IGHD2-8*01 IGHJ3*02 GGSISSSSYY 29 IYHSGST 20 CARDATYYSNGFDIWGQG 53 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01CDCLP Out-of-frame IGHV3-23*01 IGHD3-3*01 IGHJ4*02 FSNYAM 16 SGSGDRTY 23 AKAD*FLEWLFRIGDGERLLGPGN 72 MOTIF_FOUND_IN_DNA + OUTPUT: HL67IUI01AHRNH N/A IGHV3-33*01 N/A N/A WIHLQ*LW 23 YGMMEVI 23 NOT_FOUND + OUTPUT: HL67IUI01DZZ1V Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 GFTFDKYA 23 ILASG 20 LYCASEGDIVASELLSTGARV 62 MOTIF_FOUND_IN_DNA + OUTPUT: HL67IUI01DTR2Y Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 LDSPLTNM 23 LYLPVV 20 TVRVRGT*WLRSF*VLGPG 59 MOTIF_FOUND_IN_DNA + OUTPUT: HL67IUI01EQL3S In-frame IGHV7-4-1*02 IGHD6-19*01 IGHJ6*02 GYTFRTFT 23 INTNTGTP 23 CAKESGTGSAHFFYGMDVWGQG 65 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01AFG46 In-frame IGLV2-34*01 N/A IGHJ4*02 NOT_FOUND + OUTPUT: HL67IUI01EFFKO In-frame IGHV3-11*01 IGHD6-6*01 IGHJ4*02 GFTFSDYY 23 ISYSGGTI 23 CARASGAARHRPLDYWGQG 56 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01B18SG In-frame IGHV3-33*01 IGHD5-12*01 IGHJ4*02 VRQA 11 KYYANSVK 23 RLGGFDYWGQGTLVTVSS 53 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01D6LER In-frame IGHV1-24*01 IGHD3-22*01 IGHJ4*02 GYSLNELS 23 PDPEDDE 23 TVQPSRITMMAVVITRIHWGASGARE 76 MOTIF_FOUND_IN_DNA + OUTPUT: HL67IUI01CYCLF N/A IGHV4-39*01 N/A N/A GGSISSSSYY 29 IYYSGST 20 NOT_FOUND + OUTPUT: HL67IUI01B4LEE In-frame IGHV7-4-1*02 IGHD6-19*01 IGHJ6*02 GYTFRTFT 23 INTNTGTP 23 CAKESGTGSAHFFYGMDVWGQG 65 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01A4KW4 Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 LDSPLTNM 23 LYLPVV 20 TVRVRGT*WLRSF*IWGQG 58 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01E05BV In-frame IGHV1-24*01 IGHD3-22*01 IGHJ2*01 GYSLNELS 23 PDPEDDE 23 NOT_FOUND + OUTPUT: HL67IUI01CVVKY In-frame IGHV1-3*01 IGHD2-15*01 IGHJ1*01 NOT_FOUND + OUTPUT: HL67IUI01CN5P2 In-frame IGHV7-4-1*02 IGHD2-21*02 IGHJ5*02 GYSITDYG 23 LNTRTGNP 23 CAVKDARDFVSWGQG 44 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01DUUJ5 In-frame IGHV3-21*01 IGHD1-7*01 IGHJ4*02 GYTFSTYS 23 ISSSSAYR 23 CARDIRLELRDWGQG 44 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01E1AIR Out-of-frame IGHV4-39*01 N/A IGHJ3*01 WGLHRRW**L 29 FVS*RAPR 23 NOT_FOUND + OUTPUT: HL67IUI01CCZ8D Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 GFTFDKYA 23 ILASGR 20 YCASEGDIVASELLSTGARE 58 MOTIF_FOUND_IN_DNA + OUTPUT: HL67IUI01BT9IR N/A IGHV3-21*02 N/A N/A NOT_FOUND + OUTPUT: HL67IUI01COTO0 Out-of-frame IGHV4-39*01 N/A IGHJ3*01 GGFIGGGDNF 29 LYHDGRPA 23 NOT_FOUND + OUTPUT: HL67IUI01D994O In-frame IGHV7-4-1*02 IGHD2-21*02 IGHJ5*02 GYSITDYG 23 LNTRTGNP 23 CAVKDARDFVSWGQG 44 MOTIF_FOUND_IN_BOTH + OUTPUT: HL67IUI01A08CJ In-frame IGHV4-39*01 IGHD6-13*01 IGHJ5*02 GGSISSSSYY 29 IYYTWEH 21 CERARRGSSWGQLVRPLGPG 62 MOTIF_FOUN + + + + OUTPUT: # 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 Found How + OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH + ...etc... + +=head4 Also, combined grep & sed: + + $ ./IGBLAST_simple.pl igBLASTOutput.txt | grep OUTPUT: | sed 's/OUTPUT:\t//' + +=cut + +=head3 CDR3 Patterns: + +We use these two variables to try to identify the end of the CDR3 region if igBLAST doesn't report it directly: + + my $DNACDR3_Pat = "TGGGG....GG."; + my $AASequenceMotifPattern = "WG.G"; + +They are treated as regex's when tested (so use "." to mean any DNA base, rather than 'N' or 'X'). + +[NB: These are original patterns used for testing, check the code for the current ones.] + +=cut + +my $DNACDR3_Pat = "TGGGG....GG."; +my $AACDR3_Pat = "WG.G"; + +use strict; +use Data::Dumper; +# Set this as to number of the result (aka "record") you want to process or 0 for all: +my $ProcessRecord =0; +if (defined $ARGV[1]) { $ProcessRecord = pop @ARGV; } #Also accept from the command line: +if ($ProcessRecord != 0) { print "D: Request to process just record '$ProcessRecord' received\n"; } + +#Adjust the record separator: +$/="Query= "; +my $Record=0; # A simple counter, that we might not use. +#Force-loaded header / version information: +my $Header = <>; +#At the moment we don't use this - so dump it immediately: +$Header = undef; +#print "D: Force-loaded header / version information: '$Header'\n"; + +#Print the Header for the output line (we need this once, at the start) +print &printOUTPUTData ({"HEADER" => 1})."\n"; + +while (<>) + { +=head4 First check - should we be processing this record at all? + +=cut + $Record++; #Increment the record counter: + #Do we process this record - or all records? + if ($ProcessRecord != $Record && $ProcessRecord != 0) + { next; } #We need to increment the record counter before we increment + +=head4 Setup the output line storage and print the header: + +We enter this initially and work to change it: + + $DomainBoundaries{"CDR3"}{"FoundHow"} = "NOT_FOUND"; + +=cut + + my %OUTPUT_Data; #To collect data for the output line in + #Assume the first and work to find better: + $OUTPUT_Data{"CDR3 Found How"} = "NOT_FOUND"; + #The whole record - one per read - is now stored in $_ + my @Lines =split (/[\r\n]+/,$_); # split on windows/linux/mac new lines + + #If you are interested enable either of the next lines depending on how curious you are as to how the splitting went: + #print "D: Record #$Record\n"; print $_; print "\n---------\n"; + print "D: ''$Lines[0]'\nD: ...etc...'\nD: ############\n"; + +=head3 Get the ID + +Quite easy: the first field on the first line: + + Query= HL67IUI01DTR2Y length=577 xy=1452_0984 region=1 + +=cut + + (my $ID) = $Lines[0]=~ m/^(\S+)/; + unless (defined $ID && $ID ne "") + { # So a near total failure...? + $OUTPUT_Data{"ID"} = "Unknown"; + print &printOUTPUTData (\%OUTPUT_Data)."\n"; + next; #No ID is terminal for this record + } + else + { + print "D: ID is: '$ID'\n"; + $OUTPUT_Data{"ID"} = $ID; + } +=head3 Declare the variables we will need here in the next few sections to store data + +=cut + + my $CurrentRegion; + my $RegionMarkup; + + #So we can sync the coordinated of the alignment up to the domains found: + my $Query_Start = -1; my $Query_End = -1; + + #Where on the Query Sequence (i.e. the 454 read) does the alignment start & stop? + my $ThisQueryStart =-1; my $ThisQueryEnd =-1; #Think $ThisQueryEnd isn't used at the moment. + my $DNAQuerySequence =""; #The actual DNA Query sequence... + my $AAQuerySequence = ""; + + #As this changes with the alleles identified: + my $CurrentAASequence; + #The main storage variables + + my %Alginments; my %Alleles; + my %DomainBoundaries; + +=head2 Stanza 1: Get the general structure of the sequence identified + +=head3 Method 1: Use the table supplied + +Technically this valid for the top hit...realistically this is the only information we have reported to us +so we use this or nothing. This is fine for the top hit which is likely what we are interested in....but for the 2nd or 3rd? Who knows! + +Targets this block: + + Alignment summary between query and top germline V gene hit (from, to, length, matches, mismatches, gaps, percent identity) + FWR1 167 240 75 72 2 1 96 + CDR1 241 264 24 20 4 0 83.3 + FWR2 265 315 51 48 3 0 94.1 + CDR2 316 336 24 15 6 3 62.5 + FWR3 337 450 114 106 8 0 93 + CDR3 (V region only) 451 454 4 4 0 0 100 + Total N/A N/A 292 265 23 4 90.8 + +Then we split out the lines inside it in a second scanning step - less optimal but easier to read: + + FWR1 167 240 75 72 2 1 96 + CDR1 241 264 24 20 4 0 83.3 + FWR2 265 315 51 48 3 0 94.1 + CDR2 316 336 24 15 6 3 62.5 + FWR3 337 450 114 106 8 0 93 + CDR3 (V region only) 451 454 4 4 0 0 100 + +into: + + (Section, from, to, length, matches, mismatches, gaps, percent identity) + +=head3 Method 2: Use the table supplied + +The other way to do this is to split the graphical markup out of the alignment. +This works for _any_ reported alignment, not just the top hits: + +In the main alignment table processing section collect the information, collect the information: + + #Is region mark-up: + if ($#InfoColumns == -1 && $#AlignmentColumns ==0) + { +# print ": Region Markup detected\n"; + $RegionMarkup = $RegionMarkup.$AlignmentPanel; #Collect the information, then re-synthesise it at the end of record + next; + } + +Then afterwards when all the region was collected, process it like this: +#Pad the CDER3 region: + + #Remove the trailing spaces: + $RegionMarkup =~ s/ *$//g; + #Calculate the length of the CDR3 region so we can add it in: + my $CDR3PaddingNeeded = ($Query_End-$Query_Start)-length ($RegionMarkup) -length ("<-CDR3>")+1; + #Build up the CDR3 region, the 'x' operator is very helpful here (implict foreach loop): + $RegionMarkup = $RegionMarkup."<-CDR3"."-" x $CDR3PaddingNeeded. ">"; + #print "D: Need to pad with:'$CDR3PaddingNeeded' characters\n"; + + #Now really process it: + my $C_Pos = 0; + my @Domains = split (/(<*-*...[123]-*>*)/,$RegionMarkup); # + foreach my $C_Domain (@Domains) + { + if (length ($C_Domain) <=0) {next;} + my $DomainStart= $C_Pos; + my $DomainEnd = $DomainStart + length ($C_Domain)-1; + my ($DomainType) = $C_Domain =~ m/(...[123])/; +# print "D: $DomainType \t($DomainStart-$DomainEnd=",$DomainEnd-$DomainStart,"):\t$C_Domain\n"; + $DomainBoundaries{$DomainType}{"Start"} = $DomainStart; + $DomainBoundaries{$DomainType}{"End"} = $DomainEnd; + $C_Pos = $DomainEnd+1; + } + +The two pieces of code are interchangable; the table version as used below, is neater, easier to understand and works nicely. +Why stress? + + +=head3 The end of the FWR3 is the start of CDR3? + +This is an assumption made. Hence the two variables: + + my $MaxDomainReported =0 ; # In nts / bps + my $FWR3_Found_Flag = 0; # Did we find the end of the FWR3 - which is the start of the CDR3. Set to 'false' initially. + + $MaxDomainBaseFound + +=cut + my $MaxDomainBaseFound =0 ; # In nts / bps + my $AlignmentStart ; # In nts /bp #Alternative name would be: '$MinDomainBaseFound'; set to null until primed +# my $FWR3_Found_Flag = 0; # Did we find the end of the FWR3 - which is the start of the CDR3. Set to 'false' initially. + + (my @StructureSummaryTable) = returnLinesBetween (\@Lines, "Alignment summary", "Total" ); +#Enable the next line if you want the raw data we are going to parse in this section: + #print Dumper @StructureSummaryTable; + foreach my $C_Section (@StructureSummaryTable) + { + my ($DomainType, $DomainStart, $DomainEnd, $SectLength, $Matches, $Mismatches, $Gaps, $PID) = split (/\t+/,$C_Section); + #print "D: Domain type: '$DomainType'\n"; + #$DomainType =~ s/ .*$//g; + $DomainBoundaries{$DomainType}{"Start"} = $DomainStart; + $DomainBoundaries{$DomainType}{"End"} = $DomainEnd; + +#So we can do a reality check on the length / start of the CDR3 if we have to go looking: + if ($MaxDomainBaseFound <= $DomainEnd) + { $MaxDomainBaseFound = $DomainEnd; } #Store the maximum base found + if ($AlignmentStart eq undef or $AlignmentStart >= $DomainStart) + { $AlignmentStart = $DomainStart; } + } +#print Dumper %DomainBoundaries; +#die "HIT BLOCK\n"; + +=head3 Did we find the CDR3 region specifically? + +If we did fine; otherwise try to find it using the FWR3 region if we found that; otherwise give up. + +=cut + print "D: Minimum base marked-up ($AlignmentStart) - aka. \$AlignmentStart; maximum: ($MaxDomainBaseFound)\n"; + +#my @WantedSections = qw (V D J); + +=head2 Second Stanza: Parse the main Alignment Table + +=head3 Get the table, then determine the character at which to split the 'Info' & 'Alignment' panels. + +As this is a little involved and comparamentalises nicely we sub-contract this to two functions"" + + (my @Table) = returnLinesBetween (\@Lines, "Alignment", "Lambda" ); + my $PanelSplitPoint = findSplitPoint (\@Table); #Why can't they just use a fixed field width or a tab as a delimiter? + +=cut + (my @Table) = returnLinesBetween (\@Lines, "Alignment", "Lambda" ); + my $PanelSplitPoint = findSplitPoint (\@Table); #Why can't they just use a fixed field width or a tab as a delimiter? +#If you are interested, enable this line: +# print "D: The info panel was detected at: '$splitPoint'\n"; + +=head3 + +=cut + + +foreach my $C_Line (0..$#Table) + { + +=head3 Call the line type we find: There are 4: + +These are distinguished by the number of fields (one or mores spacer is a field separator) in the Info & Alignment Panels (see values in brackets) + + | <- This split is ~40 chars. from the start of the line + * InfoPanel * | * Alignment Panel * + : is a "Blank" line (-1,-1) + <----FWR1--><----------CDR1--------><-----------------------FWR2------ : is "Region Markup" (-1,0) + W A A S G F T F N T Y A V H W V R Q A P G K G : is "AA Sequence" (-1, >=0) + Query_1 27 TGGGCAGCCTCTGGATTCACCTTCAATACCTATGCTGTGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGC 96 : is "DNA Sequence" (2,1) + V 95.3% (222/233) IGHV3-30*04 64 ..T......................G..G.......A................................. 133 : is "" " + +So we split 40 chars in and then the two parts on spaces. + + +=cut + +# print "D: (sub) Line in parsed table: '$C_Line': \n"; + + my ($InfoPanel, $AlignmentPanel) = $Table[$C_Line] =~ /^(.{$PanelSplitPoint})(.*)$/; + + my @InfoColumns = split (/\s+/,$InfoPanel); + my @AlignmentColumns = split (/\s+/,$AlignmentPanel); + +#If you want to see how the line is being split enable either of these next two lines; the 2nd is more detailed than the first +# print "D: Line: $C_Line/t Number of Columns (Info, Alignment): \t$#InfoColumns \t $#AlignmentColumns\n"; +# print "D: For '$C_Line' \t line in the table there are parts: '$InfoPanel' [$#InfoColumns], '$AlignmentPanel [$#AlignmentColumns]'\n"; + +#Populate this so we can step through it + +=head4 Is a blank line: +=cut + if ($#InfoColumns == -1 && $#AlignmentColumns == -1) + { +# print ": Blank\n"; + next; + } #For now I think we just skip - is not needed (though might be implict mark-up) + +=head4 Is region mark-up: +=cut + if ($#InfoColumns == -1 && $#AlignmentColumns ==0) + { +# print ": Region Markup detected\n"; + $RegionMarkup = $RegionMarkup.$AlignmentPanel; #Collect the information, then re-synthesise it at the end of record + next; + } +=head4 Is query DNA Sequence: +=cut + if ($#InfoColumns == 2 && $#AlignmentColumns ==1) + { +# print ": DNA Query Sequence\n"; + #Detect the two coordinatates of alignment against the query sequence: (last two numbers of the two 'panels') + ($ThisQueryStart) = $InfoPanel =~ / (\d+) *$/; + ($ThisQueryEnd) = $AlignmentPanel =~ / (\d+) *$/; + my ($ThisDNASeq) = $AlignmentPanel =~ /^(.*?) /; + #If you want to know what we just found: + #print "D: This DNA Sequence: '$ThisDNASeq'\n"; + $DNAQuerySequence = $DNAQuerySequence. $ThisDNASeq; #Add it on to whatever we already have. + #Move the needle if there are smaller / greater; otherwise prime the 'needles': + if ($ThisQueryStart < $Query_Start or $Query_Start == -1) + { $Query_Start = $ThisQueryStart; } + if ($ThisQueryEnd > $Query_End or $Query_End == -1) + { $Query_End = $ThisQueryEnd; } +# print ": Query DNA Sequence detected This line: ($ThisQueryStart, $ThisQueryEnd) & Maximally: ($Query_Start, $Query_End)\n"; + next; + } +=head4 Is AA Sequence: + +This is complicated as it Need to decide whether this is the sequence of the read or that of the original V / D / J regions: + --------------> + A V Y Y C T R D M S P I M T T F A G N Y W G Q << Want this + Query_1 237 GCTGTTTATTACTGTACGAGAGATATGAGCCCCATCATGACAACGTTTGCCGGAAACTACTGGGGCCAGG 306 + V 95.3% (222/233) IGHV3-30*04 274 .....G.........G.......----------------------------------------------- 296 + A V Y Y C A R + V 95.7% (221/231) IGHV3-30-3*01 274 .....G.........G.....------------------------------------------------- 294 + + ...etc... + G T L V T V S S << Want this + Query_1 307 GAACCCTGGTCACCGTCTCCTCAG 330 + +To solve this we peak at the next line that it has the tag "Query" in it (we assume the line exists...) + +=cut + + if ($#InfoColumns == -1 && $#AlignmentColumns >=-1) + { + unless ($Table[$C_Line+1] =~ /Query/) { next; } #Is the next line the DNA sequence ? + # +# print ": AA sequence\n"; + + + $CurrentAASequence = $AlignmentPanel; + #print "D: Panel Split Point = $PanelSplitPoint, '$AlignmentPanel'\n"; + $CurrentAASequence =~ s/^ {$PanelSplitPoint}//; + #print "D: '$AAQuerySequence'\n"; +# print "D: Current AA Sequence: \t'$CurrentAASequence'\n"; + $AAQuerySequence = $AAQuerySequence.$CurrentAASequence; #Store the elongating AA Sequence as well + next; + } +=head4 Is Alignment: +=cut + if ($#InfoColumns == 4 && $#AlignmentColumns ==1) + { + #Not acutally interesting to us for this version of the parser. Delete ultimately? + next; + } + +#Is weird! Don't recognise it! + + warn "Weird! Don't recongnise this: '$ID' [$#InfoColumns,$#AlignmentColumns]// '",$Lines[$C_Line],15,"...'\n"; + } #End main iteration loop for alignment parsing. + + +=head2 The CDR3 is noted as problematic. Can we identify it? + +=cut + print "D: Starting Search for CDR3\n"; + #Do have the end of the FWR3 but not the CDR3? If so then it is worth trying to find the CDR3, otherwise...nothing we can do at this point + if (exists ($DomainBoundaries{"FWR3"}{"End"}) + && $AlignmentStart !=0 + && not (exists $DomainBoundaries{"CDR3"}{"End"}) ) #Guess we need to go looking for the end then... + { + #print "D: Placing call to markUpCDR3\n"; + my ($CDR3_Start, my $CDR3_End, my $CDR3_Found_Tag) = markUpCDR3 ($DNAQuerySequence, $AAQuerySequence, + $DomainBoundaries{"FWR3"}{"End"}, $AlignmentStart, + $DNACDR3_Pat, $AACDR3_Pat); + if ($CDR3_Start !=0 && $CDR3_End !=0) + { + $DomainBoundaries{"CDR3"}{"Start"} = $CDR3_Start; + $DomainBoundaries{"CDR3"}{"End"} = $CDR3_End ; + $DomainBoundaries{"CDR3"}{"FoundHow"} = $CDR3_Found_Tag; + print "D: CDR3 was found by pattern matching: '$CDR3_Found_Tag' ($CDR3_Start, $CDR3_End)\n"; + } + else + { print "D: CDR3 was not found [either by igBLAST or by pattern matching]\n"; + $DomainBoundaries{"CDR3"}{"FoundHow"} = "NOT_FOUND"; + } + } + else + { #Was reported by igBLAST + print "D: Found the FWR3 from the Domain Boundary Table\n"; + $DomainBoundaries{"CDR3"}{"FoundHow"} = "IGBLAST_NATIVE"; + } + +#print Dumper %DomainBoundaries; + +=head2 Get the top VDJ regions: + +=cut + +=head2 Extract General Features: + +=cut + (my $TopHit) = $_ =~ m/V-J Frame, Strand\):\n(.*?)\n/s; + print "D: Top Hits (raw)= '$TopHit' \n"; + my ($Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match, $Chain, $VJFrame, $Strand) = split (/\t/,$TopHit); + print "D: Top Hits (parsed)= '$Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match, $Chain, $VJFrame, $Strand'\n"; + +=head2 Store the V / D / J Genes used + +=cut + + if (defined $Top_V_gene_match && $Top_V_gene_match ne "") + { $OUTPUT_Data{"Top V Gene"} = $Top_V_gene_match; } + + if (defined $Top_D_gene_match && $Top_D_gene_match ne "") + { $OUTPUT_Data{"Top D Gene"} = $Top_D_gene_match; } + + if (defined $Top_J_gene_match && $Top_J_gene_match ne "") + { $OUTPUT_Data{"Top J Gene"} = $Top_J_gene_match; } + + if (defined $Strand && $Strand ne "") + { $OUTPUT_Data{"Strand"} = $Strand;} + +=head4 Preamble: ID, Frame, and V / D / J used: + +=cut + #Do a reality check: if we didn't get an ID, then skip: + unless (defined (defined $ID) && $ID ne "" && + defined $VJFrame && $VJFrame ne "") + { + print &printOUTPUTData (\%OUTPUT_Data)."\n"; + next; + } + +#Ok, so we have data...most likely: + #print "OUTPUT:\t",join ("\t", $ID, $VJFrame, $Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match); + + if (defined $VJFrame && defined $ID && $VJFrame ne "" && $ID ne "") + { $OUTPUT_Data{"VDJ Frame"} = $VJFrame;} + else + { + print &printOUTPUTData (\%OUTPUT_Data)."\n"; + next; + }#REALLY? We didn't find anything? Oh well, move to next record + +=head4 CDR1 + +=cut + #Remember that the alignment starts at the FWR1 start, not nt =0 on the read, hence we substract this off all future AA (& DNA coordinates) + + my $AlignmentOffset = $DomainBoundaries{"FWR1"}{"Start"}; + +# print "D: AA Seqeunce is: '$AAQuerySequence'\n"; + if (exists $DomainBoundaries{"CDR1"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start + { +# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... + my $CDR1Start = $DomainBoundaries{"CDR1"}{"Start"}; + my $CDR1End = $DomainBoundaries{"CDR1"}{"End"}; + my $CDR1_Length = $CDR1End - $CDR1Start; +# print "D: CDR1 $CDR1Start $CDR1End = $CDR1_Length\n"; + #Remember that the alignment starts at the FWR1 start, not nt =0 on the read + my $CDR1_Seq_AA = substr ($AAQuerySequence, $CDR1Start - $AlignmentOffset, $CDR1_Length); +# print "D: '$CDR1_Seq_AA'\n"; + $CDR1_Seq_AA =~ s/ //g; + my $CDR1_Seq_AA_Length = length ($CDR1_Seq_AA); + #Add this data to the output store specifically: + $OUTPUT_Data{"CDR1 Seq"} = $CDR1_Seq_AA; + $OUTPUT_Data{"CDR1 Length"} = $CDR1_Length; + } + #What happens if there is no CDR1 found? Leave blank - the output routine can handle this + +=head4 CDR2 + +=cut + + if (exists $DomainBoundaries{"CDR2"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start + { +# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... + my $CDR2Start = $DomainBoundaries{"CDR2"}{"Start"}; + my $CDR2End = $DomainBoundaries{"CDR2"}{"End"}; + my $CDR2_Length = $CDR2End - $CDR2Start; + my $CDR2_Seq_AA = substr ($AAQuerySequence, $CDR2Start - $AlignmentOffset , $CDR2_Length); + $CDR2_Seq_AA =~ s/ //g; + my $CDR2_Seq_AA_Length = length ($CDR2_Seq_AA); + #Add this data to the output store specifically: + $OUTPUT_Data{"CDR2 Seq"} = $CDR2_Seq_AA; + $OUTPUT_Data{"CDR2 Length"} = $CDR2_Length; + } + #What happens if there is no CDR2 found? Leave blank - the output routine can handle this. + +=head4 CDR3 + +=cut + if (exists $DomainBoundaries{"CDR3"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start + { +# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... + my $CDR3Start = $DomainBoundaries{"CDR3"}{"Start"}; + my $CDR3End = $DomainBoundaries{"CDR3"}{"End"}; + my $CDR3_Length = $CDR3End - $CDR3Start; # This variable isn't used - delete it when safe to do so + my $CDR3_Seq_AA = substr ($AAQuerySequence, $CDR3Start - $AlignmentOffset, $CDR3_Length); + my $CDR3_Seq_DNA = substr ($DNAQuerySequence, $CDR3Start - $AlignmentOffset, $CDR3_Length); + $CDR3_Seq_AA =~ s/ //g; + $CDR3_Seq_DNA =~ s/ //g; + my $CDR3_Seq_AA_Length = length ($CDR3_Seq_AA); + my $CDR3_Seq_DNA_Length = length ($CDR3_Seq_DNA); + #Add this data to the output store specifically: + $OUTPUT_Data{"CDR3 Seq"} = $CDR3_Seq_AA; + $OUTPUT_Data{"CDR3 Length"} = $CDR3_Seq_AA_Length; + $OUTPUT_Data{"CDR3 Seq DNA"} = $CDR3_Seq_DNA; + $OUTPUT_Data{"CDR3 Length DNA"} = $CDR3_Seq_DNA_Length; + #And in the case of the CDR3 how we found it: + $OUTPUT_Data{"CDR3 Found How"} = $DomainBoundaries{"CDR3"}{"FoundHow"}; + } + #What happens if there is no CDR3 found? Leave blank - the output routine can handle this. +#die "HIT BLOCK\n"; +#End of the record; output the data we have collected and move on. +print &printOUTPUTData (\%OUTPUT_Data)."\n"; +} + + + +############ +sub returnLinesBetween { +=head3 SUB: returnLinesBetween ({reference to array Index array}, {regex for top of section}, {regex for bottom of section}) + +When passed a reference to an array and two strings - interpreted as REGEX's - will return the lines of the Array +that are bounded by these tags. + +If either of the tags are not found - or are found in the wrong order - then a null list is returned. + +=cut + +my ($Text_ref, $TopTag, $BotTag) = @_; + +my @Table; +#The two boundary conditions at which we will cut the table: +#print "D: [returnLinesBetween]: '$TopTag, $BotTag'\n"; +#How we record these: +my $AlignmentLine_Top=0; my $AlignmentLine_Bot=0; + +my $LineIndex=-1; #-1 As the loop increments this line counter first, then does its checks. +#If you care: +#print "D: Lines of text passed: $$#Lines\n"; + +#Iterate through until we find what we are looking for or run out of text to search: +while (($AlignmentLine_Bot ==0 or $AlignmentLine_Top==0) && $LineIndex <=$#{$Text_ref}) + { + $LineIndex++; + #Enable if you need to care: +# print "D: Line Index = $LineIndex\n"; + + if ($$Text_ref[$LineIndex] =~ m/$TopTag/) + { + $AlignmentLine_Top = $LineIndex; +# print "D: [returnLinesBetween]: TopTag found in Line: '$$Text_ref[$LineIndex]'\n"; #Enable if you are interested + } + if ($$Text_ref[$LineIndex] =~ m/$BotTag/) + { + $AlignmentLine_Bot = $LineIndex; +# print "D: [returnLinesBetween]: Bottom Tag found in Line: '$$Text_ref[$LineIndex]'\n"; #Enable if you are interested + } + } +#Reality check: did we find anything? If not then we return null. +if ($AlignmentLine_Top ==0 && $AlignmentLine_Bot ==0) + { return; } +#Again, enable if you care: +#print "D: [returnLinesBetween] Lines for section table: '$AlignmentLine_Top to $AlignmentLine_Bot'\n"; + +#We want the lines one down and one up - so polish these. +$AlignmentLine_Top++; $AlignmentLine_Bot--; + +#Return as an array slice: +return (@$Text_ref[$AlignmentLine_Top .. $AlignmentLine_Bot]); +} +############ + +sub findSplitPoint +{ +=head2 sub: $PanelBoundaryCahracter = findSplitPoint (\@Table) + +When passed a table with the alignment in it makes an educated guess as to the precise split point to +spearate the 'info' and 'alignment' panels. +This is a right olde faff because the field / panel boundaries change. + + ' Query_6 167 GAGGTGCAGTTGTTGGAGTCTGGGGGAGGCTTGGCACAGCC-GGGGGGTCCCTGAGACTCTCCTGTGCAG 235' + ' Query_6 236 CCTCTGGATTCACCTTTGACAAATATGCCATGACCTGGGTCCGCCAGGCTCCAGGGAAGGGTCTGGAGTG 305' + ' Query_6 306 GGTCTCAACTATACTTGCCAGTGGTCG---CACAGACGACGCAGACTCCGTGAAGGGCCGGTTTGCCATC 372' + ' Query_6 373 TCCAGAGACAATTCCAAGAACACTCTGTATCTGCAAATGAACAGCCTGAGAGTCGAGGACACGGCCCTTT 442' + ' Query_6 443 ATTACTGTGCGAGTGAGGGGGACATAGTGGCTTCGGAGCTTTTGAGTACTGGGGCCAGGGAAACCTGGTC 512' +MOTIF_FOUND_IN_AA +i.e to contain just ATGC + "X" bases & the gap "-" character but not the "." character (found in the alingment proper) and have 4 fields in total + +Returns either -1 or the location of the panel boundary, issues a warning and returns -1 if is the most frequent boundary +because the pattern match has been failing more often that it suceeded. + +=cut +#A rough guess is 38 for normal sequences, 48 for reversed ones: + +my $SplitPos = 0; + +(my $Table_ref) = @_; #Get the reference to the table +my @DNALines; #We populate this for mining in the next section +foreach my $C_Line (@{$Table_ref}) + { + #print "D: $C_Line\n"; +# (my $SplitLine) = $C_Line; + #Split on consecutive tabs or spaces: + my @LineFields = split (/[\t\s]+/,$C_Line); + #print "D: Split Line: '",join (",",@LineFields),"' : $#LineFields\n"; + unless ( $LineFields[3] =~ m/[^\.]/ + && $LineFields[3] =~ m/[ATGCX]{20,}/ + && $#LineFields==4) + { next; } +#Enable if you want to know the lines we think are the DNA Query strings: + #print "D: DNA Line: '$C_Line'\n"; + push @DNALines, $C_Line; #Note it down + } + +my %PanelBounds; #Will contain the positions of the panel boundaries + +foreach my $C_DNALine (@DNALines) + { + #print "D: '$C_DNALine'\n"; + $C_DNALine =~ m/[ATGC-]+ \d+$/; #Match the DNA string and the indexingMOTIF_FOUND_IN_AA numbers afterwards, allow gap characters. + my $MatchPos = $-[0]; #This is the position of the start of the last match because we can't get the index() function to work + #(my $MatchPos) = index ($C_DNALine, / [ATGCX-]{20}/,0); + #print "D: '$C_DNALine' DNA panel starts at:'$MatchPos'\n"; + $PanelBounds{$MatchPos}++; + } +#Sort the hash values in order and then return the most frequent (will offer some resistance to the occasion pattern failure) +#The brackets around "($SplitPos)" are really necessary it seems. +($SplitPos) = (sort { $a <=> $b } keys %PanelBounds); +#If you want +#print Dumper %PanelBounds; +#Tell people if we are having difficultlty: +if ($SplitPos == -1) { warn "Couldn't identify the panel boundaries\n"; } +#print "D: $SplitPos: Returning the split position of: '$SplitPos'\n"; +return $SplitPos; +} + + +## +# +# +### + + + + + +##### +# +# +##### +sub markUpCDR3 +{ +=head3 Sub: (Start, End, Found How) = markUpCDR3 (DNASeq, AASeq, FWR3 End, FWR1 Offset, DNA Regex, AA Regex) + +Tries to identify the end of the CDR3 using the DNA and RNA Sequence patterns MOTIF_FOUND_IN_AAsupplied. The CDR3 is assumed to start +at the end of the FWR3. +To reduce FP matches only the sequences (DNA & AA) after the FWR3 are tested with the pattern. +The position of the first matching pattern is reported. + +=head4 Fuller Usage: + +my ($CDR3_Start, my $CDR3_End) = markUpCDR3 ($DNAQuerySequence, $AAQuerySequence, + $DomainBoundaries{"FWR3"}{"End"}, $DomainBoundaries{"FWR1"}{"Start"}, + $DNACDR3_Pat, $AACDR3_Pat); + + + +=head4 Returned Values + +If the CDR3 was found then we we signal like this: + + $MotifFound ==0 : Nope, didn't find either motif + $MotifFound ==1 : Found at the DNA level, not the AA level + $MotifFound ==2 : Found at the the AA level, not the DNA level + $MotifFound ==3 : Found at the the AA level & the DNA level + +(Also remember that if the FWR3 region couldn't be identified in the sequence there is a 4th option: not tested; this routine isn't called therefore) + +The Start and Ends returned are from the first sucessful match (MotifFound==3): though hopefully they are the same. +Formally the test order is: + + 1) DNA + 2) AA + +i.e. DNA bp locations have priority. + +Technically the locations are determined by a regex match then the $+[0] array (i.e. the end of the pattern match). +See pages like this: http://stackoverflow.com/questions/87380/how-can-i-find-the-location-of-a-regex-match-in-perl for an explanation. + +=head3 Manipulation of AA patternsMOTIF_FOUND_IN_AA + +Note that patterns are assumed to require white space inserting in them between the letters. +This could be a serious limitation + + +=cut + +#Get the parameters passed: +my ($DNA, $AA, $FWR3_End, $FWR1_Start, $DNAPat, $AAPat) = @_; +print "D: markUpCDR3: Passed Parameters '$FWR3_End, $FWR1_Start, $DNAPat, $AAPat' (& AA & DNA sequence)\n"; + + +#Setup our return values: +my $Start = 0; my $End =0; my $MotifFound = 0; +my $How; #Literally How the motif was found (or not if blank) + + +=head4 Prepare the sequences and the patterns for use + +Specifically: trim off the start of the AA & DNA string already allocated to other CDRs or FWRs + +Add in spaces into the AA regex pattern because we can't get regex-ex freespacing mode i.e. "$Var =~ m/$AAPat/x" working. + + +We take the "-1" as the CropPoint position to include the previous 3 nucleotides / AAs; remember to add this back on +in position calculations. + + +=cut + +#Because igBLAST doesn't always report from the start of the read (primers and things are upstream): + +my $CropPoint = $FWR3_End - $FWR1_Start - 1 ; +#print "D: markUpCDR3: Crop point is: '$CropPoint'\n"; + +#print "D: markUpCDR3: Cropping point is: '$CropPoint' characters from start\n"; +#We trim off the parts we expect to find the CDR3 motifs in leaving at extra 3nts on to allow for base miss-calling: + +my $AA_Trimmed = substr ($AA, $CropPoint); +my $DNA_Trimmed = substr ($DNA ,$CropPoint); +#print "D: markUpCDR3: AA = '$AA' (untrimmed)\nD: markUpCDR3: TR = '$AA_Trimmed' (Trimmed) ", length ($AA_Trimmed)," nts long\n"; +#print "D: markUpCDR3: Testing: AA = '$AA_Trimmed', DNA = '$DNA_Trimmed'\n"; + +#This lovely hack is to account for the spaces in the AA sequence and we can't get the "$Var =~ m/$AAPat/x" working +my $AAPat_Spaced; +foreach my $C_Char (0..length($AAPat)-1) #The -1 is because we don't want trailing spaces until the next nt -> AA translation. + { $AAPat_Spaced = $AAPat_Spaced.'\s+'.substr ($AAPat,$C_Char,1); } +#And write this back into the main pattern we were passed: +$AAPat = $AAPat_Spaced; + +#temp hack: +#$AA_Trimmed = $AA; +my $MotifFound=0; #So we can record which patterns we found +my $MotifPositionDNA =-1; +my $MotifPositionAA =-1; + +#print "D: markUpCDR3: Pattern: '$AAPat_Spaced'\n"; +=head4 At DNA level: "TGG GGx xxx GGx" [+1] + +=cut + +#print "D: markUpCDR3: '$DNA_Trimmed' (Trimmed DNA string)\n"; + +if ($DNA_Trimmed =~ m/$DNAPat/) + { + $MotifPositionDNA = $+[0]; #Just the easiest way to do this in Perl +# print "D: markUpCDR3:: Found Motif match on DNA at bp: '$MotifPositionDNA'\n"; + $MotifFound = $MotifFound + 1; + #Any more matches further on? + my $LaterString = substr ($DNA_Trimmed, $MotifPositionDNA); +# print "D: markUpCDR3: '$AA_Trimmed' (AA Trimmed string)\n"; +# print "D: markUpCDR3: '", substr ($DNA_Trimmed,0, $MotifPositionDNA)," (DNA until pattern match string)\n"; +# print "D: markUpCDR3: '$DNA_Trimmed' (Trimmed DNA string)\n"; +# print "D: markUpCDR3: '$LaterString' (Later part of DNA string)\n"; + if ($LaterString =~ m/$DNAPat/) + { print "D: markUPCDR3: Also got a match further down the DNA String: at ", $-[0] ," to ", $+ [0], " - which might be worrying\n"; } + } + +=head4 At AA level: "WGxG" [+2] + +=cut + +if ($AA_Trimmed=~ m/$AAPat/) + { + $MotifPositionAA = $+[0]; #Just the easiest way to do this in Perl + $MotifFound = $MotifFound + 2; +# print "D: markUpCDR3: Found Motif match on AA at position (on DNA remember): '$MotifPositionAA' (ie.)\n"; + (my $CDR3_seq) = substr ($AA_Trimmed, 0, $MotifPositionAA); +# print "D: markUpCDR3: Seq ='$CDR3_seq' - as detected\n"; + + } + +=head4 Assess the results of motif position finding + +=cut + +#print "D: markUpCDR3: MotifFound = '$MotifFound'\n"; + +if ($MotifFound ==0) + { return ($Start, $End, $MotifFound); } #The easy one really: return we didn't find the CDR3 + +# +$Start = $FWR3_End; #We assume the end of the FWR3 is the start of CDR3: +#Just found in DNA: +if ($MotifFound ==1) + { + $Start = $FWR3_End; #We assume the end of the FWR3 is the start of CDR3: + $End = $MotifPositionDNA; + $How = "MOTIF_FOUND_IN_DNA"; + } +#Just found in AA: +if ($MotifFound ==2) + { + $End = $MotifPositionAA; + $How = "MOTIF_FOUND_IN_AA"; + } + +#Found in both, DNA has priority: +if ($MotifFound ==3) + { + $Start = $FWR3_End ; #We assume the end of the FWR3 is the start of CDR3: + $End = $MotifPositionDNA; + $How = "MOTIF_FOUND_IN_BOTH"; + } + +#print "D: markUpCDR3: Motif found = $MotifFound\n"; + +=head4 These next few lines are for testing / diagnostics only - disable for general use + +If you are interested in getting the CDR3 directly then remember the main coordinate system is defined such that +the start of FWR1 is unlikely to be at nt 1. + +=cut + +$Start = $FWR3_End - $FWR1_Start -1; +$End = $End + $CropPoint; +my $CDR3_RegionLength = $End - $Start; +#print "D: markUpCDR3: CDR3 Length= $Start - $End = '$CDR3_RegionLength'\n"; +(my $CDR3_seq) = substr ($AA, $Start, $CDR3_RegionLength); + +#Add onto the coordinates what we trimmed off: + + +#print "D: markUpCDR3: Seq ='$CDR3_seq'\n"; + +print "D: markUpCDR3: returning: $Start, $End, $How, ($MotifFound) [NB: offset of :'+ $FWR1_Start'\n"; +#die "HIT BLOCK\n"; +return ($Start + $FWR1_Start, $End + $FWR1_Start, $How); +} + + +sub printOUTPUTData { +=head2 sub: $OutputDataString = printOUTPUTData {\%OutputData} + +When passed an array containing the appropriate CDR, Top V / D/ J genes and the seqeunce ID. +This prepared and then returned as a text string that can then be printed to STDOUT: + + print (printOUTPUTData (\%OutputData)); + +Any missing data in the Hash array it polietly ignored and a null string printed in place. +The text field is tab delimited; there are no extra trailing tabs or carriage returns in place. + +Actually the fields printed out are stored in an index array. + +=head3 Header output + +If the routine is passed a key 'HEADER' then the header columns are returned as that string. +This is tested first - so don't add this unless you mean to. + +=cut + +my @HeaderFields = ("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"); + +my $OutputString = "OUTPUT:"; #What we are going to build the output into. + +=head4 Print Header & Exit? + +=cut + +my ($Data_ref) = @_; +#print "D: printOUTPUTData: Running\n"; + +if (exists $$Data_ref {"HEADER"}) + { + $OutputString .= "\t"; + for(my $n = 0; $n <= $#HeaderFields; $n++) + { + $OutputString .= $HeaderFields[$n]; + $OutputString .= "\t" if($n < $#HeaderFields); + } + + # foreach my $C_Header (@HeaderFields) + # { $OutputString .= "$C_Header"; } # + + print "D: printOUTPUTData: HEADER Printout requested '@HeaderFields'\n"; + return ($OutputString); + } + +=head3 Assemble whatever data we have - and tab delimit the null fields + +=cut +#print "D: printOUTPUTData: Will pretty print this:\n", Dumper $Data_ref; +foreach my $C_Header (@HeaderFields) + { + + if (exists ($$Data_ref {$C_Header})) + { $OutputString .= "\t". $$Data_ref{$C_Header}; } #We have data to print out + else + { $OutputString .="\t"; } #Add a trailing space + } # + +return ($OutputString); +} + + +######################################### Code Junk ######################## + + +=head2 Code Junk Attic + +=head3 Demonstrates how to reverse translate an amino acid sequence into DNA: + +use Bio::Tools::CodonTable; +use Bio::Seq; + +# print possible codon tables + my $tables = Bio::Tools::CodonTable->tables; + while ( (my $id, my $name) = each %{$tables} ) { + print "$id = $name\n"; + } + my $CodonTable = Bio::Tools::CodonTable->new(); + + my $ExampleSeq = Bio::PrimarySeq->new(-seq=>"WGxG", -alphabet => 'protein') or die "Cannot create sequence object\n"; + + +my $rvSeq = $CodonTable->reverse_translate_all($ExampleSeq); +print "D: '$rvSeq'\n"; +die "TEST OVER\n"; + +=cut + + +=head3 For processing the 'Alignment lines' section of the alginment table + + #If we are ever interested; then enable the code below: +# print ": Alignment\n"; +# $InfoPanel =~ s/^ +//; $InfoPanel =~ s/ +$//; #Clean off trailing spaces +# my ($Germclass, $PID, $PID_Counts, $Allele) = split (/\s+/,$InfoPanel); #Split on spaces +##Enable if you need to know what we just found: +# #print "D: Fields are (Germclass, PID, PID_Counts, Allele) \t$Germclass, $PID, $PID_Counts, $Allele\n"; +# #A reality check: we should have an Allele - or some text here. +# unless (defined $Allele && $Allele ne "") +# { warn "Cannot get Allele for Line '$C_Line' - implies improper parsing: '",substr ($Lines[$C_Line],0,15),"...'\n"; } +# if (exists ($Alginments {$Germclass}{$Allele})) +# { $Alginments {$Germclass}{$Allele} = $Alginments {$Germclass}{$Allele}.$CurrentAASequence; } #Carry on adding +# else #more work needed as we need to 'pad' the sequence with fake gap characters) +# { +##Do we still need this padding? I don't think so +# +# +# my $PaddingChars = ($ThisQueryStart-$Query_Start); +# print "D: New gene found: need to pad it with ($ThisQueryStart-$Query_Start) i.e. '$PaddingChars' characters\n"; +# #To help testing, calculate this first: +# my $PaddingString = " "x $PaddingChars; +# $Alginments {$Germclass}{$Allele} = $CurrentAASequence; +# } +# next + +=head3 Demonstration of Pattern match positions + +my $Text = "12345TTT TTAAAAA"; +my $TestPat = "TTT\\s+TT"; +(my $Result)= $Text =~ m/$TestPat/; +print "D: Two vars are: - = ",$-[0], " & + =", $+[0]," for test pattern '$TestPat'\n"; + +sub printCDR3 { + +=head3 Subroutine: printCDR3 ($CDR3_Start, $CDR3_End, "SUMMARY_TABLE", $AAQuerySequence, $DNAQuerySequence); + +???? IS THIS FUNCTION IN USE ????? + +Handles the printing of the output when passed information about the CDR3 region. + + +The result is sent returned as a text string in this version hence use it like this if you want to send it to STDOUT: + + print printCDR3 ($CDR3_Start, $CDR3_End, "SUMMARY_TABLE", $AAQuerySequence, $DNAQuerySequence), "\n"; + +#=cut + +#Despite the similarity in names, these are all local copies passed to us: + +my ($Start, $End, $Tag, $FullAAQuerySequence, $FullDNAQuerySequence) = @_; + +#For DNA: +my ($CDR_DNA_Seq) = substr ($FullDNAQuerySequence, $Start, $Start+$End); +my ($CDR_DNA_Length) = length ($CDR_DNA_Seq); + +#For AA: +my ($CDR_AA_Seq) = substr ($FullAAQuerySequence, $Start, $Start+$End); +my ($CDR_AA_Length) = length ($CDR_AA_Seq); + +my $ReturnString = join ("\t", $CDR_DNA_Seq, $CDR_DNA_Length, $CDR_AA_Seq, $CDR_AA_Length, $Tag); #Create here so we can inspect it / post process it if needed: +print "D: SUB: printCDR3: As returned: '$ReturnString'\n"; +return ($ReturnString); + +} + +=cut + + + +=head2 Change Log + +=head3 Version 1.2 + + 1) Fixed the 'Process recrod request' feature' [was failed increment in $Record] + 2) Deleted / Deactivated the function 'printCDR3' [wasn't in used; kept if useful for parts]. + This function is replaced by the more general printOUTPUTData() + 3) A tag for the CDR3 status is now output for every record / read. + Initially this is set to "NOT_FOUND" and changed if evidence for the CDR3 is found. + +=head4 Version 1.3 + + 1) The tophit line was split on whitespace, however sometimes the VJFrame is something like “In-frame with stop codon”, + which means the line is also split on the spaces therein. It now splits on tabs only, and this seems to work properly. + - found by Bas Horsman. + +=head4 Version 1.3a + + 1) "MOTIF_FOUND_IN_AA" reported correctly (was impossible previously due to addition error to the $MotifFound var (never could == 3) + +=cut + +=head4 Version 1.4 + + 1) Now processes files using Mac/Unix/MS-DOS newline characters: + + $_ =~ s/\r\n/\n/g; #In case line ends are MS-DOS + $_ =~ s/\r/\n/g; #In case line ends are Mac + #The whole record - one per read - is now stored in $_ + my @Lines =split (/\R/,$_); #Split on new lines + +=head4 Version 1.4a + +1) Fixed the length of the CDR3 AA string being reported correctly: + + $OUTPUT_Data{"CDR3 Length"} = $CDR3_Length; + to: + $OUTPUT_Data{"CDR3 Length"} = $CDR3_Seq_AA_Length; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imgtconvert.py Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,139 @@ +import pandas as pd +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("--input", help="Input folder with files") +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'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'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"] + +inputFolder = args.input + +dirContents = os.listdir(inputFolder) +if len(dirContents) == 1: + inputFolder = os.path.join(inputFolder, dirContents[0]) + if os.path.isdir(inputFolder): + print "is dir" + dirContents = os.listdir(inputFolder) +files = sorted([os.path.join(inputFolder, f) for f in dirContents]) + +if len(files) % 3 is not 0: + stop_err("Files in zip not a multiple of 3, it should contain the all the 1_, 5_ and 6_ files for a sample") + import sys + sys.exit() + + +triplets = [] +step = len(files) / 3 +for i in range(0, step): + triplets.append((files[i], files[i + step], files[i + step + step])) + +outFile = args.output + +fSummary = pd.read_csv(triplets[0][0], sep="\t") +fSequence = pd.read_csv(triplets[0][1], sep="\t") +fJunction = pd.read_csv(triplets[0][2], sep="\t") +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 + +for triple in triplets[1:]: + fSummary = pd.read_csv(triple[0], sep="\t") + fSequence = pd.read_csv(triple[1], sep="\t") + fJunction = pd.read_csv(triple[2], sep="\t") + + 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 = outFrame.append(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', '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'] + +vPattern = re.compile(r"IGHV[1-9]-[0-9ab]+-?[1-9]?") +dPattern = re.compile(r"IGHD[1-9]-[0-9ab]+") +jPattern = re.compile(r"IGHJ[1-9]") + +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"] = 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 D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? +outFrame.to_csv(outFile, sep="\t", index=False) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imgtconvert.sh Mon Dec 09 06:08:52 2013 -0500 @@ -0,0 +1,57 @@ +#!/bin/bash +dir="$(cd "$(dirname "$0")" && pwd)" +mkdir $PWD/$2_$3 + + +#!/bin/bash +f=$(file $1) +zip7Type="7-zip archive" +tarType="tar archive" +bzip2Type="bzip2 compressed" +gzipType="gzip compressed" +zipType="Zip archive" +rarType="RAR archive" + +if [[ "$f" == *"$zip7Type"* ]]; then + echo "7-zip" + echo "Trying: 7za e $1 -o$PWD/$2_$3/" + 7za e $1 -o$PWD/$2_$3/ +fi + +if [[ "$f" == *"$tarType"* ]] +then + echo "tar archive" + echo "Trying: tar xvf $1 -C $PWD/$2_$3/" + tar xvf $1 -C $PWD/$2_$3/ +fi + +if [[ "$f" == *"$bzip2Type"* ]] +then + echo "bzip2 compressed data" + echo "Trying: tar jxf $1 -C $PWD/$2_$3/" + tar jxf $1 -C $PWD/$2_$3/ +fi + +if [[ "$f" == *"$gzipType"* ]] +then + echo "gzip compressed data" + echo "Trying: tar xvzf $1 -C $PWD/$2_$3/" + tar xvzf $1 -C $PWD/$2_$3/ +fi + +if [[ "$f" == *"$zipType"* ]] +then + echo "Zip archive" + echo "Trying: unzip $1 -d $PWD/$2_$3/" + unzip $1 -d $PWD/$2_$3/ +fi + +if [[ "$f" == *"$rarType"* ]] +then + echo "RAR archive" + echo "Trying: unrar e $1 $PWD/$2_$3/" + unrar e $1 $PWD/$2_$3/ +fi +find $PWD/$2_$3/ -type f | grep -v "1_Summary_\|5_AA-sequences_\|6_Junction_" | xargs rm -f +python $dir/imgtconvert.py --input $PWD/$2_$3 --output $4 +
--- a/immunerepertoirecombined_imgt/RScript.r Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,295 +0,0 @@ -#options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) - -args <- commandArgs(trailingOnly = TRUE) - -inFile = args[1] -outFile = args[2] -outDir = args[3] -clonalType = args[4] - -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/") -} -require(ggplot2) -if (!("plyr" %in% rownames(installed.packages()))) { - install.packages("plyr", repos="http://cran.xl-mirror.nl/") -} -require(plyr) - -if (!("data.table" %in% rownames(installed.packages()))) { - install.packages("data.table", repos="http://cran.xl-mirror.nl/") -} -library(data.table) - - -test = read.table(inFile, sep="\t", header=TRUE, fill=T) - -test = test[test$Sample != "",] - -test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene) -test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene) -test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene) - -#test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":")) -test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":")) - -PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ] - -NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ] - -#PRODF = PROD[ -1] - -PRODF = PROD - -#PRODF = unique(PRODF) -PRODF = PRODF[!duplicated(PRODF$VDJCDR3), ] - -PRODFV = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Sample", "Top.V.Gene")]) -PRODFV$Length = as.numeric(PRODFV$Length) -Total = 0 -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=.N), by=c("Sample", "Top.D.Gene")]) -PRODFD$Length = as.numeric(PRODFD$Length) -Total = 0 -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=.N), by=c("Sample", "Top.J.Gene")]) -PRODFJ$Length = as.numeric(PRODFJ$Length) -Total = 0 -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)) - -V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53") -tcV = textConnection(V) -Vchain = read.table(tcV, sep="\t", header=TRUE) -PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) -close(tcV) - -D = c("v.name\tchr.orderD\nIGHD1-1\t1\nIGHD2-2\t2\nIGHD3-3\t3\nIGHD6-6\t4\nIGHD1-7\t5\nIGHD2-8\t6\nIGHD3-9\t7\nIGHD3-10\t8\nIGHD4-11\t9\nIGHD5-12\t10\nIGHD6-13\t11\nIGHD1-14\t12\nIGHD2-15\t13\nIGHD3-16\t14\nIGHD4-17\t15\nIGHD5-18\t16\nIGHD6-19\t17\nIGHD1-20\t18\nIGHD2-21\t19\nIGHD3-22\t20\nIGHD4-23\t21\nIGHD5-24\t22\nIGHD6-25\t23\nIGHD1-26\t24\nIGHD7-27\t25") -tcD = textConnection(D) -Dchain = read.table(tcD, sep="\t", header=TRUE) -PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE) -close(tcD) - - -J = c("v.name\tchr.orderJ\nIGHJ1\t1\nIGHJ2\t2\nIGHJ3\t3\nIGHJ4\t4\nIGHJ5\t5\nIGHJ6\t6") -tcJ = textConnection(J) -Jchain = read.table(tcJ, sep="\t", header=TRUE) -PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) -close(tcJ) - -setwd(outDir) - -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") - -png("VPlot.png",width = 1280, height = 720) -pV -dev.off(); - -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") - -png("DPlot.png",width = 800, height = 600) -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") - -png("JPlot.png",width = 800, height = 600) -pJ -dev.off(); - -revVchain = Vchain -revDchain = Dchain -revVchain$chr.orderV = rev(revVchain$chr.orderV) -revDchain$chr.orderD = rev(revDchain$chr.orderD) - -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() -} - -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(test$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() - } - 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() -} - -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(test$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) - -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() -} - -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(test$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) - - -sampleFile <- file("samples.txt") -un = unique(test$Sample) -un = paste(un, sep="\n") -writeLines(un, sampleFile) -close(sampleFile) - - -if("Replicate" %in% colnames(test)) -{ - clonalityFrame = PROD - clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":")) - clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] - clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")]) - 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", "VDJCDR3")]) - ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) - 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(Reads), ReadsSquaredSum=sum(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$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) -}
--- a/immunerepertoirecombined_imgt/combined.sh Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,75 +0,0 @@ -#!/bin/bash - -export IGDATA=/home/david/tmp/ncbi-igblast-1.0.0; - -clonalType=${@:(-3):1} -html=${@:(-2):1} -imageDir=${@:(-1):1} -dataCount=`expr $# - 3` -inputData=${@:(1):dataCount} -dir="$(cd "$(dirname "$0")" && pwd)" -array=("$@") -echo "<html><h3>Progress</h3><table><tr><td>info</td></tr>" > $html -echo "<tr><td>-----------------------------------</td></tr>" >> $html - -id=${inputData[0]} -forwardSlash="/" -mergerInput=() -count=0 -for current in $inputData; do - if [[ "$current" != *"$forwardSlash"* ]]; then - id=$current - count=0 - mergerInput+=($id) - continue - fi - fileName=$(basename $current) - convertedFileName="${fileName%.*}" - convertedFileName="$PWD/$convertedFileName.converted" - bash $dir/imgtconvert.sh $current $id $count $convertedFileName - mergerInput+=($convertedFileName) - count=$((count+1)) -done -wait - - - -echo "<tr><td>-----------------------------------</td></tr>" >> $html -echo "<tr><td>merging</td></tr>" >> $html - -python $dir/igblastmerge.py ${mergerInput[*]} --output $PWD/merged.txt - -echo "<tr><td>done</td></tr>" >> $html -echo "<tr><td>-----------------------------------</td></tr>" >> $html -echo "<tr><td>plotting</td></tr>" >> $html - - -inputFile=$PWD/merged.txt -outputFile=$html -outputDir=$imageDir -mkdir $outputDir -Rscript --verbose $dir/RScript.r $inputFile $outputDir $outputDir $clonalType 2>&1 -echo "<html>" > $outputFile -echo "<img src='VPlot.png'/>" >> $outputFile -echo "<img src='DPlot.png'/>" >> $outputFile -echo "<img src='JPlot.png'/>" >> $outputFile - -samples=`cat $outputDir/samples.txt` -count=1 -echo "<table border='1'>" >> $outputFile -for sample in $samples; do - echo "<tr><td colspan='3' height='100'></td>" >> $outputFile - echo "<tr><td colspan='3'><h1>$sample</h1></td></tr>" >> $outputFile - echo "<tr><td><h2>V-D Heatmap:</h2></td><td><h2>V-J Heatmap:</h2></td><td><h2>D-J Heatmap:</h2></td></tr><tr>" >> $outputFile - mv "$outputDir/HeatmapVD_$sample.png" "$outputDir/VD_$sample.png" - echo "<td><img src='VD_$sample.png'/></td>" >> $outputFile - mv "$outputDir/HeatmapVJ_$sample.png" "$outputDir/VJ_$sample.png" - echo "<td><img src='VJ_$sample.png'/></td>" >> $outputFile - mv "$outputDir/HeatmapDJ_$sample.png" "$outputDir/DJ_$sample.png" - echo "<td><img src='DJ_$sample.png'/></td></tr>" >> $outputFile - count=$((count+1)) -done -echo "</table>" >> $outputFile -echo "<h1>$clonalType</h1>" >> $outputFile - -echo "</html>" >> $outputFile
--- a/immunerepertoirecombined_imgt/combined_imgt.xml Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -<tool id="immunerepertoirecombined_imgt" name="Immunerepertoire combined" version="1.0"> - <description>pipeline for IMGT</description> - <command interpreter="bash"> - combined.sh - #for $i, $f in enumerate($patients) - ${f.id} - #for $j, $g in enumerate($f.samples) - ${g.sample} - #end for - #end for - "$clonaltype_select" $out_file $out_file.files_path - </command> - <inputs> - <repeat name="patients" title="Patients" min="1" default="1"> - <repeat name="samples" title="Sample" min="1" default="1"> - <param name="sample" type="data" label="Sample to Process" /> - </repeat> - <param name="id" type="text" label="ID" /> - </repeat> - <param name="clonaltype_select" type="select" label="Clonal Type Definition"> - <option value="Top.V.Gene,CDR3.Seq">Top.V.Gene, CDR3.Seq</option> - <option value="Top.V.Gene,CDR3.Seq.DNA">Top.V.Gene, CDR3.Seq.DNA</option> - <option value="Top.V.Gene,Top.J.Gene,CDR3.Seq">Top.V.Gene, Top.J.Gene, CDR3.Seq</option> - <option value="Top.V.Gene,Top.J.Gene,CDR3.Seq.DNA">Top.V.Gene, Top.J.Gene, CDR3.Seq.DNA</option> - <option value="Top.V.Gene,Top.D.Gene,Top.J.Gene,CDR3.Seq.DNA">Top.V.Gene, Top.D.Gene, Top.J.Gene, CDR3.Seq.DNA</option> - </param> - </inputs> - <outputs> - <data format="html" name="out_file" /> - </outputs> - <help> - The entire Immune Repertoire pipeline as a single tool, input several FASTA files, give them an ID and it will BLAST, parse, merge and plot them. - </help> - <requirements> - <requirement type="package" version="1.0.0">igBlastn</requirement> - </requirements> -</tool> -
--- a/immunerepertoirecombined_imgt/igblastmerge.py Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,52 +0,0 @@ -import sys -# error -def stop_err( msg ): - sys.stderr.write( "%s\n" % msg ) - sys.exit() - -# main -def main(): - args = sys.argv[1:-2] - print args - try: - o = open(sys.argv[-1], 'w') - i = open(args[1], 'r') - separator = "\t" - newline = "\n" - line = i.readline() - #write the header - o.write(line[:line.rfind(newline)] + separator + "Sample" + separator + "Replicate" + newline) - i.close() - - current = 1 - sampleID = args[0] - count = 1 - - while True: - print str(o) - f = open(args[current], 'r') - line = f.readline() - line = f.readline() - while line: - o.write(line[:line.rfind(newline)] + separator + sampleID + separator + str(count) + newline) - line = f.readline() - f.close() - - if current >= (len(args) - 1): - break - if args[current + 1].find("/") is -1: - sampleID = args[current + 1] - current += 1 - count = 1 - else: - count += 1 - current += 1 - o.close() - - except Exception, ex: - stop_err('Error running new_column.py\n' + str(ex)) - sys.exit(0) - -if __name__ == "__main__": - print sys.argv - main()
--- a/immunerepertoirecombined_imgt/igparse.pl Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1252 +0,0 @@ -#!/usr/bin/perl -=head1 IGBLAST_simple.pl - -This version (1.4) has been heavily adapted since the original program was first created back in October 2012. -Bas Horsman (EMC, Rotterdam, The Netherlands) has contributed with minor - though important - code changes. - -From V 1.2 onwards a 'Change Log' is included at the end of the program - -=head2 Usage - -Requires no modules in general use; the Data::Dumper (supplied as part of the Perl Core module set) might be useful for debugging/adjustment -as it allows inspection of the data stores. - -The program takes a text file of the - - ./IGBLAST_simple.pl igBLASTOutput.txt <-optional: index of record to process-> - -Supply the text version of the igBLAST report in the format as in the example below. -The extra command line arugment is the record number (aka. BLAST report) to process. -If 0 or absent all are processed, if supplied that record (base 1) is processed and the program dies afterwards. - -=head2 Example Input - -A standard igBLAST record or set of them in a file; this being typical: - - BLASTN 2.2.27+ - - -Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A. -Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. -Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of -protein database search programs", Nucleic Acids Res. 25:3389-3402. - - - -Database: human_gl_V; human_gl_D; human_gl_J - 674 sequences; 179,480 total letters - - - -Query= HL67IUI01D26LR length=433 xy=1559_1437 region=1 -run=R_2012_04_10_11_57_56_ - -Length=433 - Score E -Sequences producing significant alignments: (Bits) Value - -lcl|IGHV3-30*04 330 2e-92 -lcl|IGHV3-30-3*01 330 2e-92 -lcl|IGHV3-30*01 327 2e-91 -lcl|IGHD3-16*01 14.4 11 -lcl|IGHD3-16*02 14.4 11 -lcl|IGHD1-14*01 12.4 43 -lcl|IGHJ4*02 78.3 1e-18 -lcl|IGHJ5*02 70.3 4e-16 -lcl|IGHJ4*01 68.3 2e-15 - - -Domain classification requested: imgt - - -V(D)J rearrangement summary for query sequence (Top V gene match, Top D gene match, Top J gene match, Chain type, V-J Frame, Strand): -IGHV3-30*04 IGHD3-16*01 IGHJ4*02 VH In-frame + - -V(D)J junction details (V end, V-D junction, D region, D-J junction, J start). Note that possible overlapping nucleotides at VDJ junction (i.e, nucleotides that could be assigned to either joining gene segment) are indicated in parentheses (i.e., (TACT)) but are not included under V, D, or J gene itself -AGAGA TATGAGCCCCATCATGACA ACGTTTG CCGGAA ACTAC - -Alignment summary between query and top germline V gene hit (from, to, length, matches, mismatches, gaps, percent identity) -FWR1 27 38 12 11 1 0 91.7 -CDR1 39 62 24 22 2 0 91.7 -FWR2 63 113 51 50 1 0 98 -CDR2 114 137 24 23 1 0 95.8 -FWR3 138 251 114 109 5 0 95.6 -CDR3 (V region only) 252 259 8 7 1 0 87.5 -Total N/A N/A 233 222 11 0 95.3 - - -Alignments - - <----FWR1--><----------CDR1--------><-----------------------FWR2------ - W A A S G F T F N T Y A V H W V R Q A P G K G - Query_1 27 TGGGCAGCCTCTGGATTCACCTTCAATACCTATGCTGTGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGC 96 -V 95.3% (222/233) IGHV3-30*04 64 ..T......................G..G.......A................................. 133 - C A A S G F T F S S Y A M H W V R Q A P G K G -V 95.7% (221/231) IGHV3-30-3*01 64 ..T......................G..G.......A................................. 133 -V 94.8% (221/233) IGHV3-30*01 64 ..T......................G..G.......A................................. 133 - - ----------------><----------CDR2--------><---------------------------- - L E W V A V I S Y D G S N K N Y A D S V K G R F - Query_1 97 TGGAGTGGGTGGCAGTTATATCATATGATGGAAGCAATAAAAACTACGCAGACTCCGTGAAGGGCCGATT 166 -V 95.3% (222/233) IGHV3-30*04 134 ..................................T......T............................ 203 - L E W V A V I S Y D G S N K Y Y A D S V K G R F -V 95.7% (221/231) IGHV3-30-3*01 134 .........................................T............................ 203 -V 94.8% (221/233) IGHV3-30*01 134 .A................................T......T............................ 203 - - ---------------------------FWR3--------------------------------------- - T I S R D N S K N T L Y L Q M N S L R V E D T - Query_1 167 CACCATCTCCAGAGACAATTCCAAGAACACGTTATATCTGCAAATGAACAGCCTGAGAGTTGAGGACACG 236 -V 95.3% (222/233) IGHV3-30*04 204 ...............................C.G.........................C.......... 273 - T I S R D N S K N T L Y L Q M N S L R A E D T -V 95.7% (221/231) IGHV3-30-3*01 204 ...............................C.G.........................C.......... 273 -V 94.8% (221/233) IGHV3-30*01 204 ...............................C.G.........................C.......... 273 - - --------------> - A V Y Y C T R D M S P I M T T F A G N Y W G Q - Query_1 237 GCTGTTTATTACTGTACGAGAGATATGAGCCCCATCATGACAACGTTTGCCGGAAACTACTGGGGCCAGG 306 -V 95.3% (222/233) IGHV3-30*04 274 .....G.........G.......----------------------------------------------- 296 - A V Y Y C A R -V 95.7% (221/231) IGHV3-30-3*01 274 .....G.........G.....------------------------------------------------- 294 -V 94.8% (221/233) IGHV3-30*01 274 .....G.........G.......----------------------------------------------- 296 -D 100.0% (7/7) IGHD3-16*01 12 ------------------------------------------.......--------------------- 18 -D 100.0% (7/7) IGHD3-16*02 12 ------------------------------------------.......--------------------- 18 -D 100.0% (6/6) IGHD1-14*01 8 -------------------------------------------------......--------------- 13 -J 100.0% (39/39) IGHJ4*02 10 -------------------------------------------------------............... 24 -J 100.0% (35/35) IGHJ5*02 17 -----------------------------------------------------------........... 27 -J 97.4% (38/39) IGHJ4*01 10 -------------------------------------------------------.............A. 24 - - - G T L V T V S S - Query_1 307 GAACCCTGGTCACCGTCTCCTCAG 330 -J 100.0% (39/39) IGHJ4*02 25 ........................ 48 -J 100.0% (35/35) IGHJ5*02 28 ........................ 51 -J 97.4% (38/39) IGHJ4*01 25 ........................ 48 - - -Lambda K H - 1.10 0.333 0.549 - -Gapped -Lambda K H - 1.08 0.280 0.540 - -Effective search space used: 64847385 - - -Query= HL67IUI01EQMLY length=609 xy=1826_1636 region=1 -run=R_2012_04_10_11_57_56_ - - -...etc... - -=head2 Example Output - - -Example output from the data above sent: - $ ./IGBLAST_simple.pl igBLASTOutput.txt 1 - D: Request to process just record '1' received - D: printOUTPUTData: Running - D: printOUTPUTData: HEADER Printout requested '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 Found How' - OUTPUT: # 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 Found How - D: ID is: 'HL67IUI01D26LR' - D: Minimum base marked-up (27) - aka. $AlignmentStart; maximum: (259) - D: Starting Search for CDR3 - D: markUpCDR3: Passed Parameters '251, 27, TGGGG....GG., WG.G' (& AA & DNA sequence) - D: markUpCDR3: returning: 223, 282, MOTIF_FOUND_IN_BOTH, (3) [NB: offset of :'+ 27' - D: CDR3 was found by pattern matching: 'MOTIF_FOUND_IN_BOTH' (250, 309) - D: Top Hits (raw)= 'IGHV3-30*04 IGHD3-16*01 IGHJ4*02 VH In-frame +' - D: Top Hits (parsed)= 'IGHV3-30*04, IGHD3-16*01, IGHJ4*02, VH, In-frame, +' - D: printOUTPUTData: Running - OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH - -=head4 Usage notes: - -Designed to be easy to "grep -v D:" or "grep OUTPUT:" for to select the parts you need: - - ./IGBLAST_simple.pl igBLASTOutput.txt 1 | grep OUTPUT: - - OUTPUT: # 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 Found How - OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01EQMLY In-frame IGHV4-39*01 IGHD2-8*01 IGHJ3*02 GGSISSSSYY 29 IYHSGST 20 CARDATYYSNGFDIWGQG 53 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01CDCLP Out-of-frame IGHV3-23*01 IGHD3-3*01 IGHJ4*02 FSNYAM 16 SGSGDRTY 23 AKAD*FLEWLFRIGDGERLLGPGN 72 MOTIF_FOUND_IN_DNA - OUTPUT: HL67IUI01AHRNH N/A IGHV3-33*01 N/A N/A WIHLQ*LW 23 YGMMEVI 23 NOT_FOUND - OUTPUT: HL67IUI01DZZ1V Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 GFTFDKYA 23 ILASG 20 LYCASEGDIVASELLSTGARV 62 MOTIF_FOUND_IN_DNA - OUTPUT: HL67IUI01DTR2Y Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 LDSPLTNM 23 LYLPVV 20 TVRVRGT*WLRSF*VLGPG 59 MOTIF_FOUND_IN_DNA - OUTPUT: HL67IUI01EQL3S In-frame IGHV7-4-1*02 IGHD6-19*01 IGHJ6*02 GYTFRTFT 23 INTNTGTP 23 CAKESGTGSAHFFYGMDVWGQG 65 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01AFG46 In-frame IGLV2-34*01 N/A IGHJ4*02 NOT_FOUND - OUTPUT: HL67IUI01EFFKO In-frame IGHV3-11*01 IGHD6-6*01 IGHJ4*02 GFTFSDYY 23 ISYSGGTI 23 CARASGAARHRPLDYWGQG 56 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01B18SG In-frame IGHV3-33*01 IGHD5-12*01 IGHJ4*02 VRQA 11 KYYANSVK 23 RLGGFDYWGQGTLVTVSS 53 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01D6LER In-frame IGHV1-24*01 IGHD3-22*01 IGHJ4*02 GYSLNELS 23 PDPEDDE 23 TVQPSRITMMAVVITRIHWGASGARE 76 MOTIF_FOUND_IN_DNA - OUTPUT: HL67IUI01CYCLF N/A IGHV4-39*01 N/A N/A GGSISSSSYY 29 IYYSGST 20 NOT_FOUND - OUTPUT: HL67IUI01B4LEE In-frame IGHV7-4-1*02 IGHD6-19*01 IGHJ6*02 GYTFRTFT 23 INTNTGTP 23 CAKESGTGSAHFFYGMDVWGQG 65 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01A4KW4 Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 LDSPLTNM 23 LYLPVV 20 TVRVRGT*WLRSF*IWGQG 58 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01E05BV In-frame IGHV1-24*01 IGHD3-22*01 IGHJ2*01 GYSLNELS 23 PDPEDDE 23 NOT_FOUND - OUTPUT: HL67IUI01CVVKY In-frame IGHV1-3*01 IGHD2-15*01 IGHJ1*01 NOT_FOUND - OUTPUT: HL67IUI01CN5P2 In-frame IGHV7-4-1*02 IGHD2-21*02 IGHJ5*02 GYSITDYG 23 LNTRTGNP 23 CAVKDARDFVSWGQG 44 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01DUUJ5 In-frame IGHV3-21*01 IGHD1-7*01 IGHJ4*02 GYTFSTYS 23 ISSSSAYR 23 CARDIRLELRDWGQG 44 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01E1AIR Out-of-frame IGHV4-39*01 N/A IGHJ3*01 WGLHRRW**L 29 FVS*RAPR 23 NOT_FOUND - OUTPUT: HL67IUI01CCZ8D Out-of-frame IGHV3-23*01 IGHD5-12*01 IGHJ4*02 GFTFDKYA 23 ILASGR 20 YCASEGDIVASELLSTGARE 58 MOTIF_FOUND_IN_DNA - OUTPUT: HL67IUI01BT9IR N/A IGHV3-21*02 N/A N/A NOT_FOUND - OUTPUT: HL67IUI01COTO0 Out-of-frame IGHV4-39*01 N/A IGHJ3*01 GGFIGGGDNF 29 LYHDGRPA 23 NOT_FOUND - OUTPUT: HL67IUI01D994O In-frame IGHV7-4-1*02 IGHD2-21*02 IGHJ5*02 GYSITDYG 23 LNTRTGNP 23 CAVKDARDFVSWGQG 44 MOTIF_FOUND_IN_BOTH - OUTPUT: HL67IUI01A08CJ In-frame IGHV4-39*01 IGHD6-13*01 IGHJ5*02 GGSISSSSYY 29 IYYTWEH 21 CERARRGSSWGQLVRPLGPG 62 MOTIF_FOUN - - - - OUTPUT: # 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 Found How - OUTPUT: HL67IUI01D26LR In-frame IGHV3-30*04 IGHD3-16*01 IGHJ4*02 GFTFNTYA 23 ISYDGSNK 23 CTRDMSPIMTTFAGNYWGQG 59 MOTIF_FOUND_IN_BOTH - ...etc... - -=head4 Also, combined grep & sed: - - $ ./IGBLAST_simple.pl igBLASTOutput.txt | grep OUTPUT: | sed 's/OUTPUT:\t//' - -=cut - -=head3 CDR3 Patterns: - -We use these two variables to try to identify the end of the CDR3 region if igBLAST doesn't report it directly: - - my $DNACDR3_Pat = "TGGGG....GG."; - my $AASequenceMotifPattern = "WG.G"; - -They are treated as regex's when tested (so use "." to mean any DNA base, rather than 'N' or 'X'). - -[NB: These are original patterns used for testing, check the code for the current ones.] - -=cut - -my $DNACDR3_Pat = "TGGGG....GG."; -my $AACDR3_Pat = "WG.G"; - -use strict; -use Data::Dumper; -# Set this as to number of the result (aka "record") you want to process or 0 for all: -my $ProcessRecord =0; -if (defined $ARGV[1]) { $ProcessRecord = pop @ARGV; } #Also accept from the command line: -if ($ProcessRecord != 0) { print "D: Request to process just record '$ProcessRecord' received\n"; } - -#Adjust the record separator: -$/="Query= "; -my $Record=0; # A simple counter, that we might not use. -#Force-loaded header / version information: -my $Header = <>; -#At the moment we don't use this - so dump it immediately: -$Header = undef; -#print "D: Force-loaded header / version information: '$Header'\n"; - -#Print the Header for the output line (we need this once, at the start) -print &printOUTPUTData ({"HEADER" => 1})."\n"; - -while (<>) - { -=head4 First check - should we be processing this record at all? - -=cut - $Record++; #Increment the record counter: - #Do we process this record - or all records? - if ($ProcessRecord != $Record && $ProcessRecord != 0) - { next; } #We need to increment the record counter before we increment - -=head4 Setup the output line storage and print the header: - -We enter this initially and work to change it: - - $DomainBoundaries{"CDR3"}{"FoundHow"} = "NOT_FOUND"; - -=cut - - my %OUTPUT_Data; #To collect data for the output line in - #Assume the first and work to find better: - $OUTPUT_Data{"CDR3 Found How"} = "NOT_FOUND"; - #The whole record - one per read - is now stored in $_ - my @Lines =split (/[\r\n]+/,$_); # split on windows/linux/mac new lines - - #If you are interested enable either of the next lines depending on how curious you are as to how the splitting went: - #print "D: Record #$Record\n"; print $_; print "\n---------\n"; - print "D: ''$Lines[0]'\nD: ...etc...'\nD: ############\n"; - -=head3 Get the ID - -Quite easy: the first field on the first line: - - Query= HL67IUI01DTR2Y length=577 xy=1452_0984 region=1 - -=cut - - (my $ID) = $Lines[0]=~ m/^(\S+)/; - unless (defined $ID && $ID ne "") - { # So a near total failure...? - $OUTPUT_Data{"ID"} = "Unknown"; - print &printOUTPUTData (\%OUTPUT_Data)."\n"; - next; #No ID is terminal for this record - } - else - { - print "D: ID is: '$ID'\n"; - $OUTPUT_Data{"ID"} = $ID; - } -=head3 Declare the variables we will need here in the next few sections to store data - -=cut - - my $CurrentRegion; - my $RegionMarkup; - - #So we can sync the coordinated of the alignment up to the domains found: - my $Query_Start = -1; my $Query_End = -1; - - #Where on the Query Sequence (i.e. the 454 read) does the alignment start & stop? - my $ThisQueryStart =-1; my $ThisQueryEnd =-1; #Think $ThisQueryEnd isn't used at the moment. - my $DNAQuerySequence =""; #The actual DNA Query sequence... - my $AAQuerySequence = ""; - - #As this changes with the alleles identified: - my $CurrentAASequence; - #The main storage variables - - my %Alginments; my %Alleles; - my %DomainBoundaries; - -=head2 Stanza 1: Get the general structure of the sequence identified - -=head3 Method 1: Use the table supplied - -Technically this valid for the top hit...realistically this is the only information we have reported to us -so we use this or nothing. This is fine for the top hit which is likely what we are interested in....but for the 2nd or 3rd? Who knows! - -Targets this block: - - Alignment summary between query and top germline V gene hit (from, to, length, matches, mismatches, gaps, percent identity) - FWR1 167 240 75 72 2 1 96 - CDR1 241 264 24 20 4 0 83.3 - FWR2 265 315 51 48 3 0 94.1 - CDR2 316 336 24 15 6 3 62.5 - FWR3 337 450 114 106 8 0 93 - CDR3 (V region only) 451 454 4 4 0 0 100 - Total N/A N/A 292 265 23 4 90.8 - -Then we split out the lines inside it in a second scanning step - less optimal but easier to read: - - FWR1 167 240 75 72 2 1 96 - CDR1 241 264 24 20 4 0 83.3 - FWR2 265 315 51 48 3 0 94.1 - CDR2 316 336 24 15 6 3 62.5 - FWR3 337 450 114 106 8 0 93 - CDR3 (V region only) 451 454 4 4 0 0 100 - -into: - - (Section, from, to, length, matches, mismatches, gaps, percent identity) - -=head3 Method 2: Use the table supplied - -The other way to do this is to split the graphical markup out of the alignment. -This works for _any_ reported alignment, not just the top hits: - -In the main alignment table processing section collect the information, collect the information: - - #Is region mark-up: - if ($#InfoColumns == -1 && $#AlignmentColumns ==0) - { -# print ": Region Markup detected\n"; - $RegionMarkup = $RegionMarkup.$AlignmentPanel; #Collect the information, then re-synthesise it at the end of record - next; - } - -Then afterwards when all the region was collected, process it like this: -#Pad the CDER3 region: - - #Remove the trailing spaces: - $RegionMarkup =~ s/ *$//g; - #Calculate the length of the CDR3 region so we can add it in: - my $CDR3PaddingNeeded = ($Query_End-$Query_Start)-length ($RegionMarkup) -length ("<-CDR3>")+1; - #Build up the CDR3 region, the 'x' operator is very helpful here (implict foreach loop): - $RegionMarkup = $RegionMarkup."<-CDR3"."-" x $CDR3PaddingNeeded. ">"; - #print "D: Need to pad with:'$CDR3PaddingNeeded' characters\n"; - - #Now really process it: - my $C_Pos = 0; - my @Domains = split (/(<*-*...[123]-*>*)/,$RegionMarkup); # - foreach my $C_Domain (@Domains) - { - if (length ($C_Domain) <=0) {next;} - my $DomainStart= $C_Pos; - my $DomainEnd = $DomainStart + length ($C_Domain)-1; - my ($DomainType) = $C_Domain =~ m/(...[123])/; -# print "D: $DomainType \t($DomainStart-$DomainEnd=",$DomainEnd-$DomainStart,"):\t$C_Domain\n"; - $DomainBoundaries{$DomainType}{"Start"} = $DomainStart; - $DomainBoundaries{$DomainType}{"End"} = $DomainEnd; - $C_Pos = $DomainEnd+1; - } - -The two pieces of code are interchangable; the table version as used below, is neater, easier to understand and works nicely. -Why stress? - - -=head3 The end of the FWR3 is the start of CDR3? - -This is an assumption made. Hence the two variables: - - my $MaxDomainReported =0 ; # In nts / bps - my $FWR3_Found_Flag = 0; # Did we find the end of the FWR3 - which is the start of the CDR3. Set to 'false' initially. - - $MaxDomainBaseFound - -=cut - my $MaxDomainBaseFound =0 ; # In nts / bps - my $AlignmentStart ; # In nts /bp #Alternative name would be: '$MinDomainBaseFound'; set to null until primed -# my $FWR3_Found_Flag = 0; # Did we find the end of the FWR3 - which is the start of the CDR3. Set to 'false' initially. - - (my @StructureSummaryTable) = returnLinesBetween (\@Lines, "Alignment summary", "Total" ); -#Enable the next line if you want the raw data we are going to parse in this section: - #print Dumper @StructureSummaryTable; - foreach my $C_Section (@StructureSummaryTable) - { - my ($DomainType, $DomainStart, $DomainEnd, $SectLength, $Matches, $Mismatches, $Gaps, $PID) = split (/\t+/,$C_Section); - #print "D: Domain type: '$DomainType'\n"; - #$DomainType =~ s/ .*$//g; - $DomainBoundaries{$DomainType}{"Start"} = $DomainStart; - $DomainBoundaries{$DomainType}{"End"} = $DomainEnd; - -#So we can do a reality check on the length / start of the CDR3 if we have to go looking: - if ($MaxDomainBaseFound <= $DomainEnd) - { $MaxDomainBaseFound = $DomainEnd; } #Store the maximum base found - if ($AlignmentStart eq undef or $AlignmentStart >= $DomainStart) - { $AlignmentStart = $DomainStart; } - } -#print Dumper %DomainBoundaries; -#die "HIT BLOCK\n"; - -=head3 Did we find the CDR3 region specifically? - -If we did fine; otherwise try to find it using the FWR3 region if we found that; otherwise give up. - -=cut - print "D: Minimum base marked-up ($AlignmentStart) - aka. \$AlignmentStart; maximum: ($MaxDomainBaseFound)\n"; - -#my @WantedSections = qw (V D J); - -=head2 Second Stanza: Parse the main Alignment Table - -=head3 Get the table, then determine the character at which to split the 'Info' & 'Alignment' panels. - -As this is a little involved and comparamentalises nicely we sub-contract this to two functions"" - - (my @Table) = returnLinesBetween (\@Lines, "Alignment", "Lambda" ); - my $PanelSplitPoint = findSplitPoint (\@Table); #Why can't they just use a fixed field width or a tab as a delimiter? - -=cut - (my @Table) = returnLinesBetween (\@Lines, "Alignment", "Lambda" ); - my $PanelSplitPoint = findSplitPoint (\@Table); #Why can't they just use a fixed field width or a tab as a delimiter? -#If you are interested, enable this line: -# print "D: The info panel was detected at: '$splitPoint'\n"; - -=head3 - -=cut - - -foreach my $C_Line (0..$#Table) - { - -=head3 Call the line type we find: There are 4: - -These are distinguished by the number of fields (one or mores spacer is a field separator) in the Info & Alignment Panels (see values in brackets) - - | <- This split is ~40 chars. from the start of the line - * InfoPanel * | * Alignment Panel * - : is a "Blank" line (-1,-1) - <----FWR1--><----------CDR1--------><-----------------------FWR2------ : is "Region Markup" (-1,0) - W A A S G F T F N T Y A V H W V R Q A P G K G : is "AA Sequence" (-1, >=0) - Query_1 27 TGGGCAGCCTCTGGATTCACCTTCAATACCTATGCTGTGCACTGGGTCCGCCAGGCTCCAGGCAAGGGGC 96 : is "DNA Sequence" (2,1) - V 95.3% (222/233) IGHV3-30*04 64 ..T......................G..G.......A................................. 133 : is "" " - -So we split 40 chars in and then the two parts on spaces. - - -=cut - -# print "D: (sub) Line in parsed table: '$C_Line': \n"; - - my ($InfoPanel, $AlignmentPanel) = $Table[$C_Line] =~ /^(.{$PanelSplitPoint})(.*)$/; - - my @InfoColumns = split (/\s+/,$InfoPanel); - my @AlignmentColumns = split (/\s+/,$AlignmentPanel); - -#If you want to see how the line is being split enable either of these next two lines; the 2nd is more detailed than the first -# print "D: Line: $C_Line/t Number of Columns (Info, Alignment): \t$#InfoColumns \t $#AlignmentColumns\n"; -# print "D: For '$C_Line' \t line in the table there are parts: '$InfoPanel' [$#InfoColumns], '$AlignmentPanel [$#AlignmentColumns]'\n"; - -#Populate this so we can step through it - -=head4 Is a blank line: -=cut - if ($#InfoColumns == -1 && $#AlignmentColumns == -1) - { -# print ": Blank\n"; - next; - } #For now I think we just skip - is not needed (though might be implict mark-up) - -=head4 Is region mark-up: -=cut - if ($#InfoColumns == -1 && $#AlignmentColumns ==0) - { -# print ": Region Markup detected\n"; - $RegionMarkup = $RegionMarkup.$AlignmentPanel; #Collect the information, then re-synthesise it at the end of record - next; - } -=head4 Is query DNA Sequence: -=cut - if ($#InfoColumns == 2 && $#AlignmentColumns ==1) - { -# print ": DNA Query Sequence\n"; - #Detect the two coordinatates of alignment against the query sequence: (last two numbers of the two 'panels') - ($ThisQueryStart) = $InfoPanel =~ / (\d+) *$/; - ($ThisQueryEnd) = $AlignmentPanel =~ / (\d+) *$/; - my ($ThisDNASeq) = $AlignmentPanel =~ /^(.*?) /; - #If you want to know what we just found: - #print "D: This DNA Sequence: '$ThisDNASeq'\n"; - $DNAQuerySequence = $DNAQuerySequence. $ThisDNASeq; #Add it on to whatever we already have. - #Move the needle if there are smaller / greater; otherwise prime the 'needles': - if ($ThisQueryStart < $Query_Start or $Query_Start == -1) - { $Query_Start = $ThisQueryStart; } - if ($ThisQueryEnd > $Query_End or $Query_End == -1) - { $Query_End = $ThisQueryEnd; } -# print ": Query DNA Sequence detected This line: ($ThisQueryStart, $ThisQueryEnd) & Maximally: ($Query_Start, $Query_End)\n"; - next; - } -=head4 Is AA Sequence: - -This is complicated as it Need to decide whether this is the sequence of the read or that of the original V / D / J regions: - --------------> - A V Y Y C T R D M S P I M T T F A G N Y W G Q << Want this - Query_1 237 GCTGTTTATTACTGTACGAGAGATATGAGCCCCATCATGACAACGTTTGCCGGAAACTACTGGGGCCAGG 306 - V 95.3% (222/233) IGHV3-30*04 274 .....G.........G.......----------------------------------------------- 296 - A V Y Y C A R - V 95.7% (221/231) IGHV3-30-3*01 274 .....G.........G.....------------------------------------------------- 294 - - ...etc... - G T L V T V S S << Want this - Query_1 307 GAACCCTGGTCACCGTCTCCTCAG 330 - -To solve this we peak at the next line that it has the tag "Query" in it (we assume the line exists...) - -=cut - - if ($#InfoColumns == -1 && $#AlignmentColumns >=-1) - { - unless ($Table[$C_Line+1] =~ /Query/) { next; } #Is the next line the DNA sequence ? - # -# print ": AA sequence\n"; - - - $CurrentAASequence = $AlignmentPanel; - #print "D: Panel Split Point = $PanelSplitPoint, '$AlignmentPanel'\n"; - $CurrentAASequence =~ s/^ {$PanelSplitPoint}//; - #print "D: '$AAQuerySequence'\n"; -# print "D: Current AA Sequence: \t'$CurrentAASequence'\n"; - $AAQuerySequence = $AAQuerySequence.$CurrentAASequence; #Store the elongating AA Sequence as well - next; - } -=head4 Is Alignment: -=cut - if ($#InfoColumns == 4 && $#AlignmentColumns ==1) - { - #Not acutally interesting to us for this version of the parser. Delete ultimately? - next; - } - -#Is weird! Don't recognise it! - - warn "Weird! Don't recongnise this: '$ID' [$#InfoColumns,$#AlignmentColumns]// '",$Lines[$C_Line],15,"...'\n"; - } #End main iteration loop for alignment parsing. - - -=head2 The CDR3 is noted as problematic. Can we identify it? - -=cut - print "D: Starting Search for CDR3\n"; - #Do have the end of the FWR3 but not the CDR3? If so then it is worth trying to find the CDR3, otherwise...nothing we can do at this point - if (exists ($DomainBoundaries{"FWR3"}{"End"}) - && $AlignmentStart !=0 - && not (exists $DomainBoundaries{"CDR3"}{"End"}) ) #Guess we need to go looking for the end then... - { - #print "D: Placing call to markUpCDR3\n"; - my ($CDR3_Start, my $CDR3_End, my $CDR3_Found_Tag) = markUpCDR3 ($DNAQuerySequence, $AAQuerySequence, - $DomainBoundaries{"FWR3"}{"End"}, $AlignmentStart, - $DNACDR3_Pat, $AACDR3_Pat); - if ($CDR3_Start !=0 && $CDR3_End !=0) - { - $DomainBoundaries{"CDR3"}{"Start"} = $CDR3_Start; - $DomainBoundaries{"CDR3"}{"End"} = $CDR3_End ; - $DomainBoundaries{"CDR3"}{"FoundHow"} = $CDR3_Found_Tag; - print "D: CDR3 was found by pattern matching: '$CDR3_Found_Tag' ($CDR3_Start, $CDR3_End)\n"; - } - else - { print "D: CDR3 was not found [either by igBLAST or by pattern matching]\n"; - $DomainBoundaries{"CDR3"}{"FoundHow"} = "NOT_FOUND"; - } - } - else - { #Was reported by igBLAST - print "D: Found the FWR3 from the Domain Boundary Table\n"; - $DomainBoundaries{"CDR3"}{"FoundHow"} = "IGBLAST_NATIVE"; - } - -#print Dumper %DomainBoundaries; - -=head2 Get the top VDJ regions: - -=cut - -=head2 Extract General Features: - -=cut - (my $TopHit) = $_ =~ m/V-J Frame, Strand\):\n(.*?)\n/s; - print "D: Top Hits (raw)= '$TopHit' \n"; - my ($Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match, $Chain, $VJFrame, $Strand) = split (/\t/,$TopHit); - print "D: Top Hits (parsed)= '$Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match, $Chain, $VJFrame, $Strand'\n"; - -=head2 Store the V / D / J Genes used - -=cut - - if (defined $Top_V_gene_match && $Top_V_gene_match ne "") - { $OUTPUT_Data{"Top V Gene"} = $Top_V_gene_match; } - - if (defined $Top_D_gene_match && $Top_D_gene_match ne "") - { $OUTPUT_Data{"Top D Gene"} = $Top_D_gene_match; } - - if (defined $Top_J_gene_match && $Top_J_gene_match ne "") - { $OUTPUT_Data{"Top J Gene"} = $Top_J_gene_match; } - - if (defined $Strand && $Strand ne "") - { $OUTPUT_Data{"Strand"} = $Strand;} - -=head4 Preamble: ID, Frame, and V / D / J used: - -=cut - #Do a reality check: if we didn't get an ID, then skip: - unless (defined (defined $ID) && $ID ne "" && - defined $VJFrame && $VJFrame ne "") - { - print &printOUTPUTData (\%OUTPUT_Data)."\n"; - next; - } - -#Ok, so we have data...most likely: - #print "OUTPUT:\t",join ("\t", $ID, $VJFrame, $Top_V_gene_match, $Top_D_gene_match, $Top_J_gene_match); - - if (defined $VJFrame && defined $ID && $VJFrame ne "" && $ID ne "") - { $OUTPUT_Data{"VDJ Frame"} = $VJFrame;} - else - { - print &printOUTPUTData (\%OUTPUT_Data)."\n"; - next; - }#REALLY? We didn't find anything? Oh well, move to next record - -=head4 CDR1 - -=cut - #Remember that the alignment starts at the FWR1 start, not nt =0 on the read, hence we substract this off all future AA (& DNA coordinates) - - my $AlignmentOffset = $DomainBoundaries{"FWR1"}{"Start"}; - -# print "D: AA Seqeunce is: '$AAQuerySequence'\n"; - if (exists $DomainBoundaries{"CDR1"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start - { -# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... - my $CDR1Start = $DomainBoundaries{"CDR1"}{"Start"}; - my $CDR1End = $DomainBoundaries{"CDR1"}{"End"}; - my $CDR1_Length = $CDR1End - $CDR1Start; -# print "D: CDR1 $CDR1Start $CDR1End = $CDR1_Length\n"; - #Remember that the alignment starts at the FWR1 start, not nt =0 on the read - my $CDR1_Seq_AA = substr ($AAQuerySequence, $CDR1Start - $AlignmentOffset, $CDR1_Length); -# print "D: '$CDR1_Seq_AA'\n"; - $CDR1_Seq_AA =~ s/ //g; - my $CDR1_Seq_AA_Length = length ($CDR1_Seq_AA); - #Add this data to the output store specifically: - $OUTPUT_Data{"CDR1 Seq"} = $CDR1_Seq_AA; - $OUTPUT_Data{"CDR1 Length"} = $CDR1_Length; - } - #What happens if there is no CDR1 found? Leave blank - the output routine can handle this - -=head4 CDR2 - -=cut - - if (exists $DomainBoundaries{"CDR2"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start - { -# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... - my $CDR2Start = $DomainBoundaries{"CDR2"}{"Start"}; - my $CDR2End = $DomainBoundaries{"CDR2"}{"End"}; - my $CDR2_Length = $CDR2End - $CDR2Start; - my $CDR2_Seq_AA = substr ($AAQuerySequence, $CDR2Start - $AlignmentOffset , $CDR2_Length); - $CDR2_Seq_AA =~ s/ //g; - my $CDR2_Seq_AA_Length = length ($CDR2_Seq_AA); - #Add this data to the output store specifically: - $OUTPUT_Data{"CDR2 Seq"} = $CDR2_Seq_AA; - $OUTPUT_Data{"CDR2 Length"} = $CDR2_Length; - } - #What happens if there is no CDR2 found? Leave blank - the output routine can handle this. - -=head4 CDR3 - -=cut - if (exists $DomainBoundaries{"CDR3"}{"Start"}) #It is very possible that it doesn't; assume the End does though if we find the Start - { -# my $VRegion = $Alginments{"V"}{$C_VRegion}; #Convenience.... - my $CDR3Start = $DomainBoundaries{"CDR3"}{"Start"}; - my $CDR3End = $DomainBoundaries{"CDR3"}{"End"}; - my $CDR3_Length = $CDR3End - $CDR3Start; # This variable isn't used - delete it when safe to do so - my $CDR3_Seq_AA = substr ($AAQuerySequence, $CDR3Start - $AlignmentOffset, $CDR3_Length); - my $CDR3_Seq_DNA = substr ($DNAQuerySequence, $CDR3Start - $AlignmentOffset, $CDR3_Length); - $CDR3_Seq_AA =~ s/ //g; - $CDR3_Seq_DNA =~ s/ //g; - my $CDR3_Seq_AA_Length = length ($CDR3_Seq_AA); - my $CDR3_Seq_DNA_Length = length ($CDR3_Seq_DNA); - #Add this data to the output store specifically: - $OUTPUT_Data{"CDR3 Seq"} = $CDR3_Seq_AA; - $OUTPUT_Data{"CDR3 Length"} = $CDR3_Seq_AA_Length; - $OUTPUT_Data{"CDR3 Seq DNA"} = $CDR3_Seq_DNA; - $OUTPUT_Data{"CDR3 Length DNA"} = $CDR3_Seq_DNA_Length; - #And in the case of the CDR3 how we found it: - $OUTPUT_Data{"CDR3 Found How"} = $DomainBoundaries{"CDR3"}{"FoundHow"}; - } - #What happens if there is no CDR3 found? Leave blank - the output routine can handle this. -#die "HIT BLOCK\n"; -#End of the record; output the data we have collected and move on. -print &printOUTPUTData (\%OUTPUT_Data)."\n"; -} - - - -############ -sub returnLinesBetween { -=head3 SUB: returnLinesBetween ({reference to array Index array}, {regex for top of section}, {regex for bottom of section}) - -When passed a reference to an array and two strings - interpreted as REGEX's - will return the lines of the Array -that are bounded by these tags. - -If either of the tags are not found - or are found in the wrong order - then a null list is returned. - -=cut - -my ($Text_ref, $TopTag, $BotTag) = @_; - -my @Table; -#The two boundary conditions at which we will cut the table: -#print "D: [returnLinesBetween]: '$TopTag, $BotTag'\n"; -#How we record these: -my $AlignmentLine_Top=0; my $AlignmentLine_Bot=0; - -my $LineIndex=-1; #-1 As the loop increments this line counter first, then does its checks. -#If you care: -#print "D: Lines of text passed: $$#Lines\n"; - -#Iterate through until we find what we are looking for or run out of text to search: -while (($AlignmentLine_Bot ==0 or $AlignmentLine_Top==0) && $LineIndex <=$#{$Text_ref}) - { - $LineIndex++; - #Enable if you need to care: -# print "D: Line Index = $LineIndex\n"; - - if ($$Text_ref[$LineIndex] =~ m/$TopTag/) - { - $AlignmentLine_Top = $LineIndex; -# print "D: [returnLinesBetween]: TopTag found in Line: '$$Text_ref[$LineIndex]'\n"; #Enable if you are interested - } - if ($$Text_ref[$LineIndex] =~ m/$BotTag/) - { - $AlignmentLine_Bot = $LineIndex; -# print "D: [returnLinesBetween]: Bottom Tag found in Line: '$$Text_ref[$LineIndex]'\n"; #Enable if you are interested - } - } -#Reality check: did we find anything? If not then we return null. -if ($AlignmentLine_Top ==0 && $AlignmentLine_Bot ==0) - { return; } -#Again, enable if you care: -#print "D: [returnLinesBetween] Lines for section table: '$AlignmentLine_Top to $AlignmentLine_Bot'\n"; - -#We want the lines one down and one up - so polish these. -$AlignmentLine_Top++; $AlignmentLine_Bot--; - -#Return as an array slice: -return (@$Text_ref[$AlignmentLine_Top .. $AlignmentLine_Bot]); -} -############ - -sub findSplitPoint -{ -=head2 sub: $PanelBoundaryCahracter = findSplitPoint (\@Table) - -When passed a table with the alignment in it makes an educated guess as to the precise split point to -spearate the 'info' and 'alignment' panels. -This is a right olde faff because the field / panel boundaries change. - - ' Query_6 167 GAGGTGCAGTTGTTGGAGTCTGGGGGAGGCTTGGCACAGCC-GGGGGGTCCCTGAGACTCTCCTGTGCAG 235' - ' Query_6 236 CCTCTGGATTCACCTTTGACAAATATGCCATGACCTGGGTCCGCCAGGCTCCAGGGAAGGGTCTGGAGTG 305' - ' Query_6 306 GGTCTCAACTATACTTGCCAGTGGTCG---CACAGACGACGCAGACTCCGTGAAGGGCCGGTTTGCCATC 372' - ' Query_6 373 TCCAGAGACAATTCCAAGAACACTCTGTATCTGCAAATGAACAGCCTGAGAGTCGAGGACACGGCCCTTT 442' - ' Query_6 443 ATTACTGTGCGAGTGAGGGGGACATAGTGGCTTCGGAGCTTTTGAGTACTGGGGCCAGGGAAACCTGGTC 512' -MOTIF_FOUND_IN_AA -i.e to contain just ATGC + "X" bases & the gap "-" character but not the "." character (found in the alingment proper) and have 4 fields in total - -Returns either -1 or the location of the panel boundary, issues a warning and returns -1 if is the most frequent boundary -because the pattern match has been failing more often that it suceeded. - -=cut -#A rough guess is 38 for normal sequences, 48 for reversed ones: - -my $SplitPos = 0; - -(my $Table_ref) = @_; #Get the reference to the table -my @DNALines; #We populate this for mining in the next section -foreach my $C_Line (@{$Table_ref}) - { - #print "D: $C_Line\n"; -# (my $SplitLine) = $C_Line; - #Split on consecutive tabs or spaces: - my @LineFields = split (/[\t\s]+/,$C_Line); - #print "D: Split Line: '",join (",",@LineFields),"' : $#LineFields\n"; - unless ( $LineFields[3] =~ m/[^\.]/ - && $LineFields[3] =~ m/[ATGCX]{20,}/ - && $#LineFields==4) - { next; } -#Enable if you want to know the lines we think are the DNA Query strings: - #print "D: DNA Line: '$C_Line'\n"; - push @DNALines, $C_Line; #Note it down - } - -my %PanelBounds; #Will contain the positions of the panel boundaries - -foreach my $C_DNALine (@DNALines) - { - #print "D: '$C_DNALine'\n"; - $C_DNALine =~ m/[ATGC-]+ \d+$/; #Match the DNA string and the indexingMOTIF_FOUND_IN_AA numbers afterwards, allow gap characters. - my $MatchPos = $-[0]; #This is the position of the start of the last match because we can't get the index() function to work - #(my $MatchPos) = index ($C_DNALine, / [ATGCX-]{20}/,0); - #print "D: '$C_DNALine' DNA panel starts at:'$MatchPos'\n"; - $PanelBounds{$MatchPos}++; - } -#Sort the hash values in order and then return the most frequent (will offer some resistance to the occasion pattern failure) -#The brackets around "($SplitPos)" are really necessary it seems. -($SplitPos) = (sort { $a <=> $b } keys %PanelBounds); -#If you want -#print Dumper %PanelBounds; -#Tell people if we are having difficultlty: -if ($SplitPos == -1) { warn "Couldn't identify the panel boundaries\n"; } -#print "D: $SplitPos: Returning the split position of: '$SplitPos'\n"; -return $SplitPos; -} - - -## -# -# -### - - - - - -##### -# -# -##### -sub markUpCDR3 -{ -=head3 Sub: (Start, End, Found How) = markUpCDR3 (DNASeq, AASeq, FWR3 End, FWR1 Offset, DNA Regex, AA Regex) - -Tries to identify the end of the CDR3 using the DNA and RNA Sequence patterns MOTIF_FOUND_IN_AAsupplied. The CDR3 is assumed to start -at the end of the FWR3. -To reduce FP matches only the sequences (DNA & AA) after the FWR3 are tested with the pattern. -The position of the first matching pattern is reported. - -=head4 Fuller Usage: - -my ($CDR3_Start, my $CDR3_End) = markUpCDR3 ($DNAQuerySequence, $AAQuerySequence, - $DomainBoundaries{"FWR3"}{"End"}, $DomainBoundaries{"FWR1"}{"Start"}, - $DNACDR3_Pat, $AACDR3_Pat); - - - -=head4 Returned Values - -If the CDR3 was found then we we signal like this: - - $MotifFound ==0 : Nope, didn't find either motif - $MotifFound ==1 : Found at the DNA level, not the AA level - $MotifFound ==2 : Found at the the AA level, not the DNA level - $MotifFound ==3 : Found at the the AA level & the DNA level - -(Also remember that if the FWR3 region couldn't be identified in the sequence there is a 4th option: not tested; this routine isn't called therefore) - -The Start and Ends returned are from the first sucessful match (MotifFound==3): though hopefully they are the same. -Formally the test order is: - - 1) DNA - 2) AA - -i.e. DNA bp locations have priority. - -Technically the locations are determined by a regex match then the $+[0] array (i.e. the end of the pattern match). -See pages like this: http://stackoverflow.com/questions/87380/how-can-i-find-the-location-of-a-regex-match-in-perl for an explanation. - -=head3 Manipulation of AA patternsMOTIF_FOUND_IN_AA - -Note that patterns are assumed to require white space inserting in them between the letters. -This could be a serious limitation - - -=cut - -#Get the parameters passed: -my ($DNA, $AA, $FWR3_End, $FWR1_Start, $DNAPat, $AAPat) = @_; -print "D: markUpCDR3: Passed Parameters '$FWR3_End, $FWR1_Start, $DNAPat, $AAPat' (& AA & DNA sequence)\n"; - - -#Setup our return values: -my $Start = 0; my $End =0; my $MotifFound = 0; -my $How; #Literally How the motif was found (or not if blank) - - -=head4 Prepare the sequences and the patterns for use - -Specifically: trim off the start of the AA & DNA string already allocated to other CDRs or FWRs - -Add in spaces into the AA regex pattern because we can't get regex-ex freespacing mode i.e. "$Var =~ m/$AAPat/x" working. - - -We take the "-1" as the CropPoint position to include the previous 3 nucleotides / AAs; remember to add this back on -in position calculations. - - -=cut - -#Because igBLAST doesn't always report from the start of the read (primers and things are upstream): - -my $CropPoint = $FWR3_End - $FWR1_Start - 1 ; -#print "D: markUpCDR3: Crop point is: '$CropPoint'\n"; - -#print "D: markUpCDR3: Cropping point is: '$CropPoint' characters from start\n"; -#We trim off the parts we expect to find the CDR3 motifs in leaving at extra 3nts on to allow for base miss-calling: - -my $AA_Trimmed = substr ($AA, $CropPoint); -my $DNA_Trimmed = substr ($DNA ,$CropPoint); -#print "D: markUpCDR3: AA = '$AA' (untrimmed)\nD: markUpCDR3: TR = '$AA_Trimmed' (Trimmed) ", length ($AA_Trimmed)," nts long\n"; -#print "D: markUpCDR3: Testing: AA = '$AA_Trimmed', DNA = '$DNA_Trimmed'\n"; - -#This lovely hack is to account for the spaces in the AA sequence and we can't get the "$Var =~ m/$AAPat/x" working -my $AAPat_Spaced; -foreach my $C_Char (0..length($AAPat)-1) #The -1 is because we don't want trailing spaces until the next nt -> AA translation. - { $AAPat_Spaced = $AAPat_Spaced.'\s+'.substr ($AAPat,$C_Char,1); } -#And write this back into the main pattern we were passed: -$AAPat = $AAPat_Spaced; - -#temp hack: -#$AA_Trimmed = $AA; -my $MotifFound=0; #So we can record which patterns we found -my $MotifPositionDNA =-1; -my $MotifPositionAA =-1; - -#print "D: markUpCDR3: Pattern: '$AAPat_Spaced'\n"; -=head4 At DNA level: "TGG GGx xxx GGx" [+1] - -=cut - -#print "D: markUpCDR3: '$DNA_Trimmed' (Trimmed DNA string)\n"; - -if ($DNA_Trimmed =~ m/$DNAPat/) - { - $MotifPositionDNA = $+[0]; #Just the easiest way to do this in Perl -# print "D: markUpCDR3:: Found Motif match on DNA at bp: '$MotifPositionDNA'\n"; - $MotifFound = $MotifFound + 1; - #Any more matches further on? - my $LaterString = substr ($DNA_Trimmed, $MotifPositionDNA); -# print "D: markUpCDR3: '$AA_Trimmed' (AA Trimmed string)\n"; -# print "D: markUpCDR3: '", substr ($DNA_Trimmed,0, $MotifPositionDNA)," (DNA until pattern match string)\n"; -# print "D: markUpCDR3: '$DNA_Trimmed' (Trimmed DNA string)\n"; -# print "D: markUpCDR3: '$LaterString' (Later part of DNA string)\n"; - if ($LaterString =~ m/$DNAPat/) - { print "D: markUPCDR3: Also got a match further down the DNA String: at ", $-[0] ," to ", $+ [0], " - which might be worrying\n"; } - } - -=head4 At AA level: "WGxG" [+2] - -=cut - -if ($AA_Trimmed=~ m/$AAPat/) - { - $MotifPositionAA = $+[0]; #Just the easiest way to do this in Perl - $MotifFound = $MotifFound + 2; -# print "D: markUpCDR3: Found Motif match on AA at position (on DNA remember): '$MotifPositionAA' (ie.)\n"; - (my $CDR3_seq) = substr ($AA_Trimmed, 0, $MotifPositionAA); -# print "D: markUpCDR3: Seq ='$CDR3_seq' - as detected\n"; - - } - -=head4 Assess the results of motif position finding - -=cut - -#print "D: markUpCDR3: MotifFound = '$MotifFound'\n"; - -if ($MotifFound ==0) - { return ($Start, $End, $MotifFound); } #The easy one really: return we didn't find the CDR3 - -# -$Start = $FWR3_End; #We assume the end of the FWR3 is the start of CDR3: -#Just found in DNA: -if ($MotifFound ==1) - { - $Start = $FWR3_End; #We assume the end of the FWR3 is the start of CDR3: - $End = $MotifPositionDNA; - $How = "MOTIF_FOUND_IN_DNA"; - } -#Just found in AA: -if ($MotifFound ==2) - { - $End = $MotifPositionAA; - $How = "MOTIF_FOUND_IN_AA"; - } - -#Found in both, DNA has priority: -if ($MotifFound ==3) - { - $Start = $FWR3_End ; #We assume the end of the FWR3 is the start of CDR3: - $End = $MotifPositionDNA; - $How = "MOTIF_FOUND_IN_BOTH"; - } - -#print "D: markUpCDR3: Motif found = $MotifFound\n"; - -=head4 These next few lines are for testing / diagnostics only - disable for general use - -If you are interested in getting the CDR3 directly then remember the main coordinate system is defined such that -the start of FWR1 is unlikely to be at nt 1. - -=cut - -$Start = $FWR3_End - $FWR1_Start -1; -$End = $End + $CropPoint; -my $CDR3_RegionLength = $End - $Start; -#print "D: markUpCDR3: CDR3 Length= $Start - $End = '$CDR3_RegionLength'\n"; -(my $CDR3_seq) = substr ($AA, $Start, $CDR3_RegionLength); - -#Add onto the coordinates what we trimmed off: - - -#print "D: markUpCDR3: Seq ='$CDR3_seq'\n"; - -print "D: markUpCDR3: returning: $Start, $End, $How, ($MotifFound) [NB: offset of :'+ $FWR1_Start'\n"; -#die "HIT BLOCK\n"; -return ($Start + $FWR1_Start, $End + $FWR1_Start, $How); -} - - -sub printOUTPUTData { -=head2 sub: $OutputDataString = printOUTPUTData {\%OutputData} - -When passed an array containing the appropriate CDR, Top V / D/ J genes and the seqeunce ID. -This prepared and then returned as a text string that can then be printed to STDOUT: - - print (printOUTPUTData (\%OutputData)); - -Any missing data in the Hash array it polietly ignored and a null string printed in place. -The text field is tab delimited; there are no extra trailing tabs or carriage returns in place. - -Actually the fields printed out are stored in an index array. - -=head3 Header output - -If the routine is passed a key 'HEADER' then the header columns are returned as that string. -This is tested first - so don't add this unless you mean to. - -=cut - -my @HeaderFields = ("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"); - -my $OutputString = "OUTPUT:"; #What we are going to build the output into. - -=head4 Print Header & Exit? - -=cut - -my ($Data_ref) = @_; -#print "D: printOUTPUTData: Running\n"; - -if (exists $$Data_ref {"HEADER"}) - { - $OutputString .= "\t"; - for(my $n = 0; $n <= $#HeaderFields; $n++) - { - $OutputString .= $HeaderFields[$n]; - $OutputString .= "\t" if($n < $#HeaderFields); - } - - # foreach my $C_Header (@HeaderFields) - # { $OutputString .= "$C_Header"; } # - - print "D: printOUTPUTData: HEADER Printout requested '@HeaderFields'\n"; - return ($OutputString); - } - -=head3 Assemble whatever data we have - and tab delimit the null fields - -=cut -#print "D: printOUTPUTData: Will pretty print this:\n", Dumper $Data_ref; -foreach my $C_Header (@HeaderFields) - { - - if (exists ($$Data_ref {$C_Header})) - { $OutputString .= "\t". $$Data_ref{$C_Header}; } #We have data to print out - else - { $OutputString .="\t"; } #Add a trailing space - } # - -return ($OutputString); -} - - -######################################### Code Junk ######################## - - -=head2 Code Junk Attic - -=head3 Demonstrates how to reverse translate an amino acid sequence into DNA: - -use Bio::Tools::CodonTable; -use Bio::Seq; - -# print possible codon tables - my $tables = Bio::Tools::CodonTable->tables; - while ( (my $id, my $name) = each %{$tables} ) { - print "$id = $name\n"; - } - my $CodonTable = Bio::Tools::CodonTable->new(); - - my $ExampleSeq = Bio::PrimarySeq->new(-seq=>"WGxG", -alphabet => 'protein') or die "Cannot create sequence object\n"; - - -my $rvSeq = $CodonTable->reverse_translate_all($ExampleSeq); -print "D: '$rvSeq'\n"; -die "TEST OVER\n"; - -=cut - - -=head3 For processing the 'Alignment lines' section of the alginment table - - #If we are ever interested; then enable the code below: -# print ": Alignment\n"; -# $InfoPanel =~ s/^ +//; $InfoPanel =~ s/ +$//; #Clean off trailing spaces -# my ($Germclass, $PID, $PID_Counts, $Allele) = split (/\s+/,$InfoPanel); #Split on spaces -##Enable if you need to know what we just found: -# #print "D: Fields are (Germclass, PID, PID_Counts, Allele) \t$Germclass, $PID, $PID_Counts, $Allele\n"; -# #A reality check: we should have an Allele - or some text here. -# unless (defined $Allele && $Allele ne "") -# { warn "Cannot get Allele for Line '$C_Line' - implies improper parsing: '",substr ($Lines[$C_Line],0,15),"...'\n"; } -# if (exists ($Alginments {$Germclass}{$Allele})) -# { $Alginments {$Germclass}{$Allele} = $Alginments {$Germclass}{$Allele}.$CurrentAASequence; } #Carry on adding -# else #more work needed as we need to 'pad' the sequence with fake gap characters) -# { -##Do we still need this padding? I don't think so -# -# -# my $PaddingChars = ($ThisQueryStart-$Query_Start); -# print "D: New gene found: need to pad it with ($ThisQueryStart-$Query_Start) i.e. '$PaddingChars' characters\n"; -# #To help testing, calculate this first: -# my $PaddingString = " "x $PaddingChars; -# $Alginments {$Germclass}{$Allele} = $CurrentAASequence; -# } -# next - -=head3 Demonstration of Pattern match positions - -my $Text = "12345TTT TTAAAAA"; -my $TestPat = "TTT\\s+TT"; -(my $Result)= $Text =~ m/$TestPat/; -print "D: Two vars are: - = ",$-[0], " & + =", $+[0]," for test pattern '$TestPat'\n"; - -sub printCDR3 { - -=head3 Subroutine: printCDR3 ($CDR3_Start, $CDR3_End, "SUMMARY_TABLE", $AAQuerySequence, $DNAQuerySequence); - -???? IS THIS FUNCTION IN USE ????? - -Handles the printing of the output when passed information about the CDR3 region. - - -The result is sent returned as a text string in this version hence use it like this if you want to send it to STDOUT: - - print printCDR3 ($CDR3_Start, $CDR3_End, "SUMMARY_TABLE", $AAQuerySequence, $DNAQuerySequence), "\n"; - -#=cut - -#Despite the similarity in names, these are all local copies passed to us: - -my ($Start, $End, $Tag, $FullAAQuerySequence, $FullDNAQuerySequence) = @_; - -#For DNA: -my ($CDR_DNA_Seq) = substr ($FullDNAQuerySequence, $Start, $Start+$End); -my ($CDR_DNA_Length) = length ($CDR_DNA_Seq); - -#For AA: -my ($CDR_AA_Seq) = substr ($FullAAQuerySequence, $Start, $Start+$End); -my ($CDR_AA_Length) = length ($CDR_AA_Seq); - -my $ReturnString = join ("\t", $CDR_DNA_Seq, $CDR_DNA_Length, $CDR_AA_Seq, $CDR_AA_Length, $Tag); #Create here so we can inspect it / post process it if needed: -print "D: SUB: printCDR3: As returned: '$ReturnString'\n"; -return ($ReturnString); - -} - -=cut - - - -=head2 Change Log - -=head3 Version 1.2 - - 1) Fixed the 'Process recrod request' feature' [was failed increment in $Record] - 2) Deleted / Deactivated the function 'printCDR3' [wasn't in used; kept if useful for parts]. - This function is replaced by the more general printOUTPUTData() - 3) A tag for the CDR3 status is now output for every record / read. - Initially this is set to "NOT_FOUND" and changed if evidence for the CDR3 is found. - -=head4 Version 1.3 - - 1) The tophit line was split on whitespace, however sometimes the VJFrame is something like “In-frame with stop codon”, - which means the line is also split on the spaces therein. It now splits on tabs only, and this seems to work properly. - - found by Bas Horsman. - -=head4 Version 1.3a - - 1) "MOTIF_FOUND_IN_AA" reported correctly (was impossible previously due to addition error to the $MotifFound var (never could == 3) - -=cut - -=head4 Version 1.4 - - 1) Now processes files using Mac/Unix/MS-DOS newline characters: - - $_ =~ s/\r\n/\n/g; #In case line ends are MS-DOS - $_ =~ s/\r/\n/g; #In case line ends are Mac - #The whole record - one per read - is now stored in $_ - my @Lines =split (/\R/,$_); #Split on new lines - -=head4 Version 1.4a - -1) Fixed the length of the CDR3 AA string being reported correctly: - - $OUTPUT_Data{"CDR3 Length"} = $CDR3_Length; - to: - $OUTPUT_Data{"CDR3 Length"} = $CDR3_Seq_AA_Length; -
--- a/immunerepertoirecombined_imgt/imgtconvert.py Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,139 +0,0 @@ -import pandas as pd -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("--input", help="Input folder with files") -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'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'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"] - -inputFolder = args.input - -dirContents = os.listdir(inputFolder) -if len(dirContents) == 1: - inputFolder = os.path.join(inputFolder, dirContents[0]) - if os.path.isdir(inputFolder): - print "is dir" - dirContents = os.listdir(inputFolder) -files = sorted([os.path.join(inputFolder, f) for f in dirContents]) - -if len(files) % 3 is not 0: - stop_err("Files in zip not a multiple of 3, it should contain the all the 1_, 5_ and 6_ files for a sample") - import sys - sys.exit() - - -triplets = [] -step = len(files) / 3 -for i in range(0, step): - triplets.append((files[i], files[i + step], files[i + step + step])) - -outFile = args.output - -fSummary = pd.read_csv(triplets[0][0], sep="\t") -fSequence = pd.read_csv(triplets[0][1], sep="\t") -fJunction = pd.read_csv(triplets[0][2], sep="\t") -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 - -for triple in triplets[1:]: - fSummary = pd.read_csv(triple[0], sep="\t") - fSequence = pd.read_csv(triple[1], sep="\t") - fJunction = pd.read_csv(triple[2], sep="\t") - - 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 = outFrame.append(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', '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'] - -vPattern = re.compile(r"IGHV[1-9]-[0-9ab]+-?[1-9]?") -dPattern = re.compile(r"IGHD[1-9]-[0-9ab]+") -jPattern = re.compile(r"IGHJ[1-9]") - -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"] = 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 D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? -outFrame.to_csv(outFile, sep="\t", index=False) \ No newline at end of file
--- a/immunerepertoirecombined_imgt/imgtconvert.sh Mon Dec 09 03:45:21 2013 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -#!/bin/bash -dir="$(cd "$(dirname "$0")" && pwd)" -mkdir $PWD/$2_$3 - - -#!/bin/bash -f=$(file $1) -zip7Type="7-zip archive" -tarType="tar archive" -bzip2Type="bzip2 compressed" -gzipType="gzip compressed" -zipType="Zip archive" -rarType="RAR archive" - -if [[ "$f" == *"$zip7Type"* ]]; then - echo "7-zip" - echo "Trying: 7za e $1 -o$PWD/$2_$3/" - 7za e $1 -o$PWD/$2_$3/ -fi - -if [[ "$f" == *"$tarType"* ]] -then - echo "tar archive" - echo "Trying: tar xvf $1 -C $PWD/$2_$3/" - tar xvf $1 -C $PWD/$2_$3/ -fi - -if [[ "$f" == *"$bzip2Type"* ]] -then - echo "bzip2 compressed data" - echo "Trying: tar jxf $1 -C $PWD/$2_$3/" - tar jxf $1 -C $PWD/$2_$3/ -fi - -if [[ "$f" == *"$gzipType"* ]] -then - echo "gzip compressed data" - echo "Trying: tar xvzf $1 -C $PWD/$2_$3/" - tar xvzf $1 -C $PWD/$2_$3/ -fi - -if [[ "$f" == *"$zipType"* ]] -then - echo "Zip archive" - echo "Trying: unzip $1 -d $PWD/$2_$3/" - unzip $1 -d $PWD/$2_$3/ -fi - -if [[ "$f" == *"$rarType"* ]] -then - echo "RAR archive" - echo "Trying: unrar e $1 $PWD/$2_$3/" - unrar e $1 $PWD/$2_$3/ -fi -find $PWD/$2_$3/ -type f | grep -v "1_Summary_\|5_AA-sequences_\|6_Junction_" | xargs rm -f -python $dir/imgtconvert.py --input $PWD/$2_$3 --output $4 -