pysam可用來處理bam檔案
安裝:
用 pip 或者 conda即可
使用:
pysam的函式有很多,主要的讀取函式有:
一般常用的是第乙個和第二個。
例子:
1bf是乙個迭代器,可以next()或者for讀取import
pysam
23 bf = pysam.alignmentfile("
in.bam
","rb
"); 其中r = read, b:binary. 二進位制檔案。 bam檔案index
1結果:for i in
bf:2
print i.reference_name,i.pos,i.mapq,i.isize
1 ctg000331_np121 144935 27 -284很多功能見下圖:2 ctg000331_np121 144940 48 291
3 ctg000331_np121 144941 48 309
4 ctg000331_np121 144944 48 255
5 ctg000331_np121 144946 27 -370
6 ctg000331_np121 144947 27 -346
** pysam中的座標位點是0開始,染色體起始位置為0,不是1
1一些功能:檢測index檔案是否存在存在即為true## sam 檔案依次對應的12列
2r.qname: reads 名
3r.flag :flag
4r.reference_name: 比對到的染色體
5 r.pos+1: 比對位置,必須得加一
6r.mapq: 比對質量
7r.cigarstring: cigar
8 r.next_reference_name:另外一條reads比對的參考基因組,若和第一條相同,則輸出=
9 r.mpos+1: 比對的位置,必須得加1
10r.isize: 插入片段長度
11r.seq:reads seq
12 r.qual: reads 質量
1用完記得關閉bf.check_index()
2 true
1 bf.close()計算目標區域內比對上的reads數目
1 bf.count(contig="計算目標區域內的覆蓋度。返回1個4維的array,代表acgt的覆蓋度,而每個維度的array長度為100,裡面的數字代表該鹼基在各個位置上的覆蓋度。ctg000331_np121
", start=1, stop=6000)
2 24
1 bf.count_coverage(contig="提取出比對到目標區域內的全部reads。返回的是乙個迭代器,可以通過for迴圈或者next函式從中取出reads,我們使用next()函式取出第一條reads,reads是用 alignedsegment物件表示,可以通過該物件的內建方法再對這條reads進行一些查詢操作。ctg000331_np121
",start=1,stop=100)
1 allreads=bf.fetch(contig="ctg000331_np121
",start=1,stop=10000)
2 是乙個迭代器,可以用for迴圈獲得
1 bf.get_index_statistics()有時候我們並不需要遍歷整乙份bam檔案,我們可能只想獲得區中的某乙個區域(比如chr1中301-310中的資訊),那麼這個時候可以用alignmen模組中的fetch函式:
bam檔案必須要index
參考1、如何使用pysam處理bam
2、使用pysam操作bam檔案
bam檔案讀取 bam格式檔案處理大全 一)
sam檔案是短序列比對生成的檔案,是二代測序中最核心的檔案。在rnaseq,變異檢測等分析中,都需要首先生成sam檔案格式。bam檔案是sam格式的二進位制格式,轉換為二進位制之後,可以減小檔案的儲存。掌握sam bam檔案的操作是處理二代測序資料的非常重要的內容,例如sam與bam的轉換,排序,建...
處理bam檔案提取資訊
一般來說,乙個bam檔案通常只包含乙個樣本的資訊,最多需要進行染色體位置的處理,samtools也提供了簡單的處理方式,比如要提取 chr1的reads,只需要 samtools view input.bam ch1 這幾天遇到了10x genomics的bam結果,發現單細胞的reads全包含在乙...
使用pysam操作VCF BCF檔案
讀取和寫出 from pysam import variantfile bcf in variantfile test in.vcf r bcf out variantfile test out.vcf w header bcf in.header for rec in bcf in.fecth b...