#!/usr/bin/env python

import logging, sys, optparse, re, os, datetime, gzip, urllib2, subprocess
from collections import defaultdict
from os.path import join, basename, dirname, isfile, abspath
from datetime import date, datetime, timedelta

# === COMMAND LINE INTERFACE, OPTIONS AND HELP ===
parser = optparse.OptionParser("""usage: %prog [options] filename - check and convert clinVar table to four bed files, split into CNV and shorter mutations, for both hg19 and hg38 and convert all to bigBed

Output goes into the current dir

input file is ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz
""") 

parser.add_option("-d", "--debug", dest="debug", action="store_true", help="show debug messages") 
parser.add_option("-a", "--auto", dest="auto", action="store_true", help="download the file from NCBI into the current dir and convert to bigBed")
parser.add_option("", "--maxDiff", dest="maxDiff", action="store", type="float", help="look for last month's download file in current dir and accept this much difference, expressed as a ratio. Can only be used with --auto.")
(options, args) = parser.parse_args()

if options.debug:
    logging.basicConfig(level=logging.DEBUG)
else:
    logging.basicConfig(level=logging.INFO)

clinvarExpectHeaders = "#AlleleID	Type	Name	GeneID	GeneSymbol	ClinicalSignificance	RS# (dbSNP)	nsv (dbVar)	RCVaccession	TestedInGTR	PhenotypeIDs	Origin	Assembly	Chromosome	Start	Stop	Cytogenetic	ReviewStatus	HGVS(c.)	HGVS(p.)	NumberSubmitters	LastEvaluated	Guidelines	OtherIDs	VariantID	ReferenceAllele	AlternateAllele	SubmitterCategories	ChromosomeAccession\n"
# ==== FUNCTIONs =====

def shortenName(name):
    "make the clinvar name shorter "
    cnvMatch = cnvRe.match(name)
    hgvsMatch = hgvsRe.match(name)
    if cnvMatch != None:
        # some names include the original assembly as a prefix
        band, copy = cnvMatch.groups()
        if copy == None:
            copy = ""
        name = band+copy
    #elif hgvsMatch != None:
        # some names include the whole hgvs desc, we strip the gene name + acc
        #name = hgvsMatch.groups()[2]
    elif name.startswith("GRCh"):
        assert False # name that starts with GRCh but does not match our regex should not happen
    elif "." in name:
        # many names look like NM_000274.3(OAT):c.550_552delGCT (p.Ala184del)
        # so we remove the gene and the protein description
        if ":" in name:
            name = name.split(":")[1]
        name = name.split(" ")[0]

        delMatch = delRe.match(name)
        if delMatch!=None:
            return name[-3:]

        posMatch = posRe.match(name)
        if posMatch!=None:
            dup, dele, ins, pos = posMatch.groups()[-4:]
            if dup==None:
                dup= ""
            if dele==None:
                dele = ""
            if ins==None:
                ins = ""
            name = dup+dele+ins+pos
    return name

def lastMonth(d,x=1):
    """ returns same day as this month, but last month, e.g. for 2013-01-01 
    return 2012-12-01 
    """

    days_of_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    newmonth = ((( d.month - 1) - x ) % 12 ) + 1
    newyear  = d.year + ((( d.month - 1) - x ) // 12 ) 
    if d.day > days_of_month[newmonth-1]:
      newday = days_of_month[newmonth-1]
    else:
      newday = d.day
    return datetime( newyear, newmonth, newday)

# ----------- MAIN --------------
if args==[] and not options.auto:
    parser.print_help()
    exit(1)

outSuffix = ""
# download from NCBI
if options.auto:
    today = date.today().isoformat()
    outSuffix = "%s." % today
    url = "http://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz"
    outFname = "archive/variant_summary.%s.txt.gz" % today
    logging.info("Downlading %s to %s" % (url, outFname))
    open(outFname, "w").write(urllib2.urlopen(url).read())
    filename = outFname
else:
    filename = args[0]

# e.g. GRCh38/hg38 1p36.33-36.31(chr1:1181847-5507243)x1
# GRCh38/hg38 9p24.3-q34.3(chr9:204193-138179445)
cnvRe = re.compile(r'GRCh[0-9]*/hg[0-9]* ([XY0-9pq.-]+)\([XYa-z_:0-9-]+\)(x[0-9]+)?')

hgvsRe = re.compile(r'(N[RPCGM]_[0-9.]+)(\([a-zA-Z0-9_-]+\))?:([+0-9c.ACTG>a-z_()-]+)')

# e.g.(13359800_13371688)_(13371779_13505931)del
delRe = re.compile(r'g.\([?0-9()_-]+(del|dup|ins)')

# e.g.
# c.80_83delGGATinsTGCTGTAAACTGTAACTGTAAA
# c.80A>T
# c.740+356C>T
# c.3839_3855del17
posRe = re.compile(r'[mgcpdrn]\.[*0-9_+-]+(dup)?(del)?(ins)?([>ACTG0-9]+)')

ifh = gzip.open(filename)

# check the header
line1 = ifh.readline()
if (line1 != clinvarExpectHeaders): # test if clinvar has changed their header line again
    logging.error("ClinVar has changed their header line again")
    logging.error("current header line: %s" % line1)
    logging.error("expected header line: %s" % clinvarExpectHeaders)
    raise Exception("code needs fixing")

# open out files
hg19Bed = open("archive/clinvarMain.hg19.%sbed" % outSuffix, "w")
hg38Bed = open("archive/clinvarMain.hg38.%sbed" % outSuffix, "w")
hg19BedCnv = open("archive/clinvarCnv.hg19.%sbed" % outSuffix, "w")
hg38BedCnv = open("archive/clinvarCnv.hg38.%sbed" % outSuffix, "w")

longCount = 0
noAssCount = 0

# convert lines
for line in ifh:
    line = line.replace(", ", ",").replace(";", ",") # replace, bigBed conventions
    line = line.rstrip("\n")

    alleleId, allType, name, geneId, geneSymbol, clinSign, snpAcc, dbVarAcc, \
      irvcAcc, inGtr, phenotypeIds, origin, assembly, chrom, start, end, cytogenetic, \
      reviewStatus, hgvsCod, hgvsProt, numberSubmitters, lastEval, guidelines, otherIds, varId, \
      refAll, varAll, submCategories, chromAcc = line.split("\t")

    if chrom=="" or assembly=="" or assembly=="NCBI36":
        noAssCount += 1
        continue

    if chrom=="Un" and assembly=="GRCh38":
        print "wrong chrUn chrom on assembly hg38. Skipping %s" % irvcAcc
        continue
    chrom = "chr"+chrom
    if chrom=="chrMT": # why does NCBI use chrMT but we use chrM ?
        chrom = "chrM"

    shortName = shortenName(name)
    if len(shortName)>17:
        shortName = shortName[:15]+"..."
        longCount+=1

    if start=="" or end=="":
        print "undefined start or end coordinate. record %s"% irvcAcc
        continue
        
    # sometimes they seem to inverse their coordinates
    if int(start)>int(end):
        start, end = end, start

    # reorder columns for bigBed
    start = str(int(start)-1) # convert to zero-based, half-open
    score = "0"
    strand = "."
    thickStart = start
    thickEnd = end
    itemRgb = "0"
    blockCount = "1"
    blockSizes = str(int(end)-int(start))
    blockStarts = "0"

    if dbVarAcc.startswith("nsv"):
        if "gain" in allType or "duplication" in allType:
            itemRgb = "0,0,255"
        if "deletion" in allType or "loss" in allType:
            itemRgb = "255,0,0"
    else:
        isPat, isBen = False, False
        if "pathogenic" in clinSign.lower():
            isPat = True
        if "benign" in clinSign.lower():
            isBen = True
        if (isPat==True and isBen==True) or (isPat==False and isBen==False):
            itemRgb = "128,128,128"
        elif isPat==True:
            itemRgb = "0,0,0"
        else: # is benign
            itemRgb = "200,200,200"

    row = [chrom, start, end, shortName, score, strand, thickStart, thickEnd, itemRgb, \
        blockCount, blockSizes, blockStarts,
        name, allType, geneId, geneSymbol, clinSign,
        snpAcc, dbVarAcc, irvcAcc, inGtr, phenotypeIds, origin, assembly, cytogenetic,
        reviewStatus, hgvsCod, hgvsProt, numberSubmitters, lastEval, guidelines, otherIds]

    # replace clinvar's placeholders with real empty fields
    newRow = []
    for x in row:
        if x in ["-1", "-"]:
            newRow.append("")
        else:
             newRow.append(x)
    row = newRow

    if assembly=="GRCh37":
        if chrom=="chrM":
            continue # we don't have the same MT chrom as NCBI, skip these
        ofh = hg19Bed
        if dbVarAcc.startswith("nsv"):
            ofh = hg19BedCnv

    elif assembly=="GRCh38":
        ofh = hg38Bed
        if dbVarAcc.startswith("nsv"):
            ofh = hg38BedCnv
    else:
        noAssCount +=1
        
    ofh.write("\t".join(row))
    ofh.write("\n")

hg19Bed.close()
hg38Bed.close()
hg19BedCnv.close()
hg38BedCnv.close()

logging.info("%d lines with feature name that was too long, shortened them" % longCount)
logging.info("%d lines without, with old assembly or no chrom, skipped" % noAssCount)

fnames = [hg19Bed.name, hg38Bed.name, hg19BedCnv.name, hg38BedCnv.name]

# sort the files
for fname in fnames:
    cmd = "bedSort %s %s" % (fname, fname) 
    logging.debug(cmd)
    assert(os.system(cmd)==0)

# check new bed files against last months's bed files
if options.maxDiff:
    for fname in fnames:
        # go back one month from current date
        previousDate = lastMonth(date.today())
        previousDate = previousDate.date().isoformat()
        previousFname = fname.replace(today, previousDate)
        todayCount = sum(1 for line in open(fname))
        previousCount = sum(1 for line in open(previousFname))
        percDiff = (todayCount-previousCount)/float(todayCount)
        if percDiff < 0 and options.maxDiff:
            raise Exception("line count change negative: %s %d, %s %d. Run without --maxDiff to ignore this error. Or run the cronjob script doUpdate.sh with -nocheck to ignore this error." % \
                (previousFname, previousCount, fname, todayCount))
        if percDiff > options.maxDiff:
            raise Exception("line count increased too much: yesterday %d , today %d. Check if the input file is OK. If yes, remove the --maxDiff option from the clinVarToBed, call it only with the --auto option. Or run the cronjob script doUpdate.sh with -nocheck to ignore this error." % (previousCount, todayCount))
        # see RM 14350 do more than just linecount
        # see http://stackoverflow.com/questions/1566461/how-to-count-differences-between-two-files-on-linux
        cmd = "diff %s %s --speed-large-files | grep '^[\>\<]'  | wc -l" % (fname, previousFname)
        diffCount = int(subprocess.check_output(cmd, shell=True, stderr=subprocess.STDOUT))
        logging.info("%s: %d lines changed" % (fname, diffCount))
        percDiff = float(diffCount)/todayCount
        if percDiff > options.maxDiff:
            logging.info(cmd)
            diffFname =  abspath(join("archive", "lastDiff.txt"))
            cmd = "diff %s %s --speed-large-files > %s" % (fname, previousFname, diffFname)
            os.system(cmd) # ret code is not 0, as there are differences, so no need to check
            raise Exception("too many differences: %d out of %d lines. Look at the differences n the file %s. If you are OK with them, run '/hive/data/outside/otto/clinvar/doUpdate.sh -nocheck'" % (diffCount, todayCount, diffFname))

# convert to bigBed without the date in the fname, like clinvarMain.hg19.bb
for fname in fnames:
    parts = fname.split(".")
    if len(parts) == 4: # clinvarMain.hg19.2011-11-11.bed
        base, db, date, ext = parts
    else:
        base, db, ext = parts
        
    outFname = "%s.%s.bb" % (base, db)
    outFname = basename(outFname) # strip the archive/ part, put into current dir
    cmd = "bedToBigBed -tab -type=bed9+ -as=clinvar.as %s /hive/data/genomes/%s/chrom.sizes %s" % (fname, db, outFname)
    logging.debug(cmd)
    assert(os.system(cmd)==0)
