Moved test data to test-data, tested out script

rpetit3 · Mar 23, 2016 · 3954fd3 · 3954fd3
1 parent 0dab6a1
commit 3954fd3
Show file tree

Hide file tree

Showing 10 changed files with 33,624 additions and 26,473 deletions.
diff --git a/bin/vcf-annotator b/bin/vcf-annotator
@@ -10,18 +10,17 @@ if __name__ == '__main__':
     parser = ap.ArgumentParser(
         prog='vcf-annotator.py',
         conflict_handler='resolve',
-        description="".join(["Annotate SNPs from a VCF file using the ",
-                             "reference genome's GenBank file."])
+        description=("Annotate SNPs from a VCF file using the reference "
+                     "genome's GenBank file.")
     )
-    group1 = parser.add_argument_group('Options', '')
-    group1.add_argument('gb', metavar="GENBANK_FILE", type=str,
-                        help='GenBank file of the reference genome.')
-    group1.add_argument('vcf', metavar="VCF_FILE", type=str,
+    parser.add_argument('vcf', metavar="VCF_FILE", type=str,
                         help='VCF file of SNPs')
-    group1.add_argument('--output', metavar="STRING", type=str,
+    parser.add_argument('gb', metavar="GENBANK_FILE", type=str,
+                        help='GenBank file of the reference genome.')
+    parser.add_argument('--output', metavar="STRING", type=str,
                         default='/dev/stdout',
                         help='File to write VCF output to (Default STDOUT).')
-    group1.add_argument('-h', '--help', action='help',
+    parser.add_argument('-h', '--help', action='help',
                         help='Show this help message and exit')
 
     if len(sys.argv) == 1:

diff --git a/vcfannotator/test/__init__.py → test-data/__init__.py b/vcfannotator/test/__init__.py → test-data/__init__.py
diff --git a/vcfannotator/test/circovirus.gb → test-data/circovirus.gb b/vcfannotator/test/circovirus.gb → test-data/circovirus.gb
diff --git a/vcfannotator/test/circovirus.vcf → test-data/circovirus.vcf b/vcfannotator/test/circovirus.vcf → test-data/circovirus.vcf
diff --git a/vcfannotator/test/test_01.gb → test-data/test_01.gb b/vcfannotator/test/test_01.gb → test-data/test_01.gb
diff --git a/test-data/test_01.vcf b/test-data/test_01.vcf
@@ -0,0 +1,45 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=HARD_TO_VALIDATE,Description="MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)">
+##FILTER=<ID=LowConf,Description="DP < 10|| QD < 2.5">
+##FILTER=<ID=SnpCluster,Description="SNPs found in clusters">
+##FILTER=<ID=StdFilter,Description="QUAL < 30.0 || DP < 5 || QD < 1.5 || MQ < 25.0">
+##FILTER=<ID=gatkFilter,Description="MQ < 40.0||FS > 60.0">
+##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
+##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
+##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
+##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
+##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
+##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
+##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
+##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
+##INFO=<ID=Dels,Number=1,Type=Float,Description="Fraction of Reads Containing Spanning Deletions">
+##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
+##INFO=<ID=HRun,Number=1,Type=Integer,Description="Largest Contiguous Homopolymer Run of Variant Allele In Either Direction">
+##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplotypes">
+##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
+##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
+##INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">
+##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
+##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
+##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
+##INFO=<ID=SB,Number=1,Type=Float,Description="Strand Bias">
+##UnifiedGenotyper="analysis_type=UnifiedGenotyper input_file=[107/107-SNP/107_sdrcsm.bam] read_buffer_size=null phone_home=STANDARD gatk_key=null read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL reference_sequence=107/107-SNP/tmp/bwa/index/ref.fa nonDeterministicRandomSeed=false downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=OFF baqGapOpenPenalty=40.0 performanceLog=null useOriginalQualities=false BQSR=null quantize_quals=-1 defaultBaseQualities=-1 validation_strictness=SILENT unsafe=null num_threads=1 num_cpu_threads=null num_io_threads=null num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false logging_level=INFO log_to_file=null help=false genotype_likelihoods_model=BOTH p_nonref_model=EXACT heterozygosity=0.001 pcr_error_rate=1.0E-4 genotyping_mode=DISCOVERY output_mode=EMIT_VARIANTS_ONLY standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 noSLOD=false annotateNDA=false alleles=(RodBinding name= source=UNBOUND) min_base_quality_score=17 max_deletion_fraction=0.05 max_alternate_alleles=3 min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indel_heterozygosity=1.25E-4 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 noBandedIndel=false indelDebug=false ignoreSNPAlleles=false dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub NO_HEADER=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_mismatching_base_and_quals=false"
+##VariantFiltration="analysis_type=VariantFiltration input_file=[] read_buffer_size=null phone_home=STANDARD gatk_key=null read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL reference_sequence=107/107-SNP/tmp/bwa/index/ref.fa nonDeterministicRandomSeed=false downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 performanceLog=null useOriginalQualities=false BQSR=null quantize_quals=-1 defaultBaseQualities=-1 validation_strictness=SILENT unsafe=null num_threads=1 num_cpu_threads=null num_io_threads=null num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false logging_level=INFO log_to_file=null help=false variant=(RodBinding name=variant source=107/107-SNP/tmp/GATK/GATK.vcf) mask=(RodBinding name= source=UNBOUND) out=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub NO_HEADER=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub filterExpression=[QUAL < 30.0 || DP < 5 || QD < 1.5 || MQ < 25.0, MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1), DP < 10|| QD < 2.5, MQ < 40.0||FS > 60.0] filterName=[StdFilter, HARD_TO_VALIDATE, LowConf, gatkFilter] genotypeFilterExpression=[] genotypeFilterName=[] clusterSize=3 clusterWindowSize=10 maskExtension=0 maskName=Mask missingValuesInExpressionsShouldEvaluateAsFailing=false invalidatePreviousFilters=false filter_mismatching_base_and_quals=false"
+##contig=<ID=gi|29165615|ref|NC_002745.2|,length=2814816>
+##reference=file:///data1/home/alam/Project_Staph/visa_project/107/107-SNP/tmp/bwa/index/ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	GATK	Remarks
+gi|29165615|ref|NC_002745.2|	12490	.	C	T	3479.51	PASS	AC=2;AF=1.00;AN=2;DP=88;Dels=0.00;FS=0.000;HRun=1;HaplotypeScore=0.0000;MQ=60.00;MQ0=0;QD=39.54;SB=-1409.92	GT:AD:DP:GQ:PL	1/1:0,88:88:99:3480,259,0	Superpass
+gi|29165615|ref|NC_002745.2|	24165	.	G	A	3405.52	PASS	AC=2;AF=1.00;AN=2;BaseQRankSum=1.294;DP=87;Dels=0.00;FS=3.713;HRun=2;HaplotypeScore=0.9665;MQ=59.79;MQ0=0;MQRankSum=0.299;QD=39.14;ReadPosRankSum=1.453;SB=-1412.88	GT:AD:DP:GQ:PL	1/1:1,86:87:99:3406,221,0	Superpass
+gi|29165615|ref|NC_002745.2|	26374	.	G	A	4459.12	PASS	AC=2;AF=1.00;AN=2;BaseQRankSum=-0.093;DP=112;Dels=0.00;FS=3.332;HRun=1;HaplotypeScore=0.0000;MQ=60.00;MQ0=0;MQRankSum=0.928;QD=39.81;ReadPosRankSum=0.186;SB=-2029.02	GT:AD:DP:GQ:PL	1/1:1,111:112:99:4459,296,0	Superpass
+gi|29165615|ref|NC_002745.2|	33455	.	G	T	3172.70	PASS	AC=2;AF=1.00;AN=2;DP=80;Dels=0.00;FS=0.000;HRun=0;HaplotypeScore=0.0000;MQ=58.35;MQ0=0;QD=39.66;SB=-1566.88	GT:AD:DP:GQ:PL	1/1:0,79:80:99:3206,235,0	Superpass
+gi|29165615|ref|NC_002745.2|	123267	.	A	T	436.79	SnpCluster	AC=1;AF=0.50;AN=2;BaseQRankSum=-0.589;DP=89;Dels=0.00;FS=29.543;HRun=3;HaplotypeScore=17.2973;MQ=48.03;MQ0=0;MQRankSum=-5.007;QD=4.91;ReadPosRankSum=-1.858;SB=-0.01	GT:AD:DP:GQ:PL	0/1:64,25:89:99:467,0,2160	Fail
+gi|29165615|ref|NC_002745.2|	123269	.	G	A	445.07	SnpCluster	AC=1;AF=0.50;AN=2;BaseQRankSum=0.352;DP=90;Dels=0.00;FS=24.306;HRun=2;HaplotypeScore=0.8667;MQ=47.79;MQ0=0;MQRankSum=-4.696;QD=4.95;ReadPosRankSum=-2.310;SB=-0.01	GT:AD:DP:GQ:PL	0/1:64,26:90:99:475,0,2020	Fail
+gi|29165615|ref|NC_002745.2|	123282	.	G	A	657.53	SnpCluster	AC=1;AF=0.50;AN=2;BaseQRankSum=-2.265;DP=111;Dels=0.00;FS=8.623;HRun=2;HaplotypeScore=0.8667;MQ=46.96;MQ0=0;MQRankSum=-4.142;QD=5.92;ReadPosRankSum=-3.564;SB=-224.82	GT:AD:DP:GQ:PL	0/1:73,38:111:99:688,0,2237	Fail
+gi|29165615|ref|NC_002745.2|	123285	.	G	A	742.53	SnpCluster	AC=1;AF=0.50;AN=2;BaseQRankSum=-1.419;DP=109;Dels=0.00;FS=6.621;HRun=0;HaplotypeScore=0.0000;MQ=46.68;MQ0=0;MQRankSum=-4.188;QD=6.81;ReadPosRankSum=-3.840;SB=-296.54	GT:AD:DP:GQ:PL	0/1:70,39:109:99:773,0,2105	Fail
+gi|29165615|ref|NC_002745.2|	123288	.	A	G	669.41	SnpCluster	AC=1;AF=0.50;AN=2;BaseQRankSum=-0.594;DP=111;Dels=0.00;FS=3.958;HRun=0;HaplotypeScore=0.0000;MQ=46.43;MQ0=0;MQRankSum=-2.576;QD=6.03;ReadPosRankSum=-2.874;SB=-332.32	GT:AD:DP:GQ:PL	0/1:73,38:111:99:699,0,2163	Fail
+gi|29165615|ref|NC_002745.2|	123300	.	G	A	2106.81	PASS	AC=1;AF=0.50;AN=2;BaseQRankSum=1.338;DP=112;Dels=0.00;FS=5.192;HRun=0;HaplotypeScore=0.0000;MQ=45.49;MQ0=0;MQRankSum=2.296;QD=18.81;ReadPosRankSum=1.661;SB=-1107.47	GT:AD:DP:GQ:PL	0/1:39,73:112:99:2137,0,744	Fail
+gi|29165615|ref|NC_002745.2|	506058	.	T	G	2080.85	PASS	AC=1;AF=0.50;AN=2;BaseQRankSum=3.066;DP=112;Dels=0.00;FS=3.883;HRun=1;HaplotypeScore=0.0000;MQ=45.76;MQ0=0;MQRankSum=2.070;QD=18.58;ReadPosRankSum=1.551;SB=-1102.28	GT:AD:DP:GQ:PL	0/1:39,73:112:99:2111,0,724	Fail
+gi|29165615|ref|NC_002745.2|	506165	.	G	A	1957.05	PASS	AC=1;AF=0.50;AN=2;BaseQRankSum=1.768;DP=121;Dels=0.00;FS=0.706;HRun=2;HaplotypeScore=0.0000;MQ=44.99;MQ0=0;MQRankSum=-1.013;QD=16.17;ReadPosRankSum=-0.029;SB=-1121.11	GT:AD:DP:GQ:PL	0/1:47,74:121:99:1987,0,1000	Fail
diff --git a/vcfannotator/test/test_02.gb → test-data/test_02.gb b/vcfannotator/test/test_02.gb → test-data/test_02.gb
diff --git a/vcfannotator/test/test_02.vcf → test-data/test_02.vcf b/vcfannotator/test/test_02.vcf → test-data/test_02.vcf
diff --git a/vcfannotator/annotator.py b/vcfannotator/annotator.py
@@ -7,9 +7,8 @@
 class Annotator(object):
     """Annotator class."""
 
-    def __init__(self, gb_file=False, vcf_file=False, length=15):
+    def __init__(self, gb_file=False, vcf_file=False):
         """Initialize variables."""
-        self.length = length
         self.__annotated_features = ["CDS", "tRNA", "rRNA", "ncRNA",
                                      "misc_feature"]
         self.__gb = genbank.GenBank(gb_file)
@@ -94,6 +93,10 @@ def annotate_vcf_records(self):
                             ).replace(
                                 ' ', '[space]'
                             )
+                            if v == 'anticodon':
+                                record.INFO[k] = 'anticodon{0}'.format(
+                                    record.INFO[k]
+                                )
 
             # Determine variant type
             if record.is_indel: