参考地址：https://github.com/lmdu/pyfastx
# fasta文件
# Fasta读取fasta文件，build_index设置是否构建索引，默认为true;full_name设置是否使用完整的标题。
for name, seq in pyfastx.Fasta('/data/test.fa.gz', build_index=False, full_name=True):
print(name)
print(seq)
# 也可以像这样从FASTA对象迭代序列对象
fa = pyfastx.Fasta('data/test.fa.gz')
for seq in fa:
print(seq.name)
print(seq.seq)
print(seq.description)
# 序列数
print(len(fa))
# 总序列长度
print(fa.size)
# 获取DNA序列的GC含量
print(fa.gc_content)
# 获取DNA序列的GC偏态
print(fa.gc_skew)
# 获取fasta文件中核苷酸的组成
print(fa.composition)
# fasta文件类型（DNA、RNA、蛋白质）
print(fa.type)
# 检查fasta文件是gzip压缩的
print(fa.is_gzip)
# 得到最长和最短的序列,和对应的id
print(fa.longest)
print(fa.longest.name)
print(len(fa.longest))
print(fa.shortest)
print(fa.shortest.name)
print(len(fa.shortest))
# 计算N50、L50
print(fa.nl(50))
# 计算N90、L90
print(fa.nl(90))
# 计算N75、L75
print(fa.nl(75))
# 获取序列均值和中位数长度
print(fa.mean)
print(fa.median)
# 获取长度>=200指定长度的序列的计数
print(fa.count(200))
# 得到子序列
# 获取起始和结束位置的子序列
interval = (1, 10)
print(fa.fetch('JZ822577.1', interval))
# 获取起始和结束位置列表的子序列
intervals = [(1, 10), (50, 60)]
print(fa.fetch('JZ822577.1', intervals))
# 获取反向链的子序列
print(fa.fetch('JZ822577.1', (1, 10), strand='-'))
# 获取给定子序列的右侧序列，返回包含左侧和右侧序列的元组
# 获取JZ822577:50-60子序列长度为20的侧翼序列
print(fa.flank('JZ822577.1', 50, 60, 20))
# 获得单个碱基或SNP位置100的侧翼序列
print(fa.flank('JZ822577.1', 100, 100, 20))
# 通过buffer cache获取右翼序列(使用了缓存（buffer）来提高处理速度)
print(fa.flank('JZ822577.1', 70, 90, flank_length=20, use_cache=True))
# 有时你的fasta会有一个很长的标题，其中包含多个标识符和描述，例如:>JZ822577.1 contig1 cDNA library of flower petals in tree peony by suppression subtractive hybridization Paeonia suffruticosa cDNA, mRNA sequence
# 在这种情况下，可以使用“JZ822577.1”或“contig1”作为标识符。通过指定key函数，可以选择一个作为标识符。
# 默认使用JZ822577.1作为标识符
# 指定key_func选择contig1作为标识符
# 如果索引文件已经存在，则应该删除以前的索引文件，然后使用key_func创建新的索引文件
fa = pyfastx.Fasta('/data/test.fa.gz', key_func=lambda x: x.split()[1])
# 得到序列和序列信息
s1 = fa['JZ822577.1']
print(s1)
s2 = fa[-1]
print(s2)
print(s2.id)
print(s2.name)
print(s2.description)
print(s2.seq)
print(s2.raw)
print(len(s2))
print(s2.gc_content)
print(s2.composition)
# 切片、互补、反转
print(fa[0][10:20].seq)
print(fa[0][10:20].reverse)
print(fa[0][10:20].complement)
print(fa[0][10:20].antisense)
# 逐行打印
for line in fa[0]:
print(line)
# 搜索子序列
# 从给定序列中搜索子序列，并获得第一次出现的起始位置
print(fa[0].search('GCTTCAATACA'))
# 是否存在
print('GCTTCAATACA' in fa[0])
# 搜索反义链子序列
print(fa[0].search('CCTCAAGT', '-'))
# FastaKeys
# 获取序列的所有名称作为类列表对象。
ids = fa.keys()
print(ids[0])
# 按长度降序排序键
for name in ids.sort(by='length', reverse=True):
print(name)
# 按名称升序排序键
for name in ids.sort(by='name'):
print(name)
# 按id降序排序键
for name in ids.sort(by='id', reverse=True):
print(name)
# 按序列长度和名称过滤键
# 获取长度> 600的键
ids.filter(ids > 600)
list(ids)
ids.filter(ids >= 500, ids <= 700)
ids.filter(500 < ids < 600)
ids.filter(ids % 'JZ8226')
ids.filter(ids % 'JZ8226', ids > 550)
ids.reset()
ids.filter(ids % 'JZ8226', ids > 730)
ids.filter(ids % 'JZ8226', ids > 730).sort('length', reverse=True)
ids.filter(ids % 'JZ8226', ids > 730).sort('name', reverse=True)
# fastq文件
fq = pyfastx.Fastq('/data/test.fq.gz', build_index=False)
for name, seq, qual in fq:
print(name)
print(seq)
print(qual)
for read in pyfastx.Fastq('/data/test.fq.gz'):
print(read.name)
print(read.seq)
print(read.qual)
print(read.quali)
# 序列数
print(len(fq))
# 总碱基
print(fq.size)
# gc含量
print(fq.gc_content)
# 碱基组成
print(fq.composition)
# 读取的平均长度
print(fq.avglen)
# 最大长度
print(fq.maxlen)
# 最小长度
print(fq.minlen)
# 最高质量分数
print(fq.maxqual)
# 最低质量分数
print(fq.minqual)
# phred值会影响质量分数的转换
print(fq.phred)
r1 = fq['A00129:183:H77K2DMXX:1:1101:4752:1047']
print(r1)
r2 = fq[-1]
print(r2)
print(r2.id)
print(r2.name)
print(r2.description)
print(len(r2))
print(r2.seq)
print(r2.raw)
print(r2.qual)
print(r2.quali)
ids = fq.keys()
print(len(ids))
print(ids[0])
# 命令行界面
# pyfastx -h
# 命令:
# index：为fasta/q文件建立索引
# Stat：显示fasta/q文件的详细统计信息
# Split：将fasta/q文件分割成多个文件
# fq2fa：转换fastq文件到fasta文件
# subseq：按区域从fasta文件中获取子序列
# sample：从fasta或fastq文件中随机取样
# extract：从fasta/q文件中提取完整序列或读取
