other/pipeline/trunk/gc_stat.py

Code
Comments
Other
Rev Date Author Line
5709 08 Nov 19 nicklas 1 #!/usr/bin/python
5709 08 Nov 19 nicklas 2
5709 08 Nov 19 nicklas 3 """
5719 12 Nov 19 nicklas 4 $Id$
5719 12 Nov 19 nicklas 5
5709 08 Nov 19 nicklas 6 Script to create a VCF-like file with GC content and DNA strength INFO fields. This script is
5709 08 Nov 19 nicklas 7 part of a pipeline that requires a specially prepared bed-file as input (reading from stdin).
5709 08 Nov 19 nicklas 8 Columns (0-based for convenience):
5709 08 Nov 19 nicklas 9
5709 08 Nov 19 nicklas 10 0: Chromosome (printed to VCF#CHROM)
5753 22 Nov 19 nicklas 11 1: Start (=original location - 50; or 0 if original location < 50)
5709 08 Nov 19 nicklas 12 2: Stop (=original location + 50; not used)
5709 08 Nov 19 nicklas 13 3: Original location (printed to VCF#POS)
5709 08 Nov 19 nicklas 14 4: REF  (printed to VCF#REF)
5709 08 Nov 19 nicklas 15 5: ALT (printed to VCF#ALT)
5709 08 Nov 19 nicklas 16 6: SEQUENCE between start and stop (used for calculating GC and DNA strength, which are printed to VCF#INFO)
5709 08 Nov 19 nicklas 17
5709 08 Nov 19 nicklas 18 The VCF file that is produced must be compressed (bgzip) and indexed (tabix) before it can be used by vcfanno 
5709 08 Nov 19 nicklas 19 to annotate the orginal VCF file (from VarDict) with the GC content and DNA strength. A .toml file like 
5709 08 Nov 19 nicklas 20 this is needed:
5709 08 Nov 19 nicklas 21
5709 08 Nov 19 nicklas 22 [[annotation]]
5709 08 Nov 19 nicklas 23 file="/path/to/vcf-file-produced-by-this-script.gz"
5709 08 Nov 19 nicklas 24 fields=["GC_CONT","GC_CONT_LOCAL","STRENGTH","STRENGTH_LOCAL"]
5709 08 Nov 19 nicklas 25 names=["GC_CONT","GC_CONT_LOCAL","STRENGTH","STRENGTH_LOCAL"]
5709 08 Nov 19 nicklas 26 ops=["self","self","self","self"]
5709 08 Nov 19 nicklas 27
5709 08 Nov 19 nicklas 28 """
5709 08 Nov 19 nicklas 29
5709 08 Nov 19 nicklas 30 import sys
5709 08 Nov 19 nicklas 31
5709 08 Nov 19 nicklas 32 NORM_LEN = 50
5709 08 Nov 19 nicklas 33 LOCAL_LEN = 6
5709 08 Nov 19 nicklas 34
5709 08 Nov 19 nicklas 35 # Column indexes from the BED file we are reading
5709 08 Nov 19 nicklas 36 COL_CHR = 0
5753 22 Nov 19 nicklas 37 COL_START = 1
5753 22 Nov 19 nicklas 38 COL_END = 2
5709 08 Nov 19 nicklas 39 COL_POS = 3
5709 08 Nov 19 nicklas 40 COL_REF = 4
5709 08 Nov 19 nicklas 41 COL_ALT = 5
5709 08 Nov 19 nicklas 42 COL_SEQ = 6
5709 08 Nov 19 nicklas 43
5709 08 Nov 19 nicklas 44 # GC content and DNA strength values
5747 21 Nov 19 nicklas 45 GC_CONT_VALUES = { "A": 0.0, "C": 1.0, "G": 1.0, "N": 0.5, "T": 0.0 }
5709 08 Nov 19 nicklas 46 DNA_STRENGTH_VALUES = {
5747 21 Nov 19 nicklas 47     "AA": 5.0, "AC": 10.0, "AG": 8.0, "AN": 7.5, "AT": 7.0,
5747 21 Nov 19 nicklas 48     "CA": 7.0, "CC": 11.0, "CG": 10.0, "CN": 9.0, "CT": 8.0, 
5747 21 Nov 19 nicklas 49     "GA": 8.0, "GC": 13.0, "GG": 11.0, "GN": 10.5, "GT": 10.0,
5747 21 Nov 19 nicklas 50     "NA": 6.0, "NC": 10.5, "NG": 9.0, "NN": 8.25, "NT": 7.5,
5747 21 Nov 19 nicklas 51     "TA": 4.0, "TC": 8.0, "TG": 7.0, "TN": 6.0, "TT": 5.0
5709 08 Nov 19 nicklas 52   }
5709 08 Nov 19 nicklas 53
5709 08 Nov 19 nicklas 54 def getSequenceInfo(seqIn):
5709 08 Nov 19 nicklas 55     """
5709 08 Nov 19 nicklas 56     Returns GC-content and DNA-strength of the seqIn
5709 08 Nov 19 nicklas 57     """
5709 08 Nov 19 nicklas 58     seqIn = seqIn.upper()
5709 08 Nov 19 nicklas 59     gcContSum = 0
5709 08 Nov 19 nicklas 60     strengthSum = 0
5709 08 Nov 19 nicklas 61     for seqPos in range(0, len(seqIn)):
5709 08 Nov 19 nicklas 62         if seqPos > 1:
5709 08 Nov 19 nicklas 63             strengthSum += DNA_STRENGTH_VALUES[seqIn[seqPos - 2:seqPos]]
5709 08 Nov 19 nicklas 64         gcContSum += GC_CONT_VALUES[seqIn[seqPos]]
5709 08 Nov 19 nicklas 65
5709 08 Nov 19 nicklas 66     gcCont = gcContSum / len(seqIn)
5709 08 Nov 19 nicklas 67     dnaStrength = strengthSum / len(seqIn)
5709 08 Nov 19 nicklas 68     return gcCont, dnaStrength
5709 08 Nov 19 nicklas 69
5709 08 Nov 19 nicklas 70
5709 08 Nov 19 nicklas 71 header = '''##fileformat=VCFv4.2
5709 08 Nov 19 nicklas 72 ##INFO=<ID=GC_CONT,Number=1,Type=Float,Description="Proportion of bases that are GC within {norm} bases up- and down-stream of position">
5709 08 Nov 19 nicklas 73 ##INFO=<ID=GC_CONT_LOCAL,Number=1,Type=Float,Description="Proportion of bases that are GC within {loc} bases up and downstream of position">
5709 08 Nov 19 nicklas 74 ##INFO=<ID=STRENGTH,Number=1,Type=Float,Description="DNA strength in the region within {norm} bases up- and down-stream of position">
5709 08 Nov 19 nicklas 75 ##INFO=<ID=STRENGTH_LOCAL,Number=1,Type=Float,Description="DNA strength in the region within {loc} bases up- and down-stream of position">
5709 08 Nov 19 nicklas 76 #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO'''.format(norm=NORM_LEN, loc=LOCAL_LEN)
5709 08 Nov 19 nicklas 77
5709 08 Nov 19 nicklas 78 print(header)
5709 08 Nov 19 nicklas 79
5709 08 Nov 19 nicklas 80 for line in sys.stdin:
5709 08 Nov 19 nicklas 81     cols = line.rstrip().split('\t')
5709 08 Nov 19 nicklas 82     seq = cols[COL_SEQ]
5753 22 Nov 19 nicklas 83     posInSeq = int(cols[COL_POS])-int(cols[COL_START])
5753 22 Nov 19 nicklas 84     seqLocal = cols[COL_SEQ][max(posInSeq-LOCAL_LEN, 0):posInSeq+LOCAL_LEN]
5709 08 Nov 19 nicklas 85     gcNorm, strengthNorm = getSequenceInfo(seq)
5709 08 Nov 19 nicklas 86     gcLocal, strengthLocal = getSequenceInfo(seqLocal)
5709 08 Nov 19 nicklas 87     cout = (cols[COL_CHR], cols[COL_POS], '.', cols[COL_REF], cols[COL_ALT], '.', '.', gcNorm, gcLocal, strengthNorm, strengthLocal)
5709 08 Nov 19 nicklas 88     print('{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\tGC_CONT={7};GC_CONT_LOCAL={8};STRENGTH={9};STRENGTH_LOCAL={10}'.format(*cout))

Generated on Tue Nov 21 03:40:14 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)