1、读取fasta文件:
(1)方法1:Bio库
from Bio import SeqIO# 读取包含单个序列 Fasta 格式文件
fa_seq = SeqIO.read("res/sequence1.fasta", "fasta")seq = str(fa_seq.seq)# 一个多序列文件中的所有序列
seqs = [fa.seq for fa in SeqIO.parse("res/multi.fasta", "fasta")]
(2)方法2:pysam库
fa = pysam.FastaFile(‘filename’)# 打印所有的references
fa.references# 打印所有的references的长度
fa.lengths# 获取指定references的[a, b)的序列ACTGN
seq = fa.fetch(reference=’chr21’, start=a, end=b) # 范围是前闭后开# 提取chr21的整条序列
seq = fa.fetch(reference=’chr21’)
2、将序列转化为字符串:str(fa.seq)
# 单个核苷酸计数
print("G Counts: ", fa.count("G"))# 获取反向序列
print("reverse: ", fa[::-1])# 获取反向互补序列
print("Reverse complement: ", fa.complement())
3、pysam读取bam文件时,报错:
File "pysam/libcalignmentfile.pyx", line 742, in pysam.libcalignmentfile.AlignmentFile.__cinit__
File "pysam/libcalignmentfile.pyx", line 952, in pysam.libcalignmentfile.AlignmentFile._open
File "pysam/libchtslib.pyx", line 365, in pysam.libchtslib.HTSFile.check_truncation
OSError: no BGZF EOF marker; file may be truncated
可以这样读取文件,即加入ignore_truncation=True:
bam = pysam.AlignmentFile(‘**.bam’, "rb", ignore_truncation=True)
4、CRAM格式的文件:具有高压缩的特性,并且比BAM有更高的压缩率,以后可能多数会使用CRAM格式的文件
5、比对数据(BAM/CRAM/SAM),变异数据:(VCF/BCF)
6、获得每一条read的情况
for read in bam:read.reference_name # 比对参考序列染色体名称;read.pos # read比对的位置read.mapq # read的比对质量值read.query_qualities # read sequence base qualitiesread.query_sequence # read sequence basesread.reference_length # 在reference genome上read比对的长度
7、读取cram文件
cf = pysam.AlignmentFile(‘*.cram’, ‘rc’)
8、读取sam文件
samfile = pysam.AlignmentFile(‘**.sam’, ‘r’)
9、获取bam文件中某一区域
for r in pysam.AlignmentFile(‘*.bam’, ‘rb’).fetch(‘chr21’, 300, 310):pass # 这样做的前提是:*.bam文件必须有索引