Skip to content
Snippets Groups Projects

adding cDNA generation script

Closed adding cDNA generation script
suvi into main
Closed Suvarnan Selliah requested to merge suvi into main
Compare
and
1 file
+ 239
− 0
Compare changes
  • Side-by-side
  • Inline
src/generate_cdna.py 0 → 100644
+ 239
− 0
#!/usr/bin/env python3
#Author: Suvarnan Selliah
class GeneratecDNA:
def generatecDNA(fasta, gtf, cp_nr, my_output_fasta="cDNA.fasta", my_output_csv="cDNA.csv", placeholder = "ph-file.csv"):
cDNA_transcript_id = 1
#READING INPUT FILES / PART I
#open files
with open(fasta, 'r') as fa, open(gtf, 'r') as gt, open(cp_nr, 'r') as cp:
#read fasta-file
for myFastaline in fa:
#search for transcript id and transcript sequence in fasta-file
fasta_id = ""
fasta_seq = ""
fasta_id_found = False
fasta_seq_found = False
currentFastaString = myFastaline
if fasta_id_found == False:
position_of_start = currentFastaString.find('>')
if position_of_start != 0:
continue
elif position_of_start == 0:
fasta_id = myFastaline
fasta_id = fasta_id.replace(">", "")
fasta_id_found = True
continue
else:
print("FASTA: Start position in fasta file not found")
break
if fasta_id_found == True and fasta_seq_found == False:
while fasta_seq_found == False:
currentFastaString = fa.readline()
zero_position = currentFastaString[0]
if zero_position == ";":
continue
elif zero_position == ">":
print("FASTA: No Sequence after headline")
break
else:
fasta_seq = currentFastaString
fasta_seq_found = True
#starting to work with gtf-file
#defining variables for gtf-file
gtf_seqname = ''
gtf_start = 0
gtf_end = 0
gtf_score = 0.0
gtf_info_found = False
gtf_entries = 0
gtf_list_of_lines = []
gtf_prob_list = []
#defining variables for csv-file
csv_trans_id = ""
csv_gene_id = ""
csv_count = 0
csv_info_found = False
if fasta_id_found == True and fasta_seq_found == True:
# search for transcript id from fasta-file in gtf-file
for myGTFline in gt:
currentGTFString = myGTFline
gtf_list = currentGTFString.split('\t')
gtf_seqname = gtf_list[0]
if gtf_seqname == fasta_id:
gtf_entries += 1
gtf_start = gtf_list[3]
gtf_end = gtf_list[4]
gtf_score = gtf_list[5]
gtf_temp_list = []
gtf_temp_list.append(gtf_start)
gtf_temp_list.append(gtf_end)
gtf_temp_list.append(gtf_score)
gtf_list_of_lines.append(gtf_temp_list)
gtf_prob_list.append(gtf_score)
else:
continue
if gtf_entries != 0:
gtf_info_found = True
fasta_id_found = False
fasta_seq_found = False
assert gtf_info_found, "Sequence ID from fasta-file not found in gtf-file"
#search copy number of transcript in 3. file/csv-file
for myCSVline in cp:
currentCSVString = myCSVline
csv_list = currentCSVString.split(',')
csv_trans_id = csv_list[0]
if csv_trans_id == fasta_id:
csv_gene_id = csv_list[1]
csv_count = csv_list[2]
csv_info_found = True
gtf_info_found = False
break
else:
continue
assert csv_info_found, "Data (TranscriptID,GeneID,Count) from csv-file/3.file not found"
#COMPUTATION & OUTPUT / PART II
#set score to 0 for primimg sites close to the end of sequence
csv_info_found = False
seq_len = len(fasta_seq)
seq_len -= 22
gtf_list_of_lines_len = len(gtf_list_of_lines)
for i in range(gtf_list_of_lines_len):
if gtf_list_of_lines[i][1] > seq_len:
gtf_list_of_lines[i][2] = 0.0
gtf_prob_list[i] = 0.0
#assign priming sites (according to score) to copy number
sum_of_score = 0.0
for i in gtf_prob_list:
sum_of_score += i
one_score = 100/sum_of_score
norm_list = []
gtf_prob_list_len = len(gtf_prob_list)
for i in range(gtf_prob_list_len):
norm_list.append((one_score * gtf_prob_list[i]))
one_norm = csv_count / 100
distr_RNA_to_prim_sites = []
norm_list_len = len(norm_list)
for i in range(norm_list_len):
distr_RNA_to_prim_sites.append((one_norm * norm_list[i]))
total_RNA_number = 0
distr_RNA_to_prim_sites_len = len(distr_RNA_to_prim_sites)
for i in range(distr_RNA_to_prim_sites_len):
total_RNA_number += distr_RNA_to_prim_sites[i]
new_distr_RNA_to_prim_sites = []
if total_RNA_number != csv_count:
for i in range(distr_RNA_to_prim_sites_len):
new_distr_RNA_to_prim_sites.append(int(distr_RNA_to_prim_sites[i]))
new_distr_RNA_to_prim_sites_len = len(new_distr_RNA_to_prim_sites)
counter = 0
while total_RNA_number != csv_count:
new_distr_RNA_to_prim_sites[counter] = round(distr_RNA_to_prim_sites[counter])
counter += 1
assert counter <= new_distr_RNA_to_prim_sites_len, "Calculated RNA transcripts (assigned to priming sites) are more than initial count"
#order the priming sites
prim_sites_ordered = []
for i in range(gtf_list_of_lines_len):
prim_sites_ordered.append(gtf_list_of_lines[i][0])
prim_sites_ordered.sort()
#searching for 2 priming sites
prim_sites_ordered_len = len(prim_sites_ordered)
ph_1 = 0
ph_2 = 0
for i in range(0, (prim_sites_ordered_len-1), 2):
if prim_sites_ordered[i] == 0.0:
continue
else:
search_for_1 = prim_sites_ordered[i]
search_for_2 = prim_sites_ordered[i+1]
for j in range(gtf_list_of_lines_len):
if gtf_list_of_lines[j][0] == search_for_1:
ph_1 = j
if gtf_list_of_lines[j][0] == search_for_2:
ph_2 = j
#making cDNA and comparing in library
start_1 = gtf_list_of_lines[ph_1][0]
start_2 = gtf_list_of_lines[ph_2][0]
start_1 -= 1
start_2 -= 1
trans_between_prim_sites = fasta_seq[start_1:start_2]
cDNA = ''
for element in range(0, len(trans_between_prim_sites)):
if trans_between_prim_sites[element] == 'A':
cDNA[element] = 'T'
elif trans_between_prim_sites[element] == 'U':
cDNA[element] = 'A'
elif trans_between_prim_sites[element] == 'G':
cDNA[element] = 'C'
elif trans_between_prim_sites[element] == 'C':
cDNA[element] = 'G'
else:
assert False, "cDNA synthesis failed, position is not A,U,G or C in transcript"
# open output files
if i == 0:
with open(my_output_fasta, 'a') as myfasta, open(my_output_csv, 'a') as mycsv:
myfasta.write(">" + string(cDNA_transcript_id))
myfasta.write("\n")
myfasta.write(cDNA)
mycsv.write(",".join([cDNA_transcript_id, csv_gene_id, new_distr_RNA_to_prim_sites[ph_1]]))
else:
found_cDNA_id = ''
found_cDNA_id_bool = False
with open(my_output_fasta, 'r') as myfasta, open(my_output_csv, 'r') as mycsv, open(placeholder, 'w') as phf:
for myline in myfasta:
pos = myline.find('>')
if pos != 0:
continue
if pos == 0:
fid = myline
fid = fid.replace(">", "")
fbool = False
while fbool == False:
myline = myfasta.readline()
zer = myline[0]
if zer == ";":
continue
elif zer == ">":
assert False, "Error in searching cDNA in output fasta-file"
else:
fseq = myline
fbool = True
if fseq ==cDNA:
found_cDNA_id = fid
found_cDNA_id_bool = True
break
if found_cDNA_id_bool == True:
for myline_csv_out in mycsv:
csvlist = myline_csv_out.split(',')
csvcDNAid = csvlist[0]
if csvcDNAid == found_cDNA_id:
csvgeneid = csvlist[1]
csvcDNAcount = csvlist[2]
csvcDNAcount += new_distr_RNA_to_prim_sites[ph_1]
phf.write(",".join([csvcDNAid, csvgeneid, csvcDNAcount]))
else:
phf.write(myline_csv_out)
if found_cDNA_id_bool == True:
with open(my_output_csv, 'w') as mycsv, open(placeholder, 'r') as phf:
for myplaceholder in phf:
mycsv.write(myplaceholder)
else:
cDNA_transcript_id += 1
with open(my_output_fasta, 'a') as myfasta, open(my_output_csv, 'a') as mycsv:
myfasta.write(">" + string(cDNA_transcript_id))
myfasta.write("\n")
myfasta.write(cDNA)
mycsv.write(",".join(
[cDNA_transcript_id, csv_gene_id, new_distr_RNA_to_prim_sites[ph_1]]))
return my_output_fasta, my_output_csv
\ No newline at end of file