|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +# Written by Lorena Pantano with subsequent reworking by Jonathan Manning. Released under the MIT license. |
| 4 | + |
| 5 | +import logging |
| 6 | +import argparse |
| 7 | +import glob |
| 8 | +import os |
| 9 | +import platform |
| 10 | +import re |
| 11 | +from collections import Counter, defaultdict, OrderedDict |
| 12 | +from collections.abc import Set |
| 13 | +from typing import Dict |
| 14 | + |
| 15 | +# Configure logging |
| 16 | +logging.basicConfig(format="%(name)s - %(asctime)s %(levelname)s: %(message)s") |
| 17 | +logger = logging.getLogger(__name__) |
| 18 | +logger.setLevel(logging.INFO) |
| 19 | + |
| 20 | +def format_yaml_like(data: dict, indent: int = 0) -> str: |
| 21 | + """Formats a dictionary to a YAML-like string. |
| 22 | +
|
| 23 | + Args: |
| 24 | + data (dict): The dictionary to format. |
| 25 | + indent (int): The current indentation level. |
| 26 | +
|
| 27 | + Returns: |
| 28 | + str: A string formatted as YAML. |
| 29 | + """ |
| 30 | + yaml_str = "" |
| 31 | + for key, value in data.items(): |
| 32 | + spaces = " " * indent |
| 33 | + if isinstance(value, dict): |
| 34 | + yaml_str += f"{spaces}{key}:\\n{format_yaml_like(value, indent + 1)}" |
| 35 | + else: |
| 36 | + yaml_str += f"{spaces}{key}: {value}\\n" |
| 37 | + return yaml_str |
| 38 | + |
| 39 | +def read_top_transcripts(quant_dir: str, file_pattern: str) -> Set[str]: |
| 40 | + """ |
| 41 | + Read the top 100 transcripts from the quantification file. |
| 42 | +
|
| 43 | + Parameters: |
| 44 | + quant_dir (str): Directory where quantification files are located. |
| 45 | + file_pattern (str): Pattern to match quantification files. |
| 46 | +
|
| 47 | + Returns: |
| 48 | + set: A set containing the top 100 transcripts. |
| 49 | + """ |
| 50 | + try: |
| 51 | + # Find the quantification file within the directory |
| 52 | + quant_file_path = glob.glob(os.path.join(quant_dir, "*", file_pattern))[0] |
| 53 | + with open(quant_file_path, "r") as file_handle: |
| 54 | + # Read the file and extract the top 100 transcripts |
| 55 | + return {line.split()[0] for i, line in enumerate(file_handle) if i > 0 and i <= 100} |
| 56 | + except IndexError: |
| 57 | + # Log an error and raise a FileNotFoundError if the quant file does not exist |
| 58 | + logger.error("No quantification files found.") |
| 59 | + raise FileNotFoundError("Quantification file not found.") |
| 60 | + |
| 61 | + |
| 62 | +def discover_transcript_attribute(gtf_file: str, transcripts: Set[str]) -> str: |
| 63 | + """ |
| 64 | + Discover the attribute in the GTF that corresponds to transcripts, prioritizing 'transcript_id'. |
| 65 | +
|
| 66 | + Parameters: |
| 67 | + gtf_file (str): Path to the GTF file. |
| 68 | + transcripts (Set[str]): A set of transcripts to match in the GTF file. |
| 69 | +
|
| 70 | + Returns: |
| 71 | + str: The attribute name that corresponds to transcripts in the GTF file. |
| 72 | + """ |
| 73 | + |
| 74 | + votes = Counter() |
| 75 | + with open(gtf_file) as inh: |
| 76 | + # Read GTF file, skipping header lines |
| 77 | + for line in filter(lambda x: not x.startswith("#"), inh): |
| 78 | + cols = line.split("\\t") |
| 79 | + |
| 80 | + # Use regular expression to correctly split the attributes string |
| 81 | + attributes_str = cols[8] |
| 82 | + attributes = dict(re.findall(r'(\\S+) "(.*?)(?<!\\\\)";', attributes_str)) |
| 83 | + |
| 84 | + votes.update(key for key, value in attributes.items() if value in transcripts) |
| 85 | + |
| 86 | + if not votes: |
| 87 | + # Error out if no matching attribute is found |
| 88 | + logger.error("No attribute in GTF matching transcripts") |
| 89 | + |
| 90 | + # Check if 'transcript_id' is among the attributes with the highest votes |
| 91 | + if "transcript_id" in votes and votes["transcript_id"] == max(votes.values()): |
| 92 | + logger.info("Attribute 'transcript_id' corresponds to transcripts.") |
| 93 | + return "transcript_id" |
| 94 | + |
| 95 | + # If 'transcript_id' isn't the highest, determine the most common attribute that matches the transcripts |
| 96 | + attribute, _ = votes.most_common(1)[0] |
| 97 | + logger.info(f"Attribute '{attribute}' corresponds to transcripts.") |
| 98 | + return attribute |
| 99 | + |
| 100 | + |
| 101 | +def parse_attributes(attributes_text: str) -> Dict[str, str]: |
| 102 | + """ |
| 103 | + Parse the attributes column of a GTF file. |
| 104 | +
|
| 105 | + :param attributes_text: The attributes column as a string. |
| 106 | + :return: A dictionary of the attributes. |
| 107 | + """ |
| 108 | + # Split the attributes string by semicolon and strip whitespace |
| 109 | + attributes = attributes_text.strip().split(";") |
| 110 | + attr_dict = OrderedDict() |
| 111 | + |
| 112 | + # Iterate over each attribute pair |
| 113 | + for attribute in attributes: |
| 114 | + # Split the attribute into key and value, ensuring there are two parts |
| 115 | + parts = attribute.strip().split(" ", 1) |
| 116 | + if len(parts) == 2: |
| 117 | + key, value = parts |
| 118 | + # Remove any double quotes from the value |
| 119 | + value = value.replace('"', "") |
| 120 | + attr_dict[key] = value |
| 121 | + |
| 122 | + return attr_dict |
| 123 | + |
| 124 | + |
| 125 | +def map_transcripts_to_gene( |
| 126 | + quant_type: str, gtf_file: str, quant_dir: str, gene_id: str, extra_id_field: str, output_file: str |
| 127 | +) -> bool: |
| 128 | + """ |
| 129 | + Map transcripts to gene names and write the output to a file. |
| 130 | +
|
| 131 | + Parameters: |
| 132 | + quant_type (str): The quantification method used (e.g., 'salmon'). |
| 133 | + gtf_file (str): Path to the GTF file. |
| 134 | + quant_dir (str): Directory where quantification files are located. |
| 135 | + gene_id (str): The gene ID attribute in the GTF file. |
| 136 | + extra_id_field (str): Additional ID field in the GTF file. |
| 137 | + output_file (str): The output file path. |
| 138 | +
|
| 139 | + Returns: |
| 140 | + bool: True if the operation was successful, False otherwise. |
| 141 | + """ |
| 142 | + # Read the top transcripts based on quantification type |
| 143 | + transcripts = read_top_transcripts(quant_dir, "quant.sf" if quant_type == "salmon" else "abundance.tsv") |
| 144 | + # Discover the attribute that corresponds to transcripts in the GTF |
| 145 | + transcript_attribute = discover_transcript_attribute(gtf_file, transcripts) |
| 146 | + |
| 147 | + # Open GTF and output file to write the mappings |
| 148 | + # Initialize the set to track seen combinations |
| 149 | + seen = set() |
| 150 | + |
| 151 | + with open(gtf_file) as inh, open(output_file, "w") as output_handle: |
| 152 | + # Parse each line of the GTF, mapping transcripts to genes |
| 153 | + for line in filter(lambda x: not x.startswith("#"), inh): |
| 154 | + cols = line.split("\\t") |
| 155 | + attr_dict = parse_attributes(cols[8]) |
| 156 | + if gene_id in attr_dict and transcript_attribute in attr_dict: |
| 157 | + # Create a unique identifier for the transcript-gene combination |
| 158 | + transcript_gene_pair = (attr_dict[transcript_attribute], attr_dict[gene_id]) |
| 159 | + |
| 160 | + # Check if the combination has already been seen |
| 161 | + if transcript_gene_pair not in seen: |
| 162 | + # If it's a new combination, write it to the output and add to the seen set |
| 163 | + extra_id = attr_dict.get(extra_id_field, attr_dict[gene_id]) |
| 164 | + output_handle.write(f"{attr_dict[transcript_attribute]}\\t{attr_dict[gene_id]}\\t{extra_id}\\n") |
| 165 | + seen.add(transcript_gene_pair) |
| 166 | + |
| 167 | + return True |
| 168 | + |
| 169 | + |
| 170 | +# Main function to parse arguments and call the mapping function |
| 171 | +if __name__ == "__main__": |
| 172 | + prefix = ( |
| 173 | + "$task.ext.prefix" |
| 174 | + if "$task.ext.prefix" != "null" |
| 175 | + else f"$meta.id" |
| 176 | + ) |
| 177 | + if not map_transcripts_to_gene('$quant_type', '$gtf', 'quants', '$id', '$extra', f"{prefix}.tx2gene.tsv"): |
| 178 | + logger.error("Failed to map transcripts to genes.") |
| 179 | + |
| 180 | + # Write the versions |
| 181 | + versions_this_module = {} |
| 182 | + versions_this_module["${task.process}"] = {"python": platform.python_version()} |
| 183 | + with open("versions.yml", "w") as f: |
| 184 | + f.write(format_yaml_like(versions_this_module)) |
0 commit comments