pysam 读取 bam 文件准备 Ribo-seq 质控数据

测试开头

点击上方关注老俊俊

1引言

前面关于 Ribo-seq 的质控我们是直接在 R 里面读取 bam 文件进行操作,我们也可以使用 pysam 处理好数据再导入 R 里面进行分析也可。

2读取 bam 文件

数据还是比对到转录组上:

$ zless 4.rmtRNA-data/ribo-ESC_rmtRNA.fq.gz |
  bowtie ../../index/trans-index/longest_cds_trans 
         -q -p 10 -m 1 -v 2 
         --best --strata - -S |
         samtools sort -@ 10 
         -o 7.map2trans-data/ribo-ESC.sorted.bam

安装 pysam:

linux 系统里面才能装,windows 用户这步就不用考虑了。

!pip install pysam

import pysam
import pandas as pd

读取文件:

samfile = pysam.AlignmentFile("ribo-ESC.sorted.bam", "rb")

info = []

for line in samfile:
    rname = str(line.reference_name)
    if rname == 'None':
        continue
    # split rname
    fileds = rname.split('|')
    cds_start = int(fileds[3])
    cds_end = int(fileds[4])
    # shift aln pos
    pos = int(line.pos + 1)
    # assign frame
    mod = abs(pos - cds_start) % 3
    if mod == 0:
        frame = 'frame 0'
    elif mod == 1:
        frame = 'frame 1'
    elif mod == 2:
        frame = 'frame 2'
    else:
        pass
    # calculate relative to cds_start/stop distance
    dist2start = pos - cds_start
    dist2stop = pos - cds_end
    # read length
    readlen = line.qlen
    # save in list
    info.append([rname,readlen,frame,dist2start,dist2stop])

这里计算了每个 read 的 frame 类型,长度,与 start codon 和 stop codon 的相对距离。

注意 pysam 读进来的比对的位置会减 1,所以位置得再加上 1 才是原来 bam 文件的位置。

转为数据框:

df = pd.DataFrame(info,columns=['rname','readlen','frame','dist2start','dist2stop'])

# save
# df.to_csv('df_res.csv',index=False,header=True)

df

导出来在 R 里面就可以分析了。

3读取 sam 文件

如果想在 windows 里操作,你可以直接将比对结果保存为 sam 文件:

$ zless 4.rmtRNA-data/ribo-ESC_rmtRNA.fq.gz |
  bowtie ../../index/trans-index/longest_cds_trans 
         -q -p 10 -m 1 -v 2 
         --best --strata - -S |
         samtools sort -@ 10 -O SAM 
         -o 7.map2trans-data/ribo-ESC.sorted.sam

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from collections import Counter

处理 sam 文件:

info = []

with open('ribo-ESC.sorted.sam','r') as sam:
    for line in sam:
        if line.startswith('@'):
            continue
        fields = line.split()
        if fields[2] == '*':
            continue
        # split line
        alnpos = int(fields[3])
        startpos = int(fields[2].split('|')[3])
        stoppos = int(fields[2].split('|')[4])
        # assign frame
        mod = abs(alnpos - startpos) % 3
        if mod == 0:
            frame = 'frame 0'
        elif mod == 1:
            frame = 'frame 1'
        elif mod == 2:
            frame = 'frame 2'
        else:
            pass
        # calculate relative to cds_start/stop distance
        dist2start = alnpos - startpos
        dist2stop = alnpos - stoppos
        # read length
        rlen = len(fields[9])
        info.append([fields[2],rlen,frame,dist2start,dist2stop])

转为数据框:

df = pd.DataFrame(info,columns=['rname','rlen','frame','dist2start','dist2stop'])

# save
# df.to_csv('df_res.csv',index=False,header=True)

df

我们可以简单统计一下长度:

tmp_freq = Counter(df.rlen).most_common()
# sort ascending
tmp_freq.sort(key= lambda item:item[0])
# to data frame and save
tmp_freq_df = pd.DataFrame(tmp_freq,columns=['readlen','numbers'])
tmp_freq_df

可视化一下:

sns.set_theme(style="white", context="talk")
# Generate some sequential data
sns.barplot(x=tmp_freq_df.readlen, y=tmp_freq_df.numbers, palette="Paired")