Mercurial > repos > davidvanzessen > imgt_loader_igg
changeset 7:04e72fc8b2c4 draft default tip
Uploaded
author | davidvanzessen |
---|---|
date | Fri, 05 Sep 2014 04:21:48 -0400 |
parents | 5b030e48b308 |
children | |
files | imgt_loader.py imgt_loader.sh imgt_loader.xml imgtconvert.py imgtconvert.sh imgtconvert.xml |
diffstat | 6 files changed, 217 insertions(+), 267 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imgt_loader.py Fri Sep 05 04:21:48 2014 -0400 @@ -0,0 +1,142 @@ +import pandas as pd +try: + pd.options.mode.chained_assignment = None # default='warn' +except: + pass +import re +import argparse +import os + +def stop_err( msg, ret=1 ): + sys.stderr.write( msg ) + sys.exit( ret ) + +#docs.python.org/dev/library/argparse.html +parser = argparse.ArgumentParser() +parser.add_argument("--summ", help="The 1_Summary file from the imgt output") +parser.add_argument("--aa", help="The 5_AA-Sequence file from the imgt output") +parser.add_argument("--junction", help="The 6_Junction file from the imgt output") +parser.add_argument("--output", help="Output file") + +args = parser.parse_args() + +old_summary_columns = [u'Sequence ID', u'JUNCTION frame', u'V-GENE and allele', u'D-GENE and allele', u'J-GENE and allele', u'CDR1-IMGT length', u'CDR2-IMGT length', u'CDR3-IMGT length', u'Orientation'] +old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT'] +old_junction_columns = [u'JUNCTION'] + +added_summary_columns = [u'Functionality', u'V-REGION identity %', u'V-REGION identity nt', u'D-REGION reading frame', u'AA JUNCTION', u'Functionality comment', u'Sequence'] +added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT'] +added_junction_columns = [u"P3'V-nt nb", u'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"] + +outFile = args.output + +#fSummary = pd.read_csv(triplets[0][0], sep="\t", low_memory=False) +fSummary = pd.read_csv(args.summ, sep="\t", dtype=object) +#fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False) +fSequence = pd.read_csv(args.aa, sep="\t", dtype=object) +#fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False) +fJunction = pd.read_csv(args.junction, sep="\t", dtype=object) +tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]] + +tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"] +tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"] + +tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"] +tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"] + +tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"] +tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"] + +tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"] +tmp["CDR3 Length DNA"] = '1' +tmp["Strand"] = fSummary["Orientation"] +tmp["CDR3 Found How"] = 'a' + +for col in added_summary_columns: + tmp[col] = fSummary[col] + +for col in added_sequence_columns: + tmp[col] = fSequence[col] + +for col in added_junction_columns: + tmp[col] = fJunction[col] + +outFrame = tmp + +outFrame.columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', u'CDR3 Seq', u'CDR3 Length', u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'Functionality', 'V-REGION identity %', 'V-REGION identity nt', 'D-REGION reading frame', 'AA JUNCTION', 'Functionality comment', 'Sequence', 'FR1-IMGT', 'FR2-IMGT', 'FR3-IMGT', 'CDR3-IMGT', 'JUNCTION', 'J-REGION', 'FR4-IMGT', 'P3V-nt nb', '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'] + +""" +IGHV[0-9]-[0-9ab]+-?[0-9]?D? +TRBV[0-9]{1,2}-?[0-9]?-?[123]? +IGKV[0-3]D?-[0-9]{1,2} +IGLV[0-9]-[0-9]{1,2} +TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])? +TRGV[234589] +TRDV[1-3] + +IGHD[0-9]-[0-9ab]+ +TRBD[12] +TRDD[1-3] + +IGHJ[1-6] +TRBJ[12]-[1-7] +IGKJ[1-5] +IGLJ[12367] +TRAJ[0-9]{1,2} +TRGJP?[12] +TRDJ[1-4] +""" + +vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)", + r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)", + r"(IGKV[0-3]D?-[0-9]{1,2})", + r"(IGLV[0-9]-[0-9]{1,2})", + r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)", + r"(TRGV[234589])", + r"(TRDV[1-3])"] + +dPattern = [r"(IGHD[0-9]-[0-9ab]+)", + r"(TRBD[12])", + r"(TRDD[1-3])"] + +jPattern = [r"(IGHJ[1-6])", + r"(TRBJ[12]-[1-7])", + r"(IGKJ[1-5])", + r"(IGLJ[12367])", + r"(TRAJ[0-9]{1,2})", + r"(TRGJP?[12])", + r"(TRDJ[1-4])"] + +vPattern = re.compile(r"|".join(vPattern)) + +dPattern = re.compile(r"|".join(dPattern)) + +jPattern = re.compile(r"|".join(jPattern)) + + +def filterGenes(s, pattern): + if type(s) is not str: + return "NA" + res = pattern.search(s) + if res: + return res.group(0) + return "NA" + + + +outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern)) +outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern)) +outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern)) + +print outFrame + +tmp = outFrame["VDJ Frame"] +tmp = tmp.replace("in-frame", "In-frame") +tmp = tmp.replace("null", "Out-of-frame") +tmp = tmp.replace("out-of-frame", "Out-of-frame") +outFrame["VDJ Frame"] = tmp +outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len) +safeLength = lambda x: len(x) if type(x) == str else 0 +outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? +#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? +outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imgt_loader.sh Fri Sep 05 04:21:48 2014 -0400 @@ -0,0 +1,59 @@ +#!/bin/bash +input=$1 +output=$2 +name=$3 +dir="$(cd "$(dirname "$0")" && pwd)" +mkdir -p $PWD/$name/files +f=$(file $input) +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 $input -o$PWD/files/" + 7za e $input -o$PWD/$name/files +fi + +if [[ "$f" == *"$tarType"* ]] +then + echo "tar archive" + echo "Trying: tar xvf $input -C $PWD/files/" + tar xvf $input -C $PWD/$name/files +fi + +if [[ "$f" == *"$bzip2Type"* ]] +then + echo "bzip2 compressed data" + echo "Trying: tar jxf $input -C $PWD/files/" + tar jxf $input -C $PWD/$name/files +fi + +if [[ "$f" == *"$gzipType"* ]] +then + echo "gzip compressed data" + echo "Trying: tar xvzf $input -C $PWD/files/" + tar xvzf $input -C $PWD/$name/files +fi + +if [[ "$f" == *"$zipType"* ]] +then + echo "Zip archive" + echo "Trying: unzip $input -d $PWD/files/" + unzip $input -d $PWD/$name/files > $PWD/unziplog.log +fi + +if [[ "$f" == *"$rarType"* ]] +then + echo "RAR archive" + echo "Trying: unrar e $input $PWD/files/" + unrar e $input $PWD/$name/files +fi +find $PWD/$name/files -iname "1_*" -exec cat {} + > $PWD/$name/summ.txt +find $PWD/$name/files -iname "5_*" -exec cat {} + > $PWD/$name/aa.txt +find $PWD/$name/files -iname "6_*" -exec cat {} + > $PWD/$name/junction.txt + +python $dir/imgt_loader.py --summ $PWD/$name/summ.txt --aa $PWD/$name/aa.txt --junction $PWD/$name/junction.txt --output $output
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imgt_loader.xml Fri Sep 05 04:21:48 2014 -0400 @@ -0,0 +1,16 @@ +<tool id="imgt_loader_igg" name="IMGT Loader" version="1.0"> + <description> </description> + <command interpreter="bash"> + imgt_loader.sh $in_file $out_file "tmp" + </command> + <inputs> + <param name="in_file" type="data" label="Archive with files" /> + </inputs> + <outputs> + <data format="tabular" name="out_file" label="IMGT Loader on ${in_file.name}"/> + </outputs> + <help> + Converting all 1_Summary, 5_AA-sequence and 6_Junction files belonging to a sample from an archive (supported: 7z, rar, tar, tar.bz2, tar.gz and zip) into a single file to be used in the Immune repertoire toolchain (Step 3) + </help> + +</tool>
--- a/imgtconvert.py Mon Jul 07 09:54:42 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,194 +0,0 @@ -import pandas as pd -try: - pd.options.mode.chained_assignment = None # default='warn' -except: - pass -import re -import argparse -import os - -def stop_err( msg, ret=1 ): - sys.stderr.write( msg ) - sys.exit( ret ) - -#docs.python.org/dev/library/argparse.html -parser = argparse.ArgumentParser() -parser.add_argument("--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'Functionality', u'V-REGION identity %', u'V-REGION identity nt', u'D-REGION reading frame', u'AA JUNCTION', u'Functionality comment', u'Sequence'] -added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT'] -added_junction_columns = [u"P3'V-nt nb", u'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): - dirContents = os.listdir(inputFolder) -files = sorted([os.path.join(inputFolder, f) for f in dirContents if os.path.isfile(os.path.join(inputFolder, f))]) - -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", low_memory=False) -fSummary = pd.read_csv(triplets[0][0], sep="\t") -#fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False) -fSequence = pd.read_csv(triplets[0][1], sep="\t") -#fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False) -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', u'Functionality', 'V-REGION identity %', 'V-REGION identity nt', 'D-REGION reading frame', 'AA JUNCTION', 'Functionality comment', 'Sequence', 'FR1-IMGT', 'FR2-IMGT', 'FR3-IMGT', 'CDR3-IMGT', 'JUNCTION', 'J-REGION', 'FR4-IMGT', 'P3V-nt nb', 'N1-REGION-nt nb', 'P5D-nt nb', 'P3D-nt nb', 'N2-REGION-nt nb', 'P5J-nt nb', '3V-REGION trimmed-nt nb', '5D-REGION trimmed-nt nb', '3D-REGION trimmed-nt nb', '5J-REGION trimmed-nt nb'] - -""" -IGHV[0-9]-[0-9ab]+-?[0-9]?D? -TRBV[0-9]{1,2}-?[0-9]?-?[123]? -IGKV[0-3]D?-[0-9]{1,2} -IGLV[0-9]-[0-9]{1,2} -TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])? -TRGV[234589] -TRDV[1-3] - -IGHD[0-9]-[0-9ab]+ -TRBD[12] -TRDD[1-3] - -IGHJ[1-6] -TRBJ[12]-[1-7] -IGKJ[1-5] -IGLJ[12367] -TRAJ[0-9]{1,2} -TRGJP?[12] -TRDJ[1-4] -""" - -vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)", - r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)", - r"(IGKV[0-3]D?-[0-9]{1,2})", - r"(IGLV[0-9]-[0-9]{1,2})", - r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)", - r"(TRGV[234589])", - r"(TRDV[1-3])"] - -dPattern = [r"(IGHD[0-9]-[0-9ab]+)", - r"(TRBD[12])", - r"(TRDD[1-3])"] - -jPattern = [r"(IGHJ[1-6])", - r"(TRBJ[12]-[1-7])", - r"(IGKJ[1-5])", - r"(IGLJ[12367])", - r"(TRAJ[0-9]{1,2})", - r"(TRGJP?[12])", - r"(TRDJ[1-4])"] - -vPattern = re.compile(r"|".join(vPattern)) - -dPattern = re.compile(r"|".join(dPattern)) - -jPattern = re.compile(r"|".join(jPattern)) - - -def filterGenes(s, pattern): - if type(s) is not str: - return "NA" - res = pattern.search(s) - if res: - return res.group(0) - return "NA" - - - -outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern)) -outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern)) -outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern)) - -print outFrame - -tmp = outFrame["VDJ Frame"] -tmp = tmp.replace("in-frame", "In-frame") -tmp = tmp.replace("null", "Out-of-frame") -tmp = tmp.replace("out-of-frame", "Out-of-frame") -outFrame["VDJ Frame"] = tmp -outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len) -safeLength = lambda x: len(x) if type(x) == str else 0 -outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? -#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows? -outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")
--- a/imgtconvert.sh Mon Jul 07 09:54:42 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -#!/bin/bash -dir="$(cd "$(dirname "$0")" && pwd)" -mkdir $PWD/files - - -#!/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/files/" - 7za e $1 -o$PWD/files/ -fi - -if [[ "$f" == *"$tarType"* ]] -then - echo "tar archive" - echo "Trying: tar xvf $1 -C $PWD/files/" - tar xvf $1 -C $PWD/files/ -fi - -if [[ "$f" == *"$bzip2Type"* ]] -then - echo "bzip2 compressed data" - echo "Trying: tar jxf $1 -C $PWD/files/" - tar jxf $1 -C $PWD/files/ -fi - -if [[ "$f" == *"$gzipType"* ]] -then - echo "gzip compressed data" - echo "Trying: tar xvzf $1 -C $PWD/files/" - tar xvzf $1 -C $PWD/files/ -fi - -if [[ "$f" == *"$zipType"* ]] -then - echo "Zip archive" - echo "Trying: unzip $1 -d $PWD/files/" - unzip $1 -d $PWD/files/ > $PWD/unziplog.log -fi - -if [[ "$f" == *"$rarType"* ]] -then - echo "RAR archive" - echo "Trying: unrar e $1 $PWD/files/" - unrar e $1 $PWD/files/ -fi -find $PWD/files/ -type f | grep -v "1_Summary_\|5_AA-sequences_\|6_Junction_" | xargs rm -f -python $dir/imgtconvert.py --input $PWD/files --output $2 -
--- a/imgtconvert.xml Mon Jul 07 09:54:42 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,16 +0,0 @@ -<tool id="imgt_loader_igg" name="IMGT Loader" version="1.0"> - <description> </description> - <command interpreter="bash"> - imgtconvert.sh $in_file $out_file - </command> - <inputs> - <param name="in_file" type="data" label="Archive with files" /> - </inputs> - <outputs> - <data format="tabular" name="out_file" label="IMGT Loader on ${in_file.name}"/> - </outputs> - <help> - Converting all 1_Summary, 5_AA-sequence and 6_Junction files belonging to a sample from an archive (supported: 7z, rar, tar, tar.bz2, tar.gz and zip) into a single file to be used in the Immune repertoire toolchain (Step 3) - </help> - -</tool>