学python:使用python的pysam模块统计bam文件中spliced alignment的reads的数量
2023-06-13 09:16:32 时间
使用igv查看bam文件里有cigar字段,这个是啥意思?
找到了一个解释
https://sites.google.com/site/bioinformaticsremarks/bioinfo/sam-bam-format/what-is-a-cigar
image.png
image.png
所以如果是spliced alignment 的reads cigar关键词中间会有N,只要统计cigar关键词就可以了
python的pysam模块能够统计一个给定区间内所有reads的数量,也可以统计每个reads的一些性质
import pysam
bamfile = 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 argparse
import pysam
import pandas as pd
#from multiprocessing import Pool
parser = 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
最终结果
相关文章
- 快速入门Python机器学习(15)
- python读取文件如何去除空格_python读取txt文件时怎么去掉空格
- Python入门系列(六)一篇学会python函数
- Python入门系列(十)一篇学会python文件处理
- python解析xps文件_xps文件的基本操作
- python判断linux中文件是否存在_Python判断文件是否存在的三种方法
- python中关于命名的例子_Python 命名规范入门实例「建议收藏」
- Python办公实战!按姓名拆分Excel为单独文件,微信自动发给相应联系人
- python输出unicode编码_Python以utf8编码读取文件
- pkl文件是什么_python pkl文件
- python django环境搭建_python的django框架
- Python pandas依列拆分为多个Excel文件
- python判断文件后缀名是否是jpg 或者png_python判断文件名是否包含某字段
- Python将数据写入txt文件_python将内容写入txt文件
- 事件驱动如何理解?什么场景下适合用?Python如何实现一个事件监听器?
- Python 实现子域名查询与爆破
- Python dict字典详解
- Python read()函数:按字节(字符)读取文件
- 用Python 实现的目录拷贝程序详解编程语言
- Linux创建Python文件的步骤(linux新建python文件)
- 使用Python操作MySQL数据库快速上手(python访问mysql数据库)
- python基础入门详解(文件输入/输出内建类型字典操作使用方法)
- python调用cmd复制文件代码分享
- 使用python统计文件行数示例分享