I am having problem in my code I am trying to read a fasta file i.e. "chr1.fa", then I have a mutation file which looks like this
chr1 822979 822980 CLL6.08_1_snv 88.2 +
chr1 1052781 1052782 CLL6.08_2_snv 388.9 +
chr1 1216196 1216197 CLL6.08_3_snv 625 +
chr1 5053847 5053848 CLL6.08_4_snv 722.2 +
chr1 5735093 5735094 CLL6.08_5_snv 138.9 +
this is a tab delimited file with chr1 as first column and + as last. I want to insert a N in the chr1.fa file the using the location from the second column.My code looks like this
#!/usr/bin/python
# Filename: mutation.py
import sys , os
import numpy as np
import re
#declaring the variables
lst = ''
chr_name = ''
first_cord = ''
second_cord = ''
lstFirstCord = []
sequence = ''
human_genome = ''
seqList = ''
# Method to read the Genome file (file contains data for only one chromosome)
def ReadgenomeCharacter():
header = ' '
seq = ' '
try:
human_genome = raw_input("Enter UCSC fasta file of human genome:")
human_seq = open(human_genome, 'rw+')
line = human_seq.readline()
except:
print 'File cannot be opened, wrong format you forgot something:', human_genome
exit()
while line:
line = line.rstrip('\n')
if '>' in line:
header = line
else:
seq = seq + line
line = human_seq.readline()
print header
print "Length of the chromosome is:",len(seq)
print "No. of N in the chromosome are:", seq.count("N")
return seq
#Method to replace the characters in sequence string
def ReplaceSequence():
seqList = list(sequence)
for index, item in enumerate(lstFirstCord):
if seqList[index] != "N":
seqList[index] = "N"
newSequence = ''.join(seqList)
return newSequence
#Method to write to genome file
def WriteToGenomeFile(newSequence):
try:
with open("chr1.fa", 'rw+') as f:
old_human_seq = f.read()
f.seek(0)
f.write(newSequence)
print "Data modified in the genome file"
print "Length of the new chromosome is:",len(newSequence)
print "No. of N in the new chromosome are:", newSequence.count("N")
except:
print 'File cannot be opened, wrong format you forgot something:', human_genome
exit()
sequence = ReadgenomeCharacter()
print "Here is my mutaiton file data"
data = np.genfromtxt("CLL608.txt",delimiter ="\t", dtype=None,skip_header=0) #Reading the mutation file CLL608.txt
#Storing the mutation file data in a dictionary
subarray = ({'Chr1' : data[data[:,0] == 'chr1'],'Chr2': data[data[:,0] == 'chr2'],'Chr3': data[data[:,0] == 'chr3'],
'Chr4': data[data[:,0] == 'chr4'], 'Chr5': data[data[:,0] == 'chr5'],'Chr6': data[data[:,0] == 'chr6'],
'Chr7': data[data[:,0] == 'chr7'], 'Chr8': data[data[:,0] == 'chr8'],'Chr9': data[data[:,0] == 'chr9'],
'Chr10': data[data[:,0] == 'chr10'] , 'Chr11': data[data[:,0] == 'chr11'],'Chr12': data[data[:,0] == 'chr12'],
'Chr13': data[data[:,0] == 'chr13'], 'Chr14': data[data[:,0] == 'chr14'],'Chr15': data[data[:,0] == 'chr15'],
'Chr16': data[data[:,0] == 'chr16'],'Chr17': data[data[:,0] == 'chr17'],'Chr18': data[data[:,0] == 'chr18'],
'Chr19': data[data[:,0] == 'chr19'], 'Chr20': data[data[:,0] == 'chr20'],'Chr21': data[data[:,0] == 'chr21'],
'Chr22': data[data[:,0] == 'chr22'], 'ChrX': data[data[:,0] == 'chrX']})
#For each element in the dictionary, fetch the first cord and pass this value to the method to replace the character on first chord with N in the genome file
for the_key, the_value in subarray.iteritems():
cnt = len(the_value)
for lst in the_value:
chr_name = lst[0]
first_cord = int(lst[1])
second_cord = int(lst[2])
lstFirstCord.append(first_cord)
#Call the method to replace the sequence
newSeq = ReplaceSequence()
print "length :", len(newSeq)
#Call the method to write new data to genome file
WriteToGenomeFile(newSeq)
`
I am getting output like this
Enter UCSC fasta file of human genome:chr1.fa
chr1
Length of the chromosome is: 249250622
No. of N in the chromosome are: 23970000
Here is my mutaiton file data
length : 249250622
File cannot be opened, wrong format you forgot something:
we can download the chr1.fa by typing the following command directly
rsync -avzP
rsync://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr1.fa.gz .
Somehow I am not able to insert the N in the sequence and also not able to wirte the new sequence. I will be glad if any valuable suggestion for improving the code :)
To make your life a bit easier, you may want to consider using Biopython to read in your fasta and do your converting.
Here's some documentation to get you started http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc16
Here is some starter code.
from Bio import SeqIO
handle = open("example.fasta", "rU")
output_handle = open("output.fasta", "w")
for record in SeqIO.parse(handle, "fasta"):
print record.seq
handle.close()
output_handle.close()