本文以下面这个fastq文件为例,将fastq文件转成fasta文件
@seq1 AAAAAAAAAAAAAAAAAGGCGGGGCGCGGCTTTTTTTT + ?????????????????????????????????? @seq2 GGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCC + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
思路
- 解析
fastq
文件,提取对应的序列名和序列本身,不熟悉fastq文件的同学可以复习这里Fasta和Fastq文件 - 根据
fasta
的格式,将序列名和序列按照fasta
文件的格式进行写入
def parse_fastq(fastq_path): """ 解析fastq文件 fastq文件以四行为一个单位 1. @seqname 2. seq 3. + 4. %%%???))) 质量值 """ seq_name = None seq = None with open(fastq_path, 'r') as f: for line in f: if line.startswith('@'): """ 序列名所在的行以@开头 """ seq_name = line[1:].strip('\n') if seq: yield seq_name, seq seq = None if line.startswith('+'): pass if 'A' in line or 'G' in line or 'C' in line: # 这里通过判断某个碱基是否在line里面 # 来确定这一行到底是seq行,还是质量值的行 seq = line.strip('\n') yield seq_name, seq
上面解析文件的方法是通过来判断每一行中的一些特征来确定每一行的内容,但是不够灵活,下面给出一种较灵活的写法。
def parse_fastq(fastq_path): """ 解析fastq文件 fastq文件以四行为一个单位 1. @seqname 2. seq 3. + 4. %%%???))) 质量值 """ seq_name = None seq = None count = 0 with open(fastq_path, 'r') as f: for line in f: # 每次都加1, # 后面通过将count重置为1来进行4行为一个单元的划分 count += 1 if line.startswith('@'): """ 序列名所在的行以@开头 """ seq_name = line[1:].strip('\n') if seq: yield seq_name, seq seq = None count = 1 # 遇到@开头的行,将count重置为1 if count == 2: seq = line.strip('\n') yield seq_name, seq
上面这种解析方法通过利用fastq文件四行为一个单位这个信息,按照特定的规律来进行解析,大大增加了程序的灵活性。现在我们已经成功解析出来了fastq文件转成fasta文件所必需的信息,下面就可以将信息以fasta格式进行写入。
def _convert(fastq_path): result = [] for fastq in parse_fastq(fastq_path): seq_name = fastq[0] seq = fastq[1] result.append('>' + seq_name + '\n' + seq) # fasta文件的格式 return '\n'.join(result) def fastq_to_fasta(fastq_path,fasta_path): with open(fasta_path,'w') as f: f.write(str(_convert(fastq_path)))
至此,我们已经完成了fastq文件向fastq文件进行转化的任务。你可以通过调用fastq_to_fasta
来进行转换任务。