#!/usr/bin/env python # Copyright (c) 2015 Tobias Neumann, Philipp Rescheneder. # # This file is part of Slamdunk. # # Slamdunk is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as # published by the Free Software Foundation, either version 3 of the # License, or (at your option) any later version. # # Slamdunk is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . from __future__ import print_function import tempfile import math import pysam import os import glob import numpy as np from slamdunk.utils.misc import removeExtension, replaceExtension, checkStep, getSampleInfo, complement , getPlotter, callR, SlamSeqInfo # @UnresolvedImport from slamdunk.slamseq.SlamSeqFile import SlamSeqBamFile, ReadDirection # @UnresolvedImport from slamdunk.utils import SNPtools # @UnresolvedImport from slamdunk.utils.BedReader import BedIterator # @UnresolvedImport from slamdunk.version import __version__ # @UnresolvedImport utrNormFactor = 250 baseNumber = 5 toBase = [ 'A', 'C', 'G', 'T', 'N' ] def sumLists(a, b): return [int(x) + int(y) for x, y in zip(a, b)] def sumCounts(countFile, column="ReadCount"): columnId = -1 sum = 0 with open(countFile) as f: for line in f: if (not line.startswith("#")): columns = line.rstrip().split("\t") # Find column if line.startswith("Chromosome") : id = 0 for col in columns: if col == column: columnId = id id += 1 # Column not found if (columnId < 0) : return 0 else : sum += int(columns[columnId]) return sum # Print rates in correct format for plotting def printRates(ratesFwd, ratesRev, f): print("\t", end='', file=f) for i in range(0, 5): print(toBase[i] + "\t" + toBase[i].lower() + "\t", end='', file=f) print(file=f) for i in range(0, 5): print(toBase[i] + "\t", end='', file=f) for j in range(0, 5): print(str(ratesFwd[i * 5 + j]) + "\t" + str(ratesRev[i * 5 + j]) + "\t", end='', file=f) print(file=f) def statsComputeOverallRates(referenceFile, bam, minBaseQual, outputCSV, outputPDF, log, printOnly=False, verbose=True, force=False): if(not checkStep([bam, referenceFile], [outputCSV], force)): print("Skipped computing overall rates for file " + bam, file=log) else: # Init totalRatesFwd = [0] * 25 totalRatesRev = [0] * 25 tcCount = [0] * 100 # Go through one chr after the other testFile = SlamSeqBamFile(bam, referenceFile, None) chromosomes = testFile.getChromosomes() for chromosome in chromosomes: readIterator = testFile.readsInChromosome(chromosome, minBaseQual) for read in readIterator: # Compute rates for current read rates = read.conversionRates # Get T -> C conversions for current read tc = read.tcCount tcCount[tc] += 1 # Add rates from read to total rates if(read.direction == ReadDirection.Reverse): totalRatesRev = sumLists(totalRatesRev, rates) else: totalRatesFwd = sumLists(totalRatesFwd, rates) # Print rates in correct format for plotting fo = open(outputCSV, "w") print("# slamdunk rates v" + __version__, file=fo) printRates(totalRatesFwd, totalRatesRev, fo) fo.close() if(not checkStep([bam, referenceFile], [outputPDF], force)): print("Skipped computing overall rate pdfs for file " + bam, file=log) else: #f = tempfile.NamedTemporaryFile(delete=False) #print(removeExtension(basename(bam)), outputCSV, sep='\t', file=f) #f.close() callR(getPlotter("compute_overall_rates") + " -f " + outputCSV + " -n " + removeExtension(os.path.basename(bam)) + " -O " + outputPDF, log, dry=printOnly, verbose=verbose) # This is for TC conversions # def statsComputeTCContext(referenceFile, bam, minQual, outputCSV, outputPDF, log, printOnly=False, verbose=True, force=False): # # if(not checkStep([bam, referenceFile], [outputCSV], force)): # print("Skipped computing overall rates for file " + bam, file=log) # else: # #Init # #combinations = ["AT","CT","GT","TT","NT","AA","CA","GA","TA","NA"] # frontCombinations = ["AT","CT","GT","TT","NT"] # backCombinations = ["TA","TC","TG","TT","TN"] # # counts = {} # counts['5prime']= {} # counts['3prime']= {} # counts['5prime']['fwd'] = {} # counts['5prime']['rev'] = {} # counts['3prime']['fwd'] = {} # counts['3prime']['rev'] = {} # # for combination in frontCombinations : # counts['5prime']['fwd'][combination] = 0 # counts['5prime']['rev'][combination] = 0 # # for combination in backCombinations: # counts['3prime']['fwd'][combination] = 0 # counts['3prime']['rev'][combination] = 0 # # #Go through one chr after the other # testFile = SlamSeqBamFile(bam, referenceFile, None) # # chromosomes = testFile.getChromosomes() # # for chromosome in chromosomes: # readIterator = testFile.readsInChromosome(chromosome) # # for read in readIterator: # # # Iterator over mismatches for T -> C conversion contexts # # Get mismatches rates for current read # for mismatch in read.mismatches: # if(mismatch.isTCMismatch(read.direction == ReadDirection.Reverse)): # frontContext = mismatch.referenceBase5PrimeContext + mismatch.referenceBase # backContext = mismatch.referenceBase + mismatch.referenceBase3PrimeContext # # if (read.direction == ReadDirection.Reverse) : # frontContext = complement(frontContext) # backContext = complement(backContext) # counts['5prime']['rev'][frontContext] += 1 # counts['3prime']['rev'][backContext] += 1 # else : # counts['5prime']['fwd'][frontContext] += 1 # counts['3prime']['fwd'][backContext] += 1 # # #Print rates in correct format for plotting # fo = open(outputCSV, "w") # # print("\t".join(frontCombinations),file=fo) # # frontFwdLine = "" # frontRevLine = "" # backFwdLine = "" # backRevLine = "" # # for combination in frontCombinations : # frontFwdLine += str(counts['5prime']['fwd'][combination]) + "\t" # frontRevLine += str(counts['5prime']['rev'][combination]) + "\t" # # print(frontFwdLine.rstrip(),file=fo) # print(frontRevLine.rstrip(),file=fo) # # print("\t".join(backCombinations),file=fo) # # for combination in backCombinations : # backFwdLine += str(counts['3prime']['fwd'][combination]) + "\t" # backRevLine += str(counts['3prime']['rev'][combination]) + "\t" # # print(backFwdLine.rstrip(),file=fo) # print(backRevLine.rstrip(),file=fo) # # fo.close() # # if(not checkStep([bam, referenceFile], [outputPDF], force)): # print("Skipped computing overall rate pdfs for file " + bam, file=log) # else: # f = tempfile.NamedTemporaryFile(delete=False) # print(removeExtension(basename(bam)), outputCSV, sep='\t', file=f) # f.close() # # run(pathComputeTCContext + " -f " + f.name + " -O " + outputPDF, log, dry=printOnly, verbose=verbose) # This is for Ts in reads def statsComputeTCContext(referenceFile, bam, minBaseQual, outputCSV, outputPDF, log, printOnly=False, verbose=True, force=False): if(not checkStep([bam, referenceFile], [outputCSV], force)): print("Skipped computing overall rates for file " + bam, file=log) else: # Init # combinations = ["AT","CT","GT","TT","NT","AA","CA","GA","TA","NA"] frontCombinations = ["AT", "CT", "GT", "TT", "NT"] backCombinations = ["TA", "TC", "TG", "TT", "TN"] counts = {} counts['5prime'] = {} counts['3prime'] = {} counts['5prime']['fwd'] = {} counts['5prime']['rev'] = {} counts['3prime']['fwd'] = {} counts['3prime']['rev'] = {} for combination in frontCombinations : counts['5prime']['fwd'][combination] = 0 counts['5prime']['rev'][combination] = 0 for combination in backCombinations: counts['3prime']['fwd'][combination] = 0 counts['3prime']['rev'][combination] = 0 bamFile = pysam.AlignmentFile(bam, "rb") # Go through one chr after the other testFile = SlamSeqBamFile(bam, referenceFile, None) chromosomes = testFile.getChromosomes() for chromosome in chromosomes: for read in bamFile.fetch(region=chromosome): i = 0 while i < len(read.query_sequence): if(read.query_sequence[i] == "T" and not read.is_reverse) : frontContext = None backContext = None if (i > 0) : frontContext = read.query_sequence[i - 1] if (i < (len(read.query_sequence) - 1)) : backContext = read.query_sequence[i + 1] if (frontContext != None) : counts['5prime']['fwd'][frontContext + "T"] += 1 if (backContext != None) : counts['3prime']['fwd']["T" + backContext] += 1 if(read.query_sequence[i] == "A" and read.is_reverse) : frontContext = None backContext = None if (i > 0) : backContext = read.query_sequence[i - 1] if (i < (len(read.query_sequence) - 1)) : frontContext = read.query_sequence[i + 1] if (frontContext != None) : counts['5prime']['rev'][complement(frontContext + "A")] += 1 if (backContext != None) : counts['3prime']['rev'][complement("A" + backContext)] += 1 i += 1 # Print rates in correct format for plotting fo = open(outputCSV, "w") print("\t".join(frontCombinations), file=fo) frontFwdLine = "" frontRevLine = "" backFwdLine = "" backRevLine = "" for combination in frontCombinations : frontFwdLine += str(counts['5prime']['fwd'][combination]) + "\t" frontRevLine += str(counts['5prime']['rev'][combination]) + "\t" print(frontFwdLine.rstrip(), file=fo) print(frontRevLine.rstrip(), file=fo) print("\t".join(backCombinations), file=fo) for combination in backCombinations : backFwdLine += str(counts['3prime']['fwd'][combination]) + "\t" backRevLine += str(counts['3prime']['rev'][combination]) + "\t" print(backFwdLine.rstrip(), file=fo) print(backRevLine.rstrip(), file=fo) fo.close() if(not checkStep([bam, referenceFile], [outputPDF], force)): print("Skipped computing overall rate pdfs for file " + bam, file=log) else: f = tempfile.NamedTemporaryFile(mode='w',delete=False) print(removeExtension(os.path.basename(bam)), outputCSV, sep='\t', file=f) f.close() callR(getPlotter("compute_context_TC_rates") + " -f " + f.name + " -O " + outputPDF, log, dry=printOnly, verbose=verbose) def statsComputeOverallRatesPerUTR(referenceFile, bam, minBaseQual, strictTCs, outputCSV, outputPDF, utrBed, maxReadLength, log, printOnly=False, verbose=True, force=False): sampleInfo = getSampleInfo(bam) slamseqInfo = SlamSeqInfo(bam) if(not checkStep([bam, referenceFile], [outputCSV], force)): print("Skipped computing overall rates for file " + bam, file=log) else: # Go through one chr after the other testFile = SlamSeqBamFile(bam, referenceFile, None) # UTR stats for MultiQC utrStats = dict() plotConversions = ['A>T', 'A>G', 'A>C', 'C>A', 'C>G', 'C>T', 'G>A', 'G>C', 'G>T', 'T>A', 'T>G', 'T>C', ] for conversion in plotConversions: utrStats[conversion] = list() f = tempfile.NamedTemporaryFile(mode='w', delete=False) for utr in BedIterator(utrBed): readIterator = testFile.readInRegion(utr.chromosome, utr.start, utr.stop, utr.strand, maxReadLength, minBaseQual) # Init totalRates = [0] * 25 readCount = 0 for read in readIterator: if (not read.isTcRead and strictTCs and read.tcCount > 0) : pass else : # Compute rates for current read rates = read.conversionRates # Add rates from read to total rates totalRates = sumLists(totalRates, rates) readCount += 1 print(utr.name, utr.chromosome, utr.start, utr.stop, utr.strand, readCount, "\t".join(str(x) for x in totalRates), sep="\t", file=f) # Process rates for MultiQC # Copied directly, too lazy to do it properly now utrDict = {} conversionSum = 0 A_A = totalRates[0] conversionSum =+ A_A A_C = totalRates[1] conversionSum =+ A_C A_G = totalRates[2] conversionSum =+ A_G A_T = totalRates[3] conversionSum =+ A_T C_A = totalRates[5] conversionSum =+ C_A C_C = totalRates[6] conversionSum =+ C_C C_G = totalRates[7] conversionSum =+ C_G C_T = totalRates[8] conversionSum =+ C_T G_A = totalRates[10] conversionSum =+ G_A G_C = totalRates[11] conversionSum =+ G_C G_G = totalRates[12] conversionSum =+ G_G G_T = totalRates[13] conversionSum =+ G_T T_A = totalRates[15] conversionSum =+ T_A T_C = totalRates[16] conversionSum =+ T_C T_G = totalRates[17] conversionSum =+ T_G T_T = totalRates[18] conversionSum =+ T_T if utr.strand == "-": A_A, T_T = T_T,A_A G_G, C_C = C_C,G_G A_C, T_G = T_G, A_C A_G, T_C = T_C, A_G A_T, T_A = T_A, A_T C_A, G_T = G_T, C_A C_G, G_C = G_C, C_G C_T, G_A = G_A, C_T if conversionSum > 0: Asum = A_A + A_C + A_G + A_T Csum = C_A + C_C + C_G + C_T Gsum = G_A + G_C + G_G + G_T Tsum = T_A + T_C + T_G + T_T if Asum > 0 : A_T = A_T / float(Asum) * 100 A_G = A_G / float(Asum) * 100 A_C = A_C / float(Asum) * 100 else : A_T = 0 A_G = 0 A_C = 0 if Csum > 0: C_A = C_A / float(Csum) * 100 C_G = C_G / float(Csum) * 100 C_T = C_T / float(Csum) * 100 else : C_A = 0 C_G = 0 C_T = 0 if Gsum > 0: G_A = G_A / float(Gsum) * 100 G_C = G_C / float(Gsum) * 100 G_T = G_T / float(Gsum) * 100 else : G_A = 0 G_C = 0 G_T = 0 if Tsum > 0: T_A = T_A / float(Tsum) * 100 T_G = T_G / float(Tsum) * 100 T_C = T_C / float(Tsum) * 100 else : T_A = 0 T_G = 0 T_C = 0 utrStats['A>T'].append(A_T) utrStats['A>G'].append(A_G) utrStats['A>C'].append(A_C) utrStats['C>A'].append(C_A) utrStats['C>G'].append(C_G) utrStats['C>T'].append(C_T) utrStats['G>A'].append(G_A) utrStats['G>T'].append(G_T) utrStats['G>C'].append(G_C) utrStats['T>A'].append(T_A) utrStats['T>G'].append(T_G) utrStats['T>C'].append(T_C) f.close() fo = open(outputCSV, "w") print("# slamdunk utrrates v" + __version__, file=fo) print("# Median-Conversions=",end="",file=fo) first = True for conversion in plotConversions: if (not first) : print(',',file=fo, end="") else : first = False print(conversion + ":" + str(np.median(utrStats[conversion])),file=fo, end="") print(file=fo) print("Name", "Chr", "Start", "End", "Strand", "ReadCount", sep="\t", end="\t", file=fo) for i in range(0, 5): for j in range(0, 5): print(toBase[i].upper() + "_" + toBase[j].upper(), end="", file=fo) if(i != 4 or j != 4): print("\t", end="", file=fo) print(file=fo) with open(f.name, "r") as valueFile: fo.write(valueFile.read()) fo.close() if(not checkStep([bam, referenceFile], [outputPDF], force)): print("Skipped computing global rate pdfs for file " + bam, file=log) else: f = tempfile.NamedTemporaryFile(mode='w',delete=False) print(sampleInfo.Name, outputCSV, sep='\t', file=f) f.close() callR(getPlotter("globalRatePlotter") + " -f " + f.name + " -O " + outputPDF, log, dry=printOnly, verbose=verbose) def readSummary(filteredFiles, countDirectory, outputFile, log, printOnly=False, verbose=True, force=False): # Print sort by ID contentDict = {} tsvFile = open(outputFile, "w") print("# slamdunk summary v" + __version__, file=tsvFile) if (countDirectory != None) : f = tempfile.NamedTemporaryFile(mode='w', delete=False) for bam in filteredFiles: slamseqInfo = SlamSeqInfo(bam) sampleInfo = getSampleInfo(bam) if (countDirectory != None) : countedReads = 0 countFile = os.path.join(countDirectory, replaceExtension(os.path.basename(bam), ".tsv", "_tcount")) if not os.path.exists(countFile): print("TCount directory does not seem to contain tcount file for:\t" + countFile) else : print(sampleInfo.Name, countFile, sep='\t', file=f) countedReads = sumCounts(countFile) if(int(sampleInfo.ID) in contentDict): ID = len(contentDict) + 1 else: ID = sampleInfo.ID contentDict[int(ID)] = "\t".join([bam, sampleInfo.Name, sampleInfo.Type, sampleInfo.Time, str(slamseqInfo.SequencedReads), str(slamseqInfo.MappedReads), str(slamseqInfo.DedupReads), str(slamseqInfo.MQFilteredReads), str(slamseqInfo.IdFilteredReads), str(slamseqInfo.NmFilteredReads), str(slamseqInfo.MultimapperReads), str(slamseqInfo.FilteredReads), str(countedReads), slamseqInfo.AnnotationName]) else : if(int(sampleInfo.ID) in contentDict): ID = len(contentDict) + 1 else: ID = sampleInfo.ID contentDict[int(ID)] = "\t".join([bam, sampleInfo.Name, sampleInfo.Type, sampleInfo.Time, str(slamseqInfo.SequencedReads), str(slamseqInfo.MappedReads), str(slamseqInfo.DedupReads), str(slamseqInfo.MQFilteredReads), str(slamseqInfo.IdFilteredReads), str(slamseqInfo.NmFilteredReads), str(slamseqInfo.MultimapperReads), str(slamseqInfo.FilteredReads), slamseqInfo.AnnotationName]) if (countDirectory != None) : f.close() callR(getPlotter("PCAPlotter") + " -f " + f.name + " -O " + replaceExtension(outputFile, ".pdf", "_PCA") + " -P " + replaceExtension(outputFile, ".txt", "_PCA"), log, dry=printOnly, verbose=verbose) print("FileName", "SampleName", "SampleType", "SampleTime", "Sequenced", "Mapped", "Deduplicated", "MQ-Filtered", "Identity-Filtered", "NM-Filtered", "Multimap-Filtered", "Retained", "Counted", "Annotation", sep="\t", file=tsvFile) else : print("FileName", "SampleName", "SampleType", "SampleTime", "Sequenced", "Mapped", "Deduplicated", "MQ-Filtered", "Identity-Filtered", "NM-Filtered", "Multimap-Filtered", "Retained", "Annotation", sep="\t", file=tsvFile) for key in sorted(contentDict): print(contentDict[key], file=tsvFile) tsvFile.close() def tcPerReadPos(referenceFile, bam, minQual, maxReadLength, outputCSV, outputPDF, snpsFile, log, printOnly=False, verbose=True, force=False): if(not checkStep([bam, referenceFile], [outputCSV], force)): print("Skipped computing T->C per reads position for file " + bam, file=log) else: totalReadCountFwd = [0] * maxReadLength totalReadCountRev = [0] * maxReadLength tcPerPosRev = [0] * maxReadLength tcPerPosFwd = [0] * maxReadLength allPerPosRev = [0] * maxReadLength allPerPosFwd = [0] * maxReadLength snps = SNPtools.SNPDictionary(snpsFile) snps.read() # Go through one chr after the other testFile = SlamSeqBamFile(bam, referenceFile, snps) chromosomes = testFile.getChromosomes() for chromosome in chromosomes: readIterator = testFile.readsInChromosome(chromosome, minQual) for read in readIterator: tcCounts = [0] * maxReadLength mutCounts = [0] * maxReadLength for mismatch in read.mismatches: if(mismatch.isTCMismatch(read.direction == ReadDirection.Reverse)): tcCounts[mismatch.readPosition] += 1 else : mutCounts[mismatch.readPosition] += 1 query_length = len(read.sequence) if(read.direction == ReadDirection.Reverse): tcPerPosRev = sumLists(tcPerPosRev, tcCounts) allPerPosRev = sumLists(allPerPosRev, mutCounts) for i in range(0, query_length): totalReadCountRev[i] += 1 else: tcPerPosFwd = sumLists(tcPerPosFwd, tcCounts) allPerPosFwd = sumLists(allPerPosFwd, mutCounts) for i in range(0, query_length): totalReadCountFwd[i] += 1 foTC = open(outputCSV, "w") print("# slamdunk tcperreadpos v" + __version__, file=foTC) for i in range(0, maxReadLength): print(allPerPosFwd[i], allPerPosRev[i], tcPerPosFwd[i], tcPerPosRev[i], totalReadCountFwd[i], totalReadCountRev[i], sep='\t', file=foTC) foTC.close() if(not checkStep([outputCSV], [outputPDF], force)): print("Skipped computing T->C per reads position plot for file " + bam, file=log) else: callR(getPlotter("conversion_per_read_position") + " -i " + outputCSV + " -o " + outputPDF, log, dry=printOnly, verbose=verbose) def tcPerUtr(referenceFile, utrBed, bam, minQual, maxReadLength, outputCSV, outputPDF, snpsFile, log, printOnly=False, verbose=True, force=False): if(not checkStep([bam, referenceFile], [outputCSV], force)): print("Skipped computing T->C per UTR position for file " + bam, file=log) else: counter = 0 totalUtrCountFwd = [0] * utrNormFactor totalUtrCountRev = [0] * utrNormFactor tcPerPosRev = [0] * utrNormFactor tcPerPosFwd = [0] * utrNormFactor allPerPosRev = [0] * utrNormFactor allPerPosFwd = [0] * utrNormFactor snps = SNPtools.SNPDictionary(snpsFile) snps.read() # Go through one utr after the other testFile = SlamSeqBamFile(bam, referenceFile, snps) for utr in BedIterator(utrBed): readIterator = testFile.readInRegion(utr.chromosome, utr.start, utr.stop, utr.strand, maxReadLength, minQual) tcForwardCounts = [0] * utrNormFactor mutForwardCounts = [0] * utrNormFactor tcReverseCounts = [0] * utrNormFactor mutReverseCounts = [0] * utrNormFactor for read in readIterator: tcCounts = [0] * utrNormFactor mutCounts = [0] * utrNormFactor for mismatch in read.mismatches: mismatchPos = mismatch.referencePosition # mismatchPos = read.startRefPos if (utr.strand == "+") : # New try for UTRs (remove + 1 if (mismatchPos >= (utr.getLength() - utrNormFactor) and mismatchPos < utr.getLength()) : # if (mismatchPos >= (utr.getLength() - utrNormFactor) and mismatchPos < utr.getLength() + 1) : mismatchPos = utrNormFactor - (utr.getLength() - mismatchPos) if(mismatch.isTCMismatch(read.direction == ReadDirection.Reverse)): tcCounts[mismatchPos] += 1 else : mutCounts[mismatchPos] += 1 else : if (mismatchPos >= 0 and mismatchPos < min(utr.getLength(), utrNormFactor)) : if(mismatch.isTCMismatch(read.direction == ReadDirection.Reverse)): tcCounts[mismatchPos] += 1 else : mutCounts[mismatchPos] += 1 if(read.direction == ReadDirection.Reverse): tcReverseCounts = sumLists(tcReverseCounts, tcCounts) mutReverseCounts = sumLists(mutReverseCounts, mutCounts) start = max(0, min(min(utr.getLength(), utrNormFactor), read.startRefPos)) end = max(0, min(min(utr.getLength(), utrNormFactor), read.endRefPos)) for i in range(start, end): totalUtrCountRev[i] += 1 else: tcForwardCounts = sumLists(tcForwardCounts, tcCounts) mutForwardCounts = sumLists(mutForwardCounts, mutCounts) start = min(utr.getLength(), max(utr.getLength() - utrNormFactor, read.startRefPos)) end = min(utr.getLength(), max(utr.getLength() - utrNormFactor, read.endRefPos)) for i in range(start, end): normPos = utrNormFactor - (utr.getLength() - i) totalUtrCountFwd[normPos] += 1 tcPerPosFwd = sumLists(tcPerPosFwd, tcForwardCounts) allPerPosFwd = sumLists(allPerPosFwd, mutForwardCounts) tcPerPosRev = sumLists(tcPerPosRev, tcReverseCounts) allPerPosRev = sumLists(allPerPosRev, mutReverseCounts) counter += 1 if (verbose and counter % 10000 == 0) : print("Handled " + str(counter) + " UTRs.", file=log) foTC = open(outputCSV, "w") print("# slamdunk tcperutr v" + __version__, file=foTC) reverseAllPerPosRev = allPerPosRev[::-1] reverseTcPerPosRev = tcPerPosRev[::-1] reverseTotalUtrCountRev = totalUtrCountRev[::-1] for i in range(0, utrNormFactor): print(allPerPosFwd[i], reverseAllPerPosRev[i], tcPerPosFwd[i], reverseTcPerPosRev[i], totalUtrCountFwd[i], reverseTotalUtrCountRev[i], sep='\t', file=foTC) foTC.close() if(not checkStep([outputCSV], [outputPDF], force)): print("Skipped computing T->C per UTR position plot for file " + bam, file=log) else: callR(getPlotter("conversion_per_read_position") + " -u -i " + outputCSV + " -o " + outputPDF, log, dry=printOnly, verbose=verbose) def computeSNPMaskedRates (ref, bed, snpsFile, bam, maxReadLength, minQual, coverageCutoff, variantFraction, outputCSV, outputPDF, strictTCs, log, printOnly=False, verbose=True, force=False): if(not checkStep([bam, ref], [outputCSV], force)): print("Skipped computing T->C per UTR with SNP masking for file " + bam, file=log) else: fileCSV = open(outputCSV,'w') snps = SNPtools.SNPDictionary(snpsFile) snps.read() #Go through one chr after the other testFile = SlamSeqBamFile(bam, ref, snps) progress = 0 for utr in BedIterator(bed): if(not utr.hasStrand()): raise RuntimeError("Input BED file does not contain stranded intervals.") if utr.start < 0: raise RuntimeError("Negativ start coordinate found. Please check the following entry in your BED file: " + utr) readIterator = testFile.readInRegion(utr.chromosome, utr.start, utr.stop, utr.strand, maxReadLength, minQual) unmaskedTCCount = 0 maskedTCCount = 0 readCount = 0 for read in readIterator: # Overwrite any conversions for non-TC reads (reads with < 2 TC conversions) if (not read.isTcRead and strictTCs) : read.tcCount = 0 read.mismatches = [] read.conversionRates = 0.0 read.tcRate = 0.0 isTC = False isTrueTC = False for mismatch in read.mismatches: if(mismatch.isTCMismatch(read.direction == ReadDirection.Reverse) and mismatch.referencePosition >= 0 and mismatch.referencePosition < utr.getLength()): isTrueTC = True unmasked = False if (read.direction == ReadDirection.Reverse and mismatch.referenceBase == "A" and mismatch.readBase == "G"): unmasked = True elif (read.direction != ReadDirection.Reverse and mismatch.referenceBase == "T" and mismatch.readBase == "C") : unmasked = True if (unmasked and mismatch.referencePosition >= 0 and mismatch.referencePosition < utr.getLength()) : isTC = True readCount += 1 if (isTC) : unmaskedTCCount += 1 if (isTrueTC) : maskedTCCount += 1 containsSNP = 0 if (unmaskedTCCount != maskedTCCount) : containsSNP = 1 print(utr.name + "\t" + str(readCount) + "\t" + str(unmaskedTCCount) + "\t" + str(maskedTCCount) + "\t" + str(containsSNP), file=fileCSV) progress += 1 fileCSV.close() if(not checkStep([outputCSV], [outputPDF], force)): print("Skipped computing T->C per UTR position plot for file " + bam, file=log) else: callR(getPlotter("SNPeval") + " -i " + outputCSV + " -c " + str(coverageCutoff) + " -v " + str(variantFraction) + " -o " + outputPDF, log, dry=printOnly, verbose=verbose) def halflifes(bams, outputCSV, timepoints, log, printOnly=False, verbose=True, force=False): callR(getPlotter("compute_halflifes") + " -f " + bams + " -t " + timepoints + " -o " + outputCSV, log, dry=printOnly, verbose=verbose) def mergeRates(bams, outputCSV, column, columnName, log, printOnly=False, verbose=True, force=False): cmd = getPlotter("merge_rate_files") + " -f " + bams + " -o " + outputCSV cmd = cmd + " -c \"" + column + "\"" cmd = cmd + " -n " + str(columnName) callR(cmd, log, dry=printOnly, verbose=verbose)