当前位置: 首页 » Python(人生苦短) » py_module » pysam读写Bam/SAM文件

pysam读写Bam/SAM文件

[文章目录] x

pysam是python中用来读写SAM/BAM文件的一个模块。底层是用C写的

一、简介

pysam官网

Pysam is a wrapper of the htslib C-API and provides facilities to read and write SAM/BAM/VCF/BCF/BED/GFF/GTF/FASTA/FASTQ files as well as access to the command line functionality of the samtools and bcftools packages. The module supports compression and random access through indexing.

安装

二. 函数详解

1.读写bam

#然后可以通过fetch方法来调取比对到某个区间对应的reads信息

得到结果:

#取的位置是交集还是并集呢?
初步理解为跟这个区间有交集的序列都提出来

也可以写bam文件

注:

1.AlignmentFile(filepath_or_object, mode=None, template=None, reference_names=None, reference_lengths=None, text=NULL, header=None, add_sq_text=False, check_header=True, check_sq=True, reference_filename=None, filename=None, duplicate_filehandle=True)

2.AlignmentFile打开得bam文件需要是已经建立索引,有.bai文件存在的bam文件

3.template将通过AlignmentFile打开得bam文件的header拷贝过来

 

2.提取某个碱基

输出结果:

实现逻辑:
1.根据给得区间"chr1", 100, 120,一次,从100到120点所对应的所有reads,通过pileup提取的reads有这几个属性:'n','nsegments', 'pileups', 'pos', 'reference_id', 'reference_pos', 'tid'

a.起点,终点?这里应该是先找出跟区间至少有一个碱基交集的reads,然后这些reads的所覆盖的区间的起点为pileup的起点

b.如果加上参数truncate=True,那起点就成了100,119,注意这里起始的变化,这里仅仅是输出做了一个限制,但是提取的reads的方式还是没有变化,就是只要跟指定的区间有一个碱基交集的reads就给提出来了

c.而pos代表该位点的坐标,n代表它的coverage,pileups代表对应的reads

2.每个点对应的每条read可以通过pileups来提,而这些reads有这些属性’alignment', 'indel', 'is_del', 'is_head', 'is_refskip', 'is_tail', 'level', 'query_position
3.alignment的属性包含: aend, alen, aligned_pairs, bin, blocks, cigar, cigarstring, cigartuples, compare, flag, get_aligned_pairs, get_blocks, get_overlap, get_reference_positions, get_tag, get_tags, has_tag, infer_query_length, inferred_length, is_duplicate, is_paired, is_proper_pair, is_qcfail, is_read1, is_read2, is_reverse, is_secondary, is_supplementary, is_unmapped, isize, mapping_quality, mapq, mate_is_reverse, mate_is_unmapped, mpos, mrnm, next_reference_id, next_reference_start, opt, overlap, pnext, pos, positions, qend, qlen, qname, qqual, qstart, qual, query, query_alignment_end, query_alignment_length, query_alignment_qualities, query_alignment_sequence, query_alignment_start, query_length, query_name, query_qualities, query_sequence, reference_end, reference_id, reference_length, reference_start, rlen, rname, rnext, seq, setTag, set_tag, set_tags, tags, template_length, tid, tlen

3.排序

跟samtools的效果一样

4.打开其他文件(这个功能暂时没有用到)

5.其他函数与功能略

提取成对reads的start和end

参考资料:https://github.com/etal/cnvkit/issues/212

 

 

三、讨论:

  1. ValueError: quality and sequence mismatch: 148 != 151
    原因: pysam读取后的seq和qual长度必须一样啊!!

2.文件的写入

注:
先整理到这里,后面用的时候,再回来梳理。

参考资料:
http://pysam.readthedocs.org/en/latest/glossary.html#term-column
http://pysam.readthedocs.org/en/latest/api.html#pysam.AlignmentFile.pileup


下一篇 :

上一篇 :

暂无评论

发表评论

电子邮件地址不会被公开。 必填项已用*标注

$(document).ready(function(){ $("#article-index").css('display','none');});