学python：使用python的pysam模块统计bam文件中spliced alignment的reads的数量

学python：使用python的pysam模块统计bam文件中spliced alignment的reads的数量

使用igv查看bam文件里有cigar字段，这个是啥意思？

找到了一个解释

​​alignment 的reads cigar关键词中间会有N，只要统计cigar关键词就可以了

python的pysam模块能够统计一个给定区间内所有reads的数量，也可以统计每个reads的一些性质

import pysambamfile = pysam.AlignmentFile("../barkeRTD/output.split.bam/B1/chr1H_part_1.bam",'rb')reads = bamfile.fetch("chr1H_part_1",102778300,102779978)

reads是一个可以迭代的对象，可以依次访问每个read的情况，read的性质有

image.png

image.png

可以探索的内容很多

结合gtf文件统计每个基因区间内的spliced alignment 的reads的数量

import argparseimport pysamimport pandas as pd#from multiprocessing import Poolparser = argparse.ArgumentParser(description="Stat read orientation")parser.add_argument('-g','--gtf',help="input gtf path")parser.add_argument('-b','--bam',help="input bam path")parser.add_argument('-o','--csv',help="output csv file")args = parser.parse_args()print("Let's go!")selected_col = [0,1,2,3,4,6,8]col_names = ['chromosome','source','feature','start','end','strand','gene_name']df = pd.read_table(args.gtf,sep="\t",comment="#",header=None,usecols=selected_col,names=col_names)df['gene_name'] = df["gene_name"].str.replace("ID=","")#chromo = set(df['chromosome'].tolist())chromo_name = args.bam.split("/")[-1].split(".")[0]Sam = args.bam.split("/")[-2]new_df = df.loc[df['chromosome'] == chromo_name]bamfile = pysam.AlignmentFile(args.bam,'rb')output_df = {'Sample':[], 'chromosome':[], 'gene_name':[], 'start':[], 'end':[], 'strand':[], 'positive':[], 'negative':[], 'positive_spliced':[], 'negative_spliced':[]}for i,j in new_df.iterrows(): positive = 0 negative = 0 negative_spliced = 0 positive_spliced = 0 for read in bamfile.fetch(chromo_name,j.start,j.end): if read.is_read1 and read.is_forward : positive += 1 if 'N' in read.cigarstring: positive_spliced += 1 if read.is_read1 and read.is_reverse: negative += 1 if 'N' in read.cigarstring: negative_spliced += 1 output_df['chromosome'].append(j.chromosome) output_df['gene_name'].append(j.gene_name) output_df['start'].append(j.start) output_df['end'].append(j.end) output_df['strand'].append(j.strand) output_df['positive'].append(positive) output_df['negative'].append(negative) output_df['negative_spliced'].append(negative_spliced) output_df['positive_spliced'].append(positive_spliced) output_df['Sample'].append(Sam)pd.DataFrame(output_df).to_csv(args.csv,index=False)print("Congratulations!")

这里只统计reads1中的spliced alignment

如果是双端测序的数据，pysam统计reads数量的时候会计算为2个分为reads1和reads2

脚本的使用方式

python stat_spliced_junction_read_orientation.py -g input.gtf -b input.bam -o output.csv

最终结果

image.png

小明的数据分析笔记本

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。