# See the NOTICE file distributed with this work for additional information
# regarding copyright ownership.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Set of discriminant functions that can recognize structural and compositional features
such as CpG islands, promoter regions and first splice-donor sites.
References
----------
:cite:`cpg`
"""
__all__ = ["run_cpg"]
import argparse
import logging
import logging.config
import multiprocessing
from os import PathLike
from pathlib import Path
import re
import subprocess
from typing import List, Union
from ensembl.tools.anno.utils._utils import (
check_exe,
create_dir,
check_gtf_content,
get_seq_region_length,
get_slice_id,
slice_output_to_gtf,
get_sequence,
)
logger = logging.getLogger(__name__)
[docs]
def run_cpg( # pylint:disable=too-many-arguments, too-many-positional-arguments, too-many-locals, too-many-branches
genome_file: PathLike,
output_dir: Path,
cpg_bin: Path = Path("cpg_lh"),
cpg_min_length: int = 400,
cpg_min_gc_content: int = 50,
cpg_min_oe: float = 0.6,
num_threads: int = 1,
bedtools_bin: str = "/hps/software/users/ensembl/ensw/C8-MAR21-sandybridge/linuxbrew/bin/bedtools", # pylint:disable=line-too-long
) -> None:
"""
Run CpG islands on genomic slices
:param genome_file: Genome file path.
:type genome_file: PathLike
:param output_dir: Working directory path
:type output_dir: Path
:param cpg_bin: CpG software path.
:type cpg_bin: Path
:param cpg_min_length: Min length of CpG islands
:type cpg_min_length: int
:param cpg_min_gc_content: Min GC frequency percentage
:type cpg_min_gc_content: int
:param cpg_min_oe: Min ratio of the observed to expected number of CpG (CpGo/e)
:type cpg_min_oe: float
:param num_threads: int, number of threads.
:type num_threads: int
:param bedtools_bin: Bedtools executable path.
:type bedtools_bin: str
:return: None
:rtype: None
"""
# Use default path if user didn't supply one
cpg_bin = cpg_bin or Path("cpg_lh")
check_exe(cpg_bin)
cpg_dir = create_dir(output_dir, "cpg_output")
output_file = cpg_dir / "annotation.gtf"
if output_file.exists():
transcript_count = check_gtf_content(output_file, "simple_feature")
if transcript_count > 0:
logger.info("Cpg gtf file exists")
return
logger.info("Creating list of genomic slices")
seq_region_to_length = get_seq_region_length(genome_file, 5000)
slice_ids_per_region = get_slice_id(
seq_region_to_length, slice_size=1000000, overlap=0, min_length=5000
) # pylint:disable=line-too-long
logger.info("Running CpG")
pool = multiprocessing.Pool(int(num_threads)) # pylint:disable=consider-using-with
for slice_id in slice_ids_per_region:
pool.apply_async(
_multiprocess_cpg,
args=(
cpg_bin,
slice_id,
genome_file,
cpg_dir,
cpg_min_length,
cpg_min_gc_content,
cpg_min_oe,
bedtools_bin,
),
)
pool.close()
pool.join()
slice_output_to_gtf(cpg_dir, "feature_id", "cpg", True, ".cpg.gtf")
for gtf_file in cpg_dir.glob("*.cpg.gtf"):
gtf_file.unlink()
def _multiprocess_cpg( # pylint:disable=too-many-arguments, too-many-locals, too-many-branches, too-many-positional-arguments
cpg_bin: Path,
slice_id: List[str],
genome_file: Path,
cpg_dir: Path,
cpg_min_length: int = 400,
cpg_min_gc_content: int = 50,
cpg_min_oe: float = 0.6,
bedtools_bin: str = "/hps/software/users/ensembl/ensw/C8-MAR21-sandybridge/linuxbrew/bin/bedtools", # pylint:disable=line-too-long
) -> None:
"""
Annotation of CpG islands on multiprocess on genomic slices
Args:
cpg_bin: CpG software path.
slice_id: Slice id to run CpG on.
genome_file : Genome file.
cpg_dir : Output dir.
cpg_min_length : Min length of CpG islands
cpg_min_gc_content : Min GC frequency percentage
cpg_min_oe : Min ratio of the observed to expected number of CpG (CpGo/e)
bedtools_bin: Bedtools executable path.
"""
region_name, start, end = slice_id
logger.info(
"Processing slice to find CpG islands with cpg_lh: %s:%s:%s",
region_name,
start,
end,
)
seq = get_sequence(region_name, int(start), int(end), 1, genome_file, cpg_dir, bedtools_bin)
slice_name = f"{region_name}.rs{start}.re{end}"
# with TemporaryDirectory(dir=cpg_dir) as tmpdirname:
slice_file = cpg_dir / f"{slice_name}.fa"
with open(slice_file, "w+", encoding="utf8") as region_out:
region_out.write(f">{region_name}\n{seq}\n")
region_results = cpg_dir / f"{slice_file}.cpg.gtf"
output_file = Path(f"{slice_file}.cpg")
cpg_cmd = [str(cpg_bin), str(slice_file)]
with open(output_file, "w+", encoding="utf8") as cpg_out:
subprocess.run(cpg_cmd, stdout=cpg_out, check=True)
_create_cpg_gtf(
output_file,
region_results,
region_name,
cpg_min_length,
cpg_min_gc_content,
cpg_min_oe,
)
slice_file.unlink()
output_file.unlink()
def _create_cpg_gtf( # pylint:disable=too-many-arguments, too-many-locals, too-many-branches, too-many-positional-arguments
output_file: Path,
region_results: Path,
region_name: str,
cpg_min_length: int = 400,
cpg_min_gc_content: int = 50,
cpg_min_oe: float = 0.6,
) -> None:
"""
Read the fasta file and save the content in gtf format
All the genomic slices are collected in a single gtf output
Args:
output_file : GTF file with final results.
region_results : GTF file with the results per region.
region_name :Coordinates of genomic slice.
cpg_dir : Output dir.
cpg_min_length : Min length of CpG islands
cpg_min_gc_content : Min GC frequency percentage
cpg_min_oe : Min ratio of the observed to expected number of CpG (CpGo/e)
"""
with (
open(output_file, "r", encoding="utf8") as cpg_in,
open(region_results, "w+", encoding="utf8") as cpg_out,
):
feature_count = 1
for line in cpg_in:
result_match = re.search(r"^" + region_name, line)
if result_match:
results = line.split()
start = int(results[1])
end = int(results[2])
length = end - start + 1
score = float(results[3])
gc_content = float(results[6])
oe_score_str = results[7]
oe_score: Union[float, int]
if oe_score_str in ("-", "inf"):
oe_score = 0
else:
oe_score = float(oe_score_str)
if (
int(length) >= int(cpg_min_length)
and gc_content >= int(cpg_min_gc_content)
and oe_score >= float(cpg_min_oe)
):
gtf_line = (
f"{region_name}\tCpG\tsimple_feature\t{start}\t"
f'{end}\t.\t+\t.\tfeature_id "{feature_count}"; score "{score}";\n'
)
cpg_out.write(gtf_line)
def parse_args():
"""Parse command line arguments."""
parser = argparse.ArgumentParser(description="CpG's arguments")
parser.add_argument("--genome_file", required=True, help="Genome file path")
parser.add_argument("--output_dir", required=True, help="Output directory path")
parser.add_argument("--cpg_bin", default="cpg_lh", help="CpG executable path")
parser.add_argument("--cpg_min_length", type=int, default=400, help="Min length of CpG islands")
parser.add_argument(
"--cpg_min_gc_content", type=int, default=50, help="Min GC frequency percentage"
) # pylint:disable=line-too-long
parser.add_argument(
"--cpg_min_oe",
type=float,
default=0.6,
help="Min ratio of the observed to expected number of CpG (CpGo/e)",
)
parser.add_argument("--num_threads", type=int, default=1, help="Number of threads")
parser.add_argument(
"--bedtools_bin",
default="/hps/software/users/ensembl/ensw/C8-MAR21-sandybridge/linuxbrew/bin/bedtools",
help="Bedtools executable path",
)
return parser.parse_args()
def main():
"""CpG's entry-point."""
args = parse_args()
log_file_path = create_dir(args.output_dir, "log") / "cpg.log"
loginipath = Path(__file__).parents[6] / "conf" / "logging.conf"
logging.config.fileConfig(
loginipath,
defaults={"logfilename": str(log_file_path)},
disable_existing_loggers=False,
)
run_cpg(
args.genome_file,
args.output_dir,
args.cpg_bin,
args.cpg_min_length,
args.cpg_min_gc_content,
args.cpg_min_oe,
args.num_threads,
args.bedtools_bin,
)
if __name__ == "__main__":
main()
