SoloTE/SoloTE_RepeatMasker_to_BED.py at main · bvaldebenitom/SoloTE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

import argparse
import pandas
import requests
import os
from colorama import Fore, init
from tqdm import tqdm

init(autoreset=True)

argparse_object = argparse.ArgumentParser(description="SoloTE RepeatMasker to BED")

argparse_object.add_argument("-g","--genome",help="Genome assembly identifier to obtain RepeatMasker annotation (available genomes can be checked with option -l).")
argparse_object.add_argument("-l","--list",help="List available genomes at the UCSC database.",action='store_true',required=False)
argparse_object.add_argument("-i","--inputfile",help="Process your own RepeatMasker file.")

commandargs = argparse_object.parse_args()
list_genomes = commandargs.list
genome_assembly = commandargs.genome
input_file = commandargs.inputfile

mode = "RepeatMasker"

def download(url: str, fname: str, chunk_size=1024):
	resp = requests.get(url, stream=True)
	total = int(resp.headers.get('content-length', 0))
	with open(fname, 'wb') as file, tqdm(
		desc=fname,
		total=total,
		unit='iB',
		unit_scale=True,
		unit_divisor=1024,
	) as bar:
		for data in resp.iter_content(chunk_size=chunk_size):
			size = file.write(data)
			bar.update(size)


if list_genomes == True:
	api_url = "https://api.genome.ucsc.edu/list/ucscGenomes"
	response = requests.get(api_url)
	response = response.json()
	ucscGenomes = response["ucscGenomes"]

	for genome in ucscGenomes:
		current = ucscGenomes[genome]
		genome_description = current['description']
		scientificname = current['scientificName']
		organism = current['organism']
		print(genome+"\t| "+organism+" ["+scientificname+", "+genome_description+"]")
	exit()


if input_file:
	rmsk_filename = input_file
	genome_assembly = rmsk_filename.replace(".gz","")
	genome_assembly = genome_assembly.replace(".fa.out","")
	rmsk_bed_filename = genome_assembly+'_rmsk.bed'
	print(f"[LOG] Beginning conversion of RepeatMasker file {rmsk_filename} to BED format")
	rmsk_out = pandas.read_csv(rmsk_filename,compression="gzip",skiprows=3,header=None,sep=" ",skipinitialspace=True)
	rmsk_out.columns = ['SW_score','percDiv','percDel','percIns','querySeq','queryStart','queryEnd','queryLeft','strand','matchingRepeat','repeatClass_Family','repeatBegin','repeatStart','repeatEnd','ID']
	repeatinfo = rmsk_out['repeatClass_Family'].str.split("/",expand=True)
	repeatinfo = repeatinfo.fillna(value="-")
	rmsk_out['repeatID'] = rmsk_out['matchingRepeat']+":"+repeatinfo[1]+":"+repeatinfo[0]
	rmsk_out['percDiv'] = rmsk_out['percDiv'].astype(str)
	rmsk_out['queryStart'] = rmsk_out['queryStart'].astype(str)
	rmsk_out['queryEnd'] = rmsk_out['queryEnd'].astype(str)
	rmsk_out['fixedStrand'] = rmsk_out['strand'].str.replace("C","-")
	rmsk_out['te_name'] = rmsk_out['querySeq']+"|"+rmsk_out['queryStart']+"|"+rmsk_out['queryEnd']+"|"+rmsk_out['repeatID']+"|"+rmsk_out['percDiv']+"|"+rmsk_out['fixedStrand']
	rmsk_out_bed = rmsk_out[['querySeq','queryStart','queryEnd','te_name','percDiv','fixedStrand']]
	rmsk_out_bed = rmsk_out_bed[rmsk_out['repeatClass_Family'].str.contains("LINE|SINE|LTR|DNA|RC")]
	rmsk_out_bed = rmsk_out_bed[~rmsk_out_bed['querySeq'].str.contains("chrna|_fix|_random|_alt|chrUn")]
	rmsk_out_bed.to_csv(rmsk_bed_filename,sep="\t",header=None,index=False)
	print(Fore.GREEN+"[OK] ",end='')
	print("Finished generating "+rmsk_bed_filename)
	exit()


api_url = "https://api.genome.ucsc.edu/list/ucscGenomes"
response = requests.get(api_url)
response = response.json()
ucscGenomes = response["ucscGenomes"]

if ucscGenomes.get(genome_assembly) is None:
	print(Fore.RED+"[ERROR] ",end='')
	print("Genome assembly "+genome_assembly+" not available. Check available genomes with option -l.")
	exit()
else:
	print(Fore.GREEN+"[OK] ",end='')
	print("Genome assembly "+genome_assembly+" found.")


if mode == 'UCSC':
	api_url = "http://api.genome.ucsc.edu/getData/track?genome="+genome_assembly+";track=rmsk;maxItemsOutput=1"
	response = requests.get(api_url)
	response = response.json()
	rmsk_url = response["dataDownloadUrl"]
	print("[LOG] URL to fetch RepeatMasker file: "+rmsk_url)

	rmsk_filename = genome_assembly+'_ucsc_rmsk.txt.gz'
	#r = requests.get(rmsk_url, allow_redirects=True)
	#open(rmsk_filename,'wb').write(r.content)
	print("[LOG] Downloading RepeatMasker file to "+rmsk_filename)
	download(url=rmsk_url,fname=rmsk_filename)


	rmsk_bed_filename = genome_assembly+'_ucsc_rmsk.bed'
	print("[LOG] Beginning conversion of RepeatMasker file to BED format")
	rmsk_table = pandas.read_csv(rmsk_filename,compression="gzip",header=None,sep="\t")
	rmsk_table.columns = ["bin","swScore","milliDiv","milliDel","milliIns","genoName","genoStart","genoEnd","genoLeft","strand","repName","repClass","repFamily","repStart","repEnd","repLeft","ID"]
	#print(rmsk_table)
	rmsk_table['fixed_milliDiv'] = rmsk_table['milliDiv'].astype(str).str.replace("([0-9]$)",".\\1",regex=True)
	rmsk_table['genoStart'] = rmsk_table['genoStart'].astype(str)
	rmsk_table['genoEnd'] = rmsk_table['genoEnd'].astype(str)
	rmsk_table['genoName']+"|"+str(rmsk_table['genoStart'])+"|"+str(rmsk_table['genoEnd'])
	#rmsk_table['te_name'] = rmsk_table['genoName']+"|"+rmsk_table['genoStart']+"|"+rmsk_table['genoEnd']+"|"+rmsk_table['repClass']+":"+rmsk_table['repFamily']+":"+rmsk_table['repName']+"|"+rmsk_table['fixed_milliDiv']+"|"+rmsk_table['strand']
	rmsk_table['te_name'] = rmsk_table['genoName']+"|"+rmsk_table['genoStart']+"|"+rmsk_table['genoEnd']+"|"+rmsk_table['repName']+":"+rmsk_table['repFamily']+":"+rmsk_table['repClass']+"|"+rmsk_table['fixed_milliDiv']+"|"+rmsk_table['strand']
	rmsk_table_bed = rmsk_table[['genoName','genoStart','genoEnd','te_name','fixed_milliDiv','strand']]
	rmsk_table_bed = rmsk_table_bed[rmsk_table_bed['te_name'].str.contains("LINE|SINE|LTR|DNA")]
	rmsk_table_bed.to_csv(rmsk_bed_filename,sep="\t",header=None,index=False)
	print(Fore.GREEN+"[OK] ",end='')
	print("Finished generating "+rmsk_bed_filename)


delete_rmsk = False
if mode == "RepeatMasker":
	rmsk_url = "https://hgdownload.soe.ucsc.edu/goldenPath/"+genome_assembly+"/bigZips/"+genome_assembly+".fa.out.gz"

	if genome_assembly == "hs1":
		rmsk_url = "https://hgdownload.soe.ucsc.edu/goldenPath/"+genome_assembly+"/bigZips/"+genome_assembly+".repeatMasker.out.gz"


	print("[LOG] URL to fetch RepeatMasker file: "+rmsk_url)
	rmsk_filename = genome_assembly+'.fa.out.gz'
	print("[LOG] Downloading RepeatMasker file to "+rmsk_filename)
	download(url=rmsk_url,fname=rmsk_filename)

	if input_file:
                rmsk_filename = input_file

	rmsk_bed_filename = genome_assembly+'_rmsk.bed'
	print("[LOG] Beginning conversion of RepeatMasker file to BED format")
	rmsk_out = pandas.read_csv(rmsk_filename,compression="gzip",skiprows=3,header=None,sep=" ",skipinitialspace=True)
	rmsk_out.columns = ['SW_score','percDiv','percDel','percIns','querySeq','queryStart','queryEnd','queryLeft','strand','matchingRepeat','repeatClass_Family','repeatBegin','repeatStart','repeatEnd','ID']
	repeatinfo = rmsk_out['repeatClass_Family'].str.split("/",expand=True)
	repeatinfo = repeatinfo.fillna(value="-")
	rmsk_out['repeatID'] = rmsk_out['matchingRepeat']+":"+repeatinfo[1]+":"+repeatinfo[0]
	rmsk_out['percDiv'] = rmsk_out['percDiv'].astype(str)
	rmsk_out['queryStart'] = rmsk_out['queryStart'].astype(str)
	rmsk_out['queryEnd'] = rmsk_out['queryEnd'].astype(str)
	rmsk_out['fixedStrand'] = rmsk_out['strand'].str.replace("C","-")
	rmsk_out['te_name'] = rmsk_out['querySeq']+"|"+rmsk_out['queryStart']+"|"+rmsk_out['queryEnd']+"|"+rmsk_out['repeatID']+"|"+rmsk_out['percDiv']+"|"+rmsk_out['fixedStrand']
	rmsk_out_bed = rmsk_out[['querySeq','queryStart','queryEnd','te_name','percDiv','fixedStrand']]
	rmsk_out_bed = rmsk_out_bed[rmsk_out['repeatClass_Family'].str.contains("LINE|SINE|LTR|DNA|RC")]
	rmsk_out_bed = rmsk_out_bed[~rmsk_out_bed['querySeq'].str.contains("chrna|_fix|_random|_alt|chrUn")]
	rmsk_out_bed.to_csv(rmsk_bed_filename,sep="\t",header=None,index=False)
	print(Fore.GREEN+"[OK] ",end='')
	print("Finished generating "+rmsk_bed_filename)


if delete_rmsk:
	os.remove(rmsk_filename)