3e05ac599db2f7a03948d96f5e24de16095fa70c ceisenhart Fri Apr 21 14:34:25 2017 -0700 Python script for converting an expression matrix into a barchart bed, refs #19155 diff --git src/utils/expMatrixToBarchartBed/expMatrixToBarchartBed src/utils/expMatrixToBarchartBed/expMatrixToBarchartBed new file mode 100755 index 0000000..bf12497 --- /dev/null +++ src/utils/expMatrixToBarchartBed/expMatrixToBarchartBed @@ -0,0 +1,236 @@ +#!/usr/bin/env python2.7 +# expMatrixToBarchartBed +""" +Generate a barChart bed6+2 file from a matrix, meta data, and coordinates. A bed file for each +sample type in the meta data file will be made in the directory where the script is run. These +bed files will be names <sample type>.bed and are bed6 format. +""" +import os +import sys +import argparse +import tempfile + +def parseArgs(args): + """ + Parse the command line arguments. + """ + parser= argparse.ArgumentParser(description = __doc__) + parser.add_argument ("metaFile", + help = " Two column no header, the first column is the samples which should match " + \ + "the matrix, the second is the grouping (cell type, tissue, etc)", + type = argparse.FileType("r")) + parser.add_argument ("matrixFile", + help = " The input matrix file, the first row should start with a a #." + \ + "the samples in the first row should exactly match the ones in " + \ + "the metaFile. The labels (ex ENST*****) in the first column should exactly match " + \ + "the ones in the coordinate file.", + type = argparse.FileType("r")) + parser.add_argument ("coordinateMap", + help = " Five column no header, maps the column labels from the matrix to coordinates. Tab " + \ + "separated; label, chr, strand, start coord, end coord. ", + action = "store") + parser.add_argument ("outputFile", + help = " The output file. ", + type =argparse.FileType("w")) + + parser.add_argument ("--verbose", + help = " Show runtime messages.", + action = "store_true") + + parser.set_defaults(verbose = False) + + if (len(sys.argv) == 1): + parser.print_help() + exit(1) + options = parser.parse_args() + return options + +def determineScore(tpmCuttoffs, tpm): + """ + Cast the tpm to a score between 0-1000. Since there are only + 9 visual blocks cast them to be in one of the 9 blocks. + tpmCuttoffs - A list of integers + tpm - An integer + """ + count = 0 + for val in tpmCuttoffs: + if (val > tpm): + return count*111 + count = count + 1 + return 999 + +def condenseMatrixIntoBed8Col(validTpms, sampleToGroup, matrix, bedLikeFile): + """ + Take an expression matrix and a dictionary that maps the samples to groups. + Go through the expression matrix and calculate the average for each group, outputting + it to an intermediate file as they are calculated. + validTpms - An empty list of integers. + sampleToGroup - A dictionary that maps string samples to string groups. + matrix - An expression matrix, samples are the x rows, transcripts the y rows. + bedLikeFile - An intermediate file, looks slightly like a bed. + """ + # Store some information on the bed file, most important is the order + # of the 8th column. + bedInfo = "" + firstLine = True + getBedInfo = True + # Use the first line of the matrix and the sampleToGroup dict to create a dictionary that maps + # the column to a group. + columnToGroup = dict() + + # Go through the matrix line by line. The first line is used to build an index mapping columns + # to group blocks, then for each line with TPM values merge the values based on group blocks. + for line in matrix: + splitLine = line.strip("\n").split() + + # The first line is the word 'transcript' followed by a list of the sample names. + if firstLine: + firstLine = False + count = 1 + firstCol = True + for col in splitLine: + if firstCol: + firstCol = False + continue + group = sampleToGroup[col] + columnToGroup.setdefault(count, group) + count += 1 + continue + + # Handle the tpm rows, calculating the average for each group per row. + groupAverages = dict() + groupCounts = dict() + firstCol = True + count = 1 + for col in splitLine: + if firstCol: + firstCol = False + continue + # First time this group is seen, add it to the groupCounts dict. + if (groupAverages.get(columnToGroup[count]) == None): + groupAverages.setdefault(columnToGroup[count], float(col)) + groupCounts.setdefault(columnToGroup[count], 1) + # This group has already been seen, update the TPM average. + else: + groupCounts[columnToGroup[count]] += 1 + normal = float(groupCounts[columnToGroup[count]]) + newTpm = (float(col) * (1/normal)) + oldTpm = (((normal - 1) / normal) * groupAverages[columnToGroup[count]]) + groupAverages[columnToGroup[count]] = newTpm + oldTpm + count += 1 + + # Store some information on the bed file. Most important is the groupOrder. + if getBedInfo: + getBedInfo = False + bedInfo += "#chr\tstart\tend\tname\tscore\tstrand\tgroupCount\tgroupOrder;" + for key, value in groupAverages.iteritems(): + bedInfo += key + "," + bedInfo = bedInfo[:-1] + "\toffset\tlineLength" + + # Write out the transcript name, this is needed to join with coordinates later. + bedLikeFile.write(splitLine[0] + "\t") + # Create a list of the average scores per group. + bedLine = "" + # The fullAverage is used to assign a tpm score representative of the entire bed row. + fullAverage = 0.0 + count = 0.0 + for key, value in groupAverages.iteritems(): + bedLine = bedLine + "," + str(value) + count += 1.0 + fullAverage += value + # Create what will be columns 5, 7 and 8 of the final bed. + bedLine = str(fullAverage/count) + "\t" + str(int(count)) + "\t" + bedLine[1:] + "\n" + # If the fullAverage tpm is greater than 0 then consider it in the validTpm list. + if (fullAverage > 0.0): + validTpms.append((fullAverage/count)) + # Write the bedLine to the intermediate bed-like file. + bedLikeFile.write(bedLine) + # Return the bedInfo so it can be printed right before the script ends. + return bedInfo + +def expMatrixToBarchartBed(options): + """ + Convert the expression matrix into a barchart bed file. + options - The command line options (file names, etc). + Use the meta data to map the sample names to their groups, then create a dict that maps + the columns to the groups. Go through the matrix line by line and get the average for each + group. Print this to an intermediate file, then use the unix 'join' command to link this with + the coordinates file, creating a bed like file. Go through this file once more to clean up a + bit, re arrange columns and calculate a new score creating a bed 6+2. Finally run + Max's bedJoinTabOffset to index the matrix, creating a bed 6+4 file. + """ + + # Create a dictionary that maps the sample names to their group. + sampleToGroup = dict() + for item in options.metaFile: + splitLine = item.strip("\n").split() + sampleToGroup.setdefault(splitLine[0], splitLine[1]) + + # Use an intermediate file to hold the average values for each group. + bedLikeFile = tempfile.NamedTemporaryFile( mode = "w+", bufsize = 1) + # Keep a list of TPM scores greater than 0. This will be used later + # to assign bed scores. + validTpms = [] + + # Go through the matrix and condense it into a bed like file. Populate + # the validTpms array and the bedInfo string. + bedInfo = condenseMatrixIntoBed8Col(validTpms, sampleToGroup, options.matrixFile, bedLikeFile) + # Find the number which divides the list of non 0 TPM scores into ten blocks. + tpmMedian = sorted(validTpms) + blockSizes = len(tpmMedian)/10 + + # Create a list of the ten TPM values at the edge of each block. + # These used to cast a TPM score to one of ten value between 0-1000. + tpmCuttoffs = [] + for i in range(1,10): + tpmCuttoffs.append(tpmMedian[blockSizes*i]) + + # Sort the bed like file to prepare it for the join. + sortedBedLikeFile = tempfile.NamedTemporaryFile( mode = "w+", bufsize = 1) + cmd = "sort " + bedLikeFile.name + " > " + sortedBedLikeFile.name + os.system(cmd) + + # Sort the coordinate file to prepare it for the join. + sortedCoords = tempfile.NamedTemporaryFile( mode = "w+", bufsize = 1) + cmd = "sort " + options.coordinateMap + " > " + sortedCoords.name + os.system(cmd) + + # Join the bed-like file and the coordinate file. + joinedFile = tempfile.NamedTemporaryFile(mode="w+", bufsize=1) + cmd = "join " + sortedCoords.name + " " + sortedBedLikeFile.name + " | awk -v " + \ + "OFS=\"\\t\" '$1=$1' > " + joinedFile.name + os.system(cmd) + + # Go through the joined file and re arrange the columns creating a bed 6+2 file. + # Also assign a scaled score 0 - 1000 to each tpm value. + bedFile = tempfile.NamedTemporaryFile(mode="w+", bufsize=1) + for line in joinedFile: + splitLine = line.strip("\n").split() + if ("_" in splitLine[0]): continue # Ignore alt sequences. + # Drop sequences where start is greater than end. + if (float(splitLine[3]) > float(splitLine[4])): + continue + newLine = (splitLine[1] + "\t" + str(splitLine[3]) + "\t" + str(splitLine[4]) + "\t" + \ + splitLine[0] + "\t" + str(determineScore(tpmCuttoffs, float(splitLine[5]))) + \ + "\t" + str(splitLine[2]) + "\t" + str(splitLine[6]) + "\t" + splitLine[7] + "\n") + bedFile.write(newLine) + + # Run Max's indexing script + cmd = "bedJoinTabOffset " + options.matrixFile.name + " " + bedFile.name + " " + options.outputFile.name + os.system(cmd) + + # Append the bed info to the end of the file. + cmd = "echo '" + bedInfo + "' >> " + options.outputFile.name + os.system(cmd) + + print ("The columns and order of the groups are; \n" + bedInfo) + +def main(args): + """ + Initialized options and calls other functions. + """ + options = parseArgs(args) + expMatrixToBarchartBed(options) + +if __name__ == "__main__" : + sys.exit(main(sys.argv))