參考地址：https://github.com/lmdu/pyfastx
# fasta文件
# Fasta讀取fasta文件，build_index設(shè)置是否構(gòu)建索引司志，默認為true;full_name設(shè)置是否使用完整的標(biāo)題震叙。
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)
# 序列數(shù)
print(len(fa))
# 總序列長度
print(fa.size)
# 獲取DNA序列的GC含量
print(fa.gc_content)
# 獲取DNA序列的GC偏態(tài)
print(fa.gc_skew)
# 獲取fasta文件中核苷酸的組成
print(fa.composition)
# fasta文件類型（DNA栈雳、RNA怯邪、蛋白質(zhì)）
print(fa.type)
# 檢查fasta文件是gzip壓縮的
print(fa.is_gzip)
# 得到最長和最短的序列,和對應(yīng)的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))
# 獲取序列均值和中位數(shù)長度
print(fa.mean)
print(fa.median)
# 獲取長度>=200指定長度的序列的計數(shù)
print(fa.count(200))
# 得到子序列
# 獲取起始和結(jié)束位置的子序列
interval = (1, 10)
print(fa.fetch('JZ822577.1', interval))
# 獲取起始和結(jié)束位置列表的子序列
intervals = [(1, 10), (50, 60)]
print(fa.fetch('JZ822577.1', intervals))
# 獲取反向鏈的子序列
print(fa.fetch('JZ822577.1', (1, 10), strand='-'))
# 獲取給定子序列的右側(cè)序列营密，返回包含左側(cè)和右側(cè)序列的元組
# 獲取JZ822577:50-60子序列長度為20的側(cè)翼序列
print(fa.flank('JZ822577.1', 50, 60, 20))
# 獲得單個堿基或SNP位置100的側(cè)翼序列
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會有一個很長的標(biāo)題械媒，其中包含多個標(biāo)識符和描述，例如:>JZ822577.1 contig1 cDNA library of flower petals in tree peony by suppression subtractive hybridization Paeonia suffruticosa cDNA, mRNA sequence
# 在這種情況下评汰，可以使用“JZ822577.1”或“contig1”作為標(biāo)識符纷捞。通過指定key函數(shù)，可以選擇一個作為標(biāo)識符被去。
# 默認使用JZ822577.1作為標(biāo)識符
# 指定key_func選擇contig1作為標(biāo)識符
# 如果索引文件已經(jīng)存在主儡，則應(yīng)該刪除以前的索引文件，然后使用key_func創(chuàng)建新的索引文件
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)
# 切片惨缆、互補糜值、反轉(zhuǎn)
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)
# 搜索子序列
# 從給定序列中搜索子序列，并獲得第一次出現(xiàn)的起始位置
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)
# 序列數(shù)
print(len(fq))
# 總堿基
print(fq.size)
# gc含量
print(fq.gc_content)
# 堿基組成
print(fq.composition)
# 讀取的平均長度
print(fq.avglen)
# 最大長度
print(fq.maxlen)
# 最小長度
print(fq.minlen)
# 最高質(zhì)量分數(shù)
print(fq.maxqual)
# 最低質(zhì)量分數(shù)
print(fq.minqual)
# phred值會影響質(zhì)量分數(shù)的轉(zhuǎn)換
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文件的詳細統(tǒng)計信息
# Split：將fasta/q文件分割成多個文件
# fq2fa：轉(zhuǎn)換fastq文件到fasta文件
# subseq：按區(qū)域從fasta文件中獲取子序列
# sample：從fasta或fastq文件中隨機取樣
# extract：從fasta/q文件中提取完整序列或讀取
