Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

supported to run ROSE2 in python3 environment #67

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

supported to run ROSE2 in python3 environment #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
344a490
Modified files to fix the issue from python2 to python3
gauravj49 Jul 16, 2019
aa76098
Modified files to fix the issue from python2 to python3
gauravj49 Jul 16, 2019
2a2fd18
added support for ucsc and ensembl files (ex. chr1 vs 1)
gauravj49 Sep 30, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 32 additions & 33 deletions ROSE2_geneMapper.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@

import os
import subprocess
from string import join
# from string import join
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please delete code instead of commenting it out.


from collections import defaultdict

Expand Down Expand Up @@ -189,14 +189,14 @@ def mapEnhancerToGene(annotFile,enhancerFile,transcribedFile='',uniqueGenes=True
if noFormatTable:

newEnhancerLine = list(line)
newEnhancerLine.append(join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes]),','))
newEnhancerLine.append(join(utils.uniquify([startDict[x]['name'] for x in proximalGenes]),','))
newEnhancerLine.append(','.join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes])))
newEnhancerLine.append(','.join(utils.uniquify([startDict[x]['name'] for x in proximalGenes])))
newEnhancerLine.append(closestGene)

else:
newEnhancerLine = line[0:9]
newEnhancerLine.append(join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes]),','))
newEnhancerLine.append(join(utils.uniquify([startDict[x]['name'] for x in proximalGenes]),','))
newEnhancerLine.append(','.join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes])))
newEnhancerLine.append(','.join(utils.uniquify([startDict[x]['name'] for x in proximalGenes])))
newEnhancerLine.append(closestGene)
newEnhancerLine += line[-2:]

Expand Down Expand Up @@ -238,9 +238,9 @@ def mapEnhancerToGene(annotFile,enhancerFile,transcribedFile='',uniqueGenes=True
proxEnhancers = geneDict['overlapping'][refID]+geneDict['proximal'][refID]

superStatus = max(superDict[refID])
enhancerRanks = join([str(x) for x in rankDict[refID]],',')
enhancerRanks = ','.join([str(x) for x in rankDict[refID]])

newLine = [geneName,refID,join(proxEnhancers,','),enhancerRanks,superStatus]
newLine = [geneName,refID,','.join(proxEnhancers),enhancerRanks,superStatus]
geneToEnhancerTable.append(newLine)

#resort enhancerToGeneTable
Expand Down Expand Up @@ -320,8 +320,7 @@ def mapEnhancerToGeneTop(rankByBamFile, controlBamFile, genome, annotFile, enhan
# matter
tssCollection = utils.LocusCollection(tssLoci, 50)

geneDict = {'overlapping': defaultdict(
list), 'proximal': defaultdict(list)}
geneDict = {'overlapping': defaultdict(list), 'proximal': defaultdict(list)}

# dictionaries to hold ranks and superstatus of gene nearby enhancers
rankDict = defaultdict(list)
Expand All @@ -347,16 +346,13 @@ def mapEnhancerToGeneTop(rankByBamFile, controlBamFile, genome, annotFile, enhan
else:
# set up the output tables
# first by enhancer
enhancerToGeneTable = [
header[0:9] + ['OVERLAP_GENES', 'PROXIMAL_GENES', 'CLOSEST_GENE'] + header[-2:]]
enhancerToGeneTable = [header[0:9] + ['OVERLAP_GENES', 'PROXIMAL_GENES', 'CLOSEST_GENE'] + header[-2:]]

# next by gene
geneToEnhancerTable = [
['GENE_NAME', 'REFSEQ_ID', 'PROXIMAL_ENHANCERS']]
geneToEnhancerTable = [['GENE_NAME', 'REFSEQ_ID', 'PROXIMAL_ENHANCERS']]

# next make the gene to enhancer table
geneToEnhancerTable = [
['GENE_NAME', 'REFSEQ_ID', 'PROXIMAL_ENHANCERS', 'ENHANCER_RANKS', 'IS_SUPER', 'ENHANCER_SIGNAL']]
geneToEnhancerTable = [['GENE_NAME', 'REFSEQ_ID', 'PROXIMAL_ENHANCERS', 'ENHANCER_RANKS', 'IS_SUPER', 'ENHANCER_SIGNAL']]

for line in enhancerTable:
if line[0][0] == '#' or line[0][0] == 'R':
Expand Down Expand Up @@ -405,35 +401,38 @@ def mapEnhancerToGeneTop(rankByBamFile, controlBamFile, genome, annotFile, enhan
distalGenes.remove(refID)

# Now find the closest gene
closestGene = ''
if len(allEnhancerGenes) == 0:
closestGene = ''
else:
# get enhancerCenter
enhancerCenter = (int(line[2]) + int(line[3])) / 2

# get absolute distance to enhancer center
distList = [abs(enhancerCenter - startDict[geneID]['start'][0])
for geneID in allEnhancerGenes]
# get the ID and convert to name
closestGene = startDict[
allEnhancerGenes[distList.index(min(distList))]]['name']
try:
# get absolute distance to enhancer center
distList = [abs(enhancerCenter - list(startDict[geneID]['start'])[0]) for geneID in allEnhancerGenes]

# get the ID and convert to name
closestGene = startDict[allEnhancerGenes[distList.index(min(distList))]]['name']
except:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you want this to fail silently? Or give some meaningful message to the user? Or catch a specific kind of error that you expect?

pass

# NOW WRITE THE ROW FOR THE ENHANCER TABLE
if noFormatTable:

newEnhancerLine = list(line)
newEnhancerLine.append(
join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes]), ','))
','.join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes])))
newEnhancerLine.append(
join(utils.uniquify([startDict[x]['name'] for x in proximalGenes]), ','))
','.join(utils.uniquify([startDict[x]['name'] for x in proximalGenes])))
newEnhancerLine.append(closestGene)

else:
newEnhancerLine = line[0:9]
newEnhancerLine.append(
join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes]), ','))
','.join(utils.uniquify([startDict[x]['name'] for x in overlappingGenes])))
newEnhancerLine.append(
join(utils.uniquify([startDict[x]['name'] for x in proximalGenes]), ','))
','.join(utils.uniquify([startDict[x]['name'] for x in proximalGenes])))
newEnhancerLine.append(closestGene)
newEnhancerLine += line[-2:]

Expand All @@ -460,14 +459,15 @@ def mapEnhancerToGeneTop(rankByBamFile, controlBamFile, genome, annotFile, enhan

#get the chromLists from the various bams here
cmd = 'samtools idxstats %s' % (rankByBamFile)
idxStats = subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
idxStats= idxStats.communicate()
bamChromList = [line.split('\t')[0] for line in idxStats[0].split('\n')[0:-2]]
idxStats = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True,
universal_newlines=True).communicate()[0]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Generally it's not good practice to use shell=True, but probably for a small scientific pipeline this would be okay. Here is an example of not using shell=True: https://github.com/HPC-buildtest/buildtest-framework/blob/devel/buildtest/utils/command.py#L55


bamChromList = [bytes(line.split('\t')[0]) for line in idxStats[0].split('\n')[0:-2]]

if len(controlBamFile) > 0:
cmd = 'samtools idxstats %s' % (controlBamFile)
idxStats = subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
idxStats= idxStats.communicate()
idxStats = idxStats.communicate()
bamChromListControl = [line.split('\t')[0] for line in idxStats[0].split('\n')[0:-2]]
bamChromList = [chrom for chrom in bamChromList if bamChromListControl.count(chrom) != 0]

Expand Down Expand Up @@ -572,13 +572,12 @@ def mapEnhancerToGeneTop(rankByBamFile, controlBamFile, genome, annotFile, enhan
refID] + geneDict['proximal'][refID]

superStatus = max(superDict[refID])
enhancerRanks = join([str(x) for x in rankDict[refID]], ',')
enhancerRanks = ','.join([str(x) for x in rankDict[refID]])

enhancerSignal = signalDict[refID]
geneNameSigDict[geneName].append(enhancerSignal)

newLine = [geneName, refID, join(
proxEnhancers, ','), enhancerRanks, superStatus, enhancerSignal]
newLine = [geneName, refID, ','.join(proxEnhancers), enhancerRanks, superStatus, enhancerSignal]
geneToEnhancerTable.append(newLine)
#utils.unParseTable(geneToEnhancerTable,'/grail/projects/newRose/geneMapper/foo.txt','\t')
print('MAKING ENHANCER TO TOP GENE TABLE')
Expand Down Expand Up @@ -711,7 +710,7 @@ def main():
if options.out:
outFolder = utils.formatFolder(options.out, True)
else:
outFolder = join(enhancerFile.split('/')[0:-1], '/') + '/'
outFolder = '/'.join(enhancerFile.split('/')[0:-1]) + '/'

# GETTING BAM INFO
rankByBamFile = options.rankby
Expand Down
Loading