from Bio import SeqIO
from Bio.Seq import Seq
import os
import uuid
from Bio import pairwise2
import math
import subprocess
from colorama import Fore #,Back, Style
from colorama import Back
from colorama import Style
from colorama import init
init(autoreset=True)
import numpy as np
import re
from itertools import takewhile,repeat
#import localInstall
if not os.name=='nt':
import pysam
### custom path config: (not used when not available)
#localInstall.addToPath('/hpc/hub_oudenaarden/bdebarbanson/bin/samtools-1.3.1/')
#localInstall.addToPath('/hpc/hub_oudenaarden/bin/software/bwa-0.7.10/')
#localInstall.addToPath('/media/sf_externalTools/')
### end
class AlignmentVisualiser:
def __init__(self, initialiser):
if type(initialiser) is pysam.AlignedSegment:
self.alignment = initialiser
def getPhredScoreColor(self, phredScore):
p = ord(phredScore)-33.0
if p>32:
return(Fore.GREEN + Style.BRIGHT)
if p>25:
return(Fore.GREEN)
if p>20:
return(Fore.GREEN + Style.DIM)
if p>18:
return(Fore.YELLOW)
return(Fore.RED)
def visualizeAlignment(self, annotateReference=None, annotateQuery=None, intronSummary=True):
reference = []
query = []
alignmentSymbols = []
qualitieSymbols = []
prevIntron = 0
for index,pair in enumerate(self.alignment.get_aligned_pairs( matches_only=False, with_seq=True) ):
#Pair structure = (queryPosition, referencePosition, referenceBase)
if pair[0] is None and pair[2] is None:
prevIntron+=1
continue
else:
if prevIntron>0:
intronString = ' N%s ' % prevIntron
alignmentSymbols.append(Style.DIM + intronString + Style.RESET_ALL)
qualitieSymbols.append(' '*len(intronString))
reference.append(Style.DIM + ('-'*len(intronString)) + Style.RESET_ALL)
query.append( Style.DIM + (' '*len(intronString)) + Style.RESET_ALL)
prevIntron = 0
queryPosition = pair[0]
queryNucleotide= "%s%s-%s"% (Back.RED, Fore.WHITE, Style.RESET_ALL)
referenceNucleotide= ' '# % (Fore.RED, Style.RESET_ALL) #if queryPosition>self.alignment.query_alignment_start and queryPosition<=self.alignment.query_alignment_end else ' '
queryPhred=' '
if queryPosition!=None:
queryNucleotide = self.alignment.query_sequence[queryPosition].upper()
queryPhred = self.alignment.qual[queryPosition]
qualitieSymbols.append(self.getPhredScoreColor(queryPhred) + queryPhred + Style.RESET_ALL)
else:
qualitieSymbols.append(' ')
#print("Q:%s S:%s" % (queryPosition, self.alignment.query_alignment_start))
if pair[2]!=None:
referenceNucleotide = pair[2].upper()
elif queryPosition is None:
referenceNucleotide = "%s%s-%s"% (Back.RED, Fore.WHITE, Style.RESET_ALL)
elif queryPosition>self.alignment.query_alignment_start and queryPosition<=self.alignment.query_alignment_end:
referenceNucleotide = "%s%s-%s"% (Back.RED, Fore.WHITE, Style.RESET_ALL)
if pair[2]==queryNucleotide:
alignmentSymbols.append('|')
else:
alignmentSymbols.append(' ')
reference.append(referenceNucleotide)
if queryPosition!=None and annotateQuery!=None:
if ((queryPosition - self.alignment.query_length) in annotateQuery and (queryPosition - self.alignment.query_length<0 )):
query.append(getattr(Back,annotateQuery[(queryPosition - self.alignment.query_length)]) + Fore.BLACK + queryNucleotide + Style.RESET_ALL)
elif queryPosition in annotateQuery:
query.append(getattr(Back,annotateQuery[queryPosition]) + Fore.WHITE + queryNucleotide + Style.RESET_ALL)
else:
#DIM soft clipped things:
if queryPosition<self.alignment.query_alignment_start or queryPosition>self.alignment.query_alignment_end:
query.append(Style.DIM + Fore.YELLOW + queryNucleotide + Style.RESET_ALL)
else:
query.append(queryNucleotide)
else:
query.append(queryNucleotide)
rStart = self.alignment.reference_end if self.alignment.is_reverse else self.alignment.reference_start
rEnd = self.alignment.reference_start if self.alignment.is_reverse else self.alignment.reference_end
refStartStr = str(rStart)
qStartStr = str(self.alignment.query_alignment_start)
labelSize = max(len(refStartStr), len(qStartStr))
annotatedCigar = []
for operation,length in self.alignment.cigartuples:
annotatedCigar.append( '%s%s%s' % (
[ '%s%sM%s'% (Fore.GREEN, Style.BRIGHT, Style.NORMAL),
'%s%sI%s' % (Fore.YELLOW,Style.BRIGHT, Style.NORMAL),
'%s%sD%s' % (Fore.RED,Style.BRIGHT, Style.NORMAL),
'%sN' % Fore.MAGENTA,
'%s%sS' % (Style.DIM,Fore.RED),
'%s%sH' % (Style.DIM,Fore.RED),
'%s%sP' % (Style.DIM,Fore.RED),
'%s%s=' % (Style.DIM,Fore.GREEN),
'%s%sX' % (Style.DIM,Fore.BLUE)]
[operation]
, length, Style.RESET_ALL))
multiMapping = ''
try:
multiMappingCount = int(self.alignment.get_tag('NH:i'))
if multiMappingCount==1:
multiMapping= Fore.GREEN + 'Unique map' + Style.RESET_ALL
elif multiMappingCount>1:
multiMapping= Fore.WHITE + Back.RED + ('Mapping to at least %s locations' % multiMappingCount) + Style.RESET_ALL
except:
pass
print( ('%sReference:%s %s %s %s MQ:%s %s' % ( Style.BRIGHT, Style.RESET_ALL, self.alignment.reference_name, "".join(annotatedCigar), ("FWD" if not self.alignment.is_reverse else "REV"), self.alignment.mapping_quality, multiMapping))+ (' MULTIMAPPING' if self.alignment.is_secondary else '' ) )
print(('%s%s %s %s (%s%s)' % (
Style.DIM+refStartStr+Style.RESET_ALL,
" "*(labelSize-len(refStartStr)),
"".join(reference),
(Style.DIM+str(rEnd)+Style.RESET_ALL),
"" if self.alignment.is_reverse else "+",
(Style.DIM+str(rEnd-rStart)+Style.RESET_ALL))))
print(('%s%s' % (" "*(labelSize+1), "".join(alignmentSymbols))))
print(('%s%s %s %s (%s)' % (
" "*(labelSize-len(qStartStr)),
Style.DIM+qStartStr+Style.RESET_ALL,
"".join(query),
(Style.DIM+str(self.alignment.query_alignment_end)+Style.RESET_ALL),
(Style.DIM+str(self.alignment.query_alignment_end-self.alignment.query_alignment_start)+Style.RESET_ALL),
)))
print(('%s%s' % (" "*(labelSize+1), "".join(qualitieSymbols))))
print(('%sQuery:%s %s' % (Style.BRIGHT, Style.RESET_ALL,self.alignment.query_name)))
def getLevenshteinDistance(source, target, maxDist=99999999999999):
if len(source) < len(target):
return getLevenshteinDistance(target, source)
# So now we have len(source) >= len(target).
if len(target) == 0:
return len(source)
# We call tuple() to force strings to be used as sequences
# ('c', 'a', 't', 's') - numpy uses them as values by default.
source = np.array(tuple(source))
target = np.array(tuple(target))
# We use a dynamic programming algorithm, but with the
# added optimization that we only need the last two rows
# of the matrix.
previous_row = np.arange(target.size + 1)
for s in source:
# Insertion (target grows longer than source):
current_row = previous_row + 1
# Substitution or matching:
# Target and source items are aligned, and either
# are different (cost of 1), or are the same (cost of 0).
current_row[1:] = np.minimum(
current_row[1:],
np.add(previous_row[:-1], target != s))
# Deletion (target grows shorter than source):
current_row[1:] = np.minimum(
current_row[1:],
current_row[0:-1] + 1)
previous_row = current_row
if( np.amin(previous_row)>maxDist):
return(None)
return previous_row[-1]
def cigarStringToDict(cigarString):
if cigarString==None:
cs = ''
else:
cs = cigarString
cigarStringDict = {}
for symbol in ['M','D','I','S']:
# Extract the integer which is found before the symbol supplied
# returns 0 when the symbol is not found as the sum of the list will be zero
matches = re.findall('[0-9]*%s' % symbol, cs)
cigarStringDict[symbol] = sum( [ int(x.replace(symbol,'')) for x in matches])
return( cigarStringDict)
def distanceLineToPoint(x1, y1, x2, y2, x0, y0):
return( float( abs( (y2-y1)*x0 - (x2 - x1 )*y0 + x2*y1 - y2*x1 ) ) / (math.sqrt( math.pow(y2-y1,2) + math.pow(x2-x1,2))) )
def fixGraphML(graphmlPath):
with open(graphmlPath) as f:
lines = []
for line in f:
lines.append( line.replace('"d3"','"fixedD3"') )
if lines!=None:
with open(graphmlPath,'w') as f:
for line in lines:
f.write(line)
def getHammingDistance(s1, s2,maxDist=999999):
d = 0
for index,base in enumerate(s1):
d+=(base!=s2[index] and s2[index]!="N" and base!="N")
if d>maxDist:
return(maxDist+1)
return(d)
def getHammingIndices( seqA, seqB, maxIndices=999999 ):
indices = []
for index,base in enumerate(seqA):
if len(seqB)<index:
break
if base!=seqB[index]:
indices.append(index)
if indices>=maxIndices:
return(indices)
return(indices)
def humanReadable(value, targetDigits=2,fp=0):
if value == 0:
return('0')
baseId = int(math.floor( math.log10(float(value))/3.0 ))
suffix = ""
if baseId==0:
sVal = str(round(value,targetDigits))
if len(sVal)>targetDigits and sVal.find('.'):
sVal = sVal.split('.')[0]
elif baseId>0:
sStrD = max(0,targetDigits-len(str( '{:.0f}'.format((value/(math.pow(10,baseId*3)))) )))
sVal = ('{:.%sf}' % min(fp, sStrD)).format((value/(math.pow(10,baseId*3))))
suffix = 'kMGTYZ'[baseId-1]
else:
sStrD = max(0,targetDigits-len(str( '{:.0f}'.format((value*(math.pow(10,-baseId*3)))) )))
sVal = ('{:.%sf}' % min(fp, sStrD)).format((value*(math.pow(10,-baseId*3))))
suffix = 'mnpf'[-baseId-1]
if len(sVal)+1>targetDigits:
# :(
sVal = str(round(value,fp))[1:]
suffix = ''
return('%s%s' % (sVal,suffix))
##Pathtools
def decodePath( path, spacers={'flag': ',', 'param':'_', 'kv':'=', 'invalid':['=','_']}):
invalidFileNameCharacters = spacers['invalid']
flagSpacer = spacers['flag']
paramSpacer = spacers['param']
keyValueSpacer = spacers['kv']
decodedPath = {
'runId':None,
'step':-1,
'name':'unknown',
'parameters':{},
'flags' : [],
'extension':''
}
parts = os.path.splitext(path)
basePath = parts[0]
decodedPath['extension'] = parts[1]
#Check if the anlysis base directory is in the path:
parts = basePath.split('/')
if 'analysis_' in parts[0]:
decodedPath['runId'] = parts[0]
toDecode = parts[1]
else:
toDecode = path
keyValuePairs = basePath.split(paramSpacer)
for pair in keyValuePairs:
keyValue = pair.split(keyValueSpacer)
if len(keyValue)!=2:
print(('Parameter decoding failed: %s, this probably means input files contain invalid characters such as %s' % (pair, ", ".join(invalidFileNameCharacters))))
else:
key = keyValue[0]
value = keyValue[1]
if key!='flags' and key!='name':
decodedPath['parameters'][key] = value
elif key=='flags':
decodedPath['flags'] = value.split(flagSpacer)
elif key=='name':
decodedPath['name'] = value
return(decodedPath)
def invBase(nucleotide):
if nucleotide=='A':
return(['T','C','G'])
if nucleotide=='C':
return(['T','A','G'])
if nucleotide=='T':
return(['C','A','G'])
if nucleotide=='G':
return(['C','A','T'])
return([''])
def encodePath(step, name, params, extension, spacers={'flag': ',', 'param':'_', 'kv':'=', 'invalid':['=','_']}):
invalidFileNameCharacters = spacers['invalid']
flagSpacer = spacers['flag']
paramSpacer = spacers['param']
keyValueSpacer = spacers['kv']
#Todo get keyValueSpacer in here
encodedPath = 'step={:0>2d}'.format(step)
encodedPath += '_name=%s' % name
for parname in params:
if parname.lower()!="flags":
if not( any(i in parname for i in invalidFileNameCharacters) or any(i in params[parname] for i in invalidFileNameCharacters) ):
encodedPath+='%s%s%s%s' % (paramSpacer, parname, keyValueSpacer, params[parname])
else:
raise ValueError('Parameter encoding failed: %s %s contains at least one invalid character' % (parname, params[parname]) )
#Todo check flags for mistakes
if 'flags' in params:
encodedPath+='%sflags%s%s' % (paramSpacer, keyValueSpacer,flagSpacer.join(params['flags']))
if extension!=None and extension!=False:
encodedPath = '%s.%s' % (encodedPath, extension)
return(encodedPath)
def locateIndel( seqA, seqB, verbose=False ):
alignment = pairwise2.align.globalmx(seqA, seqB, 2, -1, one_alignment_only=True)[0]
#pairwise2.format_alignment(*alignment)
indelLocA = -1
indelLocB = -1
if alignment:
align1 = alignment[0]
align2 = alignment[1]
localAlign1 = align1.strip('-')
localAlign2 = align2.strip('-')
dels = min( localAlign1.count('-'), localAlign2.count('-') )
ins = max( localAlign1.count('-'), localAlign2.count('-') )
if (ins==1) and dels==0:
seqNHasIndel = localAlign2.count('-')==0 #bool: 0: seqA has indel, 1: seqB has indel
#Find indel location:
if seqNHasIndel:
indelLocA = align1.find('-')
if verbose:
print(('seqA:%s\nseqB:%s'% (align1,align2)))
print((' '*(indelLocA+5)+'^'))
else:
indelLocB = align2.find('-')
if verbose:
print(('seqA:%s\nseqB:%s'% (align1,align2)))
print((' '*(indelLocB+5)+'^'))
#print(alignment[2]) is score
return(indelLocA,indelLocB)
class Histogram():
def __init__(self):
self.data = {}
def addCount(self,idx, count=1):
if not idx in self.data:
self.data[idx] = 0
self.data[idx]+=count
def write(self, path):
with open(path,'w') as f:
for index in self.data:
f.write('%s;%s\n' % (index, self.data[index]))
globalWritingUpdates = False
class InternalTool():
def __init__(self, screenName = "unnamed"):
self.verbose = True
self.screenName = screenName
def setVerbosity(self, verbosity):
self.verbose = verbosity
def __priorRunningPrint(self):
if self.isUpdating():
print('')
self.setUpdating(False)
def isUpdating(self):
global globalWritingUpdates
if globalWritingUpdates:
return(True)
def setUpdating(self, boolean):
global globalWritingUpdates
globalWritingUpdates = boolean
def log(self, message):
if self.verbose:
self.__priorRunningPrint()
print((Fore.GREEN + self.screenName + ', INFO: ' + Style.RESET_ALL + message))
def warn(self, message):
self.__priorRunningPrint()
print((Fore.RED + self.screenName + ', WARN: ' + Style.RESET_ALL + message))
def softWarn(self, message):
self.__priorRunningPrint()
print((Fore.YELLOW + self.screenName + ', SOFT WARN: ' + Style.RESET_ALL + message))
def update(self, message):
print('\r' + Fore.GREEN + self.screenName + ', UPDATE: ' + Style.RESET_ALL +message, end=' ')
self.setUpdating(True)
def info(self, message):
self.log(message)
def fatal(self, message):
self.__priorRunningPrint()
print((Fore.RED + self.screenName + ', FATAL ERROR: ' + Style.RESET_ALL + message))
exit()
class ExternalTool(InternalTool):
def __init__(self, screenName = "unnamed"):
InternalTool.__init__(self)
#Locations to look for binaries when they are not found on the PATH
def setExecutablePath(self, path):
self.executable = path
#Only works in unix and cygwin
def isAvailable(self):
available = subprocess.call("type " + self.executable.split()[0], shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0
return(available or os.path.isfile(self.executable))
class FastqToFasta(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('fastq_to_fasta')
self.screenName = "FASTQ_TO_FASTA"
def execute(self, fqPath, fastaPath):
os.system('%s -n -i %s -o %s' % (self.executable, fqPath, fastaPath ))
class FastQQualityFilter(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('fastq_quality_filter')
self.minBaseQuality = 36
self.sequenceCountSuffix = 'fastqQualityResultCount'
self.screenName = "FASTQ_QUALITY_FILTER"
def setMinBaseQuality(self, minBaseQuality):
self.minBaseQuality = minBaseQuality
def getCountPath(self,filteredTargetPath):
return( '%s_%s' % (filteredTargetPath,self.sequenceCountSuffix) )
def execute(self, fqPath, filteredTargetPath):
os.system('%s -q %s -i %s -o %s -v -sam > %s' % (self.executable, self.minBaseQuality, fqPath, filteredTargetPath, self.getCountPath()))
r = 0
with open(self.getCountPath(), 'w') as f:
r = int(f.readline())
return(r)
class ART(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('./art_bin_ChocolateCherryCake/art_illumina')
self.screenName = "ART"
def simulateAmpliconReads(self,sequenceAbundanceCounter, prefix, readSize=76, coverage=10, reverseComplementAmplicon=False):
sequences = []
print(( sequenceAbundanceCounter.most_common(1)))
(v,hCount) = sequenceAbundanceCounter.most_common(1)[0]
for iteration in range(0,hCount):
for sequence,count in sequenceAbundanceCounter.most_common():
if count>=hCount:
sequences.append(sequence)
else:
break
hCount-=1
bpythonSeqs=[]
for index,sequence in enumerate(sequences):
if reverseComplementAmplicon :
bpythonSeqs.append(SeqIO.SeqRecord(Seq(sequence).reverse_complement(), '%s-%s' % (str(index),str(sequence))))
else:
bpythonSeqs.append(SeqIO.SeqRecord(Seq(sequence), '%s-%s' % (str(index),str(sequence))))
fastaPath = getTempFileName('art')+'.fa'
SeqIO.write(bpythonSeqs, fastaPath, "fasta")
os.system('%s -amp -i %s -l %s -f %s -o %s -sam -ss HS25' %(self.executable, fastaPath, readSize, coverage, prefix))
def setMinBaseQuality(self, minBaseQuality):
self.minBaseQuality = minBaseQuality
def getCountPath(self,filteredTargetPath):
return( '%s_%s' % (filteredTargetPath,self.sequenceCountSuffix) )
def execute(self, fqPath, filteredTargetPath):
os.system('%s -q %s -i %s -o %s -v > %s' % (self.executable, self.minBaseQuality, fqPath, filteredTargetPath, self.getCountPath()))
r = 0
with open(self.getCountPath(), 'w') as f:
r = int(f.readline())
return(r)
class NeedleAll(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('needleall')
self.screenName = "NEEDLE_ALL"
def execute(self, fqPathA, fqPathB, outputFilePath):
command = '%s -auto -stdout -asequence "%s" -bsequence "%s" > %s' % (self.executable, fqPathA, fqPathB, outputFilePath)
print(command)
os.system(command)
r = 0
#This is an inverted distance matrix
invScores = DistanceMatrix()
#find max score:
maxScore = 0
# To perform T-SNE a normalised matrix is required: the total sum should be 1
totalScore = 0
cells = 0
with open(outputFilePath, 'r') as f:
for line in f:
#print(line)
parts = line.strip().replace('(','').replace(')','').split(' ')
if len(parts)==4:
value = float(parts[3])
if value>maxScore:
maxScore=value
totalScore+=value
cells+=1
#Normalise the scores
matrixSum = 0
with open(outputFilePath, 'r') as f:
for line in f:
#print(line)
parts = line.strip().replace('(','').replace(')','').split(' ')
#print('%s %s %s %s' % (parts[0], parts[1],parts[2],parts[3]))
if len(parts)==4:
if float(parts[3])>0:
value = (maxScore-float(parts[3]))
matrixSum+= value
invScores.setDistance(str(parts[0]), str(parts[1]), value)
print(('Sum of distance matrix is %s' % matrixSum))
return(invScores)
class ClustalOmega(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('clustalo')
self.screenName = "CLUSTAL-OMEGA"
def alignFasta(self, fastaPath,outputFastaFilePath):
self.log('aligning %s to %s' % (fastaPath, outputFastaFilePath))
os.system('%s -i %s -o %s --force' % (self.executable, fastaPath, outputFastaFilePath))
def alignFastqs(self, fqPath, outputFastaFilePath):
#self.log('converting %s to fasta' % fqpath)
f = FastqToFasta()
tmpFasta = getTempFileName()
f.execute(fqPath, tmpFasta)
self.alignFasta(tmpFasta,outputFastaFilePath)
os.remove(tmpFasta)
class SamTools(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('samtools')
self.screenName = "SAMTOOLS"
def samToSortedBam(self, samPath, bamPath, deleteSam=False):
self.log('Converting SAM to sorted BAM')
os.system('%s view -bS %s | %s sort - -o %s.bam' % (self.executable, samPath, self.executable, bamPath))
self.log('Completed conversion')
if deleteSam:
os.remove(samPath)
def bamToSortedBam(self, bamPath, sortedBamPath):
self.log('Converting BAM to sorted BAM')
os.system('%s sort -o %s.bam %s' % (self.executable, sortedBamPath, bamPath))
self.log('Completed conversion')
def index(self, bamPath):
self.log('Indexing BAM file')
os.system('%s index %s' % (self.executable, bamPath))
self.log('Completed indexing')
def readBamToDictOld(self,bamFilePath, columns={'cigar':5, 'left':3,'scarSequence':9}, index='scarSequence', appendTo={}): # columns:{ targetName:columnIndex, ... }
retDict = appendTo
#os.system("samtools view %s | cut -f1,6 | grep -oEi '^.*([0-9]+D){1,1}.*' ")
#tempfile = './temp/temp_' + str(uuid.uuid1())
tempfile = getTempFileName()
self.log("samtools view %s > %s" %(bamFilePath,tempfile))
os.system("samtools view %s > %s" %(bamFilePath,tempfile))
with open(tempfile) as f:
#Go over all lines in the BAM file
for line in f:
#Extract all columns:
parts = line.replace('\n','').split()
#Extract the primary key value:
indexKey = parts[columns[index]]
#Initialise the existence of the primary key value in the return dictionary
if not indexKey in retDict:
retDict[indexKey] = {}
#Iterate over the keysm except the index
for key in columns:
if key is not index:
#Get the value of the key on this line in the BAM file
v = parts[columns[key]]
if not v in retDict[indexKey]:
retDict[indexKey][v] = 1
else:
retDict[indexKey][v] += 1
os.remove(tempfile)
return(retDict)
def readBamToDict(self,bamFilePath, columns={'cigar':5, 'id':0,'left':3,'scarSequence':9}, index='scarSequence', appendTo={}): # columns:{ targetName:columnIndex, ... }
retDict = appendTo
#os.system("samtools view %s | cut -f1,6 | grep -oEi '^.*([0-9]+D){1,1}.*' ")
tempfile = './temp/temp_' + str(uuid.uuid1())
self.log("samtools view %s > %s" %(bamFilePath,tempfile))
os.system("samtools view %s > %s" %(bamFilePath,tempfile))
with open(tempfile) as f:
#Go over all lines in the BAM file
for line in f:
#Extract all columns:
parts = line.replace('\n','').split()
#Extract the primary key value:
indexKey = parts[columns[index]]
#Initialise the existence of the primary key value in the return dictionary
if not indexKey in retDict:
retDict[indexKey] = {}
#Iterate over the keysm except the index
for key in columns:
if key is not index:
#Get the value of the key on this line in the BAM file
if columns[key]<len(parts):
v = parts[columns[key]]
retDict[indexKey][key] = v
try:
os.remove(tempfile)
except:
print(('Temporary file could not be removed %s' % tempfile))
pass
return(retDict)
# Make sure to index the reference genome!
class BWA_MEM(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('bwa mem')
self.screenName = "BWA MEM"
self.fixedFlags = []
def setReferencePath(self, referenceGenomePath):
self.referenceGenomePath = referenceGenomePath
def execute(self, fqPath, samPath, readGroupInfo=None, secondSplitHits=False, getSecondaryAlignments=False):
if not isinstance(fqPath, str):
fqPath = " ".join(fqPath)
if self.referenceGenomePath!=None:
self.log('Aligning %s to %s' %(fqPath, self.referenceGenomePath))
flags=[]+self.fixedFlags
if secondSplitHits:
flags.append("-M")
if readGroupInfo!=None:
flags.append("-R '%s'" % readGroupInfo)
if getSecondaryAlignments:
flags.append("-a")
self.log("Flags: %s" % (" ".join(flags)))
cmd = '%s %s %s %s > %s' % (self.executable, " ".join(flags), self.referenceGenomePath, fqPath, samPath)
self.log(cmd)
os.system(cmd)
else:
self.warn('Warning; no reference genome specified for BWA MEM, canceled alignment.')
# Make sure to index the reference genome!
class STAR(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('STAR')
self.screenName = "STAR"
self.threads = 4
def index(self):
if not os.path.exists(self.refDirectory):
try:
os.mkdir(self.refDirectory )
except:
self.fatal('Could not create index directory on %s' % self.refDirectory)
os.system('%s --runMode genomeGenerate --genomeDir %s --genomeFastaFiles %s --runThreadN %s' %
(self.executable, self.refDirectory, self.referenceGenomePath, self.threads))
def setReferencePath(self, referenceGenomePath):
self.referenceGenomePath = referenceGenomePath
self.refDirectory = '%s_STAR_INDEX' % os.path.basename(self.referenceGenomePath)
if not os.path.exists(self.refDirectory):
self.info("Reference index %s does not exist." % referenceGenomePath)
self.index()
def execute(self, fqPath, outputPath):
os.system('%s --genomeDir %s --readFilesIn %s --runThreadN %s --outFileNamePrefix %s' %
(self.executable, self.refDirectory, fqPath, self.threads, outputPath))
#!!!sjdbOverhang should match the length of the reads minus one!!!
def pass2alignment(self, SJouttabPath, sjdbOverhang=74):
base = os.path.basename(SJouttabPath)
self.p2refDirectory = '%s_P2_STAR_INDEX' % os.path.basename(SJouttabPath)
if not os.path.exists(self.p2refDirectory):
try:
os.mkdir(self.p2refDirectory )
except:
self.fatal('Could not create index directory on %s' % self.p2refDirectory)
def mapSingleEnd(self, fqPath, samPath):
if not isinstance(fqPath, str):
fqPath = " ".join(fqPath)
if self.referenceGenomePath!=None:
self.log('Aligning %s to %s' %(fqPath, self.referenceGenomePath))
os.system('%s --genomeDir %s --readFilesIn %s --runThreadN %s' %
(self.executable, self.refDirectory, fqPath, self.threads))
else:
self.warn('Warning; no reference genome specified for STAR, canceled alignment.')
class GSNAP(ExternalTool):
def __init__(self):
ExternalTool.__init__(self)
self.setExecutablePath('gsnap')
self.screenName = "GSNAP"
#self.cmd= 'sudo gsnap -d lintrace analysis_bdefa5bc-6b5c-11e5-a487-0800279bd60a_fullScarList_Protocol1-DNA-seq.fq -t 4 -A sam > analysis_bdefa5bc-6b5c-11e5-a487-0800279bd60a_fullScarList_Protocol1-DNA-seq.sam'
#Reference genome path should be the name of the reference build folder
def setReferencePath(self, referenceGenomePath):
self.referenceGenomePath = referenceGenomePath
def execute(self, fqPath, samPath):
if self.referenceGenomePath!=None:
self.log('Aligning %s to %s' %(fqPath, self.referenceGenomePath))
os.system('%s -d %s %s --format=sam > %s' % (self.executable, self.referenceGenomePath, fqPath, samPath))
else:
self.warn('Warning; no reference genome specified for GSNAP, canceled alignment.')
def getTempFileName(centerFix=""):
if not os.path.exists('./temp'):
os.mkdir('./temp')
return('./temp/temp_' + centerFix + str(uuid.uuid1()))
def mapReads(readsFile, baseName, mapper, samTools=None,readGroupInfo=None, secondSplitHits=False):
if samTools==None:
samTools = SamTools()
basePath = '%s/%s' % ( '.', baseName)
samPath = '%s.sam' %( basePath)
bamPath = '%s.bam' %( basePath)
if readGroupInfo==None and secondSplitHits==False:
mapper.execute(readsFile, samPath)
else:
mapper.execute(readsFile, samPath, readGroupInfo=None, secondSplitHits=False)
samTools.samToSortedBam(samPath, basePath, True)
samTools.index(bamPath)
return({'samPath':samPath, 'bamPath':bamPath})
def mapSequencesToDict(mapper, sequences):
bpythonSeqs=[]
for index,sequence in enumerate(sequences):
bpythonSeqs.append(SeqIO.SeqRecord(Seq(sequence), str(index)))
fastaPath = getTempFileName('bwa')+'.fa'
SeqIO.write(bpythonSeqs, fastaPath, "fasta")
bPath = getTempFileName('samtools')
samPath = bPath+'.sam'
bamPath = bPath+'.bam'
mapper.execute(fastaPath, samPath)
samtools = SamTools()
samtools.samToSortedBam(samPath, bPath, True)
samtools.index(bamPath)
d = samtools.readBamToDict(bamPath)
return(d)
def mapDictOfSequences(mapper, sequenceDict, bPath=None):
bpythonSeqs=[]
for seqId in sequenceDict:
bpythonSeqs.append(SeqIO.SeqRecord(Seq(sequenceDict[seqId]), str(seqId)))
fastaPath = getTempFileName('bwa')+'.fa'
SeqIO.write(bpythonSeqs, fastaPath, "fasta")
if bPath==None:
bPath = getTempFileName('samtools')
samPath = bPath+'.sam'
bamPath = bPath+'.bam'
mapper.execute(fastaPath, samPath)
samtools = SamTools()
samtools.samToSortedBam(samPath, bPath, False)
samtools.index(bamPath)
return(bamPath)
def mapReadsToDict(mapper, fqPath, index='id',basePath=None, readGroupInfo=None, secondSplitHits=False):
if basePath==None:
bPath = getTempFileName('samtools')
else:
bPath=basePath
samPath = bPath+'.sam'
bamPath = bPath+'.bam'
mapper.execute(fqPath, samPath, readGroupInfo, secondSplitHits)
samtools = SamTools()
samtools.samToSortedBam(samPath, bPath, True)
samtools.index(bamPath)
r = {}
d = samtools.readBamToDict(bamPath, index=index, appendTo=r)
return(d)
# Reference to this function: http://stackoverflow.com/questions/845058/how-to-get-line-count-cheaply-in-python
# This method is the quickest way to count lines in big files. (Quicker than wc..)
def fileLineCount(filename):
f = open(filename, 'rb')
bufgen = takewhile(lambda x: x, (f.read(1024*1024) for _ in repeat(None)))
return sum( buf.count(b'\n') for buf in bufgen )
# Deprecated: should use Numpy sparse matrix
class DistanceMatrix():
def __init__(self, maxDist=1000000000):
self.data = {}
self.keys = set()
self.maxDist = maxDist
self.finalKeying = True # Do not store keys on every update
def setDistances(self,fromList, toList, distanceList ):
# Add keys:
if not self.finalKeying:
self.keys.update(fromList+toList)
for index,a in enumerate(fromList):
self._update(a,toList[index], distanceList[index])
def getNumpyArray(self):
#Make sure the keys are properly indexed
self._performKeying()
#Initialise empty matrix:
narray = np.zeros((len(self.keys), len(self.keys)))
for aIndex,keyA in enumerate(self.keys):
for bIndex,keyB in enumerate(self.keys):
narray[aIndex,bIndex] = self.getDistance(keyA,keyB)
return(narray)
def _update(self, a,b,distance):
if distance<=self.maxDist:
if (not a in self.data) and (not b in self.data):
self.data[a] = {}
self.data[a][b] = distance
return(True)
if a in self.data and b in self.data[a]:
self.data[a][b] = distance
return(True)
if b in self.data and a in self.data[b]:
self.data[b][a] = distance
return(True)
# A already exists but B does not: (And b also does not exists as main entry)
if a in self.data and not b in self.data and not b in self.data[a]:
self.data[a][b] = distance
return(True)
elif b in self.data:
if not a in self.data[b]:
self.data[b][a] = distance
return(True)
return(False)
def setDistance(self, a,b,distance ):
if not self.finalKeying:
if a not in self.keys and b not in self.keys:
self.keys.update([a,b])
else:
if a not in self.keys:
self.keys.add(a)
if b not in self.keys:
self.keys.add(b)
return( self._update(a,b, distance ))
def getDistance(self, a,b ):
if a in self.data and b in self.data[a]:
return(self.data[a][b])
if b in self.data and a in self.data[b]:
return(self.data[b][a])
return(self.maxDist)
def _performKeying(self):
#Retrieve all keys
self.keys = set()
addKeys = []
self.keys.update(list(self.data.keys()))
for a in self.data:
for b in self.data[a]:
if b not in self.keys:
addKeys.append(b)
self.keys.update(addKeys)
def writeToFile(self, path):
if self.finalKeying:
self._performKeying()
separator = ';'
with open(path, 'w') as f:
values= ['name']
for keyA in self.keys:
values.append(keyA)
f.write('%s\n' % (separator.join(values)))
for keyA in self.keys:
values = []
for keyB in self.keys:
values.append(str(self.getDistance(keyA, keyB)))
f.write('%s%s%s\n' % (keyA, separator,separator.join(values)))
def loadFromFile(self, path,separator=None):
with open(path) as f:
idx = 0
header = {}
for line in f:
parts = line.rstrip().split(separator)
if idx==0:
for keyId,keyName in enumerate(parts):
header[keyId] = keyName
else:
for partIndex, part in parts:
if partIndex==0:
keyA = part
else:
keyB=header[partIndex]
self.setDistance(keyA,keyB, float(part))
idx+=1
Datasets
Models
