2021/03/22
这是一个小任务
读取某个目录下的文件
import os
aim_dir = '/home/user/fastq_dir/'
fastq_files = os.listdir(aim_dir)
print(fastq_files)
结果是一个列表，包括此目录下所有的文件的全称 及 文件扩展名（.txt,.fasta,.gz等）
for fastq_file in fastq_files:
print(fastq_file)
循环输出每一个文件名
若文件名格式一致，如为 sample_name.fastq.gz ,可使用 正则表达式 获取想要的信息
若目标目录 下不是文件而是文件夹，可循环上一步获取所需要文件夹内的内容
# 例如 文件范例为 /home/user/fastq_dir/sample_name/sample_name.fastq.gz
samples = os.listdir(aim_dir)
for sample in samples:
fastq_path = f'aim_dir{sample}/{sample}.fastq.gz'
随后可以处理 fastq.gz 文件
f'{sample}'可以填入变量，是形成一段固定格式的字符串，与 '{0}'.format(sample) 相同，但更简略一点
在Linux中：查看 .gz 压缩文件
查看 .gz 压缩文件
在Linux中：
less -S *.gz | grep 'keyword' | less -S # 查看
less -S *.gz | grep 'keyword' > *.txt # 写入到txt文件中
gunzip -c *.gz | grep 'keyword' # 用gunzip不解压的查看
在python中处理.gz文件
如 fastq.gz 文件
from Bio import SeqIO # 导入biopython
from Bio.SeqIO import parse
import gzip
根据上面得到的 .fastq.gz 的绝对路径，读成字典
fastq_dt = SeqIO.to_dict(parse(gzip.open(fastq, 'rt'), 'fastq')) # 解压构建字典
seq = fastq_dt[id].seq # Seq(' ') 类型序列 无法写入xlsx
seq = str(fastq_dt[id].seq) # 普通字符串类型 可以写入xlsx
获取注释
print(fastq_dt[id].letter_annotations)
{'phred_quality': [37, 37, 37, 37, 37, 37, 37, 37]} 这是字典类型的质量值
print(fastq_dt[id].letter_annotations['phred_quality'])
[37, 37, 37, 37, 37, 37, 37, 37] 这是列表类型的质量值，方便之后过滤质量值不好的reads
可以根据关键的 id 获取序列
或者 构建一个字典 统计一下 fastq.gz 文件中 序列的 数量
from collections import defaultdict
from Bio import SeqIO
from Bio.SeqIO import parse
import gzip
seq_dict = defaultdict(int)
for record in SeqIO.parse(gzip.open(fastq_path, 'rt'), 'fastq'):
seq = str(record.seq)
seq_dict[seq] += 1
若是多个样本的 fastq.gz 文件
seq_dict = defaultdict(lambda: defauludict(int))
for sample in samples:
for record in SeqIO.parse(gzip.open(fastq_path, 'rt'), 'fastq'):
seq = str(record.seq)
seq_dict[sample][seq] += 1
这样就可以多个样本一起统计啦
可是，这样构建的字典里面是无序的，怎样把序列按 数量 排个序呢？？
对于单个样本：
seq_dict = {'aaaaaa':5, 'bbbbbb': 3, 'cccc': 6} # 这种格式的
print(seq_dict.items())
# dict_items([('aaaaaa', 5), ('bbbbbb', 3), ('cccc', 6)])
实际上变成了 由 一对键值对 组成的元组 组成的 列表
可以根据 每个 元组的 第二位，也就是 数量进行 排序
sort_seq = sorted(seq_dict.items(), key=lambda x: x[1], reverse=True)
print(sort_seq)
# [('cccc', 6), ('aaaaaa', 5), ('bbbbbb', 3)]
排列完还是个列表
可以输出你想要的前多少个
例如：前500个
for seq in sort_seq[:500]:
print(seq[0], seq[1])
对于多个样本，也差不多
seq_dict = {
'sample1': {'aaaaaa':5, 'bbbbbb': 3, 'cccc': 6}
'sample2': {'aaaaaa':5, 'bbbbbb': 3, 'cccc': 6}
} # 这种格式的
for sample in seq_dict:
sam_seq_dict = seq_dict[sample]
sort_sam_seq = sorted(sam_seq_dict.items(), key=lambda x: x[1], reverse=True)
for seq in sort_seq[:500]:
print(seq[0], seq[1])
简直一毛一样有没有
