python：手动比对序列并绘制测序饱和度图片

最近因为工作需要，有一组RNA探针测序数据要求检查其测序饱和度的情况，来评估测序的冗余度。

测序饱和度的评估参考RNA-seq的定义，并非10X定义的根据UMI计算的测序饱和度。

3）饱和度评估 得到sam文件后其实我们可以计算出每一条reads落在哪里，基因间区或者基因区（基因区的话在哪个基因上也可以列出来）。假如我们得到3,200,000条uniq mapping的reads。我们可以采用梯度为10,000，分别为100,000， 200,000， 300,000， 。。。3,100,000, 3,200,000，依次随即抽样，看抽出来的这些reads分别检测到多少基因。然后把reads数和检测的基因数画一个曲线，看看这条曲线是否饱和。若，当reads数到一定程度检测到的基因几乎不增加或者很少增加，则测序包和，否者就是测序量不够没有达到饱和。

由于测序数据是探针数据，并且数量也不是太多，考虑使用python的正则进行序列匹配，实际结果看其比对效率还是挺低的。具体如下：

只以一个测序文件为例，gzip包用于对fastq.gz格式的测序文件进行读取处理。

import gzip
import os
import re
import random
# fastq.gz文件路径，只以一个测序文件为例
os.chdir("F:\\python\\测序饱和度")
fastq_file = [i for i in os.listdir() if i.endswith('gz')]
fastq_file = fastq_file[0]

处理模板序列

模板序列都是12bp的短序列，共有8684个唯一探针序列。每两行是一个探针信息，第一行是以">"开头的探针名称，第二行是具体序列。探针的两行信息以"|分隔合并为一个字符串。为了提高正则匹配的效率，将所有的模板探针序列以逗号分隔并成一个字符串。

# 示例探针
# >RNA0048103
# AAGTGCTCCTGA
ref_seq = []
r = ""
with open('rna.fa') as f:
    for line in f:
        if line.startswith('>'):
            r = line.strip()
        else:
            ref_seq.append("{}|{}".format(r, line.strip()))
            r = ""
ref_seq_str = ",".join(ref_seq)
# 模板序列字符串
ref_seq_str[:50]
#'>RNA0048103|AAGTGCTCCTGA,>RNA0048199|TGCGGAGTTAGA,'

提取测序序列

使用循环处理fastq.gz文件，提取序列至一个列表中做后续处理。

使用n来控制提取序列，遇到@开头的行，则将n标记为1，下一次循环时则提取整行数据，将其置于预先定义的列表中。

# 示例
# @TPNB500160AR:324:H5VCKBGXJ:1:11101:10065:1093 2:N:0:CGTGGCTT+TTAATGAT
# GTGTCACTGTATCGGGCAATCACAAAA
# AAAAAAEEEAEEEEEEEEEEEEEEEE/
fastq_seq = []
n = 0
if(os.path.exists(fastq_file)):
    with gzip.open(fastq_file) as gf:
        for line in gf:
            if line.startswith(b'@'): 
                n = 1
                continue
            if n == 1:
                fastq_seq.append(line.strip().decode()[14:26])
                n = 0
                continue
else:
    print("File {} not exitst!".format(fastq_file))

正则进行序列比对

使用正则进行序列匹配，如果匹配，则返回探针序号，如果没有匹配，则返回字符串“None”。

共有86完条read，比对共运行接近8min，效率比较低，使用常规字符串操作进行序列匹配还是只适用于数据量比较少的情况。

def align(fastq_seq, ref_seq):
    s = []
    for seq in fastq_seq:
        pat = "(RNA[0-9]+)\|{}".format(seq)
        rr = re.search(pat, ref_seq)
        if rr:
            s.append(rr.group(1))
        else:
            s.append("None")
    return s
res_align = align(fastq_seq, ref_seq_str)
res_align[:3]
# ['RNA42588',
# 'RNA46325',
# 'None']

随机抽样获取饱和度数据

饱和度数据其实就是重抽样数据，筛选到比对结果后，去重并计数。

dat = []
# 抽样数量1000 - 10000，1000为步长
for i in range(100000//1000 + 1):
    res = random.sample(res_align, i * 1000)
    l = list(set(res)) # 去重
    l = [i for i in l if i != "None"] # 去除比对失败的探针
    dat.append([i * 1000, len(l)]) # 
# 抽样数量10,0000 - 总reads数，100000为步长
fastq_seq_len = len(fastq_seq) // 100000
for i in range(fastq_seq_len + 1):
    res = random.sample(res_align, i * 100000)
    l = list(set(res))
    l = [i for i in l if i != "None"]
    dat.append([i * 100000, len(l)])
# 查看重抽样的结果
dat[:3]
# [[0, 0], [1000, 784], [2000, 1399]]

绘制散点图

import pandas as pd
import matplotlib.pyplot as plt
pd_dat = pd.DataFrame(dat, columns=['Reads number', 'Detected probe'])