3683a482f97f26e6c08290f6c8acccb764d16eac lrnassar Thu Mar 2 17:58:16 2023 -0800 Removing the use of refseq select because the annotations are too sparse. Refs #30320 diff --git src/hg/makeDb/doc/evaSnp4.txt src/hg/makeDb/doc/evaSnp4.txt index aaae0a9..0bdde85 100644 --- src/hg/makeDb/doc/evaSnp4.txt +++ src/hg/makeDb/doc/evaSnp4.txt @@ -1,575 +1,575 @@ # Track for EVA snp release 4 - https://www.ebi.ac.uk/eva/?RS-Release&releaseVersion=4 # Tracks built by Lou on 3/2/2023 # Track was built for the following 36 assemblies: panTro5, melGal5, equCab3, rn5, canFam3, oreNil2, # calJac3, bosTau7, taeGut2, bosTau6, galGal4, rn7, danRer10, danRer11, gorGor3, oviAri3, ponAbe3, # rheMac10, galGal6, danRer7, bosTau9, galGal5, oviAri4, mm39, mm10, anoGam3, rheMac8, gorGor4, felCat9, # macFas5, susScr3, chlSab2, ci3, dm6, susScr11, felCat3 # The GCA accession on the eva release by accession list (https://ftp.ebi.ac.uk/pub/databases/eva/rs_releases/release_3/by_assembly/) # were compared to all native assemblies we have. In total there are 894617246 variants, # All assemblies were passed by the python pipeline described below. There were 3 assemblies # that requires post-processing, bosTau9 and galGal5 I had to remove their chrMT as we had none, # and rheMac8 I had to change the name of the chrMT to one that matches us. See RM #30320 ##################################### #### Beginning of python3 code ###### ##################################### #NEW import re, os, subprocess, sys def findAllAccFromDescPages(): """Search all gbdb description pages and pull out GCA/GCF accessions""" listOfHits = [] assemblyAlias = {} allDescPages = bash("ls /gbdb/*/html/description.html").rstrip().split("\n") regex = re.compile('GC[AF]_[0-9]*\.[0-9]*') for page in allDescPages: try: hits = bash("cat "+page+" | grep -i 'GCA\|GCF'").rstrip().split("\n") except: hits = [] if hits != []: ucscName = page.split("/")[2] hits = str(hits) assemblyAlias[ucscName] = {} try: accession = regex.findall(hits)[0] if accession.startswith("GCA"): assemblyAlias[ucscName]['GCA'] = accession elif accession.startswith("GCF"): assemblyAlias[ucscName]['GCF'] = accession except: print("The following page could not be matched with the accession format",page,hits) return(assemblyAlias) def buildGCFdicFromGenbankFile(): genbankSummaryDic = {} regexGenbankGCA = re.compile('GCA_[0-9]*\.[0-9]*') regexGenbankGCF = re.compile('GCF_[0-9]*\.[0-9]*') summaryFile = open('/hive/users/lrnassar/evaSnp4/assembly_summary_genbank.txt','r') for line in summaryFile: if "GCF_" in line: GCAacc = regexGenbankGCA.findall(line)[0] GCFacc = regexGenbankGCF.findall(line)[0] genbankSummaryDic[GCFacc] = GCAacc summaryFile.close() return(genbankSummaryDic) def findGCAaccForUCSCgcfDbs(assemblyAlias,genbankSummaryDic): """Try to find matching GCA accession to UCSC dbs based on genbank file or curling web page""" n=0 regexGenbankGCA = re.compile('GCA_[0-9]*\.[0-9]*') for key, value in assemblyAlias.items(): if 'GCF' in value: if value['GCF'] in genbankSummaryDic.keys(): assemblyAlias[key]['GCA'] = genbankSummaryDic[value['GCF']] n+=1 else: #Try to find it online at NCBI ncbiGCF = value['GCF'] ncbiWgetPage = bash("curl https://www.ncbi.nlm.nih.gov/assembly/"+ncbiGCF+"/").rstrip().split("\n") for line in ncbiWgetPage: if "GenBank assembly accession" in line and "GCA_" in line: for subLine in line.split("/dd"): if "GenBank assembly accession" in subLine: GCAacc = regexGenbankGCA.findall(subLine)[0] assemblyAlias[key]['GCA'] = GCAacc n+=1 print("Total number of UCSC GCFs matched with GCAs: ",n) return(assemblyAlias) def sanitizeListToOnlyUCSCdbsWithGCA(assemblyAlias): n = 0 m=0 ucscGCAlist = {} for key, value in assemblyAlias.items(): if 'GCA' in value: if value['GCA'] != 'GCA_000': ucscGCAlist[key] = value['GCA'] n+=1 else: m+=1 print("Total number of UCSC databases already with a GCA:",n) print("Total number of UCSC databases with GCF that will need GCA lookup:",m) return(ucscGCAlist) def buildFinalDicOfUCSCdbsWithEvaSnp4Data(ucscGCAlist): ucscEvaSnp4Dbs = {} n=0 regexEvaSnp4GCA = re.compile('GCA_[0-9]*\.[0-9]*') evaRelease4page = bash("curl https://ftp.ebi.ac.uk/pub/databases/eva/rs_releases/release_4/by_assembly/").rstrip().split("\n") activeRRassemblies = bash("""hgsql -h genome-centdb -e "select name from dbDb where active='1'" hgcentral""") for key in ucscGCAlist: if ucscGCAlist[key] in str(evaRelease4page): for line in evaRelease4page: if ucscGCAlist[key] in line: if key in activeRRassemblies: fullGCA = regexEvaSnp4GCA.findall(line)[0] ucscEvaSnp4Dbs[key] = fullGCA n+=1 ucscEvaSnp4Dbs['mm10']='GCA_000001635.6' #Manually add mm10 which is weird due to patches print("Total number of assemblies with EvaSnp4 data: ",n) return(ucscEvaSnp4Dbs) def checkDirMakeDir(workDir): if not os.path.isdir(workDir): os.mkdir(workDir) def bash(cmd): """Run the cmd in bash subprocess""" try: rawBashOutput = subprocess.run(cmd, check=True, shell=True,\ stdout=subprocess.PIPE, universal_newlines=True, stderr=subprocess.STDOUT) bashStdoutt = rawBashOutput.stdout except subprocess.CalledProcessError as e: raise RuntimeError("command '{}' return with error (code {}): {}".format(e.cmd, e.returncode, e.output)) return(bashStdoutt) def wgetFile(url, workDir): bash("wget -N %s -P %s" % (url, workDir)) fileName = workDir+url.split("/")[-1] return(fileName) def checkDupsAndUnzip(fileName,workDir): """Check for duplicate lines (common in EVA), and unzip file""" lineCount = bash("zcat %s | wc -l" % (fileName)).rstrip() uniqCount = bash("zcat %s | uniq | wc -l" % (fileName)).rstrip() currentFileName = workDir+"evaSnps.vcf" if lineCount != uniqCount: print("There are duplicated entries.\nOriginal file:%s\nUniq file:%s\nCreating new file." % (lineCount, uniqCount)) bash("zcat %s | uniq > %s" % (fileName,currentFileName)) else: bash("zcat %s > %s" % (fileName,currentFileName)) return(currentFileName) def sanitizeChromNames(currentFileName,dbs,workDir): """Swaps the chrom IDs with the accession names which are indexed in chromAlias""" inputFile = open(currentFileName,'r') outputFile = open(currentFileName+".fixed","w") IDtoAccessionDic = {} for line in inputFile: if line.startswith("##contig"): ID = line.split('=')[2].split(',')[0] accession = line.split('=')[3].split('"')[1] IDtoAccessionDic[ID] = accession outputFile.write(line) elif line.startswith("#"): outputFile.write(line) else: line = line.split("\t") line[0] = IDtoAccessionDic[line[0]] line = "\t".join(line) outputFile.write(line) inputFile.close() outputFile.close() #Replacing the same name because chromAlias should be updated soon and this function will be removed bash("mv %s.fixed %s" % (currentFileName,currentFileName)) currentFileName = convertChromToUcsc(currentFileName,dbs,workDir) print("Chroms converted") def addChrToChromNumbers(fileName,currentFileName,workDir): bash("zcat %s | uniq > %s" % (fileName,currentFileName)) with open(workDir+"evaSnps.ucscChroms.vcf",'w') as newChromVcf: with open(currentFileName,"r") as oldChromVcf: for line in oldChromVcf: if line.startswith("#"): newChromVcf.write(line) elif line.startswith("AC024175.3"): #Exception for zebrafish newChromVcf.write(line.replace("AC024175.3","chrM")) else: newChromVcf.write("chr"+line) currentFileName = workDir+"evaSnps.ucscChroms.vcf" return(currentFileName) def convertChromToUcsc(currentFileName,dbs,workDir): """Convert to UCSC-style chroms""" bash("chromToUcsc --get %s" % (dbs)) bash("chromToUcsc -i %s -o %sevaSnps.ucscChroms.vcf -a %s" % (currentFileName,workDir,dbs+".chromAlias.tsv")) bash("rm %s" % (dbs+".chromAlias.tsv")) currentFileName = workDir+"evaSnps.ucscChroms.vcf" return(currentFileName) def convertVcfToBed(currentFileName,workDir): """Convert to bed keeping only VC and SID fields""" bash("bgzip %s" % (currentFileName)) bash("tabix -p vcf %s" % (currentFileName+".gz")) bash("vcfToBed %s %sevaSnp -fields=VC,SID" % (currentFileName+".gz",workDir)) currentFileName = workDir+"evaSnp.bed" return(currentFileName) def buildChromSizesDic(workDir,dbs): """Create a dictionary the chromSizes""" chromSizes = bash("fetchChromSizes %s" % (dbs)).split("\n") chromDic = {} for chrom in chromSizes[1:-1]: chrom = chrom.split("\t") chromDic[chrom[0]] = chrom[1] return(chromDic) def splitChromsAndRunHgVai(workDir,dbs): """Split all the chroms in tenths in order to be able to run hgVai without running out of memory""" chromSizes = bash("fetchChromSizes %s" % (dbs)).split("\n") inputFile = workDir+"evaSnps.ucscChroms.vcf.gz" outputFile = workDir+"evaSnps"+dbs+"VaiResults.vep" n=0 allTables = bash("\hgsql -e \"show tables\" "+dbs+"") if dbs in ['galGal6','oviAri4','ponAbe3']: #Special exception to use refGene because ncbiRefSeq has incorrect protein sequence in Link table see #29262 geneTableToUse = "refGene" - elif "ncbiRefSeqSelect" in allTables: - geneTableToUse = "ncbiRefSeqSelect" +# elif "ncbiRefSeqSelect" in allTables: #Removing refseq select as the annotations are too sparse +# geneTableToUse = "ncbiRefSeqSelect" elif "ncbiRefSeqCurated" in allTables: geneTableToUse = "ncbiRefSeqCurated" elif "ensGene" in allTables: geneTableToUse = "ensGene" elif "refGene" in allTables: geneTableToUse = "refGene" elif "ncbiGene" in allTables: geneTableToUse = "ncbiGene" else: print(dbs) sys.exit("Could not find any tables to use for the following database: "+dbs) chromDic = buildChromSizesDic(workDir,dbs) #For function below, only bother with the chromosomes in the VCF chromsInVcf = bash("zcat "+inputFile+" | grep -v '^#' | cut -f1 | uniq").rstrip().split("\n") for chrom in chromsInVcf: chromName = chrom chromSize = chromDic[chrom] n+=1 if (n % 100) == 0: print("Have currently run "+str(n)+"/"+str(len(chromsInVcf))+" chroms. Currently running "+chromName) if int(chromSize) < 1000000: #Don't worry about splitting if chrom is less than 1M bp bash("vai.pl --variantLimit=200000000 --hgVai=/usr/local/apache/cgi-bin/hgVai --position=%s:0-%s --geneTrack=%s --hgvsG=off --hgvsCN=off \ --hgvsP=on %s %s >> %s" % (chromName,chromSize,geneTableToUse,dbs,inputFile,outputFile)) else: breakPoints = [.2,.4,.6,.8,1] for part in breakPoints: startCoord = round(int(chromSize)*(part-.2)) endCoord = round(int(chromSize)*part) if part == breakPoints[0]: #Special exception for start bash("vai.pl --variantLimit=200000000 --hgVai=/usr/local/apache/cgi-bin/hgVai --position=%s:0-%s --geneTrack=%s --hgvsG=off --hgvsCN=off \ --hgvsP=on %s %s >> %s" % (chromName,endCoord,geneTableToUse,dbs,inputFile,outputFile)) else: bash("vai.pl --variantLimit=200000000 --hgVai=/usr/local/apache/cgi-bin/hgVai --position=%s:%s-%s --geneTrack=%s --hgvsG=off --hgvsCN=off \ --hgvsP=on %s %s >> %s" % (chromName,startCoord,endCoord,geneTableToUse,dbs,inputFile,outputFile)) return(outputFile) def splitVepFileByChrom(hgVaiResultsFile,workDir): """Split the vep file to one file per chrom in order to index entire chrom into a dic and quickly lookup the rsIDs""" splitChromDir = workDir+"splitChroms/" checkDirMakeDir(splitChromDir) vepFile = open(hgVaiResultsFile,"r") for line in vepFile: if line.startswith("#") or line.startswith("Uploaded"): continue else: #Need to compute itemRGB, aaChange, and ucscClass lineContents = line.rstrip().split("\t") chrom = lineContents[1].split(":")[0] outFile = open(splitChromDir+"evaSnpVaiResults."+chrom+".vep","a") outFile.write(line) outFile.close() vepFile.close() def fetchUcscClassAndAaChange(name,alt,nameAlt,chrom,currentChrom,currentVepDic,tryAttempts,workDir): """Build dictionary out of vep results with rsID as index, populate aaChange and ucscClassList""" prevName = nameAlt #Check if chrom has changed and index current chrom if chrom != currentChrom: currentVepFile = open(workDir+"splitChroms/evaSnpVaiResults."+chrom+".vep","r") currentVepDic = {} for entry in currentVepFile: entry = entry.rstrip().split("\t") rsID = entry[0] rsIDAlt = entry[0]+entry[2] if rsIDAlt not in currentVepDic.keys(): currentVepDic[rsIDAlt] = {} currentVepDic[rsIDAlt]['ucscClassList'] = [] currentVepDic[rsIDAlt]['ucscClassList'].append(entry[6]) currentVepDic[rsIDAlt]['aaChange'] = [] if entry[13] != '-' and 'p' in entry[13] and '?' not in entry[13]: currentVepDic[rsIDAlt]['aaChange'].append(entry[13].split(";")[0].split("=")[1].split(":")[1]) else: if entry[6] not in currentVepDic[rsIDAlt]['ucscClassList']: currentVepDic[rsIDAlt]['ucscClassList'].append(entry[6]) if entry[13] != '-' and 'p' in entry[13] and '?' not in entry[13]: if entry[13].split(";")[0].split("=")[1] not in currentVepDic[rsIDAlt]['aaChange']: currentVepDic[rsIDAlt]['aaChange'].append(entry[13].split(";")[0].split("=")[1].split(":")[1]) if rsID not in currentVepDic.keys(): currentVepDic[rsID] = {} currentVepDic[rsID]['ucscClassList'] = [] currentVepDic[rsID]['ucscClassList'].append(entry[6]) currentVepDic[rsID]['aaChange'] = [] if entry[13] != '-' and 'p' in entry[13] and '?' not in entry[13]: currentVepDic[rsID]['aaChange'].append(entry[13].split(";")[0].split("=")[1].split(":")[1]) else: if entry[6] not in currentVepDic[rsID]['ucscClassList']: currentVepDic[rsID]['ucscClassList'].append(entry[6]) if entry[13] != '-' and 'p' in entry[13] and '?' not in entry[13]: if entry[13].split(";")[0].split("=")[1] not in currentVepDic[rsID]['aaChange']: currentVepDic[rsID]['aaChange'].append(entry[13].split(";")[0].split("=")[1].split(":")[1]) currentChrom=chrom currentVepFile.close() try: ucscClass = ",".join(currentVepDic[nameAlt]['ucscClassList']) aaChange = " ".join(currentVepDic[nameAlt]['aaChange']) except: tryAttempts+=1 ucscClass = ",".join(currentVepDic[name]['ucscClassList']) aaChange = " ".join(currentVepDic[name]['aaChange']) return(prevName,currentChrom,currentVepDic,ucscClass,aaChange,tryAttempts) def assignRGBcolorByConsequence(currentVepDic,name,alt,nameAlt): """Assign color based on most serious hgVai predicted consequence""" redVars = ['exon_loss_variant','frameshift_variant','inframe_deletion','inframe_insertion','initiator_codon_variant', 'missense_variant','splice_acceptor_variant','splice_donor_variant','splice_region_variant','stop_gained', 'stop_lost','coding_sequence_variant','transcript_ablation'] greenVars = ['synonymous_variant','stop_retained_variant'] blueVars = ['5_prime_UTR_variant','3_prime_UTR_variant','complex_transcript_variant', 'non_coding_transcript_exon_variant'] blackVars = ['upstream_gene_variant','downstream_gene_variant','intron_variant','intergenic_variant', 'NMD_transcript_variant','no_sequence_alteration'] try: if any(item in currentVepDic[nameAlt]['ucscClassList'] for item in redVars): itemRgb = '255,0,0' elif any(item in currentVepDic[nameAlt]['ucscClassList'] for item in greenVars): itemRgb = '0,128,0' elif any(item in currentVepDic[nameAlt]['ucscClassList'] for item in blueVars): itemRgb = '0,0,255' elif any(item in currentVepDic[nameAlt]['ucscClassList'] for item in blackVars): itemRgb = '0,0,0' else: sys.exit("One of the following consequences not found in color list, please add: "+str(currentVepDic[name]['ucscClassList'])) except: if any(item in currentVepDic[name]['ucscClassList'] for item in redVars): itemRgb = '255,0,0' elif any(item in currentVepDic[name]['ucscClassList'] for item in greenVars): itemRgb = '0,128,0' elif any(item in currentVepDic[name]['ucscClassList'] for item in blueVars): itemRgb = '0,0,255' elif any(item in currentVepDic[name]['ucscClassList'] for item in blackVars): itemRgb = '0,0,0' else: sys.exit("One of the following consequences not found in color list, please add: "+str(currentVepDic[name]['ucscClassList'])) return(itemRgb) def convertSOnumberToName(varClass): soDic={"SO:0001483": 'substitution', "SO:0000667": 'insertion', "SO:0000159": 'deletion', "SO:0001059": 'sequence alteration',"SO:1000032": 'delins',"SO:0002007": 'MNV', "SO:0000705": 'tandem_repeat'} varClass = soDic[varClass] return(varClass) def buildDuplicateRsIDset(currentFileName): """Look for multi-allelic variants to then try and more accurately match those hgVai results""" duplicateSet = set() inputFile = open(currentFileName,"r") prevName = False for line in inputFile: if line.startswith("#") or line.startswith("Uploaded"): continue else: lineContents = line.rstrip().split("\t") name = lineContents[3] if prevName == name: duplicateSet.add(name) prevName = name print("Total number of multi-allelic variants: "+str(len(duplicateSet))) return(duplicateSet) def createFinalBedWithVepAnnotations(currentFileName,workDir): """Parse bed file and output final file without FILTER field, and with SO name, item RGB, and aaChange/ucscClass from hgVai""" inputFile = open(currentFileName,"r") outputFile = open(workDir+"evaSnp.final.bed","w") prevName = "" currentChrom = "" currentVepDic = {} n=0 tryAttempts=0 duplicateSet = buildDuplicateRsIDset(currentFileName) for line in inputFile: n+=1 if n%5000000 == 0: print("Current line count: "+str(n)) if line.startswith("#") or line.startswith("Uploaded"): continue else: #Desired schema chrom chromStart chromEnd name score strand thickStart thickEnd itemRgb ref alt varClass submitters ucscClass aaChange #Need to compute itemRGB, aaChange, and ucscClass. Also need to change SO terms lineContents = line.rstrip().split("\t") chrom = lineContents[0] chromStart = lineContents[1] chromEnd = lineContents[2] name = lineContents[3] ref = lineContents[9] alt = lineContents[10] varClass = lineContents[12] submitters = lineContents[13] #Convert from SO term to name. e.g. SO:0001483 to substitution varClass = convertSOnumberToName(varClass) #Check if is the same rsID, they are repeated in ~5% of file if prevName == name: print("Repeated prevName: "+prevName) outputFile.write(chrom+"\t"+chromStart+"\t"+chromEnd+"\t"+name+"\t0\t.\t"+chromStart+"\t"+chromEnd+"\t"+ itemRgb+"\t"+ref+"\t"+alt+"\t"+varClass+"\t"+submitters+"\t"+ucscClass+"\t"+aaChange+"\n") elif name in duplicateSet: #Find the aaChange and ucscClass #Try and match the allele with what hgVai outputs if varClass == "insertion": if ',' in alt: nameAlt = name+alt.split(",")[0][1:] elif len(alt)-1 >= 59: nameAlt = name+alt[len(ref):25]+"..."+alt[len(alt)-24:]+"("+str(len(alt)-1)+" nt)" else: nameAlt = name+alt[1:] elif varClass == "deletion": nameAlt = name+"-" elif varClass == "tandem_repeat": if len(alt) > len(ref): nameAlt = name+alt[len(ref):] else: nameAlt = name+"-" else: nameAlt = name+alt prevName,currentChrom,currentVepDic,ucscClass,aaChange,tryAttempts = fetchUcscClassAndAaChange(name,alt,nameAlt,chrom,currentChrom,currentVepDic,tryAttempts,workDir) itemRgb = assignRGBcolorByConsequence(currentVepDic,name,alt,nameAlt) outputFile.write(chrom+"\t"+chromStart+"\t"+chromEnd+"\t"+name+"\t0\t.\t"+chromStart+"\t"+chromEnd+"\t"+ itemRgb+"\t"+ref+"\t"+alt+"\t"+varClass+"\t"+submitters+"\t"+ucscClass+"\t"+aaChange+"\n") else: nameAlt = name prevName,currentChrom,currentVepDic,ucscClass,aaChange,tryAttempts = fetchUcscClassAndAaChange(name,alt,nameAlt,chrom,currentChrom,currentVepDic,tryAttempts,workDir) itemRgb = assignRGBcolorByConsequence(currentVepDic,name,alt,nameAlt) outputFile.write(chrom+"\t"+chromStart+"\t"+chromEnd+"\t"+name+"\t0\t.\t"+chromStart+"\t"+chromEnd+"\t"+ itemRgb+"\t"+ref+"\t"+alt+"\t"+varClass+"\t"+submitters+"\t"+ucscClass+"\t"+aaChange+"\n") print("Total number of try attempts = "+str(tryAttempts)) print("Total number of lines in output file = "+str(n)) inputFile.close() outputFile.close() def validateFinalFile(workDir,url): """Compare and make sure that the final bed rsID# matches the input file""" originalDownloadFileName = workDir+url.split("/")[-1] bedFile = workDir+"evaSnp.final.bed" originalRsIDCount = bash("zcat %s | grep -v '^#' | cut -f3 | sort | uniq | wc -l" % (originalDownloadFileName)) finalBedRsIDCount = bash("cut -f4 %s | sort | uniq | wc -l" % (bedFile)) if originalRsIDCount != finalBedRsIDCount: sys.exit("There was an error in the pipeline, the initial numer of uniq RS IDs does not match the final bed file.\n\n"\ "%s - %s%s - %s" % (originalDownloadFileName,str(originalRsIDCount),bedFile,str(finalBedRsIDCount))) finalBedLineCount = bash("cat %s | wc -l" % (bedFile)) finalBedUniqueLineCount = bash("cat %s | sort | uniq | wc -l" % (bedFile)) if finalBedLineCount != finalBedUniqueLineCount: sys.exit("There are duplicate entries in the file: "+str(finalBedLineCount)+" vs "+str(finalBedUniqueLineCount)) def createBigBed(workDir,dbs): bedFile = workDir+"evaSnp.final.bed" bash("bedSort %s %s" % (bedFile,bedFile)) bash("bedToBigBed -tab -as=$HOME/kent/src/hg/lib/evaSnp.as -type=bed9+6 -extraIndex=name \ %s http://hgdownload.soe.ucsc.edu/goldenPath/%s/bigZips/%s.chrom.sizes %sevaSnp4.bb" % (bedFile,dbs,dbs,workDir)) def chechIfDbsDone(key): if os.path.exists("/hive/data/outside/eva4/"+key+"/evaSnp4.bb"): print("DBS already completed: "+key) return(True) else: if os.path.exists("/hive/data/outside/eva4/"+key+"/evaSnps.ucscChroms.vcf") or os.path.exists("/hive/data/outside/eva4/"+key+"/evaSnps.ucscChroms.vcf.gz"): print("Cleaning up directory for previous attempt of "+key) bash("rm /hive/data/outside/eva4/"+key+"/evaSnps*") if os.path.exists("/hive/data/outside/eva4/"+key+"/splitChroms"): try: bash("rm /hive/data/outside/eva4/"+key+"/splitChroms/*") except: print("Directly cleaned") return(False) def buildEvaSnp4Track(ucscDatabaseName,ucscEvaSnp4Dbs): """Build evaSnp4 track from start to finish""" dbs = ucscDatabaseName GCAacc = ucscEvaSnp4Dbs[dbs] workDir = "/hive/data/outside/eva4/"+dbs+"/" urlToFile = "https://ftp.ebi.ac.uk/pub/databases/eva/rs_releases/release_4/by_assembly/"+GCAacc+"/"+GCAacc+"_current_ids.vcf.gz" checkDirMakeDir(workDir) fileName = wgetFile(urlToFile, workDir) fileName = workDir+urlToFile.split("/")[-1] print("Got file: ",dbs) currentFileName = checkDupsAndUnzip(fileName,workDir) print("File unzipped") try: currentFileName = convertChromToUcsc(currentFileName,dbs,workDir) print("Chroms converted") except: try: sanitizeChromNames(currentFileName,dbs,workDir) print("Chroms sanitized") currentFileName = convertChromToUcsc(currentFileName,dbs,workDir) print("Chroms converted") except: currentFileName = addChrToChromNumbers(fileName,currentFileName,workDir) print("Chroms completed after adding '1' to starts.") currentFileName = convertVcfToBed(currentFileName,workDir) print("File converted to bed") hgVaiResultsFile = splitChromsAndRunHgVai(workDir,dbs) print("hgVai done") splitVepFileByChrom(hgVaiResultsFile,workDir) print("hgVai file split") currentFileName = workDir+"evaSnp.bed" createFinalBedWithVepAnnotations(currentFileName,workDir) print("Validating file") validateFinalFile(workDir,urlToFile) print("Created final file") createBigBed(workDir,dbs) print("bigBed made") def main(mode): #Run main with either "all" or a specific UCSC database name #Find all the assemblies we have with matching evaSnp4 data assemblyAlias = findAllAccFromDescPages() genbankSummaryDic = buildGCFdicFromGenbankFile() assemblyAlias = findGCAaccForUCSCgcfDbs(assemblyAlias,genbankSummaryDic) ucscGCAlist = sanitizeListToOnlyUCSCdbsWithGCA(assemblyAlias) ucscEvaSnp4Dbs = buildFinalDicOfUCSCdbsWithEvaSnp4Data(ucscGCAlist) #Build evaSNP4 data for either one specific assembly, or all assemblies if mode == 'all': for key in ucscEvaSnp4Dbs: if not chechIfDbsDone(key): try: buildEvaSnp4Track(key,ucscEvaSnp4Dbs) except: print("Problem with DBS: "+key) continue else: chechIfDbsDone("danRer11") buildEvaSnp4Track(mode,ucscEvaSnp4Dbs) main("all") ##################################### ######### End of python3 code ###### ##################################### # Lastly links were made with a file containing the 36 assemblies: for dbs in $(cat /hive/data/outside/eva4/assemblyReleaseList.txt); do ln -s /hive/data/outside/eva4/$dbs/evaSnp4.bb /gbdb/$dbs/bbi/; done