[Python知识库]python读取bam文件

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 qualities

??? read.query_sequence # read sequence bases

????read.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文件必须有索引