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

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

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!")

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